(potential) axial shift in scatter simulations

[KrisThielemans]

While reviewing #44, I discovered a problem in the original scatter code, existing probably since 2005.

There are 3 images passed to the scatter simulation code: the activity image, the attenuation image and the "zoomed" attenuation image used for finding the scatter points. The issue occurs when any of these images do not have the same "size" in z-direction as the original data (see below). However, this is in practice always the case, especially for the "scatter-point" image.

Coordinate systems

Projectors

For historical reasons, STIR projectors currently "centre" the image with respect to the scanner, i.e. it puts 0+(max_index.z()+min_index.z())/2.F * voxel_size.z()as the centre (i.e. the image is centred if origin.z()==0), see e.g.

STIR/src/recon_buildblock/ProjMatrixByBinUsingRayTracing.cxx

Line 655 in 29133a3

+(max_index.z()+min_index.z())/2.F * voxel_size.z();

This convention causes/d major headaches to @ashgillman in https://github.com/UCL/STIR/projects/1. It is dangerous, e.g. when using generate_image, as the location w.r.t. the scanner will change depending on the number of planes. Hence the recommendation to always use the "standard" number of planes.

Scatter

This attempts to follow the same convention. However, the scatter simulation code uses downsampled images...

In the scatter code, relevant lines are

STIR/src/scatter_buildblock/sample_scatter_points.cxx

Lines 65 to 67 in 29133a3

	const float z_to_middle =
	(image.get_max_index() + image.get_min_index())*voxel_size.z()/2.F;
	origin.z() -= z_to_middle;

and

STIR/src/scatter_buildblock/single_scatter_integrals.cxx

Lines 100 to 102 in 29133a3

	const float z_to_middle =
	(image.get_max_index() + image.get_min_index())*voxel_size.z()/2.F;
	origin.z() -= z_to_middle;

Side note: the scatter code also uses coordinates for the detectors (stored in detector_points_vector) which are in the "scanner-centred" coordinate system (used by ProjDataInfo::get_m() etc). (By the way, it still does this by using the out-dated find_cartesian_coordinates_of_detection but shifting to the centre, see here and here for the shift. See #105. This is however not a problem here).

The problem

The downsampling is generally computed using zoom_image, which makes sure that the "physical" coordinate (origin + index*voxel_size) is consistent before and after zooming, see here.

However, the z_to_middle shift depends on the number of planes and voxel-size.
These 2 conventions are incompatible.

Example

• original image has 2N-1 planes of voxel-size Z, 0 offset. This is the usual arrangement for a scanner with N rings and ring-spacing 2Z (i.e. length 2NZ).
• zoomed image has 7 planes, "covering" the whole axial FOV 2NZ, so new zoomed voxel-size ZZ=2NZ/7 (corresponding to zoom_z=7/(2N)).

In the original image the middle plane corresponds to the middle of the scanner, which has the "physical coordinate" (N-1)Z (as 0 is in centre of first plane of the original image). Indeed, the z_to_middle=(max_index.z()+min_index.z())/2.F * voxel_size.z() is in the original image (2N-2)Z/2=(N-1)Z.However, in the zoomed image 6ZZ=6*2NZ/7/2. As zoom_image will have made the physical_coordinate consistent between the 2 images, we would get the following calculations for the physical coordinate (N-1)Z for the centred coordinate system

• original image: (N-1)Z - z_to_middle = 0
• zoomed image: (N-1)Z - 6N/7Z = (N/7-1)Z

Conclusion

The scatter-code effectively shifts each of the 3 images (with potentially different amount). From the above, I believe that the shift is

 z_to_middle_org - z_to_middle_zoomed = 
    ((org_size_z-1)/2 - (new_size_z-1)/2/zoom_z) * org_voxel_size_z

Note that this conclusion seems independent of whatever offset we ask zoom_image to use.
Our current zooming strategy assumes a zero new offset

STIR/scripts/zoom_att_image.sh

Line 28 in 29133a3

offset_z=${5:-0}

In #44, I've made a test-case with a centred object, testing if the output is symmetric. I had to work somewhat hard to make the test succeed. The current code tests scatter-point image zoom_z = 1/2, new_size_z = N, which indeed gives a 0 shift according to the formula above.

[KrisThielemans]

@NikEfth possibly you saw something on this when you said that down-sampled images had to use the same downsampling. After writing this up, I think however that the problem is worse than that.

[KrisThielemans]

See some background reading at http://stir.sourceforge.net/wiki/index.php/STIR_FAQ#How_does_the_STIR_coordinate_system_work_.28e.g._for_generate_image.29

[NikEfth]

@KrisThielemans
Very well described.
This is the reason why in earlier version of the scatter estimation I was passing the zoomed images as image templates in the reconstruction process.
I was not sure how deep I was supposed to go with this. I always thought that it was a problem in zoom and should work around it.

[KrisThielemans]

Requiring the shift to be zero leads to

(org_size_z-1) = (new_size_z-1)/zoom_z

This works ok in the "symmetric" test. I've tried the following

  shared_ptr<VoxelsOnCartesianGrid<float> > tmpl_density( new VoxelsOnCartesianGrid<float>(*original_projdata_info));
  ...
  {
      const int new_size_z = 14; 
      const float zoom_z = static_cast<float>(new_size_z-1)/(water_density->size()-1);
      sss->downsample_density_image_for_scatter_points(.2,zoom_z,-1, new_size_z);
    }

with a few different values of new_size_z. All good.

It's now a question on how to handle this. Seems we can throw some errors if this condition isn't satisfied. However, we need to check on all 3 images. Needs a bit of thought.

[KrisThielemans]

I think we need to do the following check for each image

const float z_to_middle = 
   (image.get_max_index() + image.get_min_index())*voxel_size.z()/2.F; 
const float z_to_middle_standard =
   (scanner.get_num_rings()-1) * scanner.get_ring_distance()/2;
if (abs(z_to_middle - z_to_middle_standard)>.1)
  error("...");

i.e. check that it is zoomed consistently with the scanner (see formula above for the "standard" image for a scanner (N-1)Z). It seems that this needs to be checked with the downsampled scanner (note that z_to_middle_standard will change if you change the scanner but keep the same axial FOV. horrible!)

[KrisThielemans]

The above suggested test was too strict. Due to the "shift to the middle" strategy, the current code gives the same result with a global shift of all 3 images (which is admittedly weird). So, instead I've introduced a check in 06e2604 that sees if all 3 images are consistent. See here

This test seems to catch anomalous cases while letting ones through that do satisfy symmetry tests. Checks add in this commit

[KrisThielemans]

This has been resolved in the above commits.

[ashgillman]

Can this function go?

STIR/src/scatter_buildblock/ScatterSimulation.cxx

Line 402 in abd5992

check_z_to_middle_consistent(const DiscretisedDensity<3,float>& _image, const std::string& name) const

after #618