(potential) axial shift in scatter simulations

KrisThielemans commented 4 years ago

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.

KrisThielemans commented 4 years ago

@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 commented 4 years ago

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 commented 4 years ago

@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 commented 4 years ago

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 commented 4 years ago

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 commented 4 years ago

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 https://github.com/UCL/STIR/pull/44/commits/06e260458f47367fdd3bd75a29ee8f9ee436e83a 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 commented 4 years ago

This has been resolved in the above commits.

ashgillman commented 3 years ago

Can this function go? https://github.com/UCL/STIR/blob/abd59922fcbe0398a46a688f43b57c0548cf39bd/src/scatter_buildblock/ScatterSimulation.cxx#L402 after #618

markus-jehl commented 1 year ago

Is this still an issue?

danieldeidda commented 1 year ago

actually not sure but it works for blocks on cylindrical

KrisThielemans commented 1 year ago

When the code was written originally, the shift could occur. As I wrote above, I added in checks such that if the user would give images with dimensions such that the result would be wrong, it would throw an error, see here. The issue remains open as #618 should remove this restriction on image dimensions, and then we can add some tests and remove these checks.

So, in the current code, as long as the scatter simulation doesn't throw, it will be fine.

UCL / STIR