STIR 6.4.0
PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData.txx
1//
2// PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData.txx
3//
4/*
5 Copyright (C) 2006 - 2011-01-14 Hammersmith Imanet Ltd
6 Copyright (C) 2011 Kris Thielemans
7 Copyright (C) 2013, 2018, University College London
8
9 This file is part of STIR.
10
11 SPDX-License-Identifier: Apache-2.0
12
13 See STIR/LICENSE.txt for details
14*/
15/*!
16 \file
17 \ingroup GeneralisedObjectiveFunction
18 \ingroup modelling
19 \brief Implementation of class stir::PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData
20
21 \author Kris Thielemans
22 \author Sanida Mustafovic
23 \author Charalampos Tsoumpas
24*/
25#include "stir/DiscretisedDensity.h"
26#include "stir/is_null_ptr.h"
27#include "stir/recon_buildblock/TrivialBinNormalisation.h"
28#include "stir/Succeeded.h"
29#include "stir/recon_buildblock/ProjectorByBinPair.h"
30#include "stir/info.h"
31#include "stir/warning.h"
32#include "stir/error.h"
33#include "stir/format.h"
34
35// include the following to set defaults
36#ifndef USE_PMRT
37# include "stir/recon_buildblock/ForwardProjectorByBinUsingRayTracing.h"
38# include "stir/recon_buildblock/BackProjectorByBinUsingInterpolation.h"
39#else
40# include "stir/recon_buildblock/ForwardProjectorByBinUsingProjMatrixByBin.h"
41# include "stir/recon_buildblock/BackProjectorByBinUsingProjMatrixByBin.h"
42# include "stir/recon_buildblock/ProjMatrixByBinUsingRayTracing.h"
43#endif
44#include "stir/recon_buildblock/ProjectorByBinPairUsingSeparateProjectors.h"
45#include "stir/recon_buildblock/BinNormalisationWithCalibration.h"
46
47#include <algorithm>
48#include <string>
49// For the Patlak Plot Modelling
50#include "stir/modelling/ModelMatrix.h"
51#include "stir/recon_buildblock/PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData.h"
52#ifndef NDEBUG
53# include "stir/IO/write_to_file.h"
54#endif
55
56START_NAMESPACE_STIR
57
58template <typename TargetT>
59const char* const PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData<TargetT>::registered_name
60 = "PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData";
61
62template <typename TargetT>
63void
64PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData<TargetT>::set_defaults()
65{
66 base_type::set_defaults();
67
68 this->_input_filename = "";
69 this->_max_segment_num_to_process = -1; // use all segments
70 // num_views_to_add=1; // KT 20/06/2001 disabled
71
72 this->_dyn_proj_data_sptr.reset();
73 this->_zero_seg0_end_planes = 0;
74
75 this->_additive_dyn_proj_data_filename = "0";
76 this->_additive_dyn_proj_data_sptr.reset();
77
78#ifndef USE_PMRT // set default for _projector_pair_ptr
79 shared_ptr<ForwardProjectorByBin> forward_projector_ptr(new ForwardProjectorByBinUsingRayTracing());
80 shared_ptr<BackProjectorByBin> back_projector_ptr(new BackProjectorByBinUsingInterpolation());
81#else
82 shared_ptr<ProjMatrixByBin> PM(new ProjMatrixByBinUsingRayTracing());
83 shared_ptr<ForwardProjectorByBin> forward_projector_ptr(new ForwardProjectorByBinUsingProjMatrixByBin(PM));
84 shared_ptr<BackProjectorByBin> back_projector_ptr(new BackProjectorByBinUsingProjMatrixByBin(PM));
85#endif
86
87 this->_projector_pair_ptr.reset(new ProjectorByBinPairUsingSeparateProjectors(forward_projector_ptr, back_projector_ptr));
88 this->_normalisation_sptr.reset(new TrivialBinNormalisation);
89
90 this->target_parameter_parser.set_defaults();
91
92 // Modelling Stuff
93 this->_patlak_plot_sptr.reset();
94}
95
96template <typename TargetT>
97void
98PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData<TargetT>::initialise_keymap()
99{
100 base_type::initialise_keymap();
101 this->parser.add_start_key("PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData Parameters");
102 this->parser.add_stop_key("End PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData Parameters");
103 this->parser.add_key("input file", &this->_input_filename);
104
105 // parser.add_key("mash x views", &num_views_to_add); // KT 20/06/2001 disabled
106 this->parser.add_key("maximum absolute segment number to process", &this->_max_segment_num_to_process);
107 this->parser.add_key("zero end planes of segment 0", &this->_zero_seg0_end_planes);
108
109 this->target_parameter_parser.add_to_keymap(this->parser);
110 this->parser.add_parsing_key("Projector pair type", &this->_projector_pair_ptr);
111
112 // Scatter correction
113 this->parser.add_key("additive sinograms", &this->_additive_dyn_proj_data_filename);
114
115 // normalisation (and attenuation correction)
116 this->parser.add_parsing_key("Bin Normalisation type", &this->_normalisation_sptr);
117
118 // Modelling Information
119 this->parser.add_parsing_key("Kinetic Model Type", &this->_patlak_plot_sptr); // Do sth with dynamic_cast to get the PatlakPlot
120
121 // Regularization Information
122 // this->parser.add_parsing_key("prior type", &this->_prior_sptr);
123}
124
125template <typename TargetT>
126bool
127PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData<TargetT>::post_processing()
128{
129 if (base_type::post_processing() == true)
130 return true;
131 if (this->_input_filename.length() == 0)
132 {
133 warning("You need to specify an input filename");
134 return true;
135 }
136
137#if 0 // KT 20/06/2001 disabled as not functional yet
138 if (num_views_to_add!=1 && (num_views_to_add<=0 || num_views_to_add%2 != 0))
139 { warning("The 'mash x views' key has an invalid value (must be 1 or even number)"); return true; }
140#endif
141
142 this->_dyn_proj_data_sptr = DynamicProjData::read_from_file(_input_filename);
143 if (is_null_ptr(this->_dyn_proj_data_sptr))
144 {
145 warning("Error reading input file %s", _input_filename.c_str());
146 return true;
147 }
148 this->target_parameter_parser.check_values();
149
150 if (this->_additive_dyn_proj_data_filename != "0")
151 {
152 info(format("Reading additive projdata data {}", this->_additive_dyn_proj_data_filename));
153 this->_additive_dyn_proj_data_sptr = DynamicProjData::read_from_file(this->_additive_dyn_proj_data_filename);
154 if (is_null_ptr(this->_additive_dyn_proj_data_sptr))
155 {
156 warning("Error reading additive input file %s", _additive_dyn_proj_data_filename.c_str());
157 return true;
158 }
159 }
160 return false;
161}
162
163template <typename TargetT>
164PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData<
165 TargetT>::PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData()
166{
167 this->set_defaults();
168}
169
170template <typename TargetT>
171TargetT*
172PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData<TargetT>::construct_target_ptr() const
173{
174 return this->target_parameter_parser.create(this->get_input_data());
175}
176
177template <typename TargetT>
178std::unique_ptr<ExamInfo>
179PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData<TargetT>::get_exam_info_uptr_for_target() const
180{
181 auto exam_info_uptr = this->get_exam_info_uptr_for_target();
182 if (auto norm_ptr = dynamic_cast<BinNormalisationWithCalibration const* const>(get_normalisation_sptr().get()))
183 {
184 exam_info_uptr->set_calibration_factor(norm_ptr->get_calibration_factor());
185 // somehow tell the image that it's calibrated (do we have a way?)
186 }
187 else
188 {
189 exam_info_uptr->set_calibration_factor(1.F);
190 // somehow tell the image that it's not calibrated (do we have a way?)
191 }
192 return exam_info_uptr;
193}
194
195/***************************************************************
196 subset balancing
197***************************************************************/
198
199template <typename TargetT>
200bool
201PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData<TargetT>::actual_subsets_are_approximately_balanced(
202 std::string& warning_message) const
203{ // call actual_subsets_are_approximately_balanced( for first single_frame_obj_func
204 if (this->_patlak_plot_sptr->get_time_frame_definitions().get_num_frames() == 0 || this->_single_frame_obj_funcs.size() == 0)
205 error("PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData:\n"
206 "actual_subsets_are_approximately_balanced called but not frames yet.\n");
207 else if (this->_single_frame_obj_funcs.size() != 0)
208 {
209 bool frames_are_balanced = true;
210 for (unsigned int frame_num = this->_patlak_plot_sptr->get_starting_frame();
211 frame_num <= this->_patlak_plot_sptr->get_ending_frame();
212 ++frame_num)
213 frames_are_balanced &= this->_single_frame_obj_funcs[frame_num].subsets_are_approximately_balanced(warning_message);
214 return frames_are_balanced;
215 }
216 else
217 error("Something strange happened in PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData:\n"
218 "actual_subsets_are_approximately_balanced called before setup()?\n");
219 return false;
220}
221
222/***************************************************************
223 get_ functions
224***************************************************************/
225template <typename TargetT>
226const DynamicProjData&
227PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData<TargetT>::get_dyn_proj_data() const
228{
229 return *this->_dyn_proj_data_sptr;
230}
231
232template <typename TargetT>
233const shared_ptr<DynamicProjData>&
234PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData<TargetT>::get_dyn_proj_data_sptr() const
235{
236 return this->_dyn_proj_data_sptr;
237}
238
239template <typename TargetT>
240const int
241PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData<TargetT>::get_max_segment_num_to_process() const
242{
243 return this->_max_segment_num_to_process;
244}
245
246template <typename TargetT>
247const bool
248PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData<TargetT>::get_zero_seg0_end_planes() const
249{
250 return this->_zero_seg0_end_planes;
251}
252
253template <typename TargetT>
254const DynamicProjData&
255PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData<TargetT>::get_additive_dyn_proj_data() const
256{
257 return *this->_additive_dyn_proj_data_sptr;
258}
259
260template <typename TargetT>
261const shared_ptr<DynamicProjData>&
262PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData<TargetT>::get_additive_dyn_proj_data_sptr() const
263{
264 return this->_additive_dyn_proj_data_sptr;
265}
266
267template <typename TargetT>
268const ProjectorByBinPair&
269PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData<TargetT>::get_projector_pair() const
270{
271 return *this->_projector_pair_ptr;
272}
273
274template <typename TargetT>
275const shared_ptr<ProjectorByBinPair>&
276PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData<TargetT>::get_projector_pair_sptr() const
277{
278 return this->_projector_pair_ptr;
279}
280
281template <typename TargetT>
282const BinNormalisation&
283PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData<TargetT>::get_normalisation() const
284{
285 return *this->_normalisation_sptr;
286}
287
288template <typename TargetT>
289const shared_ptr<BinNormalisation>&
290PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData<TargetT>::get_normalisation_sptr() const
291{
292 return this->_normalisation_sptr;
293}
294
295/***************************************************************
296 set_ functions
297***************************************************************/
298template <typename TargetT>
299int
300PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData<TargetT>::set_num_subsets(const int num_subsets)
301{
302 for (unsigned int frame_num = this->_patlak_plot_sptr->get_starting_frame();
303 frame_num <= this->_patlak_plot_sptr->get_ending_frame();
304 ++frame_num)
305 {
306 if (this->_single_frame_obj_funcs.size() != 0)
307 if (this->_single_frame_obj_funcs[frame_num].set_num_subsets(num_subsets) != num_subsets)
308 error("set_num_subsets didn't work");
309 }
310 this->num_subsets = num_subsets;
311 return this->num_subsets;
312}
313
314/***************************************************************
315 set_up()
316***************************************************************/
317template <typename TargetT>
318Succeeded
319PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData<TargetT>::set_up_before_sensitivity(
320 shared_ptr<const TargetT> const& target_sptr)
321{
322 if (this->_max_segment_num_to_process == -1)
323 this->_max_segment_num_to_process = (this->_dyn_proj_data_sptr)->get_proj_data_sptr(1)->get_max_segment_num();
324
325 if (this->_max_segment_num_to_process > (this->_dyn_proj_data_sptr)->get_proj_data_sptr(1)->get_max_segment_num())
326 {
327 warning("_max_segment_num_to_process (%d) is too large", this->_max_segment_num_to_process);
328 return Succeeded::no;
329 }
330
331 const shared_ptr<ProjDataInfo> proj_data_info_sptr(
332 (this->_dyn_proj_data_sptr->get_proj_data_sptr(1))->get_proj_data_info_sptr()->clone());
333 proj_data_info_sptr->reduce_segment_range(-this->_max_segment_num_to_process, +this->_max_segment_num_to_process);
334
335 if (is_null_ptr(this->_projector_pair_ptr))
336 {
337 warning("You need to specify a projector pair");
338 return Succeeded::no;
339 }
340
341 if (this->num_subsets <= 0)
342 {
343 warning("Number of subsets %d should be larger than 0.", this->num_subsets);
344 return Succeeded::no;
345 }
346
347 if (is_null_ptr(this->_normalisation_sptr))
348 {
349 warning("Invalid normalisation object");
350 return Succeeded::no;
351 }
352
353 if (this->_normalisation_sptr->set_up(this->_dyn_proj_data_sptr->get_exam_info_sptr(), proj_data_info_sptr) == Succeeded::no)
354 return Succeeded::no;
355
356 if (this->_patlak_plot_sptr->set_up() == Succeeded::no)
357 return Succeeded::no;
358
359 if (this->_patlak_plot_sptr->get_starting_frame() <= 0
360 || this->_patlak_plot_sptr->get_starting_frame() > this->_patlak_plot_sptr->get_ending_frame())
361 {
362 warning("Starting frame is %d. Generally, it should be a late frame,\nbut in any case it should be less than the number of "
363 "frames %d\nand at least 1.",
364 this->_patlak_plot_sptr->get_starting_frame(),
365 this->_patlak_plot_sptr->get_ending_frame());
366 return Succeeded::no;
367 }
368 {
369 const shared_ptr<DiscretisedDensity<3, float>> density_template_sptr(
370 (target_sptr->construct_single_density(1)).get_empty_copy());
371 const shared_ptr<Scanner> scanner_sptr(new Scanner(*proj_data_info_sptr->get_scanner_ptr()));
372 this->_dyn_image_template = DynamicDiscretisedDensity(this->_patlak_plot_sptr->get_time_frame_definitions(),
373 this->_dyn_proj_data_sptr->get_start_time_in_secs_since_1970(),
374 scanner_sptr,
375 density_template_sptr);
376
377 // construct _single_frame_obj_funcs
378 this->_single_frame_obj_funcs.resize(this->_patlak_plot_sptr->get_starting_frame(),
379 this->_patlak_plot_sptr->get_ending_frame());
380
381 for (unsigned int frame_num = this->_patlak_plot_sptr->get_starting_frame();
382 frame_num <= this->_patlak_plot_sptr->get_ending_frame();
383 ++frame_num)
384 {
385 this->_single_frame_obj_funcs[frame_num].set_projector_pair_sptr(this->_projector_pair_ptr);
386 this->_single_frame_obj_funcs[frame_num].set_proj_data_sptr(this->_dyn_proj_data_sptr->get_proj_data_sptr(frame_num));
387 this->_single_frame_obj_funcs[frame_num].set_max_segment_num_to_process(this->_max_segment_num_to_process);
388 this->_single_frame_obj_funcs[frame_num].set_zero_seg0_end_planes(this->_zero_seg0_end_planes != 0);
389 if (!is_null_ptr(this->_additive_dyn_proj_data_sptr))
390 this->_single_frame_obj_funcs[frame_num].set_additive_proj_data_sptr(
391 this->_additive_dyn_proj_data_sptr->get_proj_data_sptr(frame_num));
392 this->_single_frame_obj_funcs[frame_num].set_num_subsets(this->num_subsets);
393 this->_single_frame_obj_funcs[frame_num].set_frame_num(frame_num);
394 this->_single_frame_obj_funcs[frame_num].set_frame_definitions(this->_patlak_plot_sptr->get_time_frame_definitions());
395 this->_single_frame_obj_funcs[frame_num].set_normalisation_sptr(this->_normalisation_sptr);
396 this->_single_frame_obj_funcs[frame_num].set_recompute_sensitivity(this->get_recompute_sensitivity());
397 this->_single_frame_obj_funcs[frame_num].set_use_subset_sensitivities(this->get_use_subset_sensitivities());
398 if (this->_single_frame_obj_funcs[frame_num].set_up(density_template_sptr) != Succeeded::yes)
399 error("Single frame objective functions is not set correctly!");
400 }
401 } //_single_frame_obj_funcs[frame_num]
402
403 return Succeeded::yes;
404}
405
406template <typename TargetT>
407void
408PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData<TargetT>::set_input_data(const shared_ptr<ExamData>& arg)
409{
410 this->_dyn_proj_data_sptr = dynamic_pointer_cast<DynamicProjData>(arg);
411}
412
413template <typename TargetT>
414const DynamicProjData&
415PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData<TargetT>::get_input_data() const
416{
417 return *this->_dyn_proj_data_sptr;
418}
419
420template <typename TargetT>
421void
422PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData<TargetT>::set_additive_proj_data_sptr(
423 const shared_ptr<ExamData>& arg)
424{
425 this->_additive_dyn_proj_data_sptr = dynamic_pointer_cast<DynamicProjData>(arg);
426}
427
428template <typename TargetT>
429void
430PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData<TargetT>::set_normalisation_sptr(
431 const shared_ptr<BinNormalisation>& arg)
432{
433 // this->normalisation_sptr = arg;
434 error("Not implemeted yet");
435}
436
437/*************************************************************************
438 functions that compute the value/gradient of the objective function etc
439*************************************************************************/
440
441template <typename TargetT>
442void
443PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData<TargetT>::actual_compute_subset_gradient_without_penalty(
444 TargetT& gradient, const TargetT& current_estimate, const int subset_num, const bool add_sensitivity)
445{
446 if (subset_num < 0 || subset_num >= this->get_num_subsets())
447 error("compute_sub_gradient_without_penalty subset_num out-of-range error");
448
449 DynamicDiscretisedDensity dyn_gradient = this->_dyn_image_template;
450 DynamicDiscretisedDensity dyn_image_estimate = this->_dyn_image_template;
451
452 for (unsigned int frame_num = this->_patlak_plot_sptr->get_starting_frame();
453 frame_num <= this->_patlak_plot_sptr->get_ending_frame();
454 ++frame_num)
455 std::fill(dyn_image_estimate[frame_num].begin_all(), dyn_image_estimate[frame_num].end_all(), 1.F);
456
457 this->_patlak_plot_sptr->get_dynamic_image_from_parametric_image(dyn_image_estimate, current_estimate);
458
459 // loop over single_frame and use model_matrix
460 for (unsigned int frame_num = this->_patlak_plot_sptr->get_starting_frame();
461 frame_num <= this->_patlak_plot_sptr->get_ending_frame();
462 ++frame_num)
463 {
464 std::fill(dyn_gradient[frame_num].begin_all(), dyn_gradient[frame_num].end_all(), 1.F);
465
466 this->_single_frame_obj_funcs[frame_num].actual_compute_subset_gradient_without_penalty(
467 dyn_gradient[frame_num], dyn_image_estimate[frame_num], subset_num, add_sensitivity);
468 }
469
470 this->_patlak_plot_sptr->multiply_dynamic_image_with_model_gradient(gradient, dyn_gradient);
471}
472
473template <typename TargetT>
474double
475PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData<TargetT>::actual_compute_objective_function_without_penalty(
476 const TargetT& current_estimate, const int subset_num)
477{
478 assert(subset_num >= 0);
479 assert(subset_num < this->num_subsets);
480
481 double result = 0.;
482 DynamicDiscretisedDensity dyn_image_estimate = this->_dyn_image_template;
483
484 // TODO why fill with 1?
485 for (unsigned int frame_num = this->_patlak_plot_sptr->get_starting_frame();
486 frame_num <= this->_patlak_plot_sptr->get_ending_frame();
487 ++frame_num)
488 std::fill(dyn_image_estimate[frame_num].begin_all(), dyn_image_estimate[frame_num].end_all(), 1.F);
489 this->_patlak_plot_sptr->get_dynamic_image_from_parametric_image(dyn_image_estimate, current_estimate);
490
491 // loop over single_frame
492 for (unsigned int frame_num = this->_patlak_plot_sptr->get_starting_frame();
493 frame_num <= this->_patlak_plot_sptr->get_ending_frame();
494 ++frame_num)
495 {
496 result += this->_single_frame_obj_funcs[frame_num].compute_objective_function_without_penalty(dyn_image_estimate[frame_num],
497 subset_num);
498 }
499 return result;
500}
501
502template <typename TargetT>
503void
504PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData<TargetT>::add_subset_sensitivity(TargetT& sensitivity,
505 const int subset_num) const
506{
507 DynamicDiscretisedDensity dyn_sensitivity = this->_dyn_image_template;
508
509 // loop over single_frame and use model_matrix
510 for (unsigned int frame_num = this->_patlak_plot_sptr->get_starting_frame();
511 frame_num <= this->_patlak_plot_sptr->get_ending_frame();
512 ++frame_num)
513 {
514 dyn_sensitivity[frame_num] = this->_single_frame_obj_funcs[frame_num].get_subset_sensitivity(subset_num);
515 // add_subset_sensitivity(dyn_sensitivity[frame_num],subset_num);
516 }
517
518 this->_patlak_plot_sptr->multiply_dynamic_image_with_model_gradient_and_add_to_input(sensitivity, dyn_sensitivity);
519}
520
521template <typename TargetT>
522Succeeded
523PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData<
524 TargetT>::actual_add_multiplication_with_approximate_sub_Hessian_without_penalty(TargetT& output,
525 const TargetT& input,
526 const int subset_num) const
527{
528 {
529 std::string explanation;
530 if (!input.has_same_characteristics(this->get_sensitivity(), explanation))
531 {
532 warning("PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData:\n"
533 "sensitivity and input for add_multiplication_with_approximate_Hessian_without_penalty\n"
534 "should have the same characteristics.\n%s",
535 explanation.c_str());
536 return Succeeded::no;
537 }
538 }
539#ifndef NDEBUG
540 info(
541 format("INPUT max: ({} , {})", input.construct_single_density(1).find_max(), input.construct_single_density(2).find_max()));
542#endif // NDEBUG
543 DynamicDiscretisedDensity dyn_input = this->_dyn_image_template;
544 DynamicDiscretisedDensity dyn_output = this->_dyn_image_template;
545 this->_patlak_plot_sptr->get_dynamic_image_from_parametric_image(dyn_input, input);
546
547 VectorWithOffset<float> scale_factor(this->_patlak_plot_sptr->get_starting_frame(),
548 this->_patlak_plot_sptr->get_ending_frame());
549 for (unsigned int frame_num = this->_patlak_plot_sptr->get_starting_frame();
550 frame_num <= this->_patlak_plot_sptr->get_ending_frame();
551 ++frame_num)
552 {
553 assert(dyn_input[frame_num].find_max() == dyn_input[frame_num].find_min());
554 if (dyn_input[frame_num].find_max() == dyn_input[frame_num].find_min() && dyn_input[frame_num].find_min() > 0.F)
555 scale_factor[frame_num] = dyn_input[frame_num].find_max();
556 else
557 error("The input image should be uniform even after multiplying with the Patlak Plot.\n");
558
559 /*! /note This is used to avoid higher values than these set in the precompute_denominator_of_conditioner_without_penalty()
560 function. /sa for more information see the recon_array_functions.cxx and the value of the max_quotient (originaly set to
561 10000.F)
562 */
563 dyn_input[frame_num] /= scale_factor[frame_num];
564#ifndef NDEBUG
565 info(format("scale factor[{}]: {}", frame_num, scale_factor[frame_num]));
566 info(format("dyn_input[{}] max after scale: {}", frame_num, dyn_input[frame_num].find_max()));
567#endif // NDEBUG
568 this->_single_frame_obj_funcs[frame_num].add_multiplication_with_approximate_sub_Hessian_without_penalty(
569 dyn_output[frame_num], dyn_input[frame_num], subset_num);
570#ifndef NDEBUG
571 info(format("dyn_output[{}] max before scale: {}", frame_num, dyn_output[frame_num].find_max()));
572#endif // NDEBUG
573 dyn_output[frame_num] *= scale_factor[frame_num];
574#ifndef NDEBUG
575 info(format("dyn_output[{}] max after scale: {}", frame_num, dyn_output[frame_num].find_max()));
576#endif // NDEBUG
577 } // end of loop over frames
578 shared_ptr<TargetT> unnormalised_temp(output.get_empty_copy());
579 shared_ptr<TargetT> temp(output.get_empty_copy());
580 this->_patlak_plot_sptr->multiply_dynamic_image_with_model_gradient(*unnormalised_temp, dyn_output);
581 // Trick to use a better step size for the two parameters.
582 (this->_patlak_plot_sptr->get_model_matrix()).normalise_parametric_image_with_model_sum(*temp, *unnormalised_temp);
583#ifndef NDEBUG
584 info(format("TEMP max: ({} , {})", temp->construct_single_density(1).find_max(), temp->construct_single_density(2).find_max()));
585 // Writing images
586 write_to_file("all_params_one_input.img", input);
587 write_to_file("temp_denominator.img", *temp);
588 dyn_input.write_to_ecat7("dynamic_input_from_all_params_one.img");
589 dyn_output.write_to_ecat7("dynamic_precomputed_denominator.img");
590 DynamicProjData temp_projdata = this->get_dyn_proj_data();
591 for (unsigned int frame_num = this->_patlak_plot_sptr->get_starting_frame();
592 frame_num <= this->_patlak_plot_sptr->get_ending_frame();
593 ++frame_num)
594 temp_projdata.set_proj_data_sptr(this->_single_frame_obj_funcs[frame_num].get_proj_data_sptr(), frame_num);
595
596 temp_projdata.write_to_ecat7("DynamicProjections.S");
597#endif // NDEBUG
598 // output += temp
599 typename TargetT::full_iterator out_iter = output.begin_all();
600 typename TargetT::full_iterator out_end = output.end_all();
601 typename TargetT::const_full_iterator temp_iter = temp->begin_all_const();
602 while (out_iter != out_end)
603 {
604 *out_iter += *temp_iter;
605 ++out_iter;
606 ++temp_iter;
607 }
608#ifndef NDEBUG
609 info(format(
610 "OUTPUT max: ({} , {})", output.construct_single_density(1).find_max(), output.construct_single_density(2).find_max()));
611#endif // NDEBUG
612
613 return Succeeded::yes;
614}
615
616template <typename TargetT>
617Succeeded
618PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData<
619 TargetT>::actual_accumulate_sub_Hessian_times_input_without_penalty(TargetT& output,
620 const TargetT& current_image_estimate,
621 const TargetT& input,
622 const int subset_num) const
623{
624 { // check argument characteristics
625 std::string explanation;
626 if (!input.has_same_characteristics(this->get_sensitivity(), explanation))
627 {
628 warning("PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData:\n"
629 "sensitivity and input for actual_accumulate_sub_Hessian_times_input_without_penalty\n"
630 "should have the same characteristics.\n%s",
631 explanation.c_str());
632 return Succeeded::no;
633 }
634
635 if (!current_image_estimate.has_same_characteristics(this->get_sensitivity(), explanation))
636 {
637 warning("PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData:\n"
638 "sensitivity and current_image_estimate for actual_accumulate_sub_Hessian_times_input_without_penalty\n"
639 "should have the same characteristics.\n%s",
640 explanation.c_str());
641 return Succeeded::no;
642 }
643 }
644#ifndef NDEBUG
645 info(
646 format("INPUT max: ({} , {})", input.construct_single_density(1).find_max(), input.construct_single_density(2).find_max()));
647#endif // NDEBUG
648 DynamicDiscretisedDensity dyn_input = this->_dyn_image_template;
649 DynamicDiscretisedDensity dyn_current_image_estimate = this->_dyn_image_template;
650 DynamicDiscretisedDensity dyn_output = this->_dyn_image_template;
651 this->_patlak_plot_sptr->get_dynamic_image_from_parametric_image(dyn_input, input);
652 this->_patlak_plot_sptr->get_dynamic_image_from_parametric_image(dyn_current_image_estimate, current_image_estimate);
653
654 for (unsigned int frame_num = this->_patlak_plot_sptr->get_starting_frame();
655 frame_num <= this->_patlak_plot_sptr->get_ending_frame();
656 ++frame_num)
657 {
658 assert(dyn_input[frame_num].find_max() == dyn_input[frame_num].find_min());
659 this->_single_frame_obj_funcs[frame_num].accumulate_sub_Hessian_times_input_without_penalty(
660 dyn_output[frame_num], dyn_current_image_estimate[frame_num], dyn_input[frame_num], subset_num);
661 } // end of loop over frames
662 shared_ptr<TargetT> unnormalised_temp(output.get_empty_copy());
663 shared_ptr<TargetT> temp(output.get_empty_copy());
664 this->_patlak_plot_sptr->multiply_dynamic_image_with_model_gradient(*unnormalised_temp, dyn_output);
665 // Trick to use a better step size for the two parameters.
666 (this->_patlak_plot_sptr->get_model_matrix()).normalise_parametric_image_with_model_sum(*temp, *unnormalised_temp);
667#ifndef NDEBUG
668 info(format("TEMP max: ({} , {})", temp->construct_single_density(1).find_max(), temp->construct_single_density(2).find_max()));
669 // Writing images
670 write_to_file("all_params_one_input.img", input);
671 write_to_file("temp_denominator.img", *temp);
672 dyn_input.write_to_ecat7("dynamic_input_from_all_params_one.img");
673 dyn_output.write_to_ecat7("dynamic_precomputed_denominator.img");
674 DynamicProjData temp_projdata = this->get_dyn_proj_data();
675 for (unsigned int frame_num = this->_patlak_plot_sptr->get_starting_frame();
676 frame_num <= this->_patlak_plot_sptr->get_ending_frame();
677 ++frame_num)
678 temp_projdata.set_proj_data_sptr(this->_single_frame_obj_funcs[frame_num].get_proj_data_sptr(), frame_num);
679
680 temp_projdata.write_to_ecat7("DynamicProjections.S");
681#endif // NDEBUG
682 // output += temp
683 typename TargetT::full_iterator out_iter = output.begin_all();
684 typename TargetT::full_iterator out_end = output.end_all();
685 typename TargetT::const_full_iterator temp_iter = temp->begin_all_const();
686 while (out_iter != out_end)
687 {
688 *out_iter += *temp_iter;
689 ++out_iter;
690 ++temp_iter;
691 }
692#ifndef NDEBUG
693 info(format(
694 "OUTPUT max: ({} , {})", output.construct_single_density(1).find_max(), output.construct_single_density(2).find_max()));
695#endif // NDEBUG
696
697 return Succeeded::yes;
698}
699
700END_NAMESPACE_STIR
Generated on Wed May 27 2026 17:04:05 for STIR by doxygen 1.13.2