Proposed Changes to MCIR Sensitivity

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 :=

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 :=

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

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<float>& viewgrams) const override { viewgrams /= this->_factor; }
+
+  void undo(RelatedViewgrams<float>& viewgrams) const override { viewgrams *= this->_factor; }
+
+private:
+  float _factor;
+};
+
+END_NAMESPACE_STIR
+
+#endif

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 <typename TargetT>
@@ -378,7 +381,23 @@ PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<TargetT>::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<TargetT>::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<std::pair<double, double>> frame_times(1, std::pair<double, double>(0, 1));
         this->_single_gate_obj_funcs[gate_num].set_frame_definitions(TimeFrameDefinitions(frame_times));
 
@@ -402,6 +424,31 @@ PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<TargetT>::s
             shared_ptr<ProjData> 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<float>(this->get_time_gate_definitions().get_gate_duration(gate_num) / total_gate_duration);
+        shared_ptr<BinNormalisation> 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<ProjData> 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<ProjDataInMemory> 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<float> vg = add_proj_data_sptr->get_viewgram(view_num, segment_num, /*make_num_tangential_poss_odd=*/false, timing_pos);
+                        vg /= static_cast<float>(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());