diff --git a/documentation/release_6.3.htm b/documentation/release_6.3.htm
index e7128ec3af..78442c8b8f 100644
--- a/documentation/release_6.3.htm
+++ b/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).
diff --git a/examples/python/recon_demo_gradient.py b/examples/python/recon_demo_gradient.py
index ed9c4beb15..35a3985783 100644
--- a/examples/python/recon_demo_gradient.py
+++ b/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()
 
diff --git a/recon_test_pack/SPECT/SPECTUB/README.txt b/recon_test_pack/SPECT/SPECTUB/README.txt
index dc9f2ec9a3..768f988edc 100644
--- a/recon_test_pack/SPECT/SPECTUB/README.txt
+++ b/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.
diff --git a/recon_test_pack/SPECT/SPECTUB/run_SPECTUB_tests.sh b/recon_test_pack/SPECT/SPECTUB/run_SPECTUB_tests.sh
index d898f9aab1..d81d757581 100755
--- a/recon_test_pack/SPECT/SPECTUB/run_SPECTUB_tests.sh
+++ b/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-`
diff --git a/recon_test_pack/forward_project_SPECTUB.par b/recon_test_pack/forward_project_SPECTUB.par
index ddcd2d6d6c..485e43bce6 100644
--- a/recon_test_pack/forward_project_SPECTUB.par
+++ b/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
diff --git a/recon_test_pack/generate_atten_cylinder_SPECT.par b/recon_test_pack/generate_atten_cylinder_SPECT.par
deleted file mode 100644
index 2f88210255..0000000000
--- a/recon_test_pack/generate_atten_cylinder_SPECT.par
+++ /dev/null
@@ -1,24 +0,0 @@
-generate_image Parameters :=
-output filename:=my_atten_image_SPECT
-X output image size (in pixels):=111
-Y output image size (in pixels):=111
-Z output image size (in pixels):=47
-X voxel size (in mm):= 3
-Y voxel size (in mm):= 3
-Z voxel size (in mm) :=3.27
-
-  Z number of samples to take per voxel := 1
-  Y number of samples to take per voxel := 1
-  X number of samples to take per voxel := 1
-
-shape type:= ellipsoidal cylinder
-Ellipsoidal Cylinder Parameters:=
-   radius-x (in mm):=100
-   radius-y (in mm):=100
-   length-z (in mm):=110
-   origin (in mm):={70,10,-20}
-;{70,10,-20}
-   END:=
-value := 0.096
-
-END:=
diff --git a/recon_test_pack/simulate_data_for_tests.sh b/recon_test_pack/simulate_data_for_tests.sh
index 18e04f42eb..a234f6ac29 100755
--- a/recon_test_pack/simulate_data_for_tests.sh
+++ b/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
diff --git a/src/recon_buildblock/SPECTUB_Weight3d.cxx b/src/recon_buildblock/SPECTUB_Weight3d.cxx
index 422bcdad9e..8bf71d1465 100644
--- a/src/recon_buildblock/SPECTUB_Weight3d.cxx
+++ b/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
                         }