Reinvert spect axis

documentation/release_6.3.htm

@@ -143,6 +143,22 @@ <h4>Changes to examples</h4>
 
   <h3>Bug fixes</h3>
   <ul>
+    <li>
+      Fixed a bug in the SPECTUB matrix (introduced in 4.0) that effectively inverted the "x" (horizontal in-plane) axis for the
+      emission image as part of the projection. This resulted in the need for an inconsistent attenuation
+      image orientation compared to the reconstructed emission image.
+      This axis inversion is now reverted, so
+      <ul>
+        <li><strong>reconstructions of existing SPECT data will have an x-flip with
+          respect to previous versions of STIR</strong> (unless older than 4.0).
+        </li>
+        <li>if your scripts are using an attenuation map, you need to <strong>remove
+          calls to <code>invert_axis x</code> on the attenuation map</strong> (or code to that effect) to get correct results.
+        </li>
+      </ul>
+      <br>
+      See <a href="https://github.com/UCL/STIR/pull/1631">PR #1631</a>.
+    </li>
     <li>
       Fixed a bug in the scatter estimation code (introduced in release 5.1.0) if input data is 3D and "cylindrical"
       (there was no bug for "blocksoncylindrical" data).

examples/python/recon_demo_gradient.py

@@ -6,7 +6,6 @@
 
 # go to directory with input files
 os.chdir('../recon_demo')
-
 # %%
 # initialise reconstruction object
 # we will do this here via a .par file
@@ -34,8 +33,6 @@
 s = recon.set_up(target)
 
 # %%
-
-# compute gradient of objective function
 # create a copy to store the gradient
 gradient=target.get_empty_copy();
 # compute gradient
@@ -51,6 +48,7 @@
 poissonobj.compute_sub_gradient_without_penalty_plus_sensitivity(gradient,target,subset_num)
 # extract to python for plotting
 npimage = stirextra.to_numpy(gradient)
+
 plt.plot(npimage[10, 30, :])
 plt.show()
 

recon_test_pack/SPECT/SPECTUB/README.txt

@@ -5,7 +5,6 @@ done by the STIR-SPECT team when developing the code. These are in the org direc
 Sadly, when the files were generated, we made a mistake such that the attenuation map is actually rotated
 over 90 degrees and flipped in z-direction (as compared to the reconstructed image). 
 This was spotted by Katherine Royston. This means that the "original" reconstructions are actually incorrect.
-In addition, the reconstructed images are inverted along x. So the script inverts the reconstructed images along x.
 
 However, as the test was not designed to check for correctness, but consistency over different STIR versions,
 Kris Thielemans decided to keep the files as they are.

recon_test_pack/SPECT/SPECTUB/run_SPECTUB_tests.sh

@@ -118,7 +118,6 @@ for reconpar in FBP2D OSEM_2DPSF OSEM_3DPSF; do
       output_filename=${output_filename}_${num_subiterations}
     fi
     output_image=${output_filename}.hv
-    invert_axis x ${output_image} ${output_image}
 
     # horrible way to replace "out" with "org" (as we don't want to rely on bash)
     org_output_image=org`echo ${output_image}|cut -c 4-`

recon_test_pack/forward_project_SPECTUB.par

@@ -20,7 +20,7 @@ type:=Matrix
 	;attenuation type := Simple
         ;Values in attenuation map in cm-1 (float file)                         
         ;attenuation map := attenuation.hv
-        attenuation map := my_atten_image_SPECT_modified.hv
+        attenuation map := my_atten_image_SPECT.hv
 
         ;Mask properties { Cylinder // Attenuation Map // Explicit Mask // No}
         ;mask type := Attenuation Map

recon_test_pack/simulate_data_for_tests.sh

@@ -90,8 +90,6 @@ if [ "$generate_images" -eq 1 ]; then
   generate_image  generate_uniform_cylinder.par
   echo "===  make attenuation image"
   generate_image  generate_atten_cylinder.par
-echo "===  make attenuation image for SPECT: to make up for the fliprl of SPECTUB"
-generate_image generate_atten_cylinder_SPECT.par
 fi
 
 # Function to comment out the specific line in a given file and write to a new file
@@ -110,8 +108,8 @@ uniform_output_file="my_uniform_cylinder_SPECT.hv"
 comment_out_line "$uniform_input_file" "$uniform_output_file"
 
 # Paths to the input and output files for attenuation image
-atten_input_file="my_atten_image_SPECT.hv"
-atten_output_file="my_atten_image_SPECT_modified.hv"   
+atten_input_file="my_atten_image.hv"
+atten_output_file="my_atten_image_SPECT.hv"
 
 # Comment out the specific line in the attenuation image file
 comment_out_line "$atten_input_file" "$atten_output_file"
@@ -178,7 +176,7 @@ fi
 : ${background_value:=10}
 if [ "$SPECT" -eq 1 ]; then
   # create SPECT sinograms
-  ./simulate_data.sh "my_uniform_cylinder_SPECT.hv" "my_atten_image_SPECT_modified.hv" "SPECT_test_Interfile_header.hs" "${background_value}" "${suffix}"
+  ./simulate_data.sh "my_uniform_cylinder_SPECT.hv" "my_atten_image_SPECT.hv" "SPECT_test_Interfile_header.hs" "${background_value}" "${suffix}"
   if [ $? -ne 0 ]; then
     echo "Error running SPECT simulation"
        exit 1

src/recon_buildblock/SPECTUB_Weight3d.cxx

@@ -17,7 +17,6 @@
 #include "stir/recon_buildblock/SPECTUB_Weight3d.h"
 #include "stir/error.h"
 #include <boost/math/constants/constants.hpp>
-#include "stir/spatial_transformation/InvertAxis.h"
 
 // system libraries
 #include <cstdlib>
@@ -299,9 +298,7 @@ wm_calculation(const int kOS,
 
                       if (wm.do_save_STIR)
                         {
-                          stir::InvertAxis invert;
-                          wm.nx[vox.iv] = (short int)invert.invert_axis_index(
-                              (vox.icol - (int)floor(vol.Ncold2)), vol.Ncold2 * 2, "x");  // centered index for STIR format
+                          wm.nx[vox.iv] = (short int)(vox.icol - (int)floor(vol.Ncold2)); // centered index for STIR format
                           wm.ny[vox.iv] = (short int)(vox.irow - (int)floor(vol.Nrowd2)); // centered index for STIR format
                           wm.nz[vox.iv] = (short int)vox.islc;                            // non-centered index for STIR format
                         }