diff --git a/.bandit b/.bandit new file mode 100644 index 000000000..78beca269 --- /dev/null +++ b/.bandit @@ -0,0 +1,2 @@ +skips: + - "B101" # warning about assert \ No newline at end of file diff --git a/.codacy.yml b/.codacy.yml index bd1402f16..cb3cb3107 100644 --- a/.codacy.yml +++ b/.codacy.yml @@ -1,2 +1,4 @@ exclude_paths: - - '**.md' \ No newline at end of file + - '**.md' +assert_used: + - skips: ['*test*.py', 'test*.py'] diff --git a/.mailmap b/.mailmap index 286f5402f..6aef47bab 100644 --- a/.mailmap +++ b/.mailmap @@ -23,5 +23,5 @@ Gemma Fardell Sam D. Porter <92305641+samdporter@users.noreply.github.com> Matthew Strugari Matthew Strugari <56315593+mastergari@users.noreply.github.com> - +Toni Reeves diff --git a/CHANGES.md b/CHANGES.md index 091c5526e..624bdae52 100644 --- a/CHANGES.md +++ b/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: diff --git a/examples/Python/PET/scatter_estimation.py b/examples/Python/PET/scatter_estimation.py index 140045531..98ec1b068 100755 --- a/examples/Python/PET/scatter_estimation.py +++ b/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() diff --git a/examples/parameter_files/scatter_template.hs b/examples/parameter_files/scatter_template.hs index 3e726e181..3698c30df 100644 --- a/examples/parameter_files/scatter_template.hs +++ b/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 := diff --git a/src/xSTIR/cSTIR/cstir_p.cpp b/src/xSTIR/cSTIR/cstir_p.cpp index 2e6ceec68..1bc0627ea 100644 --- a/src/xSTIR/cSTIR/cstir_p.cpp +++ b/src/xSTIR/cSTIR/cstir_p.cpp @@ -849,6 +849,21 @@ sirf::cSTIR_setScatterEstimatorParameter int value = dataFromHandle(hv); obj.set_num_iterations(value); } + + else if (sirf::iequals(name, "set_OSEM_num_subiterations")) + { + int value = dataFromHandle(hv); + obj.set_OSEM_num_subiterations(value); + } + + else if (sirf::iequals(name, "set_OSEM_num_subsets")) + { + int value = dataFromHandle(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(processor.get_num_iterations()); + if (sirf::iequals(name, "OSEM_num_subiterations")) + return dataHandle(processor.get_OSEM_num_subiterations()); + if (sirf::iequals(name, "OSEM_num_subsets")) + return dataHandle(processor.get_OSEM_num_subsets()); return parameterNotFound(name, __FILE__, __LINE__); } diff --git a/src/xSTIR/cSTIR/include/sirf/STIR/stir_x.h b/src/xSTIR/cSTIR/include/sirf/STIR/stir_x.h index 1ae732545..62a654384 100644 --- a/src/xSTIR/cSTIR/include/sirf/STIR/stir_x.h +++ b/src/xSTIR/cSTIR/include/sirf/STIR/stir_x.h @@ -620,7 +620,7 @@ 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. @@ -628,14 +628,18 @@ The actual algorithm is described in */ class PETScatterEstimator : private stir::ScatterEstimation { + using recon_type = stir::OSMAPOSLReconstruction>; 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 > > - obj_sptr(new PoissonLogLikelihoodWithLinearModelForMeanAndProjData >); - stir::shared_ptr > > - recon_sptr(new stir::OSMAPOSLReconstruction >); + auto obj_sptr = std::make_shared>>(); + auto recon_sptr = std::make_shared(); 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 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(*this->reconstruction_template_sptr); + } }; /*! diff --git a/src/xSTIR/pSTIR/STIR.py b/src/xSTIR/pSTIR/STIR.py index 0e53d435c..75bd0bf07 100644 --- a/src/xSTIR/pSTIR/STIR.py +++ b/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) diff --git a/src/xSTIR/pSTIR/tests/tests_scatter.py b/src/xSTIR/pSTIR/tests/tests_scatter.py new file mode 100644 index 000000000..bcdcd4f05 --- /dev/null +++ b/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__)