STIR 6.4.0
PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion.txx
1/*
2 Copyright (C) 2006- 2009, Hammersmith Imanet Ltd
3 Copyright (C) 2011 - 2013, King's College London
4 Copyright (C) 2018, University College London
5 This file is part of STIR.
6
7 SPDX-License-Identifier: Apache-2.0
8
9 See STIR/LICENSE.txt for details
10 */
11/*!
12 \file
13 \ingroup GeneralisedObjectiveFunction
14 \brief Implementation of class stir::PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion
15
16 \author Charalampos Tsoumpas
17 \author Kris Thielemans
18 \author Sanida Mustafovic
19
20*/
21#include "stir/DiscretisedDensity.h"
22#include "stir/is_null_ptr.h"
23#include "stir/recon_buildblock/TrivialBinNormalisation.h"
24#include "stir/Succeeded.h"
25#include "stir/recon_buildblock/ProjectorByBinPair.h"
26#include "stir/info.h"
27#include "stir/warning.h"
28#include "stir/error.h"
29#include "stir/format.h"
30#include "stir/is_null_ptr.h"
31
32// include the following to set defaults
33#ifndef USE_PMRT
34# include "stir/recon_buildblock/ForwardProjectorByBinUsingRayTracing.h"
35# include "stir/recon_buildblock/BackProjectorByBinUsingInterpolation.h"
36#else
37# include "stir/recon_buildblock/ForwardProjectorByBinUsingProjMatrixByBin.h"
38# include "stir/recon_buildblock/BackProjectorByBinUsingProjMatrixByBin.h"
39# include "stir/recon_buildblock/ProjMatrixByBinUsingRayTracing.h"
40#endif
41#include "stir/recon_buildblock/ProjectorByBinPairUsingSeparateProjectors.h"
42
43#include <algorithm>
44#include <string>
45// For Motion
46#include "stir/spatial_transformation/GatedSpatialTransformation.h"
47#include "stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion.h"
48#include "stir/recon_buildblock/BinNormalisationFromProjData.h"
49
50START_NAMESPACE_STIR
51
52template <typename TargetT>
53const char* const PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<TargetT>::registered_name
54 = "PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion";
55
56template <typename TargetT>
57void
58PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<TargetT>::set_defaults()
59{
60 base_type::set_defaults();
61
62 this->_input_filename = "";
63 this->_max_segment_num_to_process = -1; // use all segments
64 // num_views_to_add=1; // KT 20/06/2001 disabled
65
66 this->_gated_proj_data_sptr.reset();
67 this->_zero_seg0_end_planes = 0;
68 this->_reverse_motion_vectors_filename_prefix = "0";
69 this->_normalisation_gated_proj_data_filename = "1";
70 this->_normalisation_gated_proj_data_sptr.reset();
71 // this->_reverse_motion_vectors_sptr=NULL;
72 this->_motion_vectors_filename_prefix = "0";
73 // this->_motion_vectors_sptr=NULL;
74 this->_gate_definitions_filename = "0";
75 // this->_time_gate_definitions_sptr=NULL;
76 this->_additive_gated_proj_data_filename = "0";
77 this->_additive_gated_proj_data_sptr.reset();
78
79#ifndef USE_PMRT // set default for _projector_pair_ptr
80 shared_ptr<ForwardProjectorByBin> forward_projector_ptr(new ForwardProjectorByBinUsingRayTracing());
81 shared_ptr<BackProjectorByBin> back_projector_ptr(new BackProjectorByBinUsingInterpolation());
82#else
83 shared_ptr<ProjMatrixByBin> PM(new ProjMatrixByBinUsingRayTracing());
84 shared_ptr<ForwardProjectorByBin> forward_projector_ptr(new ForwardProjectorByBinUsingProjMatrixByBin(PM));
85 shared_ptr<BackProjectorByBin> back_projector_ptr(new BackProjectorByBinUsingProjMatrixByBin(PM));
86#endif
87
88 this->_projector_pair_ptr.reset(new ProjectorByBinPairUsingSeparateProjectors(forward_projector_ptr, back_projector_ptr));
89
90 this->target_parameter_parser.set_defaults();
91}
92
93template <typename TargetT>
94void
95PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<TargetT>::initialise_keymap()
96{
97 base_type::initialise_keymap();
98 this->parser.add_start_key("PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion Parameters");
99 this->parser.add_stop_key("End PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion Parameters");
100 this->parser.add_key("input filename", &this->_input_filename);
101
102 // parser.add_key("mash x views", &num_views_to_add); // KT 20/06/2001 disabled
103 this->parser.add_key("maximum absolute segment number to process", &this->_max_segment_num_to_process);
104 this->parser.add_key("zero end planes of segment 0", &this->_zero_seg0_end_planes);
105
106 this->target_parameter_parser.add_to_keymap(this->parser);
107 this->parser.add_parsing_key("Projector pair type", &this->_projector_pair_ptr);
108
109 // Scatter correction
110 this->parser.add_key("additive sinograms", &this->_additive_gated_proj_data_filename);
111
112 // normalisation (and attenuation correction)
113 this->parser.add_key("normalisation sinograms", &this->_normalisation_gated_proj_data_filename);
114
115 // Motion Information
116 this->parser.add_key("Gate Definitions filename", &this->_gate_definitions_filename);
117 this->parser.add_key("Motion Vectors filename prefix", &this->_motion_vectors_filename_prefix);
118 this->parser.add_key("Reverse Motion Vectors filename prefix", &this->_reverse_motion_vectors_filename_prefix);
119}
120
121template <typename TargetT>
122bool
123PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<TargetT>::post_processing()
124{
125 if (base_type::post_processing() == true)
126 return true;
127 if (this->_input_filename.length() == 0)
128 {
129 warning("You need to specify an input filename");
130 return true;
131 }
132
133 this->_gated_proj_data_sptr = GatedProjData::read_from_file(this->_input_filename);
134
135 // image stuff
136 this->target_parameter_parser.check_values();
137
138 if (this->_additive_gated_proj_data_filename != "0")
139 {
140 info(format("Reading additive projdata data {}", this->_additive_gated_proj_data_filename));
141 this->_additive_gated_proj_data_sptr = GatedProjData::read_from_file(this->_additive_gated_proj_data_filename);
142 }
143 if (this->_normalisation_gated_proj_data_filename != "1")
144 {
145 info(format("Reading normalisation projdata data {}", this->_normalisation_gated_proj_data_filename));
146 this->_normalisation_gated_proj_data_sptr = GatedProjData::read_from_file(this->_normalisation_gated_proj_data_filename);
147 }
148
149 this->_time_gate_definitions.read_gdef_file(this->_gate_definitions_filename);
150
151 if (this->_reverse_motion_vectors_filename_prefix != "0")
152 this->_reverse_motion_vectors.read_from_files(this->_reverse_motion_vectors_filename_prefix);
153 if (this->_motion_vectors_filename_prefix != "0")
154 this->_motion_vectors.read_from_files(this->_motion_vectors_filename_prefix);
155 return false;
156}
157
158template <typename TargetT>
159PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<
160 TargetT>::PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion()
161{
162 this->set_defaults();
163}
164
165template <typename TargetT>
166TargetT*
167PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<TargetT>::construct_target_ptr() const
168{
169 return this->target_parameter_parser.create(this->get_input_data());
170}
171/***************************************************************
172 subset balancing
173***************************************************************/
174
175template <typename TargetT>
176bool
177PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<TargetT>::actual_subsets_are_approximately_balanced(
178 std::string& warning_message) const
179{ // call actual_subsets_are_approximately_balanced() for first single_gate_obj_func
180 if (this->get_time_gate_definitions().get_num_gates() == 0 || this->_single_gate_obj_funcs.size() == 0)
181 error("PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion:\n"
182 "actual_subsets_are_approximately_balanced called but no gates yet.");
183 else if (this->_single_gate_obj_funcs.size() != 0)
184 {
185 bool gates_are_balanced = true;
186 for (unsigned int gate_num = 1; gate_num <= this->get_time_gate_definitions().get_num_gates(); ++gate_num)
187 gates_are_balanced &= this->_single_gate_obj_funcs[gate_num].subsets_are_approximately_balanced(warning_message);
188 return gates_are_balanced;
189 }
190 else
191 error("Something strange happened in PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion:\n"
192 "actual_subsets_are_approximately_balanced called before setup()?");
193 return false;
194}
195
196/***************************************************************
197 get_ functions
198***************************************************************/
199template <typename TargetT>
200const TimeGateDefinitions&
201PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<TargetT>::get_time_gate_definitions() const
202{
203 return this->_time_gate_definitions;
204}
205
206template <typename TargetT>
207const GatedProjData&
208PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<TargetT>::get_gated_proj_data() const
209{
210 return *this->_gated_proj_data_sptr;
211}
212
213template <typename TargetT>
214const shared_ptr<GatedProjData>&
215PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<TargetT>::get_gated_proj_data_sptr() const
216{
217 return this->_gated_proj_data_sptr;
218}
219
220template <typename TargetT>
221const int
222PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<TargetT>::get_max_segment_num_to_process() const
223{
224 return this->_max_segment_num_to_process;
225}
226
227template <typename TargetT>
228const bool
229PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<TargetT>::get_zero_seg0_end_planes() const
230{
231 return this->_zero_seg0_end_planes;
232}
233
234template <typename TargetT>
235const GatedProjData&
236PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<TargetT>::get_additive_gated_proj_data() const
237{
238 return *this->_additive_gated_proj_data_sptr;
239}
240
241template <typename TargetT>
242const shared_ptr<GatedProjData>&
243PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<TargetT>::get_additive_gated_proj_data_sptr() const
244{
245 return this->_additive_gated_proj_data_sptr;
246}
247
248template <typename TargetT>
249const GatedProjData&
250PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<TargetT>::get_normalisation_gated_proj_data() const
251{
252 return *this->_normalisation_gated_proj_data_sptr;
253}
254
255template <typename TargetT>
256const shared_ptr<GatedProjData>&
257PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<TargetT>::get_normalisation_gated_proj_data_sptr() const
258{
259 return this->_normalisation_gated_proj_data_sptr;
260}
261
262template <typename TargetT>
263const ProjectorByBinPair&
264PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<TargetT>::get_projector_pair() const
265{
266 return *this->_projector_pair_ptr;
267}
268
269template <typename TargetT>
270const shared_ptr<ProjectorByBinPair>&
271PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<TargetT>::get_projector_pair_sptr() const
272{
273 return this->_projector_pair_ptr;
274}
275
276/***************************************************************
277 set_ functions
278***************************************************************/
279template <typename TargetT>
280int
281PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<TargetT>::set_num_subsets(const int new_num_subsets)
282{
283 this->already_set_up = this->already_set_up && (this->num_subsets == new_num_subsets);
284 for (unsigned int gate_num = 1; gate_num <= this->get_time_gate_definitions().get_num_gates(); ++gate_num)
285 {
286 if (this->_single_gate_obj_funcs.size() != 0)
287 if (this->_single_gate_obj_funcs[gate_num].set_num_subsets(new_num_subsets) != new_num_subsets)
288 error("set_num_subsets didn't work");
289 }
290 this->num_subsets = new_num_subsets;
291 return this->num_subsets;
292}
293
294template <typename TargetT>
295void
296PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<TargetT>::set_time_gate_definitions(
297 const TimeGateDefinitions& time_gate_definitions)
298{
299 this->already_set_up = false;
300 this->_time_gate_definitions = time_gate_definitions;
301}
302
303template <typename TargetT>
304void
305PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<TargetT>::set_input_data(const shared_ptr<ExamData>& arg)
306{
307 this->already_set_up = false;
308 this->_gated_proj_data_sptr = dynamic_pointer_cast<GatedProjData>(arg);
309}
310
311template <typename TargetT>
312const GatedProjData&
313PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<TargetT>::get_input_data() const
314{
315 return *this->_gated_proj_data_sptr;
316}
317
318template <typename TargetT>
319void
320PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<TargetT>::set_additive_proj_data_sptr(
321 const shared_ptr<ExamData>& arg)
322{
323 this->already_set_up = false;
324 this->_additive_gated_proj_data_sptr = dynamic_pointer_cast<GatedProjData>(arg);
325}
326
327template <typename TargetT>
328void
329PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<TargetT>::set_normalisation_sptr(
330 const shared_ptr<BinNormalisation>& arg)
331{
332 this->already_set_up = false;
333 // this->normalisation_sptr = arg;
334 error("Not implemeted yet");
335}
336
337/***************************************************************
338 set_up()
339***************************************************************/
340template <typename TargetT>
341Succeeded
342PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<TargetT>::set_up_before_sensitivity(
343 shared_ptr<const TargetT> const& target_sptr)
344{
345 /*!todo define in the PoissonLogLikelihoodWithLinearModelForMean class to return Succeeded::yes
346 if (base_type::set_up_before_sensitivity(target_sptr) != Succeeded::yes)
347 return Succeeded::no;
348 */
349 if (this->_max_segment_num_to_process == -1)
350 this->_max_segment_num_to_process = (this->_gated_proj_data_sptr)->get_proj_data_sptr(1)->get_max_segment_num();
351
352 if (this->_max_segment_num_to_process > (this->_gated_proj_data_sptr)->get_proj_data_sptr(1)->get_max_segment_num())
353 {
354 warning("_max_segment_num_to_process (%d) is too large", this->_max_segment_num_to_process);
355 return Succeeded::no;
356 }
357
358 shared_ptr<ProjDataInfo> proj_data_info_sptr(
359 (this->_gated_proj_data_sptr->get_proj_data_sptr(1))->get_proj_data_info_sptr()->clone());
360 proj_data_info_sptr->reduce_segment_range(-this->_max_segment_num_to_process, +this->_max_segment_num_to_process);
361
362 if (is_null_ptr(this->_projector_pair_ptr))
363 {
364 warning("You need to specify a projector pair");
365 return Succeeded::no;
366 }
367
368 if (this->num_subsets <= 0)
369 {
370 warning("Number of subsets %d should be larger than 0.", this->num_subsets);
371 return Succeeded::no;
372 }
373 {
374 const shared_ptr<DiscretisedDensity<3, float>> density_template_sptr(
375 target_sptr->get_empty_copy()); // target_sptr appears not to be set up correctly
376 const shared_ptr<Scanner> scanner_sptr(new Scanner(*proj_data_info_sptr->get_scanner_ptr()));
377 this->_gated_image_template = GatedDiscretisedDensity(this->get_time_gate_definitions(), density_template_sptr);
378
379 // construct _single_gate_obj_funcs
380 this->_single_gate_obj_funcs.resize(1, this->get_time_gate_definitions().get_num_gates());
381
382 for (unsigned int gate_num = 1; gate_num <= this->get_time_gate_definitions().get_num_gates(); ++gate_num)
383 {
384 info(format("Objective Function for Gate Number: {}", gate_num));
385 this->_single_gate_obj_funcs[gate_num].set_projector_pair_sptr(this->_projector_pair_ptr);
386 this->_single_gate_obj_funcs[gate_num].set_proj_data_sptr(this->_gated_proj_data_sptr->get_proj_data_sptr(gate_num));
387 this->_single_gate_obj_funcs[gate_num].set_max_segment_num_to_process(this->_max_segment_num_to_process);
388 this->_single_gate_obj_funcs[gate_num].set_zero_seg0_end_planes(this->_zero_seg0_end_planes != 0);
389 if (!is_null_ptr(this->_additive_gated_proj_data_sptr))
390 this->_single_gate_obj_funcs[gate_num].set_additive_proj_data_sptr(
391 this->_additive_gated_proj_data_sptr->get_proj_data_sptr(gate_num));
392 this->_single_gate_obj_funcs[gate_num].set_num_subsets(this->num_subsets);
393 this->_single_gate_obj_funcs[gate_num].set_frame_num(1); // This should be gate...
394 std::vector<std::pair<double, double>> frame_times(1, std::pair<double, double>(0, 1));
395 this->_single_gate_obj_funcs[gate_num].set_frame_definitions(TimeFrameDefinitions(frame_times));
396
397 shared_ptr<BinNormalisation> current_gate_norm_factors_sptr;
398 if (is_null_ptr(this->_normalisation_gated_proj_data_sptr))
399 current_gate_norm_factors_sptr.reset(new TrivialBinNormalisation);
400 else
401 {
402 shared_ptr<ProjData> norm_data_sptr(this->_normalisation_gated_proj_data_sptr->get_proj_data_sptr(gate_num));
403 current_gate_norm_factors_sptr.reset(new BinNormalisationFromProjData(norm_data_sptr));
404 }
405 this->_single_gate_obj_funcs[gate_num].set_normalisation_sptr(current_gate_norm_factors_sptr);
406 this->_single_gate_obj_funcs[gate_num].set_recompute_sensitivity(this->get_recompute_sensitivity());
407 this->_single_gate_obj_funcs[gate_num].set_use_subset_sensitivities(this->get_use_subset_sensitivities());
408
409 if (this->_single_gate_obj_funcs[gate_num].set_up(density_template_sptr) != Succeeded::yes)
410 error("Single gate objective functions is not set correctly!");
411 }
412 } //_single_gate_obj_funcs[gate_num]
413 return Succeeded::yes;
414}
415
416/*************************************************************************
417 functions that compute the value/gradient of the objective function etc
418*************************************************************************/
419
420template <typename TargetT>
421void
422PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<TargetT>::actual_compute_subset_gradient_without_penalty(
423 TargetT& gradient, const TargetT& current_estimate, const int subset_num, const bool add_sensitivity)
424{
425 assert(subset_num >= 0);
426 assert(subset_num < this->num_subsets);
427
428 GatedDiscretisedDensity gated_gradient = this->_gated_image_template;
429 GatedDiscretisedDensity gated_image_estimate = this->_gated_image_template;
430 // The following initialization doesn't stabilize reconstruction.
431 for (unsigned int gate_num = 1; gate_num <= this->get_time_gate_definitions().get_num_gates(); ++gate_num)
432 std::fill(gated_image_estimate[gate_num].begin_all(), gated_image_estimate[gate_num].end_all(), 0.F);
433 this->_motion_vectors.warp_image(gated_image_estimate, current_estimate);
434 for (unsigned int gate_num = 1; gate_num <= this->get_time_gate_definitions().get_num_gates(); ++gate_num)
435 {
436 std::fill(gated_gradient[gate_num].begin_all(), gated_gradient[gate_num].end_all(), 0.F);
437 this->_single_gate_obj_funcs[gate_num].actual_compute_subset_gradient_without_penalty(
438 gated_gradient[gate_num], gated_image_estimate[gate_num], subset_num, add_sensitivity);
439 }
440 // if(this->_motion_correction_type==-1)
441 this->_reverse_motion_vectors.warp_image(gradient, gated_gradient);
442 // else
443 // this->_motion_vectors.warp_image(gradient,gated_gradient) ;
444}
445
446template <typename TargetT>
447double
448PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<TargetT>::actual_compute_objective_function_without_penalty(
449 const TargetT& current_estimate, const int subset_num)
450{
451 assert(subset_num >= 0);
452 assert(subset_num < this->num_subsets);
453
454 double result = 0.;
455 GatedDiscretisedDensity gated_image_estimate = this->_gated_image_template;
456 // The following initialization doesn't stabilize reconstruction.
457 for (unsigned int gate_num = 1; gate_num <= this->get_time_gate_definitions().get_num_gates(); ++gate_num)
458 std::fill(gated_image_estimate[gate_num].begin_all(), gated_image_estimate[gate_num].end_all(), 0.F);
459 this->_motion_vectors.warp_image(gated_image_estimate, current_estimate);
460 // loop over single_gate
461 for (unsigned int gate_num = 1; gate_num <= this->get_time_gate_definitions().get_num_gates(); ++gate_num)
462 {
463 result += this->_single_gate_obj_funcs[gate_num].compute_objective_function_without_penalty(gated_image_estimate[gate_num],
464 subset_num);
465 }
466 return result;
467}
468
469template <typename TargetT>
470void
471PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<TargetT>::add_subset_sensitivity(TargetT& sensitivity,
472 const int subset_num) const
473{
474 GatedDiscretisedDensity gated_subset_sensitivity = this->_gated_image_template;
475
476 // loop over single_gate to get original subset sensitivity
477 for (unsigned int gate_num = 1; gate_num <= this->get_time_gate_definitions().get_num_gates(); ++gate_num)
478 {
479 const shared_ptr<DiscretisedDensity<3, float>> single_gate_subsensitivity_sptr(
480 this->_single_gate_obj_funcs[gate_num].get_subset_sensitivity(subset_num).clone());
481 gated_subset_sensitivity.set_density_sptr(single_gate_subsensitivity_sptr, gate_num);
482 }
483
484 // perform warp
485 this->_reverse_motion_vectors.accumulate_warp_image(sensitivity, gated_subset_sensitivity);
486}
487//! /todo The
488//! PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<TargetT>::actual_add_multiplication_with_approximate_sub_Hessian_without_penalty
489//! is not validated and at the moment OSSPS does not converge with motion correction.
490template <typename TargetT>
491Succeeded
492PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<
493 TargetT>::actual_add_multiplication_with_approximate_sub_Hessian_without_penalty(TargetT& output,
494 const TargetT& input,
495 const int subset_num) const
496{
497 // TODO this does not add but replace
498 {
499 std::string explanation;
500 if (!input.has_same_characteristics(this->get_subset_sensitivity(0), ////////////////////
501 explanation))
502 {
503 warning("PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion:\n"
504 "sensitivity and input for add_multiplication_with_approximate_Hessian_without_penalty\n"
505 "should have the same characteristics.\n%s",
506 explanation.c_str());
507 return Succeeded::no;
508 }
509 }
510 GatedDiscretisedDensity gated_input = this->_gated_image_template;
511 GatedDiscretisedDensity gated_output = this->_gated_image_template;
512 this->_motion_vectors.warp_image(gated_input, input);
513
514 VectorWithOffset<float> scale_factor(1, this->get_time_gate_definitions().get_num_gates());
515 for (unsigned int gate_num = 1; gate_num <= this->get_time_gate_definitions().get_num_gates(); ++gate_num)
516 {
517 scale_factor[gate_num] = gated_input[gate_num].find_max();
518 /*! /note This is used to avoid higher values than these set in the precompute_denominator_of_conditioner_without_penalty()
519 function. /sa for more information see the recon_array_functions.cxx and the value of the max_quotient (originaly set to
520 10000.F) */
521 gated_input[gate_num] /= scale_factor[gate_num];
522 this->_single_gate_obj_funcs[gate_num].add_multiplication_with_approximate_sub_Hessian_without_penalty(
523 gated_output[gate_num], gated_input[gate_num], subset_num);
524 gated_output[gate_num] *= scale_factor[gate_num];
525 } // end of loop over gates
526 this->_reverse_motion_vectors.warp_image(output, gated_output);
527 output /= static_cast<float>(
528 this->get_time_gate_definitions().get_num_gates()); // Normalizing to get the average value to test if OSSPS works.
529 return Succeeded::yes;
530}
531
532template <typename TargetT>
533Succeeded
534PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<
535 TargetT>::actual_accumulate_sub_Hessian_times_input_without_penalty(TargetT& output,
536 const TargetT& current_image_estimate,
537 const TargetT& input,
538 const int subset_num) const
539{
540 { // check argument characteristics
541 std::string explanation;
542 if (!input.has_same_characteristics(this->get_subset_sensitivity(0), explanation))
543 {
544 warning("PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion:\n"
545 "sensitivity and input for actual_accumulate_sub_Hessian_times_input_without_penalty\n"
546 "should have the same characteristics.\n%s",
547 explanation.c_str());
548 return Succeeded::no;
549 }
550
551 if (!current_image_estimate.has_same_characteristics(this->get_subset_sensitivity(0), explanation))
552 {
553 warning("PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion:\n"
554 "sensitivity and current_image_estimate for actual_accumulate_sub_Hessian_times_input_without_penalty\n"
555 "should have the same characteristics.\n%s",
556 explanation.c_str());
557 return Succeeded::no;
558 }
559 }
560
561 GatedDiscretisedDensity gated_input = this->_gated_image_template;
562 GatedDiscretisedDensity gated_current_image_estimate = this->_gated_image_template;
563 GatedDiscretisedDensity gated_output = this->_gated_image_template;
564 this->_motion_vectors.warp_image(gated_input, input);
565 this->_motion_vectors.warp_image(gated_current_image_estimate, current_image_estimate);
566
567 for (unsigned int gate_num = 1; gate_num <= this->get_time_gate_definitions().get_num_gates(); ++gate_num)
568 {
569 this->_single_gate_obj_funcs[gate_num].accumulate_sub_Hessian_times_input_without_penalty(
570 gated_output[gate_num], gated_current_image_estimate[gate_num], gated_input[gate_num], subset_num);
571 } // end of loop over gates
572 this->_reverse_motion_vectors.warp_image(output, gated_output);
573 output /= static_cast<float>(
574 this->get_time_gate_definitions().get_num_gates()); // Normalizing to get the average value to test if OSSPS works.
575 return Succeeded::yes;
576}
577
578END_NAMESPACE_STIR
Generated on Wed May 27 2026 17:04:05 for STIR by doxygen 1.13.2