estimation of normalisation factors will fail for D690

[KrisThielemans]

Found by @francescaleek. Discussion is at UCL/STIR#838 and UCL/STIR#837.

This needs to be resolved at the STIR level.

[robbietuk]

For the record, the offset=-8 is correct for the D690, or best guess based upon the difference between the unlisted GATE simulation of a 3 point source phantom and the STIR forward projection. Difference between sinograms of unlisted GATE data and STIR forward projection with offset=-8. Blue (negitive) is GATE data and red (positive) is STIR forward projection

This is not to say the ML norm code will work as that is a subtly different issue relating to blocks and buckets.

[KrisThielemans]

indeed. the offset is there to get the rotations sorted for the projections, but it messes up the ML norm sadly (which I hadn't thought about...)

[francescaleek]

Results from GATE sim of D690FOVCylinder, unlisted with offset set to -8 in the first instance, then 0; bucket efficiency factors calculated taking only geometric factors into account.

offset (num of detectors) := -8


offset (num of detectors) := 0

Strangely, I get a warning
WARNING: CListModeDataROOT: I've set the scanner from STIR settings and ignored values in the hroot header.
but it did pick up the change...

[KrisThielemans]

alright! that worked.

Strangely, I get a warning
WARNING: CListModeDataROOT: I've set the scanner from STIR settings and ignored values in the hroot header.
but it did pick up the change...

as you probably didn't change the .hroot aside from the offset, I guess you had this before. Or you changed the template but not the hroot (as you cannot).

Question is now how to move forward with this, aside from fixing STIR... Thinking about that'll be for tomorrow!

[KrisThielemans]

can you also show GATE/eff/model. should be quite flat

[francescaleek]

It is:

[KrisThielemans]

excellent! can you enable the eff_iter now? It shouldn't change anything really.

[francescaleek]

Yes, overall change is small. Some effects are assigned to crystal efficiency, resulting in a lower overall efficiency.

[KrisThielemans]

thanks for checking. would you mind showing the quotient of the total efficiency factors for both cases (as opposed to the difference, as they are multiplicative factors) ?

[francescaleek]

Oops, of course...

[KrisThielemans]

ok. wrong title for the picture I guess.

The overall bias surprises me. Hard to understand that... The rest is likely just noise. (There are more crystal efficiencies to estimate, so therefore the result is noisier).

[robbietuk]

A temporary D690 sinogram header that will allow of the normalisation/eff_factors to be compted as per above:

!INTERFILE  :=
!imaging modality := PT
name of data file := D690FOV_f1g1d0b0.s
!version of keys := STIR4.0
!GENERAL DATA :=
!GENERAL IMAGE DATA :=
!type of data := PET
imagedata byte order := LITTLEENDIAN
!PET STUDY (General) :=
!PET data type := Emission
applied corrections := {None}
!number format := float
!number of bytes per pixel := 4
number of dimensions := 4
matrix axis label [4] := segment
!matrix size [4] := 47
matrix axis label [3] := view
!matrix size [3] := 288
matrix axis label [2] := axial coordinate
!matrix size [2] := { 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,23,22,21,20,19,18,17,16,15,14,13,12,11,10,9,8,7,6,5,4,3,2,1}
matrix axis label [1] := tangential coordinate
!matrix size [1] := 381
minimum ring difference per segment := { -23,-22,-21,-20,-19,-18,-17,-16,-15,-14,-13,-12,-11,-10,-9,-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23}
maximum ring difference per segment := { -23,-22,-21,-20,-19,-18,-17,-16,-15,-14,-13,-12,-11,-10,-9,-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23}
Scanner parameters:= 
Number of rings                          := 24
Number of detectors per ring             := 576
Inner ring diameter (cm)                 := 81.02
Average depth of interaction (cm)        := 0.94
Distance between rings (cm)              := 0.654
Default bin size (cm)                    := 0.21306
View offset (degrees)                    := -5.021
Maximum number of non-arc-corrected bins := 381
Default number of arc-corrected bins     := 331
Energy resolution := 0.124
Reference energy (in keV) := 511
Number of blocks per bucket in transaxial direction         := 1
Number of blocks per bucket in axial direction              := 4
Number of crystals per block in axial direction             := 6
Number of crystals per block in transaxial direction        := 18
Number of detector layers                                   := 1
Number of crystals per singles unit in axial direction      := 1
Number of crystals per singles unit in transaxial direction := 1
end scanner parameters:=
effective central bin size (cm) := 0.226074
number of time frames := 1
number of energy windows := 1
energy window lower level[1] := 425
energy window upper level[1] :=  650
start vertical bed position (mm) := 0
start horizontal bed position (mm) := 0
!END OF INTERFILE :=

N.B.1 This does not work on unlisting due to scanner checks in the STIR unlisting code. Therefore, this temporary modification must be made after the GATE data is unlisted.
N.B.2 The .hroot was modified to include offset (num of detectors) := 0 prior to unlisting.
N.B.3 The originating system := and Scanner type := varibales have been removed and Number of blocks per bucket in transaxial direction := 1 and Number of crystals per block in transaxial direction := 18

[robbietuk]

Back projection of unlisted GATE simulation of a ThreePointActivity phantom with various View offset (degrees) values. N.B. offset (num of detectors) := 0 in the hroot.

View offset (degrees) := 0 . Observed clockwise rotation.

View offset (degrees) := -5.021. Points lay inline (within observable rotation).

View offset (degrees) := +5.021. Observed significantly clockwise roation.

[KrisThielemans]

I guess a fused image with the ground truth could have helped, but I take it that you say that -5.021 fits very well?

Please also check forward projections of the ground truth, ideally also in a very oblique segment.

[robbietuk]

Image based rotation check

See the single blue dots from the true image in the center of each backprojected region. Image values have been rescaled.

[...]but I take it that you say that -5.021 fits very well?

Yes, the two offset voxels with value should be offset in only x and y from the central voxel, generated using https://github.com/UCL/STIR-GATE-Connection/blob/master/ExamplePhantoms/STIRparFiles/ThreePointActivity.par

It might be that it is over rotated in a clockwise direction? Very subtle and difficult to tell on such a low resolution image.

Projections

Low segment number:

Oblique segment:

(Sinogram values have been rescaled)

[KrisThielemans]

Looks pretty good to me!

[robbietuk]

Note to other users. Until ML norm code in STIR is fixed to handle buckets instead of only block (UCL/STIR#838) we cannot compute the norm using the standard D690 header.

However, the work around for the normalisation computation still exist: #59 (comment)
The norm can be computed using this modified D690 header and then re-subsitute the standard D690 header file.

[robbietuk]

I do believe that the randoms estimation (from delayed events) uses the same code as the normalisation sinogram computation.

STIR-GATE-Connection/VoxelisedSimulation/SubScripts/EstimateRandomsFromDelayed.sh

Lines 36 to 43 in ea4a62c

	## find ML normfactors
	echo "find_ML_normfactors3D"
	find_ML_normfactors3D $factors $DelayedData ones.hs 0 $eff_iters
	## mutiply ones with the norm factors to get a sino
	echo "apply_normfactors3D"
	apply_normfactors3D $OutputFilename $factors ones.hs 1 1 $eff_iters

I believe this is evident by the RandomsEstimate sinogram (shown here after SSRB computation) for XCAT data with gaps between detectors:

As I do not understand the ML code enough, I will test the randoms estimation with the modified D690 header and compare.

[robbietuk]

There is a subtle difference between using the modified and "standard" SITR scanners, probably due to noise. Nothing to worry about IMO:

Difference:

Quotient:

[KrisThielemans]

The randoms estimation uses the "crystal efficiency" part of the ML norm code only. (randoms-rate = 2tau singles-rate singles-rate). It should not care about blocks and buckets, or indeed view offset. If it does, that'd be rather weird.

[robbietuk]

Everything looks good. I was concerned it may use some of the symmetries, but it appears it doesnt.

[KrisThielemans]

but results should be identical...

[robbietuk]

The output of:
find_ML_normfactors3D $factors $DelayedData ones.hs 0 $eff_iters
is

${factors}_geo_1.out
${factors}_block_1.out
${factors}_eff_1_1.out
${factors}_eff_1_2.out
...
${factors}_eff_1_10.out

There is no diff between eff factors for the default and modified D690.

However, even though it specifies num_iterations=0 the geo and block infomation is created. Removing the geo and block files before apply_normfactors3D $OutputFilename $factors ones.hs 1 1 $eff_iters is computed results in an error. These files are different between the default and modified D690, snipets below.

Default D690: randoms_factors_block_1.txt

Modified D690: randoms_factors_mod_block_1.txt

[robbietuk]

I am unsure as the the current status of this issue or whether it is still an issue...

[robbietuk]

Closing because I feel we fixed most of the issues here with STIR changes and rotations