add already_set_up to GeneralisedObjectiveFunction

documentation/release_6.0.htm

@@ -252,6 +252,11 @@ <h3>Backward incompatibities</h3>
         <code>ProjDataInfo*NoArcCorr::get_bin_for_det_pair</code> is now private.
         Use <code>get_bin_for_det_pos_pair</code> instead.
       </li>
+      <li>
+        The <code>GeneralisedObjectiveFunction</code> hierarchy now has a <code>already_set_up</code>
+        member variable that needs to be set to <code>false</code> by <code>set_*</code>
+        functions and checked by callers.
+      </li>
     </ul>
 
     <h3>New functionality</h3>

src/buildblock/ProjDataInfoGenericNoArcCorr.cxx

@@ -74,18 +74,8 @@ ProjDataInfoGenericNoArcCorr(const shared_ptr<Scanner> scanner_sptr,
   this->initialise_det1det2_to_uncompressed_view_tangpos();
 #endif
 
-  CartesianCoordinate3D< float> b1,b2;
-  Bin bin;
-  bin.segment_num() = 0;
-  bin.axial_pos_num() = 0;
-  bin.view_num() = 0;
-  bin.tangential_pos_num() = 0;
-// setting shift_z to 0 before it is actually estimated. Otherwise the next function will use it
-  this->z_shift.z()=0;
-  find_cartesian_coordinates_of_detection(b1,b2,bin);
-  float shift=b2.z();
-
-  this->z_shift.z()=shift;
+  // find shift between "new" centre-of-scanner and "old" centre-of-first-ring coordinate system
+  this->z_shift.z()=this->get_scanner_ptr()->get_coordinate_for_det_pos(DetectionPosition<>(0,0,0)).z();
   this->z_shift.y()=0;
   this->z_shift.x()=0;
 }

src/experimental/motion/PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion.cxx

@@ -285,6 +285,7 @@ void
 PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<TargetT>::
 set_proj_data_sptr(const shared_ptr<GatedProjData>& arg)
 {
+  this->already_set_up = false;
   this->_gated_proj_data_sptr = arg;
 }
 
@@ -293,15 +294,16 @@ void
 PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<TargetT>::
 set_max_segment_num_to_process(const int arg)
 {
+  this->already_set_up = this->already_set_up && (this->max_timing_pos_num_to_process == arg);
   this->max_segment_num_to_process = arg;
-
 }
 
 template<typename TargetT>
 void
 PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<TargetT>::
 set_zero_seg0_end_planes(const bool arg)
 {
+  this->already_set_up = this->already_set_up && (this->zero_seg0_end_planes == arg);
   this->zero_seg0_end_planes = arg;
 }
 
@@ -310,7 +312,7 @@ void
 PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<TargetT>::
 set_additive_proj_data_sptr(const shared_ptr<GatedProjData>& arg)
 {
-
+  this->already_set_up = false;
   this->_gated_additive_proj_data_sptr = arg;
 }
 
@@ -319,6 +321,7 @@ void
 PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<TargetT>::
 set_projector_pair_sptr(const shared_ptr<ProjectorByBinPair>& arg) 
 {
+  this->already_set_up = false;  
   this->projector_pair_ptr = arg;
 
 }
@@ -328,6 +331,7 @@ void
 PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<TargetT>::
 set_frame_num(const int arg)
 {
+  this->already_set_up = this->already_set_up && (this->frame_num == arg);
   this->frame_num = arg;
 }
 
@@ -336,6 +340,7 @@ void
 PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<TargetT>::
 set_frame_definitions(const TimeFrameDefinitions& arg)
 {
+  this->already_set_up = this->already_set_up && (this->frame_defs == arg);
   this->frame_defs = arg;
 }
 

src/include/stir/recon_buildblock/GeneralisedObjectiveFunction.h

@@ -86,12 +86,19 @@ class GeneralisedObjectiveFunction:
 {
 public:
 
-  //GeneralisedObjectiveFunction(); 
+  GeneralisedObjectiveFunction()
+    : already_set_up(false)
+  {}
+
 
   virtual ~GeneralisedObjectiveFunction(); 
 
 
   //! Creates a suitable target as determined by the parameters
+  /*!
+    \warning This should <b>not</b> check \c already_set_up (unfortunately),
+    as it is currently called in Reconstruction::reconstruct() before calling set_up().
+  */
   virtual TargetT *
     construct_target_ptr() const = 0; 
 
@@ -331,6 +338,7 @@ class GeneralisedObjectiveFunction:
 
 protected:
   int num_subsets;
+  bool already_set_up;
 
   shared_ptr<GeneralisedPrior<TargetT> > prior_sptr;
 

src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion.txx

@@ -175,7 +175,7 @@ template <typename TargetT>
 TargetT *
 PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<TargetT>::
 construct_target_ptr() const
-{  
+{
   return 
     this->target_parameter_parser.create(this->get_input_data());
 }
@@ -280,29 +280,34 @@ get_projector_pair_sptr() const
 template<typename TargetT>
 int
 PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<TargetT>::
-set_num_subsets(const int num_subsets)
+set_num_subsets(const int new_num_subsets)
 {
+  this->already_set_up = this->already_set_up && (this->num_subsets == new_num_subsets);
   for(unsigned int gate_num=1;gate_num<=this->get_time_gate_definitions().get_num_gates();++gate_num)
     {
       if(this->_single_gate_obj_funcs.size() != 0)
-		if(this->_single_gate_obj_funcs[gate_num].set_num_subsets(num_subsets) != num_subsets)
+		if(this->_single_gate_obj_funcs[gate_num].set_num_subsets(new_num_subsets) != new_num_subsets)
 			error("set_num_subsets didn't work");
     }
-  this->num_subsets=num_subsets;
+  this->num_subsets=new_num_subsets;
   return this->num_subsets;
 }
 
 template<typename TargetT>
 void
 PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<TargetT>::
 set_time_gate_definitions(const TimeGateDefinitions & time_gate_definitions)
-{ this->_time_gate_definitions=time_gate_definitions; }
+{
+  this->already_set_up = false;
+  this->_time_gate_definitions=time_gate_definitions;
+}
 
 template<typename TargetT>
 void
 PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<TargetT>::
 set_input_data(const shared_ptr<ExamData> & arg)
 {
+  this->already_set_up = false;
     this->_gated_proj_data_sptr = dynamic_pointer_cast<GatedProjData>(arg);
 }
 
@@ -319,6 +324,7 @@ void
 PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<TargetT>::
 set_additive_proj_data_sptr(const shared_ptr<ExamData> &arg)
 {
+  this->already_set_up = false;
     this->_additive_gated_proj_data_sptr = dynamic_pointer_cast<GatedProjData>(arg);
 }
 
@@ -327,6 +333,7 @@ void
 PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<TargetT>::
 set_normalisation_sptr(const shared_ptr<BinNormalisation>& arg)
 {
+  this->already_set_up = false;
 //  this->normalisation_sptr = arg;
     error("Not implemeted yet");
 }

src/recon_buildblock/GeneralisedObjectiveFunction.cxx

@@ -142,6 +142,8 @@ compute_sub_gradient(TargetT& gradient,
 		     const TargetT &current_estimate, 
 		     const int subset_num)
 {
+  if (!this->already_set_up)
+    error("Need to call set_up() for objective function first");
   assert(gradient.get_index_range() == current_estimate.get_index_range());
 
   if (subset_num<0 || subset_num>=this->get_num_subsets())
@@ -194,6 +196,8 @@ GeneralisedObjectiveFunction<TargetT>::
 compute_objective_function_without_penalty(const TargetT& current_estimate,
 					   const int subset_num)
 {
+  if (!this->already_set_up)
+    error("Need to call set_up() for objective function first");
   if (subset_num<0 || subset_num>=this->get_num_subsets())
     error("compute_objective_function_without_penalty subset_num out-of-range error");
 
@@ -231,6 +235,8 @@ add_multiplication_with_approximate_sub_Hessian_without_penalty(TargetT& output,
 								const TargetT& input,
 								const int subset_num) const
 {
+  if (!this->already_set_up)
+    error("Need to call set_up() for objective function first");
   if (subset_num<0 || subset_num>=this->get_num_subsets())
     error("add_multiplication_with_approximate_sub_Hessian_without_penalty subset_num out-of-range error");
 
@@ -259,6 +265,8 @@ add_multiplication_with_approximate_sub_Hessian(TargetT& output,
 						const TargetT& input,
 						const int subset_num) const
 {
+  if (!this->already_set_up)
+    error("Need to call set_up() for objective function first");
   if (this->add_multiplication_with_approximate_sub_Hessian_without_penalty(output, input, subset_num) ==
       Succeeded::no)
     return Succeeded::no;
@@ -406,6 +414,8 @@ accumulate_sub_Hessian_times_input_without_penalty(TargetT& output,
                                                    const TargetT& input,
                                                    const int subset_num) const
 {
+  if (!this->already_set_up)
+    error("Need to call set_up() for objective function first");
   if (subset_num<0 || subset_num>=this->get_num_subsets())
     error("accumulate_sub_Hessian_times_input_without_penalty subset_num out-of-range error");
 
@@ -487,6 +497,12 @@ bool
 GeneralisedObjectiveFunction<TargetT>::
 subsets_are_approximately_balanced(std::string& warning_message) const
 {
+#if 0
+  // TODO cannot do this yet, as this function is called in
+  // `PoissonLogLikelihoodWithLinearModelForMean::compute_sensitivities` during set_up()
+  if (!this->already_set_up)
+    error("Need to call set_up() for objective function first");
+#endif
   return this->actual_subsets_are_approximately_balanced(warning_message);
 }
 

src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMean.cxx

@@ -95,6 +95,7 @@ void
 PoissonLogLikelihoodWithLinearModelForMean<TargetT>::
 set_sensitivity_filename(const std::string& filename)
 {
+  this->already_set_up = false;
   this->sensitivity_filename = filename;
 }
 
@@ -103,6 +104,7 @@ void
 PoissonLogLikelihoodWithLinearModelForMean<TargetT>::
 set_subsensitivity_filenames(const std::string& filenames)
 {
+  this->already_set_up = false;
   this->subsensitivity_filenames = filenames;
   try
     {
@@ -171,6 +173,7 @@ void
 PoissonLogLikelihoodWithLinearModelForMean<TargetT>::
 set_use_subset_sensitivities(const bool arg)
 {
+  this->already_set_up = this->already_set_up && (this->use_subset_sensitivities == arg);
   this->use_subset_sensitivities = arg;
 }
 
@@ -179,6 +182,7 @@ void
 PoissonLogLikelihoodWithLinearModelForMean<TargetT>::
 set_subset_sensitivity_sptr(const shared_ptr<TargetT>& arg, const int subset_num)
 {
+  this->already_set_up = false;
   this->subsensitivity_sptrs[subset_num] = arg;
 }
 
@@ -348,7 +352,7 @@ set_up(shared_ptr<TargetT> const& target_sptr)
           return Succeeded::no;
         }
     }
-      
+  this->already_set_up = true;
   return Succeeded::yes;
 }
 
@@ -359,6 +363,8 @@ compute_sub_gradient_without_penalty(TargetT& gradient,
                                      const TargetT &current_estimate, 
                                      const int subset_num)
 {
+  if (!this->already_set_up)
+    error("Need to call set_up() for objective function first");
   this->actual_compute_subset_gradient_without_penalty(gradient, current_estimate, subset_num, false);
 }
 
@@ -369,6 +375,8 @@ compute_sub_gradient_without_penalty_plus_sensitivity(TargetT& gradient,
                                                       const TargetT &current_estimate,
                                                       const int subset_num)
 {
+  if (!this->already_set_up)
+    error("Need to call set_up() for objective function first");
   actual_compute_subset_gradient_without_penalty(gradient, current_estimate, subset_num, true);
 }
 

src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndListModeData.cxx

@@ -100,10 +100,10 @@ PoissonLogLikelihoodWithLinearModelForMeanAndListModeData<TargetT>::post_process
   if (this->list_mode_filename.length() == 0 && !this->skip_lm_input_file)
   { warning("You need to specify an input file\n"); return true; }
   if (!this->skip_lm_input_file)
-  {
+    {
       this->list_mode_data_sptr=
           read_from_file<ListModeData>(this->list_mode_filename);
-  }
+    }
 
   if (this->additive_projection_data_filename != "0")
     {
@@ -123,6 +123,7 @@ PoissonLogLikelihoodWithLinearModelForMeanAndListModeData<TargetT>::post_process
     } 
   target_parameter_parser.check_values();
 
+  this->already_set_up = false;
   return false;
 } 
 
@@ -131,6 +132,7 @@ void
 PoissonLogLikelihoodWithLinearModelForMeanAndListModeData<TargetT>::
 set_input_data(const shared_ptr<ExamData> & arg)
 {
+  this->already_set_up = false;
   try
     {
       this->list_mode_data_sptr = dynamic_pointer_cast<ListModeData>(arg);
@@ -146,6 +148,7 @@ void
 PoissonLogLikelihoodWithLinearModelForMeanAndListModeData<TargetT>::
 set_max_segment_num_to_process(const int arg)
 {
+  this->already_set_up = this->already_set_up && (this->max_segment_num_to_process == arg);
   this->max_segment_num_to_process = arg;
 }
 
@@ -258,6 +261,7 @@ void
 PoissonLogLikelihoodWithLinearModelForMeanAndListModeData<TargetT>::
 set_additive_proj_data_sptr(const shared_ptr<ExamData> &arg)
 {
+  this->already_set_up = false;
   try
     {
       this->additive_proj_data_sptr = dynamic_pointer_cast<ProjData>(arg);
@@ -281,6 +285,7 @@ void
 PoissonLogLikelihoodWithLinearModelForMeanAndListModeData<TargetT>::
 set_normalisation_sptr(const shared_ptr<BinNormalisation>& arg)
 {
+  this->already_set_up = false;
   this->normalisation_sptr = arg;
 }
 

src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin.cxx

@@ -116,6 +116,7 @@ int
 PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin<TargetT>::
 set_num_subsets(const int new_num_subsets)
 {
+  this->already_set_up = this->already_set_up && (this->num_subsets == new_num_subsets);
   this->num_subsets = new_num_subsets;
   return this->num_subsets;
 }
@@ -125,7 +126,8 @@ void
 PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin<TargetT>::
 set_proj_matrix(const shared_ptr<ProjMatrixByBin>& arg)
 {
-    this->PM_sptr = arg;
+  this->already_set_up = false;
+  this->PM_sptr = arg;
 }
 
 #if 0
@@ -134,6 +136,7 @@ void
 PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin<TargetT>::
 set_proj_data_info(const ProjData& arg)
 {
+  this->already_set_up = false;
   // this will be broken now. Why did we need this?
     this->proj_data_info_sptr = arg.get_proj_data_info_sptr()->create_shared_clone();
     if(this->skip_lm_input_file)
@@ -724,8 +727,7 @@ TargetT *
 PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin<TargetT>:: 
 construct_target_ptr() const 
 { 
-
- return 
+  return 
    this->target_parameter_parser.create(this->get_input_data());
 }