Skip to content

Commit

Permalink
Merge pull request #631 from UCL/ECAT8axial_effects
Browse files Browse the repository at this point in the history
use axial_effects for ECAT8 normalisation
  • Loading branch information
KrisThielemans authored Sep 30, 2023
2 parents a6efef6 + b7d69e6 commit 131216e
Show file tree
Hide file tree
Showing 10 changed files with 430 additions and 157 deletions.
15 changes: 14 additions & 1 deletion documentation/release_5.2.htm
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@
<body>
<h1>Summary of changes in STIR release 5.2</h1>

<p>This version is 100% backwards compatible with STIR 5.0. However, there is a <b> change in the output of scatter estimation<b>, see below for more information.
<p>This version is 100% backwards compatible with STIR 5.0 as far as usage goes. However, there are <b> changes in the output of scatter estimation and ECAT8 normalisation<b>, see below for more information.
</p>
<h2>Overall summary</h2>

Expand All @@ -30,6 +30,19 @@ <h3>Bug fixes</h3>
<li>Setting SPECTUB resolution model with STIR python or SIRF divided slope by 10 in error. The problem did not occur when set using parameter file</li>
</ul>

<h3>Changed functionality</h3>
<ul>
<li> The ECAT8 normalisation (used for the Siemens mMR) code now takes the 4th component <i>axial effects</i> into account.
These normalisation factors are therefore different (even up to ~10%). This gives improved axial uniformity in the images.
The use of the axial effects can be switched off by adding setting <tt>use_axial_effects_factors:=0</tt> to the
parameter file (see an example in <tt>examples/Siemens-mMR/correct_projdata_no_axial_effects.par</tt>), or the class member
of the same name.<br />
In addition, the Siemens normalisation header is now read (using a new class <code>InterfileNormHeaderSiemens</code>) such that
hard-coded variables for the Siemens mMR have been removed. Further testing of this functionality is still required however.
<br /><a href="https://github.com/UCL/STIR/pull/1182/">PR #1182</a>.
</li>
</ul>

<h3>New functionality</h3>
<ul>
<li>The <tt>Discretised Shape3D</tt> shape/ROI has now an extra value <tt>label index</tt>. For ROIs, this allows
Expand Down
3 changes: 3 additions & 0 deletions examples/Siemens-mMR/correct_projdata.par
Original file line number Diff line number Diff line change
Expand Up @@ -40,6 +40,9 @@ correct_projdata Parameters :=
Bin Normalisation type := from ecat8
Bin Normalisation From ecat8 :=
normalisation filename:= ${ECATNORM}
; keyword that can be used to write the components to a separate text files for debugging
; files are written in the current directory and are called geom_out.txt etc.
; write_components_to_file := 0
End Bin Normalisation From ecat8:=

; attenuation image, will be forward projected to get attenuation factors
Expand Down
62 changes: 62 additions & 0 deletions examples/Siemens-mMR/correct_projdata_no_axial_effects.par
Original file line number Diff line number Diff line change
@@ -0,0 +1,62 @@
correct_projdata Parameters :=
; a sample file that switches the use of the "axial effects"
; (the 4th component in the ECAT8 norm file) off
input file := ${INPUT}

; Current way of specifying time frames, pending modifications to
; STIR to read time info from the headers
; see class documentation for TimeFrameDefinitions for the format of this file
; time frame definition filename := frames.fdef

; if a frame definition file is specified, you can say that the input data
; corresponds to a specific time frame
; the number should be between 1 and num_frames and defaults to 1
; this is currently only used to pass the relevant time to the normalisation
; time frame number := 1

; output file
; for future compatibility, do not use an extension in the name of the
; output file. It will be added automatically
output filename := ${OUTPUT}

; default value for next is -1, meaning 'all segments'
; maximum absolute segment number to process :=


; use data in the input file, or substitute data with all 1's
; (useful to get correction factors only)
; default is '1'
use data (1) or set to one (0) := 0

; precorrect data, or undo precorrection
; default is '1'
; apply (1) or undo (0) correction :=

; parameters specifying correction factors
; if no value is given, the corresponding correction will not be performed

; random coincidences estimate, subtracted before anything else is done
;randoms projdata filename := random.hs
; normalisation (or binwise multiplication, so can contain attenuation factors as well)
Bin Normalisation type := from ecat8
Bin Normalisation From ecat8 :=
normalisation filename:= ${ECATNORM}
use_axial_effects_factors:=0
End Bin Normalisation From ecat8:=

; attenuation image, will be forward projected to get attenuation factors
; OBSOLETE
;attenuation image filename := attenuation_image.hv

; forward projector used to estimate attenuation factors, defaults to Ray Tracing
; OBSOLETE
;forward_projector type := Ray Tracing

; scatter term to be subtracted AFTER norm+atten correction
; defaults to 0
;scatter projdata filename := scatter.hs

; to interpolate to uniform sampling in 's', set value to 1
; arc correction := 0

END:=
5 changes: 4 additions & 1 deletion examples/Siemens-mMR/correct_projdata_only_counts.par
Original file line number Diff line number Diff line change
@@ -1,4 +1,6 @@
correct_projdata Parameters :=
; a par file to investigate what BinNormalisationFromECAT8 does when only
; taking span/mashing into account (i.e. actually ignoring the norm factors themselves)

input file := ${INPUT}

Expand Down Expand Up @@ -44,7 +46,8 @@ correct_projdata Parameters :=
use_detector_efficiencies:=0
use_geometric_factors:=0
use_crystal_interference_factors:=0
End Bin Normalisation From ecat8:=
use_axial_effects_factors:=0
End Bin Normalisation From ecat8:=

; attenuation image, will be forward projected to get attenuation factors
; OBSOLETE
Expand Down
6 changes: 0 additions & 6 deletions src/IO/InterfileHeader.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -306,12 +306,6 @@ bool InterfileHeader::post_processing()
warning("Interfile error: no matrix size keywords present\n");
return true;
}
if (matrix_size[matrix_size.size()-1].size()!=1)
{
warning("Interfile error: last dimension (%d) of 'matrix size' cannot be a list of numbers\n",
matrix_size[matrix_size.size()-1].size());
return true;
}
for (unsigned int dim=0; dim != matrix_size.size(); ++dim)
{
if (matrix_size[dim].size()==0)
Expand Down
142 changes: 130 additions & 12 deletions src/IO/InterfileHeaderSiemens.cxx
Original file line number Diff line number Diff line change
@@ -1,12 +1,9 @@
/*
Copyright (C) 2000 PARAPET partners
Copyright (C) 2000 - 2009-04-30, Hammersmith Imanet Ltd
Copyright (C) 2011-07-01 - 2012, Kris Thielemans
Copyright (C) 2013, 2016, 2018, 2020 University College London
Copyright (C) 2018, 2020, 2021, 2023 University College London
Copyright (C) 2018 STFC
This file is part of STIR.
SPDX-License-Identifier: Apache-2.0 AND License-ref-PARAPET-license
SPDX-License-Identifier: Apache-2.0
See STIR/LICENSE.txt for details
*/
Expand All @@ -29,6 +26,7 @@
#include "stir/IO/stir_ecat_common.h"
#include <numeric>
#include <functional>
#include "stir/stream.h"
#include "stir/warning.h"
#include "stir/error.h"

Expand All @@ -48,15 +46,17 @@ InterfileHeaderSiemens::InterfileHeaderSiemens()
{
// always PET
exam_info_sptr->imaging_modality = ImagingModality::PT;


byte_order_values.clear();
byte_order_values.push_back("LITTLEENDIAN");
byte_order_values.push_back("BIGENDIAN");


PET_data_type_values.clear();
PET_data_type_values.push_back("Emission");
PET_data_type_values.push_back("Transmission");
PET_data_type_values.push_back("Blank");
PET_data_type_values.push_back("AttenuationCorrection");
PET_data_type_values.push_back("Normalisation");
PET_data_type_values.push_back("Normalization");
PET_data_type_values.push_back("Image");

for (int patient_position_idx = 0; patient_position_idx <= PatientPosition::unknown_position; ++patient_position_idx)
Expand Down Expand Up @@ -169,7 +169,7 @@ InterfileRawDataHeaderSiemens::InterfileRawDataHeaderSiemens()
: InterfileHeaderSiemens()
{
// first set to some crazy values
num_segments = -1;
num_segments = 0;
num_rings = -1;
maximum_ring_difference = -1;
axial_compression = -1;
Expand Down Expand Up @@ -247,8 +247,8 @@ bool InterfileRawDataHeaderSiemens::post_processing()

const std::string PET_data_type =
standardise_interfile_keyword(PET_data_type_values[PET_data_type_index]);
if (PET_data_type != "emission" && PET_data_type != "transmission")
{ error("Interfile error: expecting emission or transmission for 'PET data type'"); }
if (PET_data_type != "emission" && PET_data_type != "transmission" && PET_data_type != "normalization")
{ error("Interfile error: expecting 'emission' or 'transmission' or 'normalization' for 'PET data type'"); }

// handle scanner

Expand Down Expand Up @@ -277,7 +277,7 @@ bool InterfileRawDataHeaderSiemens::post_processing()
error("Interfile warning: 'number of segments' and length of 'segment table' are not consistent");
}
segment_sequence = ecat::find_segment_sequence(*data_info_ptr);
//XXX check if order here and segment_table are consistent
//TODO check if order here and segment_table are consistent
}

// Set the bed position
Expand Down Expand Up @@ -530,5 +530,123 @@ bool InterfileListmodeHeaderSiemens::post_processing()
return false;
}

InterfileNormHeaderSiemens::InterfileNormHeaderSiemens()
: InterfileRawDataHeaderSiemens()
{
// some defaults
calib_factor = 1.F;
cross_calib_factor = 1.F;
num_buckets = 0; // should be set normally
num_components = 0; // should be set to 8 normally
axial_compression = 11; // should be set normally but seems to be this always
is_arccorrected = false; // norm data is never arc-corrected

ignore_key("data description");
ignore_key("%expiration date (yyyy:mm:dd)");
ignore_key("%expiration time (hh:mm:ss GMT-05:00)");
// currently keywords are truncated at :
ignore_key("%expiration time (hh");
ignore_key("%expiration date (yyyy");
ignore_key("%raw normalization scans description");

// remove some standard keys, which Siemens has replaced with similar names
remove_key("matrix size");
remove_key("matrix axis label");
remove_key("scaling factor (mm/pixel)");

// keywords for the components
add_key("%number of normalization components",
KeyArgument::INT, (KeywordProcessor)&InterfileNormHeaderSiemens::read_num_components, &num_components);
add_vectorised_key("%matrix size", &matrix_size);
add_vectorised_key("%matrix axis label", &matrix_labels);
ignore_key("%matrix axis unit");
ignore_key("%normalization component");
ignore_key("%normalization components description");
add_vectorised_key("data offset in bytes", &data_offset_each_dataset);
remove_key("number of dimensions");
add_vectorised_key("number of dimensions", &number_of_dimensions);
ignore_key("%scale factor");

// other things
add_key("%number of buckets", &num_buckets);

ignore_key("%global scanner calibration factor");
add_key("%scanner quantification factor (Bq*s/ECAT counts)",& calib_factor);
add_key("%cross calibration factor",& cross_calib_factor);
ignore_key("%calibration date (yyyy:mm:dd)");
ignore_key("%calibration time (hh:mm:ss GMT+00:00)");
// currently keywords are truncated at :
ignore_key("%calibration time (hh");
ignore_key("%calibration date (yyyy");

// isotope things are vectorised in norm files and not in other raw data, so we could
// fix that, but as we are not interested in it anyway (tends to be Ge-68), let's just ignore it.
remove_key("isotope name");
ignore_key("isotope name");

ignore_key("%number of normalization scans");
ignore_key("%normalization scan");
remove_key("image duration (sec)");
ignore_key("image duration (sec)");
#if 0
// change to vectorised key
// would need to set image_durations length from "number of normalization scans"
add_vectorised_key("image duration (sec)", &image_durations);
#endif
ignore_key("%data format");
ignore_key("%data set description");
ignore_key("total number of data sets");
ignore_key("%data set");
}

void InterfileNormHeaderSiemens::read_num_components()
{
set_variable();

matrix_labels.resize(num_components);
matrix_size.resize(num_components);
// matrix_axis_units.resize(num_components);
// matrix_axis_units.resize(num_components);
// pixel_sizes.resize(num_components);
// normalization_components.resize(num_components);
data_offset_each_dataset.resize(num_components);
number_of_dimensions.resize(num_components);
}

bool InterfileNormHeaderSiemens::post_processing()
{
if (matrix_size.size() < 4)
error("Error parsing ECAT8 norm file header: '%number of normalization components='" +
std::to_string(matrix_size.size())+
"' but should be at least 4");
// %normalization component [1]:=geometric effects
if (matrix_size[0].size() != 2)
error("Error parsing ECAT8 norm file header: '%matrix size[1]' should have length 2");
// %normalization component [3]:=crystal efficiencies
if (matrix_size[2].size() != 2)
error("Error parsing ECAT8 norm file header: '%matrix size[3]' should have length 2");

// TODO should do far more checks...

// remove trailing \r (or other white space) occuring in mMR norm files (they sometimes have \r\r at end of line)
std::string s=exam_info_sptr->originating_system;
s.erase( std::remove_if( s.begin(), s.end(), isspace ), s.end() );
exam_info_sptr->originating_system=s;
s=data_file_name;
s.erase( std::remove_if( s.begin(), s.end(), isspace ), s.end() );
data_file_name=s;

// norm headers don't seem to have "number of views". We need to get it from elsewhere...
// The crystal efficiencies have as first dimension the number of crystals, so let's use that.
const int num_detectors_per_ring = matrix_size[2][0];
this->num_views = num_detectors_per_ring/2;
// find num_bins from geometric effects
this->num_bins = matrix_size[0][0];

if (InterfileRawDataHeaderSiemens::post_processing() == true)
return true;

return false;
}

END_NAMESPACE_STIR
Loading

0 comments on commit 131216e

Please sign in to comment.