From 6aee189d830e16787330162452ff877864cec300 Mon Sep 17 00:00:00 2001 From: richard Date: Thu, 18 Jun 2020 10:33:16 +0100 Subject: [PATCH 01/14] poisson gated w motion: c++ changes --- src/CMakeLists.txt | 17 +++ src/Synergistic/tests/CMakeLists.txt | 7 ++ .../tests/test_PoissonGatedMotion.cpp | 113 ++++++++++++++++++ src/xSTIR/cSTIR/cstir.cpp | 35 ++++++ src/xSTIR/cSTIR/include/sirf/STIR/cstir.h | 8 ++ .../cSTIR/include/sirf/STIR/stir_types.h | 11 ++ src/xSTIR/cSTIR/include/sirf/STIR/stir_x.h | 47 +++++++- src/xSTIR/cSTIR/stir_x.cpp | 111 ++++++++++++++++- 8 files changed, 342 insertions(+), 7 deletions(-) create mode 100644 src/Synergistic/tests/test_PoissonGatedMotion.cpp diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index bb162df38..ca5a91097 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -79,6 +79,7 @@ if (DISABLE_STIR) message(STATUS "STIR support disabled.") else() find_package(STIR 4.0.0 REQUIRED) + ## NiftyPET if (STIR_WITH_NiftyPET_PROJECTOR) set(NiftyPET_BOOL_STR "1") message(STATUS "STIR was built with NiftyPET, GPU projectors will be enabled.") @@ -91,6 +92,22 @@ else() MESSAGE(STATUS "STIR not built with NiftyPET.") set(NiftyPET_BOOL_STR "0") endif() + + ## Gated recon w/ motion + if (STIR_WITH_EXPERIMENTAL AND STIR_BUILT_WITH_ITK) + set(STIR_GATED_MOTION "1") + message(STATUS "STIR built with EXPERIMENTAL and ITK, gated motion recon enabled.") + if(${CMAKE_VERSION} VERSION_LESS "3.12.0") + add_definitions(-DSTIR_GATED_MOTION) + else() + add_compile_definitions(STIR_GATED_MOTION) + endif() + else() + MESSAGE(STATUS "STIR not built with EXPERIMENTAL or ITK, gated motion recon disabled.") + set(STIR_GATED_MOTION "0") + remove_definitions(-DSTIR_GATED_MOTION) + endif() + ADD_SUBDIRECTORY(xSTIR) endif() diff --git a/src/Synergistic/tests/CMakeLists.txt b/src/Synergistic/tests/CMakeLists.txt index a9fb52aa6..749e59e5d 100644 --- a/src/Synergistic/tests/CMakeLists.txt +++ b/src/Synergistic/tests/CMakeLists.txt @@ -56,6 +56,13 @@ ADD_TEST(NAME SYN_TEST_CPLUSPLUS COMMAND SYN_TEST_CPLUSPLUS ${test_data} WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) +# Poisson likelihood with gated proj data and motion +if (STIR_WITH_EXPERIMENTAL AND STIR_BUILT_WITH_ITK) + ADD_EXECUTABLE (test_PoissonGatedMotion test_PoissonGatedMotion.cpp ${STIR_REGISTRIES}) + TARGET_LINK_LIBRARIES(test_PoissonGatedMotion PUBLIC cstir Reg) + ADD_TEST(NAME SYN_TEST_CPLUSPLUS_POISSONWMOTION COMMAND test_PoissonGatedMotion WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) +endif() + # If python if (BUILD_PYTHON) # Make into test diff --git a/src/Synergistic/tests/test_PoissonGatedMotion.cpp b/src/Synergistic/tests/test_PoissonGatedMotion.cpp new file mode 100644 index 000000000..c9f5ac172 --- /dev/null +++ b/src/Synergistic/tests/test_PoissonGatedMotion.cpp @@ -0,0 +1,113 @@ +/* +SyneRBI Synergistic Image Reconstruction Framework (SIRF) +Copyright 2020 University College London + +This is software developed for the Collaborative Computational +Project in Synergistic Reconstruction for Biomedical Imaging (formerly CCP PETMR) +(http://www.ccpsynerbi.ac.uk/). + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at +http://www.apache.org/licenses/LICENSE-2.0 +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. + +*/ + +/*! +\file +\ingroup STIR Tests + +\author Richard Brown +\author SyneRBI +*/ + +// STIR stuff +#include "sirf/STIR/stir_x.h" +#include "sirf/common/getenv.h" + +// Reg stuff +#include "sirf/Reg/NiftiImageData3DDisplacement.h" + +using namespace sirf; +using namespace stir; + +int test4() +{ + try { + + std::string SIRF_path = sirf::getenv("SIRF_PATH"); + if (SIRF_path.length() < 1) { + std::cout << "SIRF_PATH not defined, cannot find data" << std::endl; + return 1; + } + + // Open acquisition data + std::string path = SIRF_path + "/data/examples/PET/"; + std::string f_acq_data = path + "Utahscat600k_ca_seg4.hs"; + auto acq_data_sptr = std::make_shared(f_acq_data.c_str()); + + // Recon image + auto im_sptr = std::make_shared(*acq_data_sptr); + im_sptr->fill(0.f); + + // Create blank displacement field and convert to STIR format + NiftiImageData3DDisplacement disp; + disp.create_from_3D_image(*im_sptr); + const std::string temp_fname = "temp_disp.nii"; + disp.write(temp_fname); + auto disp_sptr = MAKE_SHARED(temp_fname, 3 /*bspline_order*/); + + // create additive term + auto a_sptr = acq_data_sptr->new_acquisition_data(); + a_sptr->fill(0.05f); + // create background term + auto b_sptr = acq_data_sptr->new_acquisition_data(); + b_sptr->fill(0.1f); + + // Create acquisition model + auto matrix_sptr = MAKE_SHARED(); + matrix_sptr->set_num_tangential_LORs(2); + auto am_sptr = MAKE_SHARED(); + am_sptr->set_matrix(matrix_sptr); + am_sptr->set_additive_term(a_sptr); + am_sptr->set_background_term(b_sptr); + am_sptr->set_up(acq_data_sptr, im_sptr); + + // Create poisson loglikelihood + PoissonLogLhLinModMeanGatedProjDataWMotion3DF obj_fn; + + // Add motion gates + const unsigned num_motion_gates = 4; + for (unsigned i=0; idata_sptr()); + + // Recon + OSMAPOSLReconstruction recon; + recon.set_up(im_sptr->data_sptr()); + recon.set_num_subiterations(2); + recon.reconstruct(im_sptr->data_sptr()); + + std::cout << "\nSuccess!\n"; + + return 0; + } + catch (...) + { + return 1; + } +} + +//int test5(); + +int main() +{ + return test4(); + //return test5(); +} diff --git a/src/xSTIR/cSTIR/cstir.cpp b/src/xSTIR/cSTIR/cstir.cpp index 172035553..499b74d90 100644 --- a/src/xSTIR/cSTIR/cstir.cpp +++ b/src/xSTIR/cSTIR/cstir.cpp @@ -89,6 +89,10 @@ void* cSTIR_newObject(const char* name) "PoissonLogLikelihoodWithLinearModelForMeanAndProjData")) return NEW_OBJECT_HANDLE (xSTIR_PoissonLogLikelihoodWithLinearModelForMeanAndProjData3DF); + if (boost::iequals(name, + "PoissonLogLhLinModMeanGatedProjDataWMotion3DF")) + return NEW_OBJECT_HANDLE + (PoissonLogLhLinModMeanGatedProjDataWMotion3DF); if (boost::iequals(name, "AcqModUsingMatrix")) return NEW_OBJECT_HANDLE(AcqModUsingMatrix3DF); #ifdef STIR_WITH_NiftyPET_PROJECTOR @@ -884,6 +888,37 @@ cSTIR_objectiveFunctionGradientNotDivided(void* ptr_f, void* ptr_i, int subset) CATCH; } +#ifdef STIR_GATED_MOTION +extern "C" +void* +cSTIR_PoissonGatedWMotion_add_gate( + void* poisson_ptr, void* ad_ptr, void* am_ptr, + const char* disp_fname, const int b_spline_order) +{ + try { + SPTR_FROM_HANDLE(PoissonLogLhLinModMeanGatedProjDataWMotion3DF, + poisson_sptr, poisson_ptr); + SPTR_FROM_HANDLE(PETAcquisitionData, ad_sptr, ad_ptr); + SPTR_FROM_HANDLE(PETAcquisitionModel, am_sptr, am_ptr); + // set it + poisson_sptr->add_gate(ad_sptr, am_sptr, disp_fname, b_spline_order); + } + CATCH; +} + +extern "C" +void* +cSTIR_PoissonGatedWMotion_clear_gates(void* poisson_ptr) +{ + try { + SPTR_FROM_HANDLE(PoissonLogLhLinModMeanGatedProjDataWMotion3DF, + poisson_sptr, poisson_ptr); + poisson_sptr->clear_gates(); + } + CATCH; +} +#endif + extern "C" void* cSTIR_setupPrior(void* ptr_p, void* ptr_i) diff --git a/src/xSTIR/cSTIR/include/sirf/STIR/cstir.h b/src/xSTIR/cSTIR/include/sirf/STIR/cstir.h index 5969ca83d..a8ef8fdce 100644 --- a/src/xSTIR/cSTIR/include/sirf/STIR/cstir.h +++ b/src/xSTIR/cSTIR/include/sirf/STIR/cstir.h @@ -119,6 +119,14 @@ extern "C" { void* cSTIR_objectiveFunctionGradientNotDivided (void* ptr_f, void* ptr_i, int subset); +#ifdef STIR_GATED_MOTION + // Poisson gated proj data w motion + void* cSTIR_PoissonGatedWMotion_add_gate( + void* poisson_ptr, void* ad_ptr, void* am_ptr, + const char* disp_fname, const int b_spline_order); + void* cSTIR_PoissonGatedWMotion_clear_gates(void* poisson_ptr); +#endif + // Prior methods void* cSTIR_setupPrior(void* ptr_p, void* ptr_i); void* cSTIR_priorGradient(void* ptr_p, void* ptr_i); diff --git a/src/xSTIR/cSTIR/include/sirf/STIR/stir_types.h b/src/xSTIR/cSTIR/include/sirf/STIR/stir_types.h index 885606221..489317086 100644 --- a/src/xSTIR/cSTIR/include/sirf/STIR/stir_types.h +++ b/src/xSTIR/cSTIR/include/sirf/STIR/stir_types.h @@ -67,6 +67,11 @@ limitations under the License. #ifdef STIR_WITH_NiftyPET_PROJECTOR #include "stir/recon_buildblock/NiftyPET_projector/ProjectorByBinPairUsingNiftyPET.h" #endif +#ifdef STIR_GATED_MOTION +#include "stir_experimental/motion/PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion.h" +#include "stir_experimental/motion/NonRigidObjectTransformationUsingBSplines.h" +#include "stir_experimental/motion/Transform3DObjectImageProcessor.h" +#endif #include "stir/StirException.h" #include "stir/TextWriter.h" @@ -98,6 +103,12 @@ namespace sirf { typedef stir::QuadraticPrior QuadPrior3DF; typedef stir::DataProcessor DataProcessor3DF; typedef stir::TruncateToCylindricalFOVImageProcessor CylindricFilter3DF; +#ifdef STIR_GATED_MOTION + typedef stir::ObjectTransformation<3,float> STIRTrans3DF; + typedef stir::NonRigidObjectTransformationUsingBSplines<3,float> STIRDisp3DF; + typedef stir::PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotionNew < Image3DF > + stirPoissonLogLhLinModMeanGatedWMotion3DF; +#endif } diff --git a/src/xSTIR/cSTIR/include/sirf/STIR/stir_x.h b/src/xSTIR/cSTIR/include/sirf/STIR/stir_x.h index 51d4225c8..324bda7ab 100644 --- a/src/xSTIR/cSTIR/include/sirf/STIR/stir_x.h +++ b/src/xSTIR/cSTIR/include/sirf/STIR/stir_x.h @@ -400,8 +400,8 @@ The actual algorithm is described in } virtual stir::Succeeded set_up( - stir::shared_ptr sptr_acq, - stir::shared_ptr sptr_image); + const stir::shared_ptr sptr_acq, + const stir::shared_ptr sptr_image); // computes and returns a subset of forward-projected data stir::shared_ptr @@ -456,8 +456,8 @@ The actual algorithm is described in get_proj_matrix_sptr(); } virtual stir::Succeeded set_up( - stir::shared_ptr sptr_acq, - stir::shared_ptr sptr_image) + const stir::shared_ptr sptr_acq, + const stir::shared_ptr sptr_image) { if (!sptr_matrix_.get()) return stir::Succeeded::no; @@ -591,6 +591,43 @@ The actual algorithm is described in typedef xSTIR_PoissonLogLikelihoodWithLinearModelForMeanAndProjData3DF PoissonLogLhLinModMeanProjData3DF; +#ifdef STIR_GATED_MOTION + class xSTIR_PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion3DF : + public stirPoissonLogLhLinModMeanGatedWMotion3DF { + private: + typedef stirPoissonLogLhLinModMeanGatedWMotion3DF base_type; + public: + + /// Add a motion gate. Required: acq data, acq model and disp field + void add_gate( + stir::shared_ptr ad_sptr, + stir::shared_ptr am_sptr, + stir::shared_ptr trans_sptr); + + /// Add a motion gate. Required: acq data, acq model, + /// disp field (as string) and the b-spline order + void add_gate( + stir::shared_ptr ad_sptr, + stir::shared_ptr am_sptr, + const std::string &disp_fname, + const int b_spline_order); + + /// Set up + virtual stir::Succeeded set_up(stir::shared_ptr const& target_sptr); + + /// Clear all gates + void clear_gates(); + + private: + std::vector > _ad_vec; + std::vector > _am_vec; + std::vector > _trans_vec; + }; + + typedef xSTIR_PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion3DF + PoissonLogLhLinModMeanGatedProjDataWMotion3DF; +#endif + class xSTIR_IterativeReconstruction3DF : public stir::IterativeReconstruction < Image3DF > { public: @@ -661,7 +698,7 @@ The actual algorithm is described in ("wrong frequency cut-off", __FILE__, __LINE__); fc_ramp = fc; } - stir::Succeeded set_up(stir::shared_ptr sptr_id) + stir::Succeeded set_up(const stir::shared_ptr sptr_id) { _sptr_image_data.reset(new STIRImageData(*sptr_id)); _is_set_up = true; diff --git a/src/xSTIR/cSTIR/stir_x.cpp b/src/xSTIR/cSTIR/stir_x.cpp index 0d3dda717..270db6631 100644 --- a/src/xSTIR/cSTIR/stir_x.cpp +++ b/src/xSTIR/cSTIR/stir_x.cpp @@ -523,8 +523,8 @@ PETAttenuationModel::normalise(PETAcquisitionData& ad) const Succeeded PETAcquisitionModel::set_up( - shared_ptr sptr_acq, - shared_ptr sptr_image) + const shared_ptr sptr_acq, + const shared_ptr sptr_image) { Succeeded s = Succeeded::no; if (sptr_projectors_.get()) { @@ -668,3 +668,110 @@ PETAcquisitionModel::backward(PETAcquisitionData& ad, return sptr_id; } + +void +PoissonLogLhLinModMeanGatedProjDataWMotion3DF:: +add_gate( + stir::shared_ptr ad_sptr, + stir::shared_ptr am_sptr, + stir::shared_ptr trans_sptr) +{ + _ad_vec.push_back(ad_sptr); + _am_vec.push_back(am_sptr); + _trans_vec.push_back(trans_sptr); +} + +void +PoissonLogLhLinModMeanGatedProjDataWMotion3DF:: +add_gate( + stir::shared_ptr ad_sptr, + stir::shared_ptr am_sptr, + const std::string &disp_fname, + const int b_spline_order) +{ + auto disp_sptr = MAKE_SHARED(disp_fname, b_spline_order); + add_gate(ad_sptr, am_sptr, disp_sptr); +} + +static stir::shared_ptr +vec_AcqData_to_GatedProjData(const std::vector > &vec_ad) +{ + auto sino_sptr = MAKE_SHARED(); + if (!vec_ad.empty()) { + sino_sptr->set_exam_info(*vec_ad.at(0)->get_exam_info_sptr()); + sino_sptr->resize(vec_ad.size()); + for (unsigned i=0; iset_proj_data_sptr(vec_ad.at(i)->data(), i+1); + } + return sino_sptr; +} + +template +void check_all_or_none(const std::vector > &vec, const std::string &explanation) +{ + const bool first_is_null_ptr = vec.at(0) == nullptr; + for (unsigned i=1; i &target_sptr) +{ + const unsigned num_gates = _ad_vec.size(); + if (num_gates==0) + throw std::runtime_error("PoissonLogLhLinModMeanGatedProjDataWMotion3DF::set_up: no gates added."); + + // Projector pair + set_projector_pair_sptr(_am_vec.at(0)->projectors_sptr()); + + // Raw data: vector -> GatedProjData + auto sino_sptr = vec_AcqData_to_GatedProjData(_ad_vec); + set_proj_data_sptr(sino_sptr); + + // Additives + std::vector > add_vec; + for (unsigned i=0; iadditive_term_sptr()); + check_all_or_none(add_vec, "additive"); + auto add_sptr = vec_AcqData_to_GatedProjData(add_vec); + set_additive_proj_data_sptr(add_sptr); + + // Norms + std::vector > norm_vec; + for (unsigned i=0; inormalisation_sptr()); + check_all_or_none(norm_vec, "norms"); + if (norm_vec.at(0)) { + stir::VectorWithOffset > norm_stir_vec(1, num_gates); + for (unsigned i=1; i<=num_gates; ++i) + norm_stir_vec.at(i) = norm_vec.at(i-1); + set_normalisations(norm_stir_vec); + } + + // Transformations + stir::VectorWithOffset > trans_stir_vec(1,num_gates); + for (unsigned i=0; i >(_trans_vec.at(i)); + set_forward_transformations(trans_stir_vec); + + // Blank TimeFrameDefinitions + std::vector > frame_times(num_gates); + for (unsigned i=0; i(i,i+1); + set_frame_definitions(frame_times); + + // Call the base type set up + return base_type::set_up(target_sptr); +} + +void +PoissonLogLhLinModMeanGatedProjDataWMotion3DF:: +clear_gates() +{ + _ad_vec.clear(); + _am_vec.clear(); + _trans_vec.clear(); +} From 4403de6a53bdbc31ffcb6343f0f8ec11affa8e92 Mon Sep 17 00:00:00 2001 From: richard Date: Thu, 18 Jun 2020 11:20:45 +0100 Subject: [PATCH 02/14] poisson gated w motion: python --- examples/Python/PET/osem_gated_w_motion.py | 336 +++++++++++++++++++++ src/xSTIR/pSTIR/STIR.py.in | 34 +++ 2 files changed, 370 insertions(+) create mode 100644 examples/Python/PET/osem_gated_w_motion.py diff --git a/examples/Python/PET/osem_gated_w_motion.py b/examples/Python/PET/osem_gated_w_motion.py new file mode 100644 index 000000000..45f771711 --- /dev/null +++ b/examples/Python/PET/osem_gated_w_motion.py @@ -0,0 +1,336 @@ +'''OSEM reconstruction demo for gated data with motion. +We actually use the OSMAPOSL reconstructor in this demo. This reconstructor +implements an Ordered Subsets (OS) version of the One Step Late algorithm (OSL) +from Green et al for Maximum a Posteriori (MAP) maximisation. Here we use it +for Maximum Likelihood (ML) in which case it is equivalent to OSEM. + +Usage: + osem_reconstruction_gated_w_motion [--help | options] + +Options: + -S , --sino= sinogram pattern, * or % wildcard + (e.g., sino_ms*.hs). Enclose in quotations. + -a , --attn= attenuation pattern, * or % wildcard + (e.g., attn_ms*.hv). Enclose in quotations. + -R , --rand= randoms pattern, * or % wildcard + (e.g., rand_ms*.hs). Enclose in quotations. + -n , --norm= ECAT8 bin normalization file + -T , --trans= transformation pattern, * or % wildcard + (e.g., tm_ms*.txt). Enclose in quotations. + -t , --trans_type= transformation type (tm, disp, def) + [default: tm] + -s , --subs= number of subsets [default: 21] + -i , --subiter= number of sub-iterations [default: 2] + -e , --engine= reconstruction engine [default: STIR] + -o , --outp= output file prefix [default: recon] + --initial= initial estimate + --nxny= image x and y dimension [default: 127] + --dxdy= image x and y spacing + (default: determined by scanner) + --b_spline_order= b-spline order for resampler [default: 1] + --verbosity= Verbosity [default: 0] + --visualisations show visualisations + --nifti save output as nifti + --gpu use gpu +''' + +# SyneRBI Synergistic Image Reconstruction Framework (SIRF) +# Copyright 2015 - 2018 Rutherford Appleton Laboratory STFC +# Copyright 2015 - 2020 University College London. +# +# This is software developed for the Collaborative Computational +# Project in Synergistic Reconstruction for Biomedical Imaging +# (formerly CCP PETMR) +# (http://www.ccpsynerbi.ac.uk/). +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# http://www.apache.org/licenses/LICENSE-2.0 +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +import os +from glob import glob +from docopt import docopt +import sirf.Reg as reg + + +__version__ = '0.1.0' +args = docopt(__doc__, version=__version__) + +# import engine module +if args['--engine'] == 'STIR': + import sirf.STIR as pet + +pet.AcquisitionData.set_storage_scheme('memory') + +# Multiple files +sino_pattern = str(args['--sino']).replace('%', '*') +attn_pattern = str(args['--attn']).replace('%', '*') +rand_pattern = str(args['--rand']).replace('%', '*') +trans_pattern = str(args['--trans']).replace('%', '*') +trans_type = str(args['--trans_type']) + +if attn_pattern is None: + attn_pattern = "" +if rand_pattern is None: + rand_pattern = "" + +# Norm +norm_file = None +if args['--norm']: + norm_file = str(args['--norm']) + if not os.path.isfile(norm_file): + raise error("Norm file not found: " + norm_file) + +# Initial estimate +initial_estimate = None +if args['--initial']: + initial_estimate = str(args['--initial']) + +# Verbosity +verbosity = int(args['--verbosity']) +pet.set_verbosity(verbosity) +if verbosity == 0: + msg_red = pet.MessageRedirector(None, None, None) + +# Rest of the input options +num_subsets = int(args['--subs']) +num_subiterations = int(args['--subiter']) +nxny = int(args['--nxny']) +outp_file = args['--outp'] +b_spline_order = args['--b_spline_order'] +visualisations = True if args['--visualisations'] else False +nifti = True if args['--nifti'] else False +use_gpu = True if args['--gpu'] else False + + +def get_resampler(image, ref=None, trans=None): + """returns a NiftyResample object for the specified transform and image""" + if ref is None: + ref = image + resampler = reg.NiftyResample() + resampler.set_reference_image(ref) + resampler.set_floating_image(image) + resampler.set_padding_value(0) + resampler.set_interpolation_type_to_linear() + if trans is not None: + resampler.add_transformation(trans) + return resampler + + +def get_asm_attn(sino, attn, acq_model): + """Get attn ASM from sino, attn image and acq model""" + asm_attn = pet.AcquisitionSensitivityModel(attn, acq_model) + # temporary fix pending attenuation offset fix in STIR: + # converting attenuation into 'bin efficiency' + asm_attn.set_up(sino) + bin_eff = pet.AcquisitionData(sino) + bin_eff.fill(1.0) + asm_attn.unnormalise(bin_eff) + asm_attn = pet.AcquisitionSensitivityModel(bin_eff) + return asm_attn + + +def main(): + + ########################################################################### + # Parse input files + ########################################################################### + + if trans_pattern is None: + raise AssertionError("--trans missing") + if sino_pattern is None: + raise AssertionError("--sino missing") + trans_files = sorted(glob(trans_pattern)) + sino_files = sorted(glob(sino_pattern)) + attn_files = sorted(glob(attn_pattern)) + rand_files = sorted(glob(rand_pattern)) + + num_ms = len(sino_files) + # Check some sinograms found + if num_ms == 0: + raise AssertionError("No sinograms found!") + # Should have as many trans as sinos + if num_ms != len(trans_files): + raise AssertionError("#trans should match #sinos. " + "#sinos = " + str(num_ms) + + ", #trans = " + str(len(trans_files))) + # If any rand, check num == num_ms + if len(rand_files) > 0 and len(rand_files) != num_ms: + raise AssertionError("#rand should match #sinos. " + "#sinos = " + str(num_ms) + + ", #rand = " + str(len(rand_files))) + + # For attn, there should be 0, 1 or num_ms images + if len(attn_files) > 1 and len(attn_files) != num_ms: + raise AssertionError("#attn should be 0, 1 or #sinos") + + ########################################################################### + # Read input + ########################################################################### + + trans = num_ms*[None] + for i in range(num_ms): + if trans_type == "disp": + trans[i] = reg.NiftiImageData3DDisplacement(trans_files[i]) + elif trans_type == "tm": + trans[i] = reg.NiftiImageData3DDisplacement( + reg.AffineTransformation(trans_files[i])) + elif trans_type == "def": + trans[i] = reg.NiftiImageData3DDisplacement( + reg.NiftiImageData3DDeformation(trans_files[i])) + else: + raise error("Unknown transformation type") + + sinos = [pet.AcquisitionData(file) for file in sino_files] + attns = [pet.ImageData(file) for file in attn_files] + rands = [pet.AcquisitionData(file) for file in rand_files] + + ########################################################################### + # Initialise recon image + ########################################################################### + + if initial_estimate: + image = pet.ImageData(initial_estimate) + else: + image = sinos[0].create_uniform_image(0.0, nxny) + # If using GPU, need to make sure that image is right size. + if use_gpu: + image.initialise(dim=(127, 320, 320), + vsize=(2.03125, 2.08626, 2.08626)) + image.fill(0.0) + + ########################################################################### + # Resample attenuation images (if necessary) + ########################################################################### + + resampled_attns = None + if len(attns) > 0: + resampled_attns = [0]*num_ms + # if using GPU, dimensions of attn and recon images have to match + ref = image if use_gpu else None + for i in range(len(attns)): + # if we only have 1 attn image, then we need to resample into + # space of each gate. However, if we have num_ms attn images, then + # assume they are already in the correct position, so use None as + # transformation. + tran = trans[i] if len(attns) == 1 else None + # If only 1 attn image, then resample that. If we have num_ms attn + # images, then use each attn image of each frame. + attn = attns[0] if len(attns) == 1 else attns[i] + resam = get_resampler(attn, ref=ref, trans=tran) + resampled_attns[i] = resam.forward(attn) + + ########################################################################### + # Set up acquisition models + ########################################################################### + + print("Setting up acquisition models...") + if not use_gpu: + acq_models = num_ms * [pet.AcquisitionModelUsingRayTracingMatrix()] + else: + acq_models = num_ms * [pet.AcquisitionModelUsingNiftyPET()] + for acq_model in acq_models: + acq_model.set_use_truncation(True) + acq_model.set_cuda_verbosity(verbosity) + + # If present, create ASM from ECAT8 normalisation data + asm_norm = None + if norm_file: + asm_norm = pet.AcquisitionSensitivityModel(norm_file) + + # Loop over each motion state + for ind in range(num_ms): + # Create attn ASM if necessary + asm_attn = None + if resampled_attns: + asm_attn = get_asm_attn(sinos[ind], resampled_attns[i], + acq_models[ind]) + + # Get ASM dependent on attn and/or norm + asm = None + if asm_norm and asm_attn: + if ind == 0: + print("ASM contains norm and attenuation...") + asm = pet.AcquisitionSensitivityModel(asm_norm, asm_attn) + elif asm_norm: + if ind == 0: + print("ASM contains norm...") + asm = asm_norm + elif asm_attn: + if ind == 0: + print("ASM contains attenuation...") + asm = asm_attn + if asm: + acq_models[ind].set_acquisition_sensitivity(asm) + + if len(rands) > 0: + acq_models[ind].set_background_term(rands[ind]) + + # Set up + acq_models[ind].set_up(sinos[ind], image) + + ########################################################################### + # Set up reconstructor + ########################################################################### + + print("Setting up reconstructor...") + + # define objective function to be maximized as + # Poisson logarithmic likelihood (with linear model for mean) + obj_fun = pet.PoissonLogLhLinModMeanGatedProjDataWMotion() + + # Loop over each motion state and set + for ms in num_ms: + obj_fun.add_gate(sinos[i], acq_models[i], trans[i]) + + # select Ordered Subsets Maximum A-Posteriori One Step Late as the + # reconstruction algorithm (since we are not using a penalty, or prior, in + # this example, we actually run OSEM); + # this algorithm does not converge to the maximum of the objective function + # but is used in practice to speed-up calculations + recon = OSMAPOSLReconstructor() + recon.set_objective_function(obj_fun) + recon.set_num_subsets(num_subsets) + recon.set_num_subiterations(num_subiterations) + recon.set_input(acq_data) + + # set up the reconstructor based on a sample image + # (checks the validity of parameters, sets up objective function + # and other objects involved in the reconstruction, which involves + # computing/reading sensitivity image etc etc.) + print('setting up, please wait...') + recon.set_up(image) + + # set the initial image estimate + recon.set_current_estimate(image) + + # reconstruct + print('reconstructing, please wait...') + recon.process() + out = recon.get_output() + if not args['--nifti']: + out.write(outp_file) + else: + sirf.Reg.NiftiImageData(out).write(outp_file) + + if visualisations: + # show reconstructed image + image_array = out.as_array() + show_2D_array('Reconstructed image', image_array[z, :, :]) + pylab.show() + + +# if anything goes wrong, an exception will be thrown +# (cf. Error Handling section in the spec) +try: + main() + print('done') +except pet.error as err: + # display error information + print('%s' % err.value) diff --git a/src/xSTIR/pSTIR/STIR.py.in b/src/xSTIR/pSTIR/STIR.py.in index 5e9d2412f..cf4de1c7d 100644 --- a/src/xSTIR/pSTIR/STIR.py.in +++ b/src/xSTIR/pSTIR/STIR.py.in @@ -1789,6 +1789,40 @@ class PoissonLogLikelihoodWithLinearModelForMeanAndProjData\ parms.set_parameter\ (self.handle, self.name, 'acquisition_data', ad.handle) + +if @STIR_GATED_MOTION@: + class PoissonLogLhLinModMeanGatedProjDataWMotion(ObjectiveFunction): + """ + Class for STIR PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion object. + There are two versions of this class in the STIR code base. The one used here lives in the + experimental part, and was chosen becase it uses the adjoint for the motion, as opposed to the + version in the main code base, which uses the inverse. + """ + def __init__(self): + """Init""" + self.handle = None + self.name = 'PoissonLogLhLinModMeanGatedProjDataWMotion3DF' + self.handle = pystir.cSTIR_newObject(self.name) + check_status(self.handle) + def __del__(self): + if self.handle is not None: + pyiutil.deleteDataHandle(self.handle) + def add_gate(self, acq_data, acq_model, disp, b_spline_order=1): + """Add a gate by giving an acquisition data, acquisition model, + displacement filename, and optionally b-spline order.""" + assert_validity(acq_data, AcquisitionData) + assert_validity(acq_model, AcquisitionModel) + assert_validity(disp, sirf.Reg.NiftiImageData3DDisplacement) + disp_fname = "temp_disp_" + uuid.uuid4().hex + ".nii" + disp.write(disp_fname) + try_calling(pystir.cSTIR_PoissonGatedWMotion_add_gate( + self.handle, acq_data.handle, acq_model.handle, disp_fname, b_spline_order)) + os.remove(disp_fname) + def add_gate(self): + """Clear all gates.""" + try_calling(pystir.cSTIR_PoissonGatedWMotion_clear_gates(self.handle)) + + class Reconstructor(object): ''' Class for a generic PET reconstructor. From b0acebd40d1005f8306a7956051b0aab05a03b99 Mon Sep 17 00:00:00 2001 From: Richard Date: Thu, 18 Jun 2020 10:53:39 +0000 Subject: [PATCH 03/14] extra ifdefs for when stir not built with experimental --- src/xSTIR/cSTIR/cstir.cpp | 2 ++ src/xSTIR/cSTIR/stir_x.cpp | 2 ++ 2 files changed, 4 insertions(+) diff --git a/src/xSTIR/cSTIR/cstir.cpp b/src/xSTIR/cSTIR/cstir.cpp index 499b74d90..7de9ebdb3 100644 --- a/src/xSTIR/cSTIR/cstir.cpp +++ b/src/xSTIR/cSTIR/cstir.cpp @@ -89,10 +89,12 @@ void* cSTIR_newObject(const char* name) "PoissonLogLikelihoodWithLinearModelForMeanAndProjData")) return NEW_OBJECT_HANDLE (xSTIR_PoissonLogLikelihoodWithLinearModelForMeanAndProjData3DF); +#ifdef STIR_GATED_MOTION if (boost::iequals(name, "PoissonLogLhLinModMeanGatedProjDataWMotion3DF")) return NEW_OBJECT_HANDLE (PoissonLogLhLinModMeanGatedProjDataWMotion3DF); +#endif if (boost::iequals(name, "AcqModUsingMatrix")) return NEW_OBJECT_HANDLE(AcqModUsingMatrix3DF); #ifdef STIR_WITH_NiftyPET_PROJECTOR diff --git a/src/xSTIR/cSTIR/stir_x.cpp b/src/xSTIR/cSTIR/stir_x.cpp index 270db6631..5d98e4044 100644 --- a/src/xSTIR/cSTIR/stir_x.cpp +++ b/src/xSTIR/cSTIR/stir_x.cpp @@ -669,6 +669,7 @@ PETAcquisitionModel::backward(PETAcquisitionData& ad, return sptr_id; } +#ifdef STIR_GATED_MOTION void PoissonLogLhLinModMeanGatedProjDataWMotion3DF:: add_gate( @@ -775,3 +776,4 @@ clear_gates() _am_vec.clear(); _trans_vec.clear(); } +#endif From 4ff86f532781c95e49fa5f7e2710eb7da9259152 Mon Sep 17 00:00:00 2001 From: Richard Date: Thu, 18 Jun 2020 11:11:59 +0000 Subject: [PATCH 04/14] codacy changes1 --- examples/Python/PET/osem_gated_w_motion.py | 22 +++++++--------------- src/xSTIR/pSTIR/STIR.py.in | 2 +- 2 files changed, 8 insertions(+), 16 deletions(-) diff --git a/examples/Python/PET/osem_gated_w_motion.py b/examples/Python/PET/osem_gated_w_motion.py index 45f771711..79fd16432 100644 --- a/examples/Python/PET/osem_gated_w_motion.py +++ b/examples/Python/PET/osem_gated_w_motion.py @@ -1,4 +1,5 @@ '''OSEM reconstruction demo for gated data with motion. + We actually use the OSMAPOSL reconstructor in this demo. This reconstructor implements an Ordered Subsets (OS) version of the One Step Late algorithm (OSL) from Green et al for Maximum a Posteriori (MAP) maximisation. Here we use it @@ -29,7 +30,6 @@ (default: determined by scanner) --b_spline_order= b-spline order for resampler [default: 1] --verbosity= Verbosity [default: 0] - --visualisations show visualisations --nifti save output as nifti --gpu use gpu ''' @@ -85,7 +85,7 @@ if args['--norm']: norm_file = str(args['--norm']) if not os.path.isfile(norm_file): - raise error("Norm file not found: " + norm_file) + raise pet.error("Norm file not found: " + norm_file) # Initial estimate initial_estimate = None @@ -104,13 +104,12 @@ nxny = int(args['--nxny']) outp_file = args['--outp'] b_spline_order = args['--b_spline_order'] -visualisations = True if args['--visualisations'] else False nifti = True if args['--nifti'] else False use_gpu = True if args['--gpu'] else False def get_resampler(image, ref=None, trans=None): - """returns a NiftyResample object for the specified transform and image""" + """returns a NiftyResample object for the specified transform and image.""" if ref is None: ref = image resampler = reg.NiftyResample() @@ -124,7 +123,7 @@ def get_resampler(image, ref=None, trans=None): def get_asm_attn(sino, attn, acq_model): - """Get attn ASM from sino, attn image and acq model""" + """Get attn ASM from sino, attn image and acq model.""" asm_attn = pet.AcquisitionSensitivityModel(attn, acq_model) # temporary fix pending attenuation offset fix in STIR: # converting attenuation into 'bin efficiency' @@ -185,7 +184,7 @@ def main(): trans[i] = reg.NiftiImageData3DDisplacement( reg.NiftiImageData3DDeformation(trans_files[i])) else: - raise error("Unknown transformation type") + raise pet.error("Unknown transformation type") sinos = [pet.AcquisitionData(file) for file in sino_files] attns = [pet.ImageData(file) for file in attn_files] @@ -286,7 +285,7 @@ def main(): obj_fun = pet.PoissonLogLhLinModMeanGatedProjDataWMotion() # Loop over each motion state and set - for ms in num_ms: + for i in range(num_ms): obj_fun.add_gate(sinos[i], acq_models[i], trans[i]) # select Ordered Subsets Maximum A-Posteriori One Step Late as the @@ -298,7 +297,6 @@ def main(): recon.set_objective_function(obj_fun) recon.set_num_subsets(num_subsets) recon.set_num_subiterations(num_subiterations) - recon.set_input(acq_data) # set up the reconstructor based on a sample image # (checks the validity of parameters, sets up objective function @@ -317,13 +315,7 @@ def main(): if not args['--nifti']: out.write(outp_file) else: - sirf.Reg.NiftiImageData(out).write(outp_file) - - if visualisations: - # show reconstructed image - image_array = out.as_array() - show_2D_array('Reconstructed image', image_array[z, :, :]) - pylab.show() + reg.NiftiImageData(out).write(outp_file) # if anything goes wrong, an exception will be thrown diff --git a/src/xSTIR/pSTIR/STIR.py.in b/src/xSTIR/pSTIR/STIR.py.in index cf4de1c7d..6e59a3e6f 100644 --- a/src/xSTIR/pSTIR/STIR.py.in +++ b/src/xSTIR/pSTIR/STIR.py.in @@ -1818,7 +1818,7 @@ if @STIR_GATED_MOTION@: try_calling(pystir.cSTIR_PoissonGatedWMotion_add_gate( self.handle, acq_data.handle, acq_model.handle, disp_fname, b_spline_order)) os.remove(disp_fname) - def add_gate(self): + def clear_gates(self): """Clear all gates.""" try_calling(pystir.cSTIR_PoissonGatedWMotion_clear_gates(self.handle)) From 5d7d775b206fc0b851bf0cee9a850d7b907fc9a6 Mon Sep 17 00:00:00 2001 From: Richard Date: Thu, 18 Jun 2020 11:36:27 +0000 Subject: [PATCH 05/14] codacy changes 2 --- examples/Python/PET/osem_gated_w_motion.py | 10 +++++----- src/xSTIR/pSTIR/STIR.py.in | 6 +++++- 2 files changed, 10 insertions(+), 6 deletions(-) diff --git a/examples/Python/PET/osem_gated_w_motion.py b/examples/Python/PET/osem_gated_w_motion.py index 79fd16432..78f68f893 100644 --- a/examples/Python/PET/osem_gated_w_motion.py +++ b/examples/Python/PET/osem_gated_w_motion.py @@ -1,4 +1,4 @@ -'''OSEM reconstruction demo for gated data with motion. +"""OSEM reconstruction demo for gated data with motion. We actually use the OSMAPOSL reconstructor in this demo. This reconstructor implements an Ordered Subsets (OS) version of the One Step Late algorithm (OSL) @@ -32,7 +32,7 @@ --verbosity= Verbosity [default: 0] --nifti save output as nifti --gpu use gpu -''' +""" # SyneRBI Synergistic Image Reconstruction Framework (SIRF) # Copyright 2015 - 2018 Rutherford Appleton Laboratory STFC @@ -109,7 +109,7 @@ def get_resampler(image, ref=None, trans=None): - """returns a NiftyResample object for the specified transform and image.""" + """Returns a NiftyResample object for the specified transform and image.""" if ref is None: ref = image resampler = reg.NiftyResample() @@ -213,7 +213,7 @@ def main(): resampled_attns = [0]*num_ms # if using GPU, dimensions of attn and recon images have to match ref = image if use_gpu else None - for i in range(len(attns)): + for i in range(num_ms): # if we only have 1 attn image, then we need to resample into # space of each gate. However, if we have num_ms attn images, then # assume they are already in the correct position, so use None as @@ -293,7 +293,7 @@ def main(): # this example, we actually run OSEM); # this algorithm does not converge to the maximum of the objective function # but is used in practice to speed-up calculations - recon = OSMAPOSLReconstructor() + recon = pet.OSMAPOSLReconstructor() recon.set_objective_function(obj_fun) recon.set_num_subsets(num_subsets) recon.set_num_subiterations(num_subiterations) diff --git a/src/xSTIR/pSTIR/STIR.py.in b/src/xSTIR/pSTIR/STIR.py.in index 6e59a3e6f..66887e03a 100644 --- a/src/xSTIR/pSTIR/STIR.py.in +++ b/src/xSTIR/pSTIR/STIR.py.in @@ -42,6 +42,10 @@ import sirf.pystir as pystir import sirf.STIR_params as parms +if @STIR_GATED_MOTION@: + import sirf.Reg as reg + import uuid + try: input = raw_input except NameError: @@ -1812,7 +1816,7 @@ if @STIR_GATED_MOTION@: displacement filename, and optionally b-spline order.""" assert_validity(acq_data, AcquisitionData) assert_validity(acq_model, AcquisitionModel) - assert_validity(disp, sirf.Reg.NiftiImageData3DDisplacement) + assert_validity(disp, reg.NiftiImageData3DDisplacement) disp_fname = "temp_disp_" + uuid.uuid4().hex + ".nii" disp.write(disp_fname) try_calling(pystir.cSTIR_PoissonGatedWMotion_add_gate( From debde8a66b0927570f3082cc2d655b50b30fe4af Mon Sep 17 00:00:00 2001 From: Richard Date: Thu, 18 Jun 2020 11:52:39 +0000 Subject: [PATCH 06/14] codacy changes 3 --- examples/Python/PET/osem_gated_w_motion.py | 115 +++++++++++++++------ 1 file changed, 84 insertions(+), 31 deletions(-) diff --git a/examples/Python/PET/osem_gated_w_motion.py b/examples/Python/PET/osem_gated_w_motion.py index 78f68f893..a08b0b95d 100644 --- a/examples/Python/PET/osem_gated_w_motion.py +++ b/examples/Python/PET/osem_gated_w_motion.py @@ -135,11 +135,8 @@ def get_asm_attn(sino, attn, acq_model): return asm_attn -def main(): - - ########################################################################### - # Parse input files - ########################################################################### +def parse_input_files(): + """Parse input files.""" if trans_pattern is None: raise AssertionError("--trans missing") @@ -169,10 +166,11 @@ def main(): if len(attn_files) > 1 and len(attn_files) != num_ms: raise AssertionError("#attn should be 0, 1 or #sinos") - ########################################################################### - # Read input - ########################################################################### + return num_ms, trans_files, sino_files, attn_files, rand_files + +def read_trans(num_ms, trans_files): + """Read transformations.""" trans = num_ms*[None] for i in range(num_ms): if trans_type == "disp": @@ -186,28 +184,36 @@ def main(): else: raise pet.error("Unknown transformation type") - sinos = [pet.AcquisitionData(file) for file in sino_files] - attns = [pet.ImageData(file) for file in attn_files] - rands = [pet.AcquisitionData(file) for file in rand_files] - - ########################################################################### - # Initialise recon image - ########################################################################### +def get_initial_estimate(): + """Get initial estimate.""" if initial_estimate: image = pet.ImageData(initial_estimate) else: - image = sinos[0].create_uniform_image(0.0, nxny) + # Create image based on ProjData + image = sinos[0].create_uniform_image(0.0, (nxny, nxny)) # If using GPU, need to make sure that image is right size. if use_gpu: - image.initialise(dim=(127, 320, 320), - vsize=(2.03125, 2.08626, 2.08626)) - image.fill(0.0) - - ########################################################################### - # Resample attenuation images (if necessary) - ########################################################################### - + dim = (127, 320, 320) + spacing = (2.03125, 2.08626, 2.08626) + # elif non-default spacing desired + elif args['--dxdy']: + dim = image.dimensions() + dxdy = float(args['--dxdy']) + spacing = (image.voxel_sizes()[0], dxdy, dxdy) + if use_gpu or args['--dxdy']: + image.initialise(dim=dim, + vsize=spacing) + return image + + +def get_attn_images(num_ms, attns, trans, image): + """Get attn images. + If none supplied, return None. + If 1 supplied, resample to each motion state. + If num_ms supplied return them. + If GPU projector, resample attn images for correct dimensions. + """ resampled_attns = None if len(attns) > 0: resampled_attns = [0]*num_ms @@ -224,11 +230,11 @@ def main(): attn = attns[0] if len(attns) == 1 else attns[i] resam = get_resampler(attn, ref=ref, trans=tran) resampled_attns[i] = resam.forward(attn) + return resampled_attns - ########################################################################### - # Set up acquisition models - ########################################################################### +def get_acq_models(num_ms, resampled_attns, sinos, rands): + """Get acquisition models.""" print("Setting up acquisition models...") if not use_gpu: acq_models = num_ms * [pet.AcquisitionModelUsingRayTracingMatrix()] @@ -273,11 +279,11 @@ def main(): # Set up acq_models[ind].set_up(sinos[ind], image) + return acq_models - ########################################################################### - # Set up reconstructor - ########################################################################### +def set_up_reconstructor(sinos, acq_models, trans, image): + """Set up reconstructor.""" print("Setting up reconstructor...") # define objective function to be maximized as @@ -307,8 +313,55 @@ def main(): # set the initial image estimate recon.set_current_estimate(image) + return recon + + +def main(): + + ########################################################################### + # Parse input files + ########################################################################### + + num_ms, trans_files, sino_files, attn_files, rand_files = \ + parse_input_files() + + ########################################################################### + # Read input + ########################################################################### + + trans = read_trans(num_ms, trans_files) + sinos = [pet.AcquisitionData(file) for file in sino_files] + attns = [pet.ImageData(file) for file in attn_files] + rands = [pet.AcquisitionData(file) for file in rand_files] + + ########################################################################### + # Initialise recon image + ########################################################################### + + image = get_initial_estimate() + + ########################################################################### + # Resample attenuation images (if necessary) + ########################################################################### + + resampled_attns = get_attn_images(num_ms, attns, trans, image) + + ########################################################################### + # Set up acquisition models + ########################################################################### + + acq_models = get_acq_models(num_ms, resampled_attns, sinos, rands) + + ########################################################################### + # Set up reconstructor + ########################################################################### + + recon = set_up_reconstructor(sinos, acq_models, trans, image) + + ########################################################################### + # Reconstruct + ########################################################################### - # reconstruct print('reconstructing, please wait...') recon.process() out = recon.get_output() From 2d334f725b7a65eea0d4574f34ce0b7c2ea3b456 Mon Sep 17 00:00:00 2001 From: Richard Date: Thu, 18 Jun 2020 12:10:53 +0000 Subject: [PATCH 07/14] codacy changes 4 --- examples/Python/PET/osem_gated_w_motion.py | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/examples/Python/PET/osem_gated_w_motion.py b/examples/Python/PET/osem_gated_w_motion.py index a08b0b95d..317a2103f 100644 --- a/examples/Python/PET/osem_gated_w_motion.py +++ b/examples/Python/PET/osem_gated_w_motion.py @@ -183,9 +183,10 @@ def read_trans(num_ms, trans_files): reg.NiftiImageData3DDeformation(trans_files[i])) else: raise pet.error("Unknown transformation type") + return trans -def get_initial_estimate(): +def get_initial_estimate(sinos): """Get initial estimate.""" if initial_estimate: image = pet.ImageData(initial_estimate) @@ -209,6 +210,7 @@ def get_initial_estimate(): def get_attn_images(num_ms, attns, trans, image): """Get attn images. + If none supplied, return None. If 1 supplied, resample to each motion state. If num_ms supplied return them. @@ -233,7 +235,7 @@ def get_attn_images(num_ms, attns, trans, image): return resampled_attns -def get_acq_models(num_ms, resampled_attns, sinos, rands): +def get_acq_models(num_ms, resampled_attns, sinos, rands, image): """Get acquisition models.""" print("Setting up acquisition models...") if not use_gpu: @@ -254,7 +256,7 @@ def get_acq_models(num_ms, resampled_attns, sinos, rands): # Create attn ASM if necessary asm_attn = None if resampled_attns: - asm_attn = get_asm_attn(sinos[ind], resampled_attns[i], + asm_attn = get_asm_attn(sinos[ind], resampled_attns[ind], acq_models[ind]) # Get ASM dependent on attn and/or norm @@ -282,7 +284,7 @@ def get_acq_models(num_ms, resampled_attns, sinos, rands): return acq_models -def set_up_reconstructor(sinos, acq_models, trans, image): +def set_up_reconstructor(num_ms, sinos, acq_models, trans, image): """Set up reconstructor.""" print("Setting up reconstructor...") @@ -338,7 +340,7 @@ def main(): # Initialise recon image ########################################################################### - image = get_initial_estimate() + image = get_initial_estimate(sinos) ########################################################################### # Resample attenuation images (if necessary) @@ -350,13 +352,13 @@ def main(): # Set up acquisition models ########################################################################### - acq_models = get_acq_models(num_ms, resampled_attns, sinos, rands) + acq_models = get_acq_models(num_ms, resampled_attns, sinos, rands, image) ########################################################################### # Set up reconstructor ########################################################################### - recon = set_up_reconstructor(sinos, acq_models, trans, image) + recon = set_up_reconstructor(num_ms, sinos, acq_models, trans, image) ########################################################################### # Reconstruct From de37b0c6117427101ce07ef5ef54c9ea721a8923 Mon Sep 17 00:00:00 2001 From: Richard Date: Thu, 18 Jun 2020 12:12:08 +0000 Subject: [PATCH 08/14] codacy changes 5 --- examples/Python/PET/osem_gated_w_motion.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/examples/Python/PET/osem_gated_w_motion.py b/examples/Python/PET/osem_gated_w_motion.py index 317a2103f..6205c1afc 100644 --- a/examples/Python/PET/osem_gated_w_motion.py +++ b/examples/Python/PET/osem_gated_w_motion.py @@ -1,4 +1,5 @@ -"""OSEM reconstruction demo for gated data with motion. +""" +OSEM reconstruction demo for gated data with motion. We actually use the OSMAPOSL reconstructor in this demo. This reconstructor implements an Ordered Subsets (OS) version of the One Step Late algorithm (OSL) @@ -137,7 +138,6 @@ def get_asm_attn(sino, attn, acq_model): def parse_input_files(): """Parse input files.""" - if trans_pattern is None: raise AssertionError("--trans missing") if sino_pattern is None: @@ -209,7 +209,8 @@ def get_initial_estimate(sinos): def get_attn_images(num_ms, attns, trans, image): - """Get attn images. + """ + Get attn images. If none supplied, return None. If 1 supplied, resample to each motion state. From 03c1550dd7821eb5e2f823d93b09e3afa3ac22dd Mon Sep 17 00:00:00 2001 From: Richard Date: Thu, 18 Jun 2020 12:32:04 +0000 Subject: [PATCH 09/14] bug fix --- src/xSTIR/cSTIR/cstir.cpp | 10 ++++++++-- src/xSTIR/cSTIR/include/sirf/STIR/cstir.h | 2 -- 2 files changed, 8 insertions(+), 4 deletions(-) diff --git a/src/xSTIR/cSTIR/cstir.cpp b/src/xSTIR/cSTIR/cstir.cpp index 7de9ebdb3..f9849dc7a 100644 --- a/src/xSTIR/cSTIR/cstir.cpp +++ b/src/xSTIR/cSTIR/cstir.cpp @@ -890,7 +890,6 @@ cSTIR_objectiveFunctionGradientNotDivided(void* ptr_f, void* ptr_i, int subset) CATCH; } -#ifdef STIR_GATED_MOTION extern "C" void* cSTIR_PoissonGatedWMotion_add_gate( @@ -898,12 +897,16 @@ cSTIR_PoissonGatedWMotion_add_gate( const char* disp_fname, const int b_spline_order) { try { +#ifdef STIR_GATED_MOTION SPTR_FROM_HANDLE(PoissonLogLhLinModMeanGatedProjDataWMotion3DF, poisson_sptr, poisson_ptr); SPTR_FROM_HANDLE(PETAcquisitionData, ad_sptr, ad_ptr); SPTR_FROM_HANDLE(PETAcquisitionModel, am_sptr, am_ptr); // set it poisson_sptr->add_gate(ad_sptr, am_sptr, disp_fname, b_spline_order); +#else + throw std::runtime_error("cSTIR_PoissonGatedWMotion_add_gate: shouldn't be here"); +#endif } CATCH; } @@ -913,13 +916,16 @@ void* cSTIR_PoissonGatedWMotion_clear_gates(void* poisson_ptr) { try { +#ifdef STIR_GATED_MOTION SPTR_FROM_HANDLE(PoissonLogLhLinModMeanGatedProjDataWMotion3DF, poisson_sptr, poisson_ptr); poisson_sptr->clear_gates(); +#else + throw std::runtime_error("cSTIR_PoissonGatedWMotion_clear_gates: shouldn't be here"); +#endif } CATCH; } -#endif extern "C" void* diff --git a/src/xSTIR/cSTIR/include/sirf/STIR/cstir.h b/src/xSTIR/cSTIR/include/sirf/STIR/cstir.h index a8ef8fdce..069903d95 100644 --- a/src/xSTIR/cSTIR/include/sirf/STIR/cstir.h +++ b/src/xSTIR/cSTIR/include/sirf/STIR/cstir.h @@ -119,13 +119,11 @@ extern "C" { void* cSTIR_objectiveFunctionGradientNotDivided (void* ptr_f, void* ptr_i, int subset); -#ifdef STIR_GATED_MOTION // Poisson gated proj data w motion void* cSTIR_PoissonGatedWMotion_add_gate( void* poisson_ptr, void* ad_ptr, void* am_ptr, const char* disp_fname, const int b_spline_order); void* cSTIR_PoissonGatedWMotion_clear_gates(void* poisson_ptr); -#endif // Prior methods void* cSTIR_setupPrior(void* ptr_p, void* ptr_i); From df7a0ceae5674d22138e9bab60fcc284cb2b10d8 Mon Sep 17 00:00:00 2001 From: Richard Date: Thu, 18 Jun 2020 12:35:37 +0000 Subject: [PATCH 10/14] codacy changes 6 --- examples/Python/PET/osem_gated_w_motion.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/examples/Python/PET/osem_gated_w_motion.py b/examples/Python/PET/osem_gated_w_motion.py index 6205c1afc..4a0d15f97 100644 --- a/examples/Python/PET/osem_gated_w_motion.py +++ b/examples/Python/PET/osem_gated_w_motion.py @@ -1,5 +1,4 @@ -""" -OSEM reconstruction demo for gated data with motion. +"""OSEM reconstruction demo for gated data with motion. We actually use the OSMAPOSL reconstructor in this demo. This reconstructor implements an Ordered Subsets (OS) version of the One Step Late algorithm (OSL) From 5c88ca54f7878169dfff07712dc81fd16215bd6d Mon Sep 17 00:00:00 2001 From: Richard Date: Thu, 18 Jun 2020 13:06:55 +0000 Subject: [PATCH 11/14] bug fix --- src/xSTIR/cSTIR/stir_x.cpp | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/src/xSTIR/cSTIR/stir_x.cpp b/src/xSTIR/cSTIR/stir_x.cpp index 5d98e4044..f512c1194 100644 --- a/src/xSTIR/cSTIR/stir_x.cpp +++ b/src/xSTIR/cSTIR/stir_x.cpp @@ -697,13 +697,14 @@ add_gate( static stir::shared_ptr vec_AcqData_to_GatedProjData(const std::vector > &vec_ad) { + if (vec_ad.empty() || vec_ad.at(0) == nullptr) + return nullptr; + auto sino_sptr = MAKE_SHARED(); - if (!vec_ad.empty()) { - sino_sptr->set_exam_info(*vec_ad.at(0)->get_exam_info_sptr()); - sino_sptr->resize(vec_ad.size()); - for (unsigned i=0; iset_proj_data_sptr(vec_ad.at(i)->data(), i+1); - } + sino_sptr->set_exam_info(*vec_ad.at(0)->get_exam_info_sptr()); + sino_sptr->resize(vec_ad.size()); + for (unsigned i=0; iset_proj_data_sptr(vec_ad.at(i)->data(), i+1); return sino_sptr; } From 77f913e8ef901c2751e601cb1a288f6ffeda301a Mon Sep 17 00:00:00 2001 From: Richard Date: Thu, 18 Jun 2020 14:25:37 +0000 Subject: [PATCH 12/14] bug fix 3 --- src/xSTIR/cSTIR/cstir.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/xSTIR/cSTIR/cstir.cpp b/src/xSTIR/cSTIR/cstir.cpp index f9849dc7a..2c29f1ce2 100644 --- a/src/xSTIR/cSTIR/cstir.cpp +++ b/src/xSTIR/cSTIR/cstir.cpp @@ -907,6 +907,7 @@ cSTIR_PoissonGatedWMotion_add_gate( #else throw std::runtime_error("cSTIR_PoissonGatedWMotion_add_gate: shouldn't be here"); #endif + return new DataHandle; } CATCH; } @@ -923,6 +924,7 @@ cSTIR_PoissonGatedWMotion_clear_gates(void* poisson_ptr) #else throw std::runtime_error("cSTIR_PoissonGatedWMotion_clear_gates: shouldn't be here"); #endif + return new DataHandle; } CATCH; } From 9375082f4c570e2dc4b7a0c1ae38e65c2fbc3b40 Mon Sep 17 00:00:00 2001 From: Richard Date: Thu, 18 Jun 2020 14:53:20 +0000 Subject: [PATCH 13/14] codacy changes 7 --- examples/Python/PET/osem_gated_w_motion.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/examples/Python/PET/osem_gated_w_motion.py b/examples/Python/PET/osem_gated_w_motion.py index 4a0d15f97..6205c1afc 100644 --- a/examples/Python/PET/osem_gated_w_motion.py +++ b/examples/Python/PET/osem_gated_w_motion.py @@ -1,4 +1,5 @@ -"""OSEM reconstruction demo for gated data with motion. +""" +OSEM reconstruction demo for gated data with motion. We actually use the OSMAPOSL reconstructor in this demo. This reconstructor implements an Ordered Subsets (OS) version of the One Step Late algorithm (OSL) From 2897c033725ff4ed018ba779156ed5149a4cd82b Mon Sep 17 00:00:00 2001 From: Richard Date: Fri, 26 Jun 2020 08:32:57 +0000 Subject: [PATCH 14/14] New -> UsingAdjoint --- src/xSTIR/cSTIR/include/sirf/STIR/stir_types.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/xSTIR/cSTIR/include/sirf/STIR/stir_types.h b/src/xSTIR/cSTIR/include/sirf/STIR/stir_types.h index 489317086..5036a4180 100644 --- a/src/xSTIR/cSTIR/include/sirf/STIR/stir_types.h +++ b/src/xSTIR/cSTIR/include/sirf/STIR/stir_types.h @@ -68,7 +68,7 @@ limitations under the License. #include "stir/recon_buildblock/NiftyPET_projector/ProjectorByBinPairUsingNiftyPET.h" #endif #ifdef STIR_GATED_MOTION -#include "stir_experimental/motion/PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion.h" +#include "stir_experimental/motion/PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotionUsingAdjoint.h" #include "stir_experimental/motion/NonRigidObjectTransformationUsingBSplines.h" #include "stir_experimental/motion/Transform3DObjectImageProcessor.h" #endif @@ -106,7 +106,7 @@ namespace sirf { #ifdef STIR_GATED_MOTION typedef stir::ObjectTransformation<3,float> STIRTrans3DF; typedef stir::NonRigidObjectTransformationUsingBSplines<3,float> STIRDisp3DF; - typedef stir::PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotionNew < Image3DF > + typedef stir::PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotionUsingAdjoint < Image3DF > stirPoissonLogLhLinModMeanGatedWMotion3DF; #endif