diff --git a/documentation/release_6.0.htm b/documentation/release_6.0.htm index b9e7ed1817..cc8be1cfa7 100644 --- a/documentation/release_6.0.htm +++ b/documentation/release_6.0.htm @@ -252,6 +252,11 @@

Backward incompatibities

ProjDataInfo*NoArcCorr::get_bin_for_det_pair is now private. Use get_bin_for_det_pos_pair instead. +
  • + The GeneralisedObjectiveFunction hierarchy now has a already_set_up + member variable that needs to be set to false by set_* + functions and checked by callers. +
  • New functionality

    diff --git a/src/buildblock/ProjDataInfoGenericNoArcCorr.cxx b/src/buildblock/ProjDataInfoGenericNoArcCorr.cxx index e900897d24..ee21c30280 100644 --- a/src/buildblock/ProjDataInfoGenericNoArcCorr.cxx +++ b/src/buildblock/ProjDataInfoGenericNoArcCorr.cxx @@ -74,18 +74,8 @@ ProjDataInfoGenericNoArcCorr(const shared_ptr 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; } diff --git a/src/experimental/motion/PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion.cxx b/src/experimental/motion/PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion.cxx index 6d2b4de604..4453c06ea8 100644 --- a/src/experimental/motion/PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion.cxx +++ b/src/experimental/motion/PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion.cxx @@ -285,6 +285,7 @@ void PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion:: set_proj_data_sptr(const shared_ptr& arg) { + this->already_set_up = false; this->_gated_proj_data_sptr = arg; } @@ -293,8 +294,8 @@ void PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion:: 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 @@ -302,6 +303,7 @@ void PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion:: 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:: set_additive_proj_data_sptr(const shared_ptr& arg) { - + this->already_set_up = false; this->_gated_additive_proj_data_sptr = arg; } @@ -319,6 +321,7 @@ void PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion:: set_projector_pair_sptr(const shared_ptr& arg) { + this->already_set_up = false; this->projector_pair_ptr = arg; } @@ -328,6 +331,7 @@ void PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion:: 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:: set_frame_definitions(const TimeFrameDefinitions& arg) { + this->already_set_up = this->already_set_up && (this->frame_defs == arg); this->frame_defs = arg; } diff --git a/src/include/stir/recon_buildblock/GeneralisedObjectiveFunction.h b/src/include/stir/recon_buildblock/GeneralisedObjectiveFunction.h index 021089257a..d479dd8940 100644 --- a/src/include/stir/recon_buildblock/GeneralisedObjectiveFunction.h +++ b/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 not 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 > prior_sptr; diff --git a/src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion.txx b/src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion.txx index 2983842152..42a4a86e37 100644 --- a/src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion.txx +++ b/src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion.txx @@ -175,7 +175,7 @@ template TargetT * PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion:: construct_target_ptr() const -{ +{ return this->target_parameter_parser.create(this->get_input_data()); } @@ -280,15 +280,16 @@ get_projector_pair_sptr() const template int PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion:: -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; } @@ -296,13 +297,17 @@ template void PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion:: 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 void PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion:: set_input_data(const shared_ptr & arg) { + this->already_set_up = false; this->_gated_proj_data_sptr = dynamic_pointer_cast(arg); } @@ -319,6 +324,7 @@ void PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion:: set_additive_proj_data_sptr(const shared_ptr &arg) { + this->already_set_up = false; this->_additive_gated_proj_data_sptr = dynamic_pointer_cast(arg); } @@ -327,6 +333,7 @@ void PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion:: set_normalisation_sptr(const shared_ptr& arg) { + this->already_set_up = false; // this->normalisation_sptr = arg; error("Not implemeted yet"); } diff --git a/src/recon_buildblock/GeneralisedObjectiveFunction.cxx b/src/recon_buildblock/GeneralisedObjectiveFunction.cxx index ee25939088..05867223cb 100644 --- a/src/recon_buildblock/GeneralisedObjectiveFunction.cxx +++ b/src/recon_buildblock/GeneralisedObjectiveFunction.cxx @@ -142,6 +142,8 @@ compute_sub_gradient(TargetT& gradient, const TargetT ¤t_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:: 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:: 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); } diff --git a/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMean.cxx b/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMean.cxx index 6af6f741d5..031a346ac6 100644 --- a/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMean.cxx +++ b/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMean.cxx @@ -95,6 +95,7 @@ void PoissonLogLikelihoodWithLinearModelForMean:: set_sensitivity_filename(const std::string& filename) { + this->already_set_up = false; this->sensitivity_filename = filename; } @@ -103,6 +104,7 @@ void PoissonLogLikelihoodWithLinearModelForMean:: set_subsensitivity_filenames(const std::string& filenames) { + this->already_set_up = false; this->subsensitivity_filenames = filenames; try { @@ -171,6 +173,7 @@ void PoissonLogLikelihoodWithLinearModelForMean:: 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:: set_subset_sensitivity_sptr(const shared_ptr& arg, const int subset_num) { + this->already_set_up = false; this->subsensitivity_sptrs[subset_num] = arg; } @@ -348,7 +352,7 @@ set_up(shared_ptr 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 ¤t_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 ¤t_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); } diff --git a/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndListModeData.cxx b/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndListModeData.cxx index bf95a274ef..6af8e55f97 100644 --- a/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndListModeData.cxx +++ b/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndListModeData.cxx @@ -100,10 +100,10 @@ PoissonLogLikelihoodWithLinearModelForMeanAndListModeData::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(this->list_mode_filename); - } + } if (this->additive_projection_data_filename != "0") { @@ -123,6 +123,7 @@ PoissonLogLikelihoodWithLinearModelForMeanAndListModeData::post_process } target_parameter_parser.check_values(); + this->already_set_up = false; return false; } @@ -131,6 +132,7 @@ void PoissonLogLikelihoodWithLinearModelForMeanAndListModeData:: set_input_data(const shared_ptr & arg) { + this->already_set_up = false; try { this->list_mode_data_sptr = dynamic_pointer_cast(arg); @@ -146,6 +148,7 @@ void PoissonLogLikelihoodWithLinearModelForMeanAndListModeData:: 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:: set_additive_proj_data_sptr(const shared_ptr &arg) { + this->already_set_up = false; try { this->additive_proj_data_sptr = dynamic_pointer_cast(arg); @@ -281,6 +285,7 @@ void PoissonLogLikelihoodWithLinearModelForMeanAndListModeData:: set_normalisation_sptr(const shared_ptr& arg) { + this->already_set_up = false; this->normalisation_sptr = arg; } diff --git a/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin.cxx b/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin.cxx index 94bd7c288e..b5467eeaa4 100644 --- a/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin.cxx +++ b/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin.cxx @@ -116,6 +116,7 @@ int PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin:: 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:: set_proj_matrix(const shared_ptr& arg) { - this->PM_sptr = arg; + this->already_set_up = false; + this->PM_sptr = arg; } #if 0 @@ -134,6 +136,7 @@ void PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin:: 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:: construct_target_ptr() const { - - return + return this->target_parameter_parser.create(this->get_input_data()); } diff --git a/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx b/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx index f8401fba05..ac686cf449 100644 --- a/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx +++ b/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx @@ -391,6 +391,7 @@ int PoissonLogLikelihoodWithLinearModelForMeanAndProjData:: 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 = std::max(new_num_subsets,1); return this->num_subsets; @@ -401,6 +402,7 @@ void PoissonLogLikelihoodWithLinearModelForMeanAndProjData:: set_proj_data_sptr(const shared_ptr& arg) { + this->already_set_up = false; this->proj_data_sptr = arg; } @@ -409,6 +411,7 @@ void PoissonLogLikelihoodWithLinearModelForMeanAndProjData:: 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; } @@ -417,6 +420,7 @@ void PoissonLogLikelihoodWithLinearModelForMeanAndProjData:: set_max_timing_pos_num_to_process(const int arg) { + this->already_set_up = this->already_set_up && (this->max_timing_pos_num_to_process == arg); this->max_timing_pos_num_to_process = arg; } @@ -425,6 +429,7 @@ void PoissonLogLikelihoodWithLinearModelForMeanAndProjData:: 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; } @@ -433,7 +438,8 @@ void PoissonLogLikelihoodWithLinearModelForMeanAndProjData:: set_additive_proj_data_sptr(const shared_ptr &arg) { - this->additive_proj_data_sptr = dynamic_pointer_cast(arg); + this->already_set_up = false; + this->additive_proj_data_sptr = dynamic_pointer_cast(arg); } template @@ -441,6 +447,7 @@ void PoissonLogLikelihoodWithLinearModelForMeanAndProjData:: set_projector_pair_sptr(const shared_ptr& arg) { + this->already_set_up = false; this->projector_pair_ptr = arg; } @@ -449,6 +456,7 @@ void PoissonLogLikelihoodWithLinearModelForMeanAndProjData:: set_frame_num(const int arg) { + this->already_set_up = this->already_set_up && (this->frame_num == arg); this->frame_num = arg; } @@ -457,6 +465,7 @@ void PoissonLogLikelihoodWithLinearModelForMeanAndProjData:: set_frame_definitions(const TimeFrameDefinitions& arg) { + this->already_set_up = this->already_set_up && (this->frame_defs == arg); this->frame_defs = arg; } @@ -465,6 +474,7 @@ void PoissonLogLikelihoodWithLinearModelForMeanAndProjData:: set_normalisation_sptr(const shared_ptr& arg) { + this->already_set_up = false; this->normalisation_sptr = arg; } @@ -473,7 +483,8 @@ void PoissonLogLikelihoodWithLinearModelForMeanAndProjData:: set_input_data(const shared_ptr & arg) { - this->proj_data_sptr = dynamic_pointer_cast(arg); + this->already_set_up = false; + this->proj_data_sptr = dynamic_pointer_cast(arg); } template diff --git a/src/recon_test/test_blocks_on_cylindrical_projectors.cxx b/src/recon_test/test_blocks_on_cylindrical_projectors.cxx index d325ff685f..c1d60b6693 100644 --- a/src/recon_test/test_blocks_on_cylindrical_projectors.cxx +++ b/src/recon_test/test_blocks_on_cylindrical_projectors.cxx @@ -217,7 +217,7 @@ BlocksTests::run_symmetry_test(ForwardProjectorByBin& forw_projector1, ForwardPr scannerBlocks_sptr->set_up(); auto proj_data_info_blocks_sptr = std::make_shared(); - proj_data_info_blocks_sptr = set_direct_projdata_info(scannerBlocks_sptr); + proj_data_info_blocks_sptr = set_direct_projdata_info(scannerBlocks_sptr, 2); // now forward-project images // info(boost::format("Test blocks on Cylindrical: Forward projector used: %1%") % forw_projector1.parameter_info());