STIR  6.2.0
ML_norm.h
Go to the documentation of this file.
1 //
2 //
3 /*
4  This file is part of STIR.
5 
6  Copyright (C) 2001- 2011, Hammersmith Imanet Ltd
7  Copyright (C) 2020, University College London
8  Copyright (C) 2016-2017, PETsys Electronics
9  Copyright (C) 2022, National Physical Laboratory
10  This file is part of STIR.
11 
12  SPDX-License-Identifier: Apache-2.0
13 
14  See STIR/LICENSE.txt for details
15 */
44 #ifndef __stir_ML_norm_H__
45 #define __stir_ML_norm_H__
46 
47 #include "stir/ProjData.h"
48 #include "stir/Array.h"
49 #include "stir/ProjDataInfoCylindricalNoArcCorr.h"
50 #include "stir/ProjDataInfoBlocksOnCylindricalNoArcCorr.h"
51 #include "stir/IndexRange2D.h"
52 #include "stir/Sinogram.h"
53 #include "stir/warning.h"
54 #include <iostream>
55 
56 START_NAMESPACE_STIR
57 
58 typedef Array<2, float> GeoData;
59 typedef Array<2, float> BlockData;
60 
61 class DetPairData : private Array<2, float>
62 {
63 public:
64  DetPairData();
65  DetPairData(const IndexRange<2>& range);
66  DetPairData& operator=(const DetPairData&);
67 
68  float& operator()(const int a, const int b);
69 
70  float operator()(const int a, const int b) const;
71 
72  void fill(const float d);
73 
74  int get_min_index() const;
75  int get_max_index() const;
76  int get_min_index(const int a) const;
77  int get_max_index(const int a) const;
78  bool is_in_data(const int a, const int b) const;
79  float sum() const;
80  float sum(const int a) const;
81  float find_max() const;
82  float find_min() const;
83  int get_num_detectors() const;
84  void grow(const IndexRange<2>&) override;
85 
86 private:
87  typedef Array<2, float> base_type;
88  int num_detectors;
89 };
90 
91 void
92 get_fan_info(int& num_rings, int& num_detectors_per_ring, int& max_ring_diff, int& fan_size, const ProjDataInfo& proj_data_info);
93 
95 void
96 make_det_pair_data(DetPairData& det_pair_data, const ProjDataInfo& proj_data_info, const int segment_num, const int ax_pos_num);
97 
98 void make_det_pair_data(DetPairData& det_pair_data, const ProjData& proj_data, const int segment_num, const int ax_pos_num);
99 
100 void set_det_pair_data(ProjData& proj_data, const DetPairData& det_pair_data, const int segment_num, const int ax_pos_num);
101 void apply_block_norm(DetPairData& det_pair_data, const BlockData& geo_data, const bool apply = true);
102 void apply_geo_norm(DetPairData& det_pair_data, const GeoData& geo_data, const bool apply = true);
103 void apply_efficiencies(DetPairData& det_pair_data, const Array<1, float>& efficiencies, const bool apply = true);
104 
105 void make_fan_sum_data(Array<1, float>& data_fan_sums, const DetPairData& det_pair_data);
106 
107 void make_geo_data(GeoData& geo_data, const DetPairData& det_pair_data);
108 
109 void make_block_data(BlockData& block_data, const DetPairData& det_pair_data);
110 
111 void iterate_efficiencies(Array<1, float>& efficiencies, const Array<1, float>& data_fan_sums, const DetPairData& model);
112 
113 void iterate_geo_norm(GeoData& norm_geo_data, const GeoData& measured_geo_data, const DetPairData& model);
114 
115 void iterate_block_norm(BlockData& norm_block_data, const BlockData& measured_block_data, const DetPairData& model);
116 
117 //******3 D
118 typedef Array<2, float> DetectorEfficiencies;
119 
120 class GeoData3D : public Array<4, float>
121 {
122 public:
123  GeoData3D();
124  GeoData3D(const int num_axial_crystals_per_block,
125  const int half_num_transaxial_crystals_per_block,
126  const int num_rings,
127  const int num_detectors_per_ring);
128  ~GeoData3D() override;
129  GeoData3D& operator=(const GeoData3D&);
130 
131  float& operator()(const int ra, const int a, const int rb, const int b);
132 
133  float operator()(const int ra, const int a, const int rb, const int b) const;
134 
135  void fill(const float d);
136 
137  int get_min_ra() const;
138  int get_max_ra() const;
139  int get_min_a() const;
140  int get_max_a() const;
141  int get_min_b(const int a) const;
142  int get_max_b(const int a) const;
143  int get_min_rb(const int ra) const;
144  int get_max_rb(const int ra) const;
145  bool is_in_data(const int ra, const int a, const int rb, const int b) const;
146  float sum() const;
147  float sum(const int ra, const int a) const;
148  float find_max() const;
149  float find_min() const;
150  int get_num_axial_crystals_per_block() const;
151  int get_half_num_transaxial_crystals_per_block() const;
152 
153 private:
154  friend std::ostream& operator<<(std::ostream&, const GeoData3D&);
155  friend std::istream& operator>>(std::istream&, GeoData3D&);
156  typedef Array<4, float> base_type;
157  int num_axial_crystals_per_block;
158  int half_num_transaxial_crystals_per_block;
159  int num_rings;
160  int num_detectors_per_ring;
161 };
162 
163 class FanProjData : private Array<4, float>
164 {
165 public:
166  FanProjData();
167  FanProjData(const int num_rings, const int num_detectors_per_ring, const int max_ring_diff, const int fan_size);
168  ~FanProjData() override;
169  FanProjData& operator=(const FanProjData&);
170 
171  float& operator()(const int ra, const int a, const int rb, const int b);
172 
173  float operator()(const int ra, const int a, const int rb, const int b) const;
174 
175  void fill(const float d);
176 
177  int get_min_ra() const;
178  int get_max_ra() const;
179  int get_min_a() const;
180  int get_max_a() const;
181  int get_min_b(const int a) const;
182  int get_max_b(const int a) const;
183  int get_min_rb(const int ra) const;
184  int get_max_rb(const int ra) const;
185  bool is_in_data(const int ra, const int a, const int rb, const int b) const;
186  float sum() const;
187  float sum(const int ra, const int a) const;
188  float find_max() const;
189  float find_min() const;
190  int get_max_delta() const;
191  int get_num_detectors_per_ring() const;
192  int get_num_rings() const;
193 
194 private:
195  friend std::ostream& operator<<(std::ostream&, const FanProjData&);
196  friend std::istream& operator>>(std::istream&, FanProjData&);
197  typedef Array<4, float> base_type;
198  // FanProjData(const IndexRange<4>& range);
199  // void grow(const IndexRange<4>&);
200  int num_rings;
201  int num_detectors_per_ring;
202  int max_ring_diff;
203  int half_fan_size;
204 };
205 
206 typedef FanProjData BlockData3D;
207 
208 void make_fan_data_remove_gaps(FanProjData& fan_data, const ProjData& proj_data);
209 
210 void set_fan_data_add_gaps(ProjData& proj_data, const FanProjData& fan_data, const float gap_value = 0.F);
211 
212 void apply_block_norm(FanProjData& fan_data, const BlockData3D& block_data, const bool apply = true);
213 
214 void apply_geo_norm(FanProjData& fan_data, const GeoData3D& geo_data, const bool apply = true);
215 
216 void apply_efficiencies(FanProjData& fan_data, const DetectorEfficiencies& efficiencies, const bool apply = true);
217 
218 void make_fan_sum_data(Array<2, float>& data_fan_sums, const FanProjData& fan_data);
219 
220 void make_fan_sum_data(Array<2, float>& data_fan_sums, const ProjData& proj_data);
221 
222 void make_fan_sum_data(Array<2, float>& data_fan_sums,
223  const DetectorEfficiencies& efficiencies,
224  const int max_ring_diff,
225  const int half_fan_size);
226 
227 void make_geo_data(GeoData3D& geo_data, const FanProjData& fan_data);
228 
229 void make_block_data(BlockData3D& block_data, const FanProjData& fan_data);
230 
231 void iterate_efficiencies(DetectorEfficiencies& efficiencies, const Array<2, float>& data_fan_sums, const FanProjData& model);
232 
233 // version without model
234 void iterate_efficiencies(DetectorEfficiencies& efficiencies,
235  const Array<2, float>& data_fan_sums,
236  const int max_ring_diff,
237  const int half_fan_size);
238 
239 void iterate_geo_norm(GeoData3D& geo_data, const GeoData3D& measured_geo_data, const FanProjData& model);
240 
241 void iterate_block_norm(BlockData3D& norm_block_data, const BlockData3D& measured_block_data, const FanProjData& model);
242 
243 inline double
244 KL(const double a, const double b, const double threshold_a = 0)
245 {
246  assert(a >= 0);
247  assert(b >= 0);
248  double res = a <= threshold_a ? b : (a * (log(a) - log(b)) + b - a);
249 #ifndef NDEBUG
250  if (res != res)
251  warning("KL nan at a=%g b=%g, threshold %g\n", a, b, threshold_a);
252  if (res > 1.E20)
253  warning("KL large at a=%g b=%g, threshold %g\n", a, b, threshold_a);
254 #endif
255  assert(res >= -1.e-4);
256  return res;
257 }
258 
259 template <int num_dimensions, typename elemT>
260 double
261 KL(const Array<num_dimensions, elemT>& a, const Array<num_dimensions, elemT>& b, const double threshold_a = 0)
262 {
263  double sum = 0;
264  typename Array<num_dimensions, elemT>::const_full_iterator iter_a = a.begin_all();
265  typename Array<num_dimensions, elemT>::const_full_iterator iter_b = b.begin_all();
266  while (iter_a != a.end_all())
267  {
268  sum += static_cast<double>(KL(*iter_a++, *iter_b++, threshold_a));
269  }
270  return static_cast<double>(sum);
271 }
272 
273 double KL(const DetPairData& d1, const DetPairData& d2, const double threshold);
274 
275 double KL(const FanProjData& d1, const FanProjData& d2, const double threshold);
276 
277 // double KL(const GeoData3D& d1, const GeoData3D& d2, const double threshold);
278 
279 END_NAMESPACE_STIR
280 #endif
Sinogram.h
Declaration of class stir::Sinogram.
ProjDataInfoCylindricalNoArcCorr.h
Declaration of class stir::ProjDataInfoCylindricalNoArcCorr.
ProjData.h
Declaration of class stir::ProjData.
Array.h
defines the Array class for multi-dimensional (numeric) arrays
stir::make_det_pair_data
void make_det_pair_data(DetPairData &det_pair_data, const ProjDataInfo &proj_data_info_general_type, const int segment_num, const int ax_pos_num)
Makes a DetPairData of appropriate dimensions and fills it with 0.
Definition: ML_norm.cxx:173
stir::warning
void warning(const char *const s,...)
Print warning with format string a la printf.
Definition: warning.cxx:41
warning.h
Declaration of stir::warning()
ProjDataInfoBlocksOnCylindricalNoArcCorr.h
Declaration of class stir::ProjDataInfoBlocksOnCylindricalNoArcCorr.
IndexRange2D.h
This file declares the class stir::IndexRange2D.
stir::sum
elemT sum(IterT start, IterT end, elemT init)
Compute the sum of a sequence using operator+=(), using an initial value.
Definition: more_algorithms.inl:52
Generated on Tue Jul 23 2024 09:13:10 for STIR by   doxygen 1.8.13