Reinvert spect axis

documentation/release_6.3.htm

@@ -131,6 +131,13 @@ <h3>Changed functionality</h3>
 
   <h3>Bug fixes</h3>
   <ul>
+    <li>
+      Fixed a bug in the SPECTUB matrix (introduced in 4.0) that inverted the SPECT reconstruction along the axis.
+      This had resulted with an inconsitent attenuation correction orientation compared to the reonctructed image. 
+      The axis inversion is now reverted, so if your scripts were using an inverted attenuation map they need to be
+      updated with this release. 
+      See <a href="https://github.com/UCL/STIR/pull/1631"">PR #1631</a>" . </li>
+    </ul>
     <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/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

@@ -299,9 +299,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
                         }