merged master, resolved conflicts

.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

@@ -2,6 +2,11 @@
 
 ## v3.6.0
 
+* PET:
+  - added extra members to ScatterEstimation to set behaviour of OSEM used during scatter estimation
+  - added test for scatter simulation and estimation
+  - added missing `set`/`get` methods for OSSPS `relaxation_parameter`, `relaxation_gamma` and `upper_bound`.
+
 * CMake/building:
   - default `DISABLE_MATLAB` to `ON` as our Matlab support is out-of-date and could
   generate conflicts with Python shared libraries.
@@ -18,10 +23,7 @@
   - `return None` in the method `Datacontainer.shape()` replaced with more Pythonesque `return (0,)`.
 
 * MR
-  - Further improvements in handling of "irregular" ISMRMRD acquisitions.
-
-* PET
-  - Added missing `set`/`get` methods for OSSPS `relaxation_parameter`, `relaxation_gamma` and `upper_bound`.
+  - Handling of "irregular" ISMRMRD acquisitions hit what appears to be a bug in recent Gadgetron (HEAD detached at 0670db84). A quick fix applied, Gadgetron issue to be raised.
 
 ## v3.5.0
 

README.md

@@ -10,12 +10,12 @@ This software is the main output of [SyneRBI](https://www.ccpsynerbi.ac.uk), the
 Platform for Synergistic Reconstruction for Biomedical Imaging (formerly CCP PETMR).
 
 Please start with our latest [User's Guide](doc/UserGuide.md).
-See [our Wiki page for installation instructions](https://github.com/SyneRBI/SIRF/wiki/Installation-instructions).
 
-## How do obtain SIRF
+## How to obtain SIRF
 There are multiple ways to obtain a binary version of SIRF, please check them out in our documentation page [how to obtain SIRF](https://github.com/SyneRBI/SIRF/wiki/How-to-obtain-SIRF)
+or [our Wiki page for installation instructions](https://github.com/SyneRBI/SIRF/wiki/Installation-instructions).
 
-# Where is everything installed?
+## Where is everything installed?
 
 SIRF largely follows the usual directory structure with some minor tweaks. When using
 the default options when building SIRF, you will get (most of the) following, depending
@@ -41,7 +41,7 @@ what was found/built:
 [gh-action-link]: https://github.com/SyneRBI/SIRF/actions/workflows/build-test.yml
 [travis-badge]: https://travis-ci.org/SyneRBI/SIRF.svg?branch=master
 [travis]: https://travis-ci.org/SyneRBI/SIRF
-[style-badge]: https://api.codacy.com/project/badge/Grade/392861b4085f4f438d12c41029f86b47
-[style-link]: https://www.codacy.com/gh/SyneRBI/SIRF?utm_source=github.com&utm_medium=referral&utm_content=SyneRBI/SIRF&utm_campaign=Badge_Grade
+[style-badge]: https://app.codacy.com/project/badge/Grade/392861b4085f4f438d12c41029f86b47
+[style-link]: https://app.codacy.com/gh/SyneRBI/SIRF/dashboard?utm_source=gh&utm_medium=referral&utm_content=&utm_campaign=Badge_grade
 [zenodo-badge]: https://zenodo.org/badge/DOI/10.5281/zenodo.2707911.svg
 [zenodo-link]: https://doi.org/10.5281/zenodo.2707911

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__)