44#ifndef __stir_ML_norm_H__
45#define __stir_ML_norm_H__
58typedef Array<2, float> GeoData;
59typedef Array<2, float> BlockData;
61class DetPairData :
private Array<2, float>
65 DetPairData(
const IndexRange<2>& range);
66 DetPairData& operator=(
const DetPairData&);
68 float& operator()(
const int a,
const int b);
70 float operator()(
const int a,
const int b)
const;
72 void fill(
const float d);
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;
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;
87 typedef Array<2, float> base_type;
92get_fan_info(
int& num_rings,
int& num_detectors_per_ring,
int& max_ring_diff,
int& fan_size,
const ProjDataInfo& proj_data_info);
96make_det_pair_data(DetPairData& det_pair_data,
const ProjDataInfo& proj_data_info,
const int segment_num,
const int ax_pos_num);
98void make_det_pair_data(DetPairData& det_pair_data,
const ProjData& proj_data,
const int segment_num,
const int ax_pos_num);
100void set_det_pair_data(ProjData& proj_data,
const DetPairData& det_pair_data,
const int segment_num,
const int ax_pos_num);
101void apply_block_norm(DetPairData& det_pair_data,
const BlockData& geo_data,
const bool apply =
true);
102void apply_geo_norm(DetPairData& det_pair_data,
const GeoData& geo_data,
const bool apply =
true);
103void apply_efficiencies(DetPairData& det_pair_data,
const Array<1, float>& efficiencies,
const bool apply =
true);
105void make_fan_sum_data(Array<1, float>& data_fan_sums,
const DetPairData& det_pair_data);
107void make_geo_data(GeoData& geo_data,
const DetPairData& det_pair_data);
109void make_block_data(BlockData& block_data,
const DetPairData& det_pair_data);
111void iterate_efficiencies(Array<1, float>& efficiencies,
const Array<1, float>& data_fan_sums,
const DetPairData& model);
113void iterate_geo_norm(GeoData& norm_geo_data,
const GeoData& measured_geo_data,
const DetPairData& model);
115void iterate_block_norm(BlockData& norm_block_data,
const BlockData& measured_block_data,
const DetPairData& model);
118typedef Array<2, float> DetectorEfficiencies;
120class GeoData3D :
public Array<4, float>
124 GeoData3D(
const int num_axial_crystals_per_block,
125 const int half_num_transaxial_crystals_per_block,
127 const int num_detectors_per_ring);
128 ~GeoData3D()
override;
129 GeoData3D& operator=(
const GeoData3D&);
131 float& operator()(
const int ra,
const int a,
const int rb,
const int b);
133 float operator()(
const int ra,
const int a,
const int rb,
const int b)
const;
135 void fill(
const float d);
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;
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;
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;
160 int num_detectors_per_ring;
163class FanProjData :
private Array<4, float>
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&);
171 float& operator()(
const int ra,
const int a,
const int rb,
const int b);
173 float operator()(
const int ra,
const int a,
const int rb,
const int b)
const;
175 void fill(
const float d);
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;
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;
195 friend std::ostream& operator<<(std::ostream&,
const FanProjData&);
196 friend std::istream& operator>>(std::istream&, FanProjData&);
197 typedef Array<4, float> base_type;
201 int num_detectors_per_ring;
206typedef FanProjData BlockData3D;
208void make_fan_data_remove_gaps(FanProjData& fan_data,
const ProjData& proj_data);
210void set_fan_data_add_gaps(ProjData& proj_data,
const FanProjData& fan_data,
const float gap_value = 0.F);
212void apply_block_norm(FanProjData& fan_data,
const BlockData3D& block_data,
const bool apply =
true);
214void apply_geo_norm(FanProjData& fan_data,
const GeoData3D& geo_data,
const bool apply =
true);
216void apply_efficiencies(FanProjData& fan_data,
const DetectorEfficiencies& efficiencies,
const bool apply =
true);
218void make_fan_sum_data(Array<2, float>& data_fan_sums,
const FanProjData& fan_data);
220void make_fan_sum_data(Array<2, float>& data_fan_sums,
const ProjData& proj_data);
222void 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);
227void make_geo_data(GeoData3D& geo_data,
const FanProjData& fan_data);
229void make_block_data(BlockData3D& block_data,
const FanProjData& fan_data);
231void iterate_efficiencies(DetectorEfficiencies& efficiencies,
const Array<2, float>& data_fan_sums,
const FanProjData& model);
234void iterate_efficiencies(DetectorEfficiencies& efficiencies,
235 const Array<2, float>& data_fan_sums,
236 const int max_ring_diff,
237 const int half_fan_size);
239void iterate_geo_norm(GeoData3D& geo_data,
const GeoData3D& measured_geo_data,
const FanProjData& model);
241void iterate_block_norm(BlockData3D& norm_block_data,
const BlockData3D& measured_block_data,
const FanProjData& model);
244KL(
const double a,
const double b,
const double threshold_a = 0)
248 double res = a <= threshold_a ? b : (a * (log(a) - log(b)) + b - a);
251 warning(
"KL nan at a=%g b=%g, threshold %g\n", a, b, threshold_a);
253 warning(
"KL large at a=%g b=%g, threshold %g\n", a, b, threshold_a);
255 assert(res >= -1.e-4);
259template <
int num_dimensions,
typename elemT>
261KL(
const Array<num_dimensions, elemT>& a,
const Array<num_dimensions, elemT>& b,
const double threshold_a = 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())
268 sum +=
static_cast<double>(KL(*iter_a++, *iter_b++, threshold_a));
270 return static_cast<double>(
sum);
273double KL(
const DetPairData& d1,
const DetPairData& d2,
const double threshold);
275double KL(
const FanProjData& d1,
const FanProjData& d2,
const double threshold);
defines the stir::Array class for multi-dimensional (numeric) arrays
This file declares the class stir::IndexRange2D.
Declaration of class stir::ProjDataInfoBlocksOnCylindricalNoArcCorr.
Declaration of class stir::ProjDataInfoCylindricalNoArcCorr.
Declaration of class stir::ProjData.
Declaration of class stir::Sinogram.
void warning(const char *const s,...)
Print warning with format string a la printf.
Definition warning.cxx:41
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
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
Declaration of stir::warning()
1.13.2