Create SPECT subsets from ProjData (#1596)

* fix SPECT subsets functionality 

* separated subset tests for pet & spect

---------
Co-authored-by: Kris Thielemans <[email protected]>

Author: samdporter

documentation/release_6.3.htm

@@ -75,6 +75,11 @@ <h4>General</h4>
       larger than the default from the vendor).<br>
       <a href=https://github.com/UCL/STIR/pull/1315>PR #1315</a>
     </li>
+    <li>
+      Adapt the SPECTUB projector to handle data with only a subset of the views (constructed using <code>ProjData::get_subsets()</code>.
+      This is useful for SIRF/CIL applications using stochastic optimisation.<br>
+      <a href=https://github.com/UCL/STIR/pull/1596>PR #1596</a>
+    </li>
   </ul>
   <h4>Python</h4>
   <ul>

examples/samples/SPECT_Interfile_header.hs

@@ -5,7 +5,10 @@
 !imaging modality := nucmed
 
 ; name of file with binary data
-name of data file := somefile.s
+; Note: we are "self-referring" to the header here such that this file can be read by STIR as a "template".
+; (This is used in one of our tests.)
+; Obviously, this would normally not be the case.
+name of data file := SPECT_Interfile_header.hs
 
 !version of keys := 3.3
 !GENERAL DATA :=

src/include/stir/recon_buildblock/SPECTUB_Tools.h

@@ -14,6 +14,7 @@
 #define _WM_SPECTUB_H
 
 #include <string>
+#include <vector>
 
 namespace SPECTUB
 {
@@ -84,10 +85,9 @@ typedef struct
   int Nsli;   // number of slices
   float thcm; // slice thickness in cm
 
-  int Nang;   // number of projection angles
-  float ang0; // initial projection angle. degrees from upper detection plane (parallel to table). Negative for CW rotacions (see
-              // manual)
-  float incr; // angle increment between two consecutive projection angles. Degrees. Negative for CW, Positive for CCW
+  int Nang;                  // number of projection angles
+  std::vector<float> angles; // projection angles. degrees from upper detection plane (parallel to table). Negative for CW
+                             // rotations (see manual)
 
   int NOS;    // number of subsets in which to split the matrix
   int NangOS; // Number of angles in each subset = Nang/NOS
@@ -339,4 +339,4 @@ void free_wm_da(wm_da_type* f); // to free weight_mat_da
 
 } // namespace SPECTUB
 
-#endif //_WM_SPECT_H
+#endif //_WM_SPECT_H

src/recon_buildblock/ProjMatrixByBinSPECTUB.cxx

@@ -24,12 +24,14 @@
 #include "stir/recon_buildblock/ProjMatrixByBinSPECTUB.h"
 #include "stir/recon_buildblock/TrivialDataSymmetriesForBins.h"
 #include "stir/ProjDataInfoCylindricalArcCorr.h"
+#include "stir/ProjDataInfoSubsetByView.h"
 #include "stir/IO/read_from_file.h"
 #include "stir/ProjDataInfo.h"
 #include "stir/VoxelsOnCartesianGrid.h"
 #include "stir/Succeeded.h"
 #include "stir/is_null_ptr.h"
 #include "stir/Coordinate3D.h"
+#include "stir/stream.h"
 #include "stir/info.h"
 #include "stir/warning.h"
 #include "stir/error.h"
@@ -254,9 +256,6 @@ ProjMatrixByBinSPECTUB::set_up(const shared_ptr<const ProjDataInfo>& proj_data_i
   this->voxel_size = image_info_ptr->get_voxel_size();
   this->origin = image_info_ptr->get_origin();
 
-  const ProjDataInfoCylindricalArcCorr* proj_Data_Info_Cylindrical
-      = dynamic_cast<const ProjDataInfoCylindricalArcCorr*>(this->proj_data_info_ptr.get());
-
   CPUTimer timer;
   timer.start();
 
@@ -299,15 +298,23 @@ ProjMatrixByBinSPECTUB::set_up(const shared_ptr<const ProjDataInfo>& proj_data_i
   vox.thcm = vol.thcm;
 
   //... projecction parameters ..........................................
-  prj.ang0 = this->proj_data_info_ptr->get_scanner_ptr()->get_intrinsic_azimuthal_tilt() * float(180 / _PI);
-  prj.incr = proj_Data_Info_Cylindrical->get_azimuthal_angle_sampling() * float(180 / _PI);
-  prj.thcm = proj_Data_Info_Cylindrical->get_axial_sampling(0) / 10;
+  prj.angles = std::vector<float>(prj.Nang);
+  {
+    Bin bin;
+    for (int i = 0; i < prj.Nang; i++)
+      {
+        bin.view_num() = i;
+        prj.angles[i] = static_cast<float>(this->proj_data_info_ptr->get_phi(bin) * 180 / _PI);
+      }
+    // all bins will have same axial sampling
+    prj.thcm = this->proj_data_info_ptr->get_sampling_in_m(bin) / 10;
+  }
 
   //.......geometrical and other derived parameters of projection structure...........
-  prj.Nsli = proj_Data_Info_Cylindrical->get_num_axial_poss(0); // number of slices
-  prj.lngcm = prj.Nbin * prj.szcm;                              // length in cm of the detection line
-  prj.Nbp = prj.Nbin * prj.Nsli;                                // number of bins for each projection angle (2D-projection)
-  prj.Nbt = prj.Nbp * prj.Nang;                                 // total number of bins considering all the projection angles
+  prj.Nsli = this->proj_data_info_ptr->get_num_axial_poss(0); // number of slices
+  prj.lngcm = prj.Nbin * prj.szcm;                            // length in cm of the detection line
+  prj.Nbp = prj.Nbin * prj.Nsli;                              // number of bins for each projection angle (2D-projection)
+  prj.Nbt = prj.Nbp * prj.Nang;                               // total number of bins considering all the projection angles
   prj.Nbind2 = (float)prj.Nbin / (float)2.; // half of the number of bins in a detection line (integer or semi-integer)
   prj.lngcmd2 = prj.lngcm / (float)2.;      // half of the length of detection line (cm)
   prj.Nslid2 = (float)prj.Nsli / (float)2.; // half of the number of slices (integer or semi-integer)
@@ -329,7 +336,35 @@ ProjMatrixByBinSPECTUB::set_up(const shared_ptr<const ProjDataInfo>& proj_data_i
                  wmh.prj.Nsli,
                  vol.Nsli));
   //....rotation radius .................................................
-  const VectorWithOffset<float> radius_all_views = proj_Data_Info_Cylindrical->get_ring_radii_for_all_views();
+  VectorWithOffset<float> radius_all_views(prj.Nang);
+  {
+    if (auto proj_Data_Info_Cylindrical = dynamic_cast<const ProjDataInfoCylindricalArcCorr*>(this->proj_data_info_ptr.get()))
+      {
+        radius_all_views = proj_Data_Info_Cylindrical->get_ring_radii_for_all_views();
+      }
+    else if (auto proj_data_info_subset_ptr = dynamic_cast<const ProjDataInfoSubsetByView*>(this->proj_data_info_ptr.get()))
+      {
+        if (auto proj_Data_Info_Cylindrical = dynamic_cast<const ProjDataInfoCylindricalArcCorr*>(
+                proj_data_info_subset_ptr->get_original_proj_data_info_sptr().get()))
+          {
+            for (int i = 0; i < prj.Nang; ++i)
+              {
+                Bin bin;
+                bin.view_num() = i;
+                const int org_view = proj_data_info_subset_ptr->get_original_bin(bin).view_num();
+                radius_all_views[i] = proj_Data_Info_Cylindrical->get_ring_radius(org_view);
+              }
+          }
+        else
+          {
+            error("ProjMatrixByBinSPECTUB: can only handle ProjDataInfoCylindricalArcCorr");
+          }
+      }
+    else
+      {
+        error("ProjMatrixByBinSPECTUB: can only handle ProjDataInfoCylindricalArcCorr");
+      }
+  }
 
   {
     const auto max_radius = *std::max_element(radius_all_views.begin(), radius_all_views.end());
@@ -493,8 +528,7 @@ ProjMatrixByBinSPECTUB::set_up(const shared_ptr<const ProjDataInfo>& proj_data_i
   info_stream << "Number of slices: " << wmh.vol.Nsli << "\tslice_thickness: " << wmh.vol.thcm << std::endl;
   info_stream << "Number of bins: " << wmh.prj.Nbin << "\tbin size: " << wmh.prj.szcm << "\taxial size: " << wmh.prj.thcm
               << std::endl;
-  info_stream << "Number of angles: " << wmh.prj.Nang << "\tAngle increment: " << wmh.prj.incr
-              << "\tFirst angle: " << wmh.prj.ang0 << std::endl;
+  info_stream << "Number of angles: " << wmh.prj.Nang << "\tangles (deg): " << wmh.prj.angles << std::endl;
   info_stream << "Number of subsets: " << wmh.prj.NOS << std::endl;
   if (wmh.do_att)
     {

src/recon_buildblock/SPECTUB_Tools.cxx

@@ -13,6 +13,7 @@
 // system libraries
 #include <stdlib.h>
 #include <math.h>
+
 using namespace std;
 
 //... user defined libraries .......................................
@@ -301,9 +302,9 @@ fill_ang(angle_type* ang, const SPECTUB::wmh_type& wmh, const float* Rrad)
 
       //... ratios calculation .......................................................
 
-      const float deg = wmh.prj.ang0 + (float)i * wmh.prj.incr; // angle in degrees
-      ang[i].cos = cos(deg * dg2rd);                            // cosinus of the angle
-      ang[i].sin = sin(deg * dg2rd);                            // sinus of the angle
+      const float deg = wmh.prj.angles[i];
+      ang[i].cos = cos(deg * dg2rd); // cosinus of the angle
+      ang[i].sin = sin(deg * dg2rd); // sinus of the angle
 
       //... first octave (0->45degrees) equivalent angle and its trigonometric ratios .......
 

src/test/CMakeLists.txt

@@ -26,7 +26,7 @@ set(${dir_SIMPLE_TEST_EXE_SOURCES}
 	test_DynamicDiscretisedDensity.cxx
 	test_ScatterSimulation.cxx
         test_ML_norm.cxx
-	test_proj_data_info_subsets.cxx
+	test_proj_data_info_subsets_pet.cxx
 )
 
 set(${dir_SIMPLE_TEST_EXE_SOURCES_NO_REGISTRIES}
@@ -43,6 +43,7 @@ Set(${dir_INVOLVED_TEST_EXE_SOURCES}
 	test_linear_regression.cxx
 	test_stir_math.cxx
 	test_time_of_flight.cxx
+	test_proj_data_info_subsets_spect.cxx
         # the next 2 are interactive, so we don't add a test for it, but only compile them
 	test_display.cxx
 	test_interpolate.cxx
@@ -98,6 +99,11 @@ endforeach()
 
 if (BUILD_TESTING)
 
+# Add test for test_proj_data_info_subsets_spect with arguments
+ADD_TEST(NAME test_proj_data_info_subsets_spect
+			COMMAND ${CMAKE_CURRENT_BINARY_DIR}/test_proj_data_info_subsets_spect ${CMAKE_SOURCE_DIR}/examples/samples/SPECT_Interfile_header.hs)
+add_STIR_CONFIG_DIR(test_proj_data_info_subsets_spect)
+
 ADD_TEST(test_linear_regression
    ${CMAKE_CURRENT_BINARY_DIR}/test_linear_regression ${CMAKE_CURRENT_SOURCE_DIR}/input/test_linear_regression.in
 )
@@ -166,6 +172,7 @@ ADD_TEST(test_IO_DynamicDiscretisedDensity_Interfile_short
 ADD_TEST(test_IO_DynamicDiscretisedDensity_Multi
 	${CMAKE_CURRENT_BINARY_DIR}/test_IO_DynamicDiscretisedDensity ${CMAKE_CURRENT_SOURCE_DIR}/input/test_MultiOutputFileFormat.in)
 
+
 endif(BUILD_TESTING)
 
 endif() # MINI_STIR