Merge pull request #1223 from rmapree/scatter_subsets

Add set get functions for OSEM reconstructor in  STIR.ScatterEstimator
Add tests for STIR.ScatterSimulator and STIR.ScatterSimulator

.bandit

@@ -0,0 +1,2 @@
+skips:
+  - "B101" # warning about assert

.codacy.yml

@@ -1,2 +1,4 @@
 exclude_paths:
-  - '**.md'
+  - '**.md'
+assert_used:
+  - skips: ['*test*.py', 'test*.py']

.mailmap

@@ -23,5 +23,5 @@ Gemma Fardell <[email protected]> <[email protected]>
 Sam D. Porter <[email protected]>
 Matthew Strugari <[email protected]>
 Matthew Strugari <[email protected]> <[email protected]>
-
+Toni Reeves <[email protected]>
 

CHANGES.md

@@ -1,5 +1,9 @@
 # ChangeLog
 
+* PET:
+  - added extra members to ScatterEstimation to set behaviour of OSEM used during scatter estimation
+  - added test for scatter simulation and estimation
+
 ## v3.5.1
 
 * CMake/building:

examples/Python/PET/scatter_estimation.py

@@ -42,7 +42,7 @@
 ##   See the License for the specific language governing permissions and
 ##   limitations under the License.
 
-__version__ = '1.0.0'
+__version__ = '1.0.1'
 from docopt import docopt
 
 args = docopt(__doc__, version=__version__)
@@ -94,9 +94,12 @@ def main():
         se.set_asm(PET.AcquisitionSensitivityModel(norm_file))
     if not(acf_file is None):
         se.set_attenuation_correction_factors(PET.AcquisitionData(acf_file))
-    # could set number of iterations if you want to
+    # Could set number of iterations if you want to (we recommend at least 3)
     se.set_num_iterations(1)
     print("number of scatter iterations that will be used: %d" % se.get_num_iterations())
+    # Could set number of subsets used by the OSEM reconstruction inside the scatter estimation loop.
+    # Here we will set it to 7 (which is in fact the default), which is appropriate for the mMR
+    set.set_OSEM_num_subsets(7)
     se.set_output_prefix(output_prefix)
     se.set_up()
     se.process()

examples/parameter_files/scatter_template.hs

@@ -5,6 +5,7 @@
 ; This is sufficient to act as a template
 name of data file := scatter_template.hs
 originating system := unknown
+!imaging modality := PT
 !GENERAL DATA :=
 !GENERAL IMAGE DATA :=
 !type of data := PET
@@ -25,19 +26,30 @@ matrix axis label [1] := tangential coordinate
 !matrix size [1] := 35
 minimum ring difference per segment := { 0}
 maximum ring difference per segment := { 0}
-Scanner parameters:= 
-Scanner type := unknown
-Energy resolution := 0.145
-Reference energy (in keV) := 511
 number of energy windows := 1
 energy window lower level[1] := 430
 energy window upper level[1] := 610
+number of time frames := 1
+image duration (sec)[1] := 60
+image relative start time (sec)[1] := 0
 
+Scanner parameters:= 
+Scanner type := unknown
+Energy resolution := 0.145
+Reference energy (in keV) := 511
 Number of rings                          := 8
 Number of detectors per ring             := 64
 Inner ring diameter (cm)                 := 65.6
 Average depth of interaction (cm)        := 0.7
 Distance between rings (cm)              := 3.25
 View offset (degrees)                    := 0
+; some block info to avoid warnings
+Number of blocks per bucket in transaxial direction         := 1
+Number of blocks per bucket in axial direction              := 1
+Number of crystals per block in axial direction             := 1
+Number of crystals per block in transaxial direction        := 1
+Number of detector layers                                   := 1
+Number of crystals per singles unit in axial direction      := 1
+Number of crystals per singles unit in transaxial direction := 1
 end scanner parameters:=
 !END OF INTERFILE :=

src/xSTIR/cSTIR/cstir_p.cpp

@@ -849,6 +849,21 @@ sirf::cSTIR_setScatterEstimatorParameter
         int value = dataFromHandle<int>(hv);
         obj.set_num_iterations(value);
     }
+
+    else if (sirf::iequals(name, "set_OSEM_num_subiterations"))
+    {
+        int value = dataFromHandle<int>(hv);
+        obj.set_OSEM_num_subiterations(value);
+    }
+
+    else if (sirf::iequals(name, "set_OSEM_num_subsets"))
+    {
+        int value = dataFromHandle<int>(hv);
+        obj.set_OSEM_num_subsets(value);
+    }
+
+
+
     else if (sirf::iequals(name, "set_output_prefix"))
     {
         obj.set_output_prefix(charDataFromHandle(hv));
@@ -867,6 +882,10 @@ sirf::cSTIR_ScatterEstimatorParameter(DataHandle* hp, const char* name)
 		return newObjectHandle(processor.get_output());
 	if (sirf::iequals(name, "num_iterations"))
           return dataHandle<int>(processor.get_num_iterations());
+	if (sirf::iequals(name, "OSEM_num_subiterations"))
+          return dataHandle<int>(processor.get_OSEM_num_subiterations());
+	if (sirf::iequals(name, "OSEM_num_subsets"))
+          return dataHandle<int>(processor.get_OSEM_num_subsets());
 	return parameterNotFound(name, __FILE__, __LINE__);
 }
 

src/xSTIR/cSTIR/include/sirf/STIR/stir_x.h

@@ -620,22 +620,26 @@ The actual algorithm is described in
       \brief Class for estimating the scatter contribution in PET projection data
 
       This class implements the SSS iterative algorithm from STIR. It
-      is an iterative loop of reconstruction, single scatter estimation,
+      is an iterative loop of reconstruction via OSEM, single scatter estimation,
       upsampling, tail-fitting.
 
       Output is an acquisition_data object with the scatter contribution.
       This can be added to the randoms to use in PETAcquisitionModel.set_background_term().
     */
     class PETScatterEstimator : private stir::ScatterEstimation
     {
+      using recon_type = stir::OSMAPOSLReconstruction<DiscretisedDensity<3,float>>;
     public:
-        //!
+        //! constructor.
+        /*! sets reconstruction to OSEM with 8 subiterations, 7 subsets, and a postfilter of (15mm)^3.
+
+         \warning The default settings might not work for your data, in particular the number of subsets.
+         Use set_OSEM_num_subsets() to change it.
+        */
         PETScatterEstimator() : stir::ScatterEstimation()
         {
-          stir::shared_ptr<stir::PoissonLogLikelihoodWithLinearModelForMeanAndProjData<DiscretisedDensity<3,float> > >
-            obj_sptr(new PoissonLogLikelihoodWithLinearModelForMeanAndProjData<DiscretisedDensity<3,float> >);
-          stir::shared_ptr<stir::OSMAPOSLReconstruction<DiscretisedDensity<3,float> > >
-            recon_sptr(new stir::OSMAPOSLReconstruction<DiscretisedDensity<3,float> >);
+          auto obj_sptr = std::make_shared<PoissonLogLikelihoodWithLinearModelForMeanAndProjData<DiscretisedDensity<3,float>>>();
+          auto recon_sptr = std::make_shared<recon_type>();
           recon_sptr->set_num_subiterations(8);
           recon_sptr->set_num_subsets(7); // this works for the mMR. TODO
           recon_sptr->set_disable_output(true);
@@ -705,6 +709,27 @@ The actual algorithm is described in
         {
           return stir::ScatterEstimation::get_num_iterations();
         }
+
+        void set_OSEM_num_subiterations(int arg)
+        {
+          this->get_reconstruction_method().set_num_subiterations(arg);
+        }
+
+        int get_OSEM_num_subiterations() const
+        {
+          return this->get_reconstruction_method().get_num_subiterations();
+        }
+
+        void set_OSEM_num_subsets(int arg)
+        {
+          this->get_reconstruction_method().set_num_subsets(arg);
+          this->_already_set_up_sirf = false;
+        }
+
+        int get_OSEM_num_subsets() const
+        {
+          return this->get_reconstruction_method().get_num_subsets();
+        }
 
         std::shared_ptr<STIRAcquisitionData> get_scatter_estimate(int est_num = -1) const
         {
@@ -751,15 +776,27 @@ The actual algorithm is described in
           stir::ScatterEstimation::set_initial_activity_image_sptr(image_sptr);
           if (stir::ScatterEstimation::set_up() == Succeeded::no)
             THROW("scatter estimation set_up failed");
+          this->_already_set_up_sirf = true;
           return Succeeded::yes;
         }
         void process()
         {
+          if (!this->_already_set_up_sirf)
+            THROW("scatter estimation needs to be set-up first");
           if (!stir::ScatterEstimation::already_setup())
             THROW("scatter estimation needs to be set-up first");
           if (stir::ScatterEstimation::process_data() == Succeeded::no)
             THROW("scatter estimation failed");
         }
+    private:
+        //! work-around for private variable in stir::ScatterEstimation
+        bool _already_set_up_sirf;
+
+        //! work-around to missing method in stir::ScatterEstimation
+        recon_type& get_reconstruction_method() const
+        {
+          return dynamic_cast<recon_type&>(*this->reconstruction_template_sptr);
+        }
     };
 
 	/*!

src/xSTIR/pSTIR/STIR.py

@@ -3354,6 +3354,14 @@ def get_output(self):
     def get_num_iterations(self):
         """Get number of iterations of the SSS algorithm to use."""
         return parms.int_par(self.handle, 'PETScatterEstimator', 'num_iterations')
+
+    def get_OSEM_num_subiterations(self):
+        """Get number of subiterations used by OSEM in the SSS algorithm."""
+        return parms.int_par(self.handle, 'PETScatterEstimator', 'OSEM_num_subiterations')
+
+    def get_OSEM_num_subsets(self):
+        """Get number of subsets used by OSEM in the SSS algorithm."""
+        return parms.int_par(self.handle, 'PETScatterEstimator', 'OSEM_num_subsets')
 
     def set_attenuation_image(self, image):
         assert_validity(image, ImageData)
@@ -3376,6 +3384,14 @@ def set_asm(self, asm):
         assert_validity(asm, AcquisitionSensitivityModel)
         parms.set_parameter(self.handle, self.name, 'setASM', asm.handle)
 
+    def set_OSEM_num_subiterations(self, v):
+        """Set number of subiterations used by OSEM in the SSS algorithm."""
+        parms.set_int_par(self.handle, 'PETScatterEstimator', 'set_OSEM_num_subiterations', v)
+
+    def set_OSEM_num_subsets(self, v):
+        """Set number of subsets used by OSEM in the SSS algorithm."""
+        parms.set_int_par(self.handle, 'PETScatterEstimator', 'set_OSEM_num_subsets', v)
+
     def set_num_iterations(self, v):
         """Set number of iterations of the SSS algorithm to use."""
         parms.set_int_par(self.handle, 'PETScatterEstimator', 'set_num_iterations', v)

src/xSTIR/pSTIR/tests/tests_scatter.py

@@ -0,0 +1,108 @@
+# -*- coding: utf-8 -*-
+"""sirf.STIR tests for ScatterEstimation
+v{version}
+
+Usage:
+  tests_scatter [--help | options]
+
+Options:
+  -r, --record   record the measurements rather than check them [actually not used for current test]
+  -v, --verbose  report each test status
+
+{author}
+
+{licence}
+"""
+import sirf.STIR as pet
+import os
+from sirf.Utilities import runner, __license__
+__version__ = "0.1.0"
+__author__ = "Kris Thielemans"
+
+
+def test_main(rec=False, verb=False, throw=True):
+    _ = pet.MessageRedirector()
+
+    #data_path = pet.examples_data_path('PET')
+    #raw_data_file = pet.existing_filepath(data_path,'Utahscat600k_ca_seg4.hs')
+    data_path = os.path.join(os.path.dirname(__file__), '..', '..', '..', '..', 'examples', 'parameter_files')
+    print(data_path)
+    acq_template_filename = pet.existing_filepath(data_path, 'scatter_template.hs')
+    #image_file = pet.existing_filepath(data_path, 'test_image_PM_QP_6.hv')
+
+    #acq_data = pet.AcquisitionData(raw_data_file)
+    acq_template = pet.AcquisitionData(acq_template_filename)
+    act_image = pet.ImageData(acq_template)
+    atten_image = act_image.get_uniform_copy(0)
+    image_size = atten_image.dimensions()
+    voxel_size = atten_image.voxel_sizes()
+
+    # create a cylindrical water shape
+    water_cyl = pet.EllipticCylinder()
+    water_cyl.set_length(image_size[0]*voxel_size[0])
+    water_cyl.set_radii((100, 100))
+    water_cyl.set_origin((image_size[0]*voxel_size[0]*0.5, 0, 0))
+
+    # add the shape to the image
+    atten_image.add_shape(water_cyl, scale = 9.687E-02) # use mu of water
+    act_image.add_shape(water_cyl, scale = 20)
+    # Need to create ACFs
+    # create acquisition model
+    am = pet.AcquisitionModelUsingRayTracingMatrix()
+    am.set_up(acq_template, atten_image)
+    # create acquisition sensitivity model from attenuation image
+    asm = pet.AcquisitionSensitivityModel(atten_image, am)
+    asm.set_up(acq_template)
+    am.set_acquisition_sensitivity(asm)
+    # apply attenuation to the uniform acquisition data to obtain ACFs
+    ACFs = acq_template.get_uniform_copy(1)
+    asm.normalise(ACFs)
+
+    sss = pet.SingleScatterSimulator()
+    sss.set_attenuation_image(atten_image)
+    sss.set_up(acq_template, act_image)
+    scatter_data = sss.forward(act_image)
+
+    acq_model = pet.AcquisitionModelUsingRayTracingMatrix()
+    acq_model.set_acquisition_sensitivity(asm)
+    acq_model.set_up(acq_template, act_image)
+    unscattered_data = acq_model.forward(act_image)
+
+    acq_data = unscattered_data + scatter_data
+
+    # I get around 21% scatter fraction for this data
+    scatter_fraction = scatter_data.norm()/acq_data.norm()
+    if scatter_fraction < .18 or scatter_fraction > .25:
+        scatter_data.write("out_scatter_data.hs")
+        unscattered_data.write("out_unscattered.hs")
+        assert False, f"Scatter fraction ({scatter_fraction}) is out of range (should be around .2 for this data)"
+
+    scat_est = pet.ScatterEstimator()
+    scat_est.set_input(acq_data)
+    scat_est.set_attenuation_image(atten_image)
+    scat_est.set_attenuation_correction_factors(ACFs)
+    scat_est.set_asm(pet.AcquisitionSensitivityModel(acq_data.get_uniform_copy(1)))
+    scat_est.set_randoms(acq_data.get_uniform_copy(0))
+    scat_est.set_OSEM_num_subsets(4)
+    assert scat_est.get_OSEM_num_subsets() == 4
+    scat_est.set_OSEM_num_subiterations(3)
+    assert scat_est.get_OSEM_num_subiterations() == 3
+    scat_est.set_num_iterations(5)
+    assert scat_est.get_num_iterations() == 5
+    scat_est.set_up()
+    scat_est.process()
+    scatter_estimate = scat_est.get_output()
+
+    rel_err = (scatter_data - scatter_estimate).norm() / scatter_estimate.norm()
+    # I get around 10% error, due to randomness etc, so set our threshold a bit larger
+    if  rel_err > .2:
+        scatter_estimate.write("out_scatter_estimate.hs")
+        scatter_data.write("out_scatter_data.hs")
+        unscattered_data.write("out_unscattered.hs")
+        assert False, f"Difference between simulated and estimated scatter is too large (rel err {rel_err}). Data written to file as out*.hs"
+
+    return 0, 1
+
+
+if __name__ == "__main__":
+    runner(test_main, __doc__, __version__, __author__)