(potential) axial shift in scatter simulations

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.
https://github.com/UCL/STIR/blob/29133a3ffa690b444894b4392aa710e8bb33a42e/src/recon_buildblock/ProjMatrixByBinUsingRayTracing.cxx#L655

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 https://github.com/UCL/STIR/blob/29133a3ffa690b444894b4392aa710e8bb33a42e/src/scatter_buildblock/sample_scatter_points.cxx#L65-L67 and https://github.com/UCL/STIR/blob/29133a3ffa690b444894b4392aa710e8bb33a42e/src/scatter_buildblock/single_scatter_integrals.cxx#L100-L102

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 https://github.com/UCL/STIR/blob/29133a3ffa690b444894b4392aa710e8bb33a42e/scripts/zoom_att_image.sh#L28

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.