first attempt to use axial_effects for ECAT8 normalisation

[KrisThielemans]

closes #624

[KrisThielemans]

@danieldeidda let me know if you cannot push to this branch

[KrisThielemans]

I've put most things in place now I believe, but haven't even checked if it compiled

[danieldeidda]

@danieldeidda let me know if you cannot push to this branch

it seems I cannot push

[KrisThielemans]

@danieldeidda let me know if you cannot push to this branch

it seems I cannot push

seems that I cannot enable write access for you to this branch only. I've therefore given you write access to the repo, but protected a few more branches. let me know. (just be careful!)

[danieldeidda]

@danieldeidda let me know if you cannot push to this branch

it seems I cannot push

seems that I cannot enable write access for you to this branch only. I've therefore given you write access to the repo, but protected a few more branches. let me know. (just be careful!)

probably better to go back once this is done ;)

[danieldeidda]

with these changes it compiles, now i need to try with the data

[KrisThielemans]

we can re-use things from InterfileRawDataHeaderSiemens. It also parses

%axial compression:=11
%maximum ring difference:=49
number of rings:=55
number of energy windows:=1
%energy window lower level (keV) [1]:=435
%energy window upper level (keV) [1]:=650

• Need to change https://github.com/UCL/STIR/blob/4c25191bfc78214de9a1e3d261080493732b3751/src/IO/InterfileHeaderSiemens.cxx#L66 to use "Normalization"
• probably change https://github.com/UCL/STIR/blob/7e0282f0593c0eef157b5e99d8f05eb9504c74dc/src/recon_buildblock/BinNormalisationFromECAT8.cxx#L261-L264 to

#if 1
  InterfileRawDataHeaderSiemens interfile_parser;
  // interfile_parser.ignore_key("data format");
  interfile_parser.parse(filename.c_str());

and remove some lines below, but use interfile_parser.post_processing() instead to do some of the work for you. will need to add the "normalization" case to https://github.com/UCL/STIR/blob/4c25191bfc78214de9a1e3d261080493732b3751/src/IO/InterfileHeaderSiemens.cxx#L258

Might need to add

interfile_parser.ignore_key("%normalization component");
interfile_parser.ignore_key("data offset in bytes");
interfile_parser.ignore_key("number of dimensions");
interfile_parser.ignore_key("%matrix size");
interfile_parser.ignore_key("%matrix axis label");
interfile_parser.ignore_key("%matrix axis unit");
interfile_parser.ignore_key("%scale factor");

If it warns about existing key, should first remove_key.

Note that we could add

interfile_parser.ignore_key("%GLOBAL SCANNER CALIBRATION FACTOR");
interfile_parser.add_key("%scanner quantification factor (Bq*s/ECAT counts)", &calibration_factor)

[KrisThielemans]

yes indeed. a bad case of copy-paste. There are really two different "z" here. What about saying

for (int Siemens_sino_index=0; Siemens_sino_index < this->num_Siemens_sinograms; ++Siemens_sino_index)
{
  int z = Siemens_sino_index;
  ....
}

[danieldeidda]

https://github.com/UCL/STIR/blob/7e0282f0593c0eef157b5e99d8f05eb9504c74dc/src/recon_buildblock/BinNormalisationFromECAT8.cxx#L676

this expects ring1 and ring2 in input. is there a way to get ring1 and ring2 from bin?

[KrisThielemans]

seems https://github.com/UCL/STIR/blob/7e0282f0593c0eef157b5e99d8f05eb9504c74dc/src/recon_buildblock/BinNormalisationFromECAT8.cxx#L676 needs to change to

        const float axial_effect_factor = find_axial_effects(pos1.axial_coord(), pos2.axial_coord());	

[danieldeidda]

using the axial effect we have the following normsino (segment 0) the bottom is without the axial effects

[KrisThielemans]

great. that's looking more realistic! Are we getting closer to the output of e7tools?

[danieldeidda]

great. that's looking more realistic! Are we getting closer to the output of e7tools?

Not exactly! there is also a big difference in scale

[KrisThielemans]

@danieldeidda I've fixed a bug here I believe (last commit), as well as made the norm header parsing much safer/cleaner by moving it to a separate class. Maybe you can try it again on some mCT data?

@ashgillman @Ede1994 this would be of interest to you as well

[KrisThielemans]

@danieldeidda https://github.com/UCL/STIR/pull/631/commits/5c2938bfe3232e152bc56765882dcd341c59876e fixes a bug in BinNormalisation that likely affects you.

[KrisThielemans]

Travis failure in OSX LLVM job due to unrelated https://github.com/UCL/STIR/issues/910.

[KrisThielemans]

Now also fixes #162

[danieldeidda]

@danieldeidda I've fixed a bug here I believe (last commit), as well as made the norm header parsing much safer/cleaner by moving it to a separate class. Maybe you can try it again on some mCT data?

@ashgillman @Ede1994 this would be of interest to you as well

Trying this on the data I used previosly. I just run into few errors: ERROR: KeyParser: the list corresponding to the keyword "data offset in bytes" has to be resized to size 2. This means you have a problem in the keyword values. terminate called after throwing an instance of 'std::__cxx11::basic_string<char, std::char_traits, std::allocator >' Aborted (core dumped) also it does not like vectorised image duration(sec)

[KrisThielemans]

Can you post the .n.hdr here please? (or even post main diffs with an mMR one, e.g. from the NEMA data). Presumably it's different from an mMR one. If you don't have a %number of normalization components keyword, maybe you could just try to add it in for now. should normally be 8 (and occur before any of the other keywords that occur 8 times)

[danieldeidda]

Can you post the .n.hdr here please? (or even post main diffs with an mMR one, e.g. from the NEMA data). Presumably it's different from an mMR one. If you don't have a %number of normalization components keyword, maybe you could just try to add it in for now. should normally be 8 (and occur before any of the other keywords that occur 8 times)

!INTERFILE:= %comment:=CPS ECAT8 sinogram common attributes !originating system:=1104 %SMS-MI header name space:=normalization header %SMS-MI version number:=3.7

!GENERAL DATA:= data description:=PET scanner normalization Coefficients !name of data file:= use_EXPORT_PET_RAW_DATA-norm.n %expiration date (yyyy:mm:dd):=2099:01:01 %expiration time (hh:mm:ss GMT-05:00):=12:00:00

!GENERAL IMAGE DATA:= %study date (yyyy:mm:dd):=2020:06:25 %study time (hh:mm:ss GMT+00:00):=06:48:13 image data byte order:=LITTLEENDIAN !PET data type:=Normalization %data format:=sinogram - expandable from norm components number format:=float number of bytes per pixel:=4

%RAW NORMALIZATION SCANS DESCRIPTION:= %number of normalization scans:=2 %normalization scan[1]:=geometric profile scan %normalization scan[2]:=axial profile scan isotope name[1]:=Ge-68 isotope name[2]:=Ge-68 radiopharmaceutical[1]:=PHARM_TEST radiopharmaceutical[2]:=PHARM_TEST total prompts[1]:=0 total prompts[2]:=606438670 %total randoms[1]:=0 %total randoms[2]:=239382307 image duration (sec) :=7200 %number of buckets:=12 %total uncorrected singles rate[1]:=0 %total uncorrected singles rate[2]:=0

%NORMALIZATION COMPONENTS DESCRIPTION:= %number of normalization components:=7

%normalization component[1]:=geometric effects %normalization component[2]:=crystal interference %normalization component[3]:=crystal efficiencies %normalization component[4]:=axial effects %normalization component[5]:=paralyzing ring DT parameters %normalization component[6]:=non-paralyzing ring DT parameters %normalization component[7]:=TX crystal DT parameter data offset in bytes[1]:=0 data offset in bytes[2]:=174400 data offset in bytes[3]:=196800 data offset in bytes[4]:=344640 data offset in bytes[5]:=347124 data offset in bytes[6]:=347344 data offset in bytes[7]:=347564 number of dimensions[1]:=2 number of dimensions[2]:=2 number of dimensions[3]:=2 number of dimensions[4]:=1 number of dimensions[5]:=1 number of dimensions[6]:=1 number of dimensions[7]:=1 %matrix size[1]:={400,109} %matrix size[2]:={14,400} %matrix size[3]:={672,55} %matrix size[4]:={621} %matrix size[5]:={55} %matrix size[6]:={55} %matrix size[7]:={14} %matrix axis label[1]:={sinogram projection bins,sinogram planes} %matrix axis label[2]:={crystal number,sinogram projection bins} %matrix axis label[3]:={crystal number,ring number} %matrix axis label[4]:={plane number} %matrix axis label[5]:={ring number} %matrix axis label[6]:={ring number} %matrix axis label[7]:={crystal number} %matrix axis unit[1]:={mm/pixel,mm/pixel} %matrix axis unit[2]:={mm/pixel,mm/pixel} %matrix axis unit[3]:={mm/pixel,mm/pixel} %matrix axis unit[4]:={mm/pixel} %matrix axis unit[5]:={mm/pixel} %matrix axis unit[6]:={mm/pixel} %matrix axis unit[7]:={mm/pixel} %scale factor[1]:={2.005,2.027} %scale factor[2]:={2.005,2.005} %scale factor[3]:={2.005,4.054} %scale factor[4]:={2.027} %scale factor[5]:={4.054} %scale factor[6]:={4.054} %scale factor[7]:={2.027}

%axial compression:=11 %maximum ring difference:=49 number of rings:=55 number of energy windows:=1 %energy window lower level (keV) [1]:=435 %energy window upper level (keV) [1]:=650

%GLOBAL SCANNER CALIBRATION FACTOR:= %scanner quantification factor (Bq*s/ECAT counts):=3.12039e+007 %calibration date (yyyy:mm:dd):=2020:06:25 %calibration time (hh:mm:ss GMT+00:00):=06:48:13 %cross calibration factor:=1

%DATA SET DESCRIPTION:= !total number of data sets:=1 %data set[1]:={0,,use_EXPORT_PET_RAW_DATA-norm.n} %number of bad blocks:=0

[danieldeidda]

And this is the diff output:

diff ../../../NEMA_IQ/20170809_NEMA_UCL.n.hdr norm_sinogram/new_norm.n.hdr 1,112c1,109 < !INTERFILE:= < %comment:=CPS ECAT8 sinogram common attributes < !originating system:=2008 < %SMS-MI header name space:=normalization header < %SMS-MI version number:=3.4 < < !GENERAL DATA:= < data description:=PET scanner normalization Coefficients %expiration date (yyyy:mm:dd):=2099:01:01 < %expiration time (hh:mm:ss GMT-05:00):=12:00:00 < < !GENERAL IMAGE DATA:= < %study date (yyyy:mm:dd):=2017:08:09 < %study time (hh:mm:ss GMT+00:00):=07:23:33 < image data byte order:=LITTLEENDIAN < !PET data type:=normalization < %data format:=sinogram - expandable from norm components < number format:=float < number of bytes per pixel:=4 < < %RAW NORMALIZATION SCANS DESCRIPTION:= < %number of normalization scans:=2 < %normalization scan [1]:=geometric profile scan < %normalization scan [2]:=axial profile scan < isotope name [1]:=Ge-68 < isotope name [2]:=Ge-68 < radiopharmaceutical [1]:=PHARM_TEST < radiopharmaceutical [2]:=PHARM_TEST < total prompts [1]:=0 < total prompts [2]:=606438670 < %total randoms [1]:=0 < %total randoms [2]:=239382307 < image duration (sec) [1]:=0 < image duration (sec) [2]:=7200 < %number of buckets:=12 < %total uncorrected singles rate [1]:=0 < %total uncorrected singles rate [2]:=0 < < %NORMALIZATION COMPONENTS DESCRIPTION:= < %number of normalization components:=8 < %normalization component [1]:=geometric effects < %normalization component [2]:=crystal interference < %normalization component [3]:=crystal efficiencies < %normalization component [4]:=axial effects < %normalization component [5]:=paralyzing ring DT parameters < %normalization component [6]:=non-paralyzing ring DT parameters < %normalization component [7]:=TX crystal DT parameter < %normalization component [8]:=additional axial effects < data offset in bytes [1]:=0 < data offset in bytes [2]:=174752 < data offset in bytes [3]:=187136 < data offset in bytes [4]:=316160 < data offset in bytes [5]:=319508 < data offset in bytes [6]:=319764 < data offset in bytes [7]:=320020 < data offset in bytes [8]:=320056 < number of dimensions [1]:=2 < number of dimensions [2]:=2 < number of dimensions [3]:=2 < number of dimensions [4]:=1 < number of dimensions [5]:=1 < number of dimensions [6]:=1 < number of dimensions [7]:=1 < number of dimensions [8]:=1 < %matrix size [1]:={344,127} < %matrix size [2]:={9,344} < %matrix size [3]:={504,64} < %matrix size [4]:={837} < %matrix size [5]:={64} < %matrix size [6]:={64} < %matrix size [7]:={9} < %matrix size [8]:={837} < %matrix axis label [1]:={sinogram projection bins,sinogram planes} < %matrix axis label [2]:={crystal number,sinogram projection bins} < %matrix axis label [3]:={crystal number,ring number} < %matrix axis label [4]:={plane number} < %matrix axis label [5]:={ring number} < %matrix axis label [6]:={ring number} < %matrix axis label [7]:={crystal number} < %matrix axis label [8]:={plane number} < %matrix axis unit [1]:={mm/pixel,mm/pixel} < %matrix axis unit [2]:={mm/pixel,mm/pixel} < %matrix axis unit [3]:={mm/pixel,mm/pixel} < %matrix axis unit [4]:={mm/pixel} < %matrix axis unit [5]:={mm/pixel} < %matrix axis unit [6]:={mm/pixel} < %matrix axis unit [7]:={mm/pixel} < %matrix axis unit [8]:={mm/pixel} < %scale factor [1]:={2.0445,2.03125} < %scale factor [2]:={2.0445,2.0445} < %scale factor [3]:={2.0445,4.0625} < %scale factor [4]:={2.03125} < %scale factor [5]:={4.0625} < %scale factor [6]:={4.0625} < %scale factor [7]:={2.03125} < %scale factor [8]:={2.03125} < %axial compression:=11 < %maximum ring difference:=60 < number of rings:=64 < number of energy windows:=1 < %energy window lower level (keV) [1]:=430 < %energy window upper level (keV) [1]:=610 < < %GLOBAL SCANNER CALIBRATION FACTOR:= < %scanner quantification factor (Bq*s/ECAT counts):=1.98459e+007 < %calibration date (yyyy:mm:dd):=2017:08:09 < %calibration time (hh:mm:ss GMT+00:00):=07:23:33 < %cross calibration factor:=0.96 < < %DATA SET DESCRIPTION:= < !total number of data sets:=1 < %data set [1]:={0,,20170809_NEMA_UCL.n}

!INTERFILE:= %comment:=CPS ECAT8 sinogram common attributes !originating system:=1104 %SMS-MI header name space:=normalization header %SMS-MI version number:=3.7

!GENERAL DATA:= data description:=PET scanner normalization Coefficients !name of data file:= use_EXPORT_PET_RAW_DATA-norm.n %expiration date (yyyy:mm:dd):=2099:01:01 %expiration time (hh:mm:ss GMT-05:00):=12:00:00

!GENERAL IMAGE DATA:= %study date (yyyy:mm:dd):=2020:06:25 %study time (hh:mm:ss GMT+00:00):=06:48:13 image data byte order:=LITTLEENDIAN !PET data type:=Normalization %data format:=sinogram - expandable from norm components number format:=float number of bytes per pixel:=4

%RAW NORMALIZATION SCANS DESCRIPTION:= %number of normalization scans:=2 %normalization scan[1]:=geometric profile scan %normalization scan[2]:=axial profile scan isotope name[1]:=Ge-68 isotope name[2]:=Ge-68 radiopharmaceutical[1]:=PHARM_TEST radiopharmaceutical[2]:=PHARM_TEST total prompts[1]:=0 total prompts[2]:=606438670 %total randoms[1]:=0 %total randoms[2]:=239382307 image duration (sec) :=7200 %number of buckets:=12 %total uncorrected singles rate[1]:=0 %total uncorrected singles rate[2]:=0

%NORMALIZATION COMPONENTS DESCRIPTION:= %number of normalization components:=7

%normalization component[1]:=geometric effects %normalization component[2]:=crystal interference %normalization component[3]:=crystal efficiencies %normalization component[4]:=axial effects %normalization component[5]:=paralyzing ring DT parameters %normalization component[6]:=non-paralyzing ring DT parameters %normalization component[7]:=TX crystal DT parameter data offset in bytes[1]:=0 data offset in bytes[2]:=174400 data offset in bytes[3]:=196800 data offset in bytes[4]:=344640 data offset in bytes[5]:=347124 data offset in bytes[6]:=347344 data offset in bytes[7]:=347564 number of dimensions[1]:=2 number of dimensions[2]:=2 number of dimensions[3]:=2 number of dimensions[4]:=1 number of dimensions[5]:=1 number of dimensions[6]:=1 number of dimensions[7]:=1 %matrix size[1]:={400,109} %matrix size[2]:={14,400} %matrix size[3]:={672,55} %matrix size[4]:={621} %matrix size[5]:={55} %matrix size[6]:={55} %matrix size[7]:={14} %matrix axis label[1]:={sinogram projection bins,sinogram planes} %matrix axis label[2]:={crystal number,sinogram projection bins} %matrix axis label[3]:={crystal number,ring number} %matrix axis label[4]:={plane number} %matrix axis label[5]:={ring number} %matrix axis label[6]:={ring number} %matrix axis label[7]:={crystal number} %matrix axis unit[1]:={mm/pixel,mm/pixel} %matrix axis unit[2]:={mm/pixel,mm/pixel} %matrix axis unit[3]:={mm/pixel,mm/pixel} %matrix axis unit[4]:={mm/pixel} %matrix axis unit[5]:={mm/pixel} %matrix axis unit[6]:={mm/pixel} %matrix axis unit[7]:={mm/pixel} %scale factor[1]:={2.005,2.027} %scale factor[2]:={2.005,2.005} %scale factor[3]:={2.005,4.054} %scale factor[4]:={2.027} %scale factor[5]:={4.054} %scale factor[6]:={4.054} %scale factor[7]:={2.027}

%axial compression:=11 %maximum ring difference:=49 number of rings:=55 number of energy windows:=1 %energy window lower level (keV) [1]:=435 %energy window upper level (keV) [1]:=650

%GLOBAL SCANNER CALIBRATION FACTOR:= %scanner quantification factor (Bq*s/ECAT counts):=3.12039e+007 %calibration date (yyyy:mm:dd):=2020:06:25 %calibration time (hh:mm:ss GMT+00:00):=06:48:13 %cross calibration factor:=1

%DATA SET DESCRIPTION:= !total number of data sets:=1 %data set[1]:={0,,use_EXPORT_PET_RAW_DATA-norm.n} %number of bad blocks:=0

[KrisThielemans]

@danieldeidda I tried to with the mCT header you've pasted above and it parses fine. I wonder if you used an old executable by mistake.

I don't have the .n data to check the output so I created a .n file which is just 1 (which should have the same effect as setting all use fields to 0). With that I get as norm factors (heavily thresholded to account for infinity in the gaps)

Of course, I guess that Siemens still uses gap filling for the mCT, so the output of e7tools norm should be compared to STIR output with use_gaps:=0 (.par file provided).

There's also a sample .par file that switches off using axial effects. When using that, the result should be identical to what we had before this PR.

[danieldeidda]

sorry I was actually using the wrong parfile. This is what I get:

top=e7tools bottom=STIR

[danieldeidda]

This is with gap set to 1: this could be potentially inverted?

[danieldeidda]

and this is with inverted e7 tools using the same color scale

[KrisThielemans]

@NicoleJurjew can you please check if this (still) works? Ideally also add a plot illustrating axial uniformity of cylinder data before and after. (something like an ROI value per plane)

@danieldeidda we found out we had to divide by axial_effects as opposed to multiply. If you have time, could you try the mCT again? (I will merge on Friday or so, but we can create a new mCT PR of course)

[NicoleJurjew]

@KrisThielemans I re-built just now and everything still works.

I've added a plot showing the average over a 20x20 pixel ROI in the transaxial center for each plane. Let me know if you'd like to see more data.

[danieldeidda]

@NicoleJurjew can you please check if this (still) works? Ideally also add a plot illustrating axial uniformity of cylinder data before and after. (something like an ROI value per plane)

@danieldeidda we found out we had to divide by axial_effects as opposed to multiply. If you have time, could you try the mCT again? (I will merge on Friday or so, but we can create a new mCT PR of course)

This is what I get. The bottom one is with the latest ECT8axial_effects

[KrisThielemans]

sorry. what the top and middles ones again? (it's a bit long ago). Also did you set "use gaps:=0"? Otherwise we cannot compare to the output of e7tools

[KrisThielemans]

Thanks @NicoleJurjew . It's flatter (less block structure) but clearly not flat. An improvement, but some way to go. Not so clear why now. Probably not purely norm related.

However, I did find a problem in the code that will make a ~1% difference in the oblique segments (which we didn't check yet I guess). I don't think this will make any appreciable difference on the "dip" but I'll push a new version soon.