diff --git a/recon_test_pack/OSMAPOSL_with_motion_test_f1.par b/recon_test_pack/OSMAPOSL_with_motion_test_f1.par new file mode 100644 index 0000000000..54adcb487f --- /dev/null +++ b/recon_test_pack/OSMAPOSL_with_motion_test_f1.par @@ -0,0 +1,40 @@ +OSMAPOSLParameters := + +objective function type:= PoissonLogLikelihoodWithLinearModelForMeanAndProjData +PoissonLogLikelihoodWithLinearModelForMeanAndProjData Parameters:= + +input file := my_att_fwd_test_object_g1.hs +; if disabled, defaults to maximum segment number in the file +maximum absolute segment number to process := -1 +; see User's Guide to see when you need this +zero end planes of segment 0 := 0 + +projector pair type := Matrix +Projector Pair Using Matrix Parameters := +Matrix type := Ray Tracing +Ray Tracing Matrix Parameters:= + restrict to cylindrical fov := 1 + ; use a slightly better approximation than simple ray tracing + number of rays in tangential direction to trace for each bin := 10 + ; remove one symmetry such that we can use more subsets + do symmetry 90degrees min phi := 0 +End Ray Tracing Matrix Parameters:= +End Projector Pair Using Matrix Parameters := + +bin normalisation type := From ProjData + bin normalisation from projdata := + normalisation projdata filename := my_ACF_test_object_g1.hs +end bin normalisation from projdata := + +zoom:= 1 +xy output image size (in pixels) := 101 + +end PoissonLogLikelihoodWithLinearModelForMeanAndProjData Parameters:= + +number of subsets:= 28 +number of subiterations:= 46 +start at subiteration number:=28 + +output filename prefix := my_rec_test_object_g1 + +END := diff --git a/recon_test_pack/OSMAPOSL_with_motion_test_f2.par b/recon_test_pack/OSMAPOSL_with_motion_test_f2.par new file mode 100644 index 0000000000..3d67307676 --- /dev/null +++ b/recon_test_pack/OSMAPOSL_with_motion_test_f2.par @@ -0,0 +1,40 @@ +OSMAPOSLParameters := + +objective function type:= PoissonLogLikelihoodWithLinearModelForMeanAndProjData +PoissonLogLikelihoodWithLinearModelForMeanAndProjData Parameters:= + +input file := my_att_fwd_test_object_g2.hs +; if disabled, defaults to maximum segment number in the file +maximum absolute segment number to process := -1 +; see User's Guide to see when you need this +zero end planes of segment 0 := 0 + +projector pair type := Matrix +Projector Pair Using Matrix Parameters := +Matrix type := Ray Tracing +Ray Tracing Matrix Parameters:= + restrict to cylindrical fov := 1 + ; use a slightly better approximation than simple ray tracing + number of rays in tangential direction to trace for each bin := 10 + ; remove one symmetry such that we can use more subsets + do symmetry 90degrees min phi := 0 +End Ray Tracing Matrix Parameters:= +End Projector Pair Using Matrix Parameters := + +bin normalisation type := From ProjData + bin normalisation from projdata := + normalisation projdata filename := my_ACF_test_object_g2.hs +end bin normalisation from projdata := + +zoom:= 1 +xy output image size (in pixels) := 101 + +end PoissonLogLikelihoodWithLinearModelForMeanAndProjData Parameters:= + +number of subsets:= 28 +number of subiterations:= 46 +start at subiteration number:=28 + +output filename prefix := my_rec_test_object_g2 + +END := diff --git a/recon_test_pack/run_test_simulate_and_recon_with_motion.sh b/recon_test_pack/run_test_simulate_and_recon_with_motion.sh index 2dec1319a7..bc97e3df05 100755 --- a/recon_test_pack/run_test_simulate_and_recon_with_motion.sh +++ b/recon_test_pack/run_test_simulate_and_recon_with_motion.sh @@ -145,6 +145,27 @@ fi error_log_files="" cp my_template.gdef my_att_fwd_test_object.gdef cp my_template.gdef my_ACF_test_object.gdef + # adjusted to run to compare with RTA instead of ground truth because of Gibbs ringing artefacts around + # the edge of object. + echo "Running OSMAPOSL OSMAPOSL_with_motion_test_f1.par for 1 iteration (28 subsets)" + ${MPIRUN} OSMAPOSL OSMAPOSL_with_motion_test_f1.par > my_rec_test_object_f1.log 2>&1 + if [ $? -ne 0 ]; then + echo "Error running reconstruction. CHECK RECONSTRUCTION LOG my_rec_test_object_f1.log" + error_log_files="${error_log_files} my_rec_test_object_f1.log" + exit 1 + fi + echo "Running OSMAPOSL OSMAPOSL_with_motion_test_f2.par for 1 iteration (28 subsets)" + ${MPIRUN} OSMAPOSL OSMAPOSL_with_motion_test_f2.par > my_rec_test_object_f2.log 2>&1 + if [ $? -ne 0 ]; then + echo "Error running reconstruction. CHECK RECONSTRUCTION LOG my_rec_test_object_f2.log" + error_log_files="${error_log_files} my_rec_test_object_f2.log" + exit 1 + fi + + mv my_rec_test_object_g1_46.hv my_rec_test_object_g1.hv + mv my_rec_test_object_g2_46.hv my_rec_test_object_g2.hv + cp my_template.gdef my_rec_test_object.gdef + warp_and_accumulate_gated_images my_twice_test_object_RTA my_rec_test_object my_reverse_translation echo "Running OSMAPOSL OSMAPOSL_with_motion_test.par for 1 iteration (28 subsets)" ${MPIRUN} OSMAPOSL OSMAPOSL_with_motion_test.par > my_rec_test_object.log 2>&1 @@ -153,7 +174,8 @@ cp my_template.gdef my_ACF_test_object.gdef error_log_files="${error_log_files} my_rec_test_object.log" exit 1 fi - compare_image -t 0.25 my_rec_test_object_46.hv my_test_object_g1.hv + + compare_image -t 0.25 my_rec_test_object_46.hv my_twice_test_object_RTA.hv if [ $? -ne 0 ]; then echo "ERROR comparison of reconstructed image fails"; exit 1; fi @@ -167,13 +189,22 @@ fi fi # compare ROI value output_voxel_size_x=`stir_print_voxel_sizes.sh ${output_image}|awk '{print $3}'` - output_ROI_mean=`awk "NR>2 {print \\$2*${input_voxel_size_x}/${output_voxel_size_x}}" ${output_image}.roistats` - echo "Input ROI mean: $input_ROI_mean" + expected_image=my_twice_test_object_RTA.hv + list_ROI_values ${expected_image}.roistats ${expected_image} ${ROI} 0 > ${expected_image}.roistats.log 2>&1 + if [ $? -ne 0 ]; then + echo "Error running list_ROI_values for reference image. CHECK LOG ${expected_image}.roistats.log" + error_log_files="${error_log_files} ${expected_image}.roistats.log" + break + fi + expected_voxel_size_x=`stir_print_voxel_sizes.sh ${expected_image}|awk '{print $3}'` + expected_ROI_mean=`awk "NR>2 {print \\$2*${expected_voxel_size_x}/${input_voxel_size_x}}" ${expected_image}.roistats` + output_ROI_mean=`awk "NR>2 {print \\$2*${output_voxel_size_x}/${input_voxel_size_x}}" ${output_image}.roistats` + echo "Expected ROI mean: $expected_ROI_mean" echo "Output ROI mean: $output_ROI_mean" - error_bigger_than_1percent=`echo $input_ROI_mean $output_ROI_mean| awk '{ print(($2/$1 - 1)*($2/$1 - 1)>0.0001) }'` + error_bigger_than_1percent=`echo $expected_ROI_mean $output_ROI_mean| awk '{ print(($2/$1 - 1)*($2/$1 - 1)>0.0001) }'` if [ ${error_bigger_than_1percent} -eq 1 ]; then - echo "DIFFERENCE IN ROI VALUES IS TOO LARGE. CHECK RECONSTRUCTION LOG ${parfile}.log" - error_log_files="${error_log_files} ${parfile}.log" + echo "DIFFERENCE IN ROI VALUES IS TOO LARGE. CHECK RECONSTRUCTION LOG ${output_image}.roistats.log" + error_log_files="${error_log_files} ${output_image}.roistats.log" else echo "This seems fine." fi diff --git a/src/include/stir/recon_buildblock/BinNormalisationFromConstantFactor.h b/src/include/stir/recon_buildblock/BinNormalisationFromConstantFactor.h new file mode 100644 index 0000000000..b261d9e321 --- /dev/null +++ b/src/include/stir/recon_buildblock/BinNormalisationFromConstantFactor.h @@ -0,0 +1,56 @@ +// +// +/*! + \file + \ingroup normalisation + + \brief Declaration of class stir::BinNormalisationFromConstantfactor + + \author Viet Dao +*/ +/* + Copyright (C) 2000- 2011, Hammersmith Imanet Ltd + Copyright (C) 2023, University College London + This file is part of STIR. + + SPDX-License-Identifier: Apache-2.0 + + See STIR/LICENSE.txt for details +*/ +#ifndef __stir_recon_buildblock_BinNormalisationFromConstantFactor_H__ +#define __stir_recon_buildblock_BinNormalisationFromConstantFactor_H__ + +#include "stir/recon_buildblock/BinNormalisation.h" +#include "stir/RelatedViewgrams.h" +START_NAMESPACE_STIR + +/*! + \ingroup normalisation + \brief A BinNormalisation class that gets the normalisation factors from + a constant factor such as time scale or decay factors. +*/ + +class BinNormalisationFromConstantFactor : public BinNormalisation +{ +public: + explicit BinNormalisationFromConstantFactor(const float factor) + : _factor(factor) + {} + + bool is_trivial() const override { return this->_factor == 1.F; } + + std::string get_registered_name() const override { return "Constant factor"; } + + float get_bin_efficiency(const Bin&) const override { return this->_factor; } + + void apply(RelatedViewgrams& viewgrams) const override { viewgrams /= this->_factor; } + + void undo(RelatedViewgrams& viewgrams) const override { viewgrams *= this->_factor; } + +private: + float _factor; +}; + +END_NAMESPACE_STIR + +#endif \ No newline at end of file diff --git a/src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion.txx b/src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion.txx index 7ab6fdb123..12375c9d93 100644 --- a/src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion.txx +++ b/src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion.txx @@ -46,7 +46,10 @@ #include "stir/spatial_transformation/GatedSpatialTransformation.h" #include "stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion.h" #include "stir/recon_buildblock/BinNormalisationFromProjData.h" - +#include "stir/recon_buildblock/BinNormalisationFromConstantFactor.h" +#include "stir/recon_buildblock/ChainedBinNormalisation.h" +#include "stir/ProjDataInMemory.h" +#include "stir/Viewgram.h" START_NAMESPACE_STIR template @@ -378,7 +381,23 @@ PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion::s // construct _single_gate_obj_funcs this->_single_gate_obj_funcs.resize(1, this->get_time_gate_definitions().get_num_gates()); - + + double total_gate_duration = 0.; + for (unsigned int gate_num = 1; gate_num <= this->get_time_gate_definitions().get_num_gates(); ++gate_num) + { + const double gate_duration = this->get_time_gate_definitions().get_gate_duration(gate_num); + if (gate_duration <= 0.) + { + warning("Gate duration for gate %u should be larger than 0.", gate_num); + return Succeeded::no; + } + total_gate_duration += gate_duration; + } + if (total_gate_duration <= 0.) + { + warning("Total gate duration should be larger than 0."); + return Succeeded::no; + } for (unsigned int gate_num = 1; gate_num <= this->get_time_gate_definitions().get_num_gates(); ++gate_num) { info(format("Objective Function for Gate Number: {}", gate_num)); @@ -391,6 +410,9 @@ PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion::s this->_additive_gated_proj_data_sptr->get_proj_data_sptr(gate_num)); this->_single_gate_obj_funcs[gate_num].set_num_subsets(this->num_subsets); this->_single_gate_obj_funcs[gate_num].set_frame_num(1); // This should be gate... + const double gate_duration = this->get_time_gate_definitions().get_gate_duration(gate_num); + // Keep gate timing available for time-dependent normalisation, but do not use it + // as an extra sensitivity scale factor below. std::vector> frame_times(1, std::pair(0, 1)); this->_single_gate_obj_funcs[gate_num].set_frame_definitions(TimeFrameDefinitions(frame_times)); @@ -402,6 +424,31 @@ PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion::s shared_ptr norm_data_sptr(this->_normalisation_gated_proj_data_sptr->get_proj_data_sptr(gate_num)); current_gate_norm_factors_sptr.reset(new BinNormalisationFromProjData(norm_data_sptr)); } + const float gate_duration_fraction = static_cast(this->get_time_gate_definitions().get_gate_duration(gate_num) / total_gate_duration); + shared_ptr gate_duration_norm_sptr(new BinNormalisationFromConstantFactor(gate_duration_fraction)); + current_gate_norm_factors_sptr.reset(new ChainedBinNormalisation(gate_duration_norm_sptr, current_gate_norm_factors_sptr)); + // Apply the inverse of the gate-duration constant factor to additive projection data for this gate + if (!is_null_ptr(this->_additive_gated_proj_data_sptr) && gate_duration_fraction != 1.F && gate_duration_fraction > 0.F) + { + shared_ptr add_proj_data_sptr(this->_additive_gated_proj_data_sptr->get_proj_data_sptr(gate_num)); + if (!is_null_ptr(add_proj_data_sptr)) + { + shared_ptr scaled_add_proj_sptr(new ProjDataInMemory(add_proj_data_sptr->get_exam_info_sptr(), + add_proj_data_sptr->get_proj_data_info_sptr(), + /*initialise_with_0=*/true)); + for (int segment_num = add_proj_data_sptr->get_min_segment_num(); + segment_num <= add_proj_data_sptr->get_max_segment_num(); ++segment_num) + for (int view_num = add_proj_data_sptr->get_min_view_num(); view_num <= add_proj_data_sptr->get_max_view_num(); ++view_num) + for (int timing_pos = add_proj_data_sptr->get_min_tof_pos_num(); + timing_pos <= add_proj_data_sptr->get_max_tof_pos_num(); ++timing_pos) + { + Viewgram vg = add_proj_data_sptr->get_viewgram(view_num, segment_num, /*make_num_tangential_poss_odd=*/false, timing_pos); + vg /= static_cast(gate_duration_fraction); + scaled_add_proj_sptr->set_viewgram(vg); + } + this->_single_gate_obj_funcs[gate_num].set_additive_proj_data_sptr(scaled_add_proj_sptr); + } + } this->_single_gate_obj_funcs[gate_num].set_normalisation_sptr(current_gate_norm_factors_sptr); this->_single_gate_obj_funcs[gate_num].set_recompute_sensitivity(this->get_recompute_sensitivity()); this->_single_gate_obj_funcs[gate_num].set_use_subset_sensitivities(this->get_use_subset_sensitivities());