44#ifndef __stir_ML_norm_H__
45#define __stir_ML_norm_H__
55typedef Array<2, float> GeoData;
56typedef Array<2, float> BlockData;
58class DetPairData :
private Array<2, float>
62 DetPairData(
const IndexRange<2>& range);
63 DetPairData& operator=(
const DetPairData&);
65 float& operator()(
const int a,
const int b);
67 float operator()(
const int a,
const int b)
const;
69 void fill(
const float d);
71 int get_min_index()
const;
72 int get_max_index()
const;
73 int get_min_index(
const int a)
const;
74 int get_max_index(
const int a)
const;
75 bool is_in_data(
const int a,
const int b)
const;
77 float sum(
const int a)
const;
78 float find_max()
const;
79 float find_min()
const;
80 int get_num_detectors()
const;
81 void grow(
const IndexRange<2>&)
override;
84 typedef Array<2, float> base_type;
89get_fan_info(
int& num_rings,
int& num_detectors_per_ring,
int& max_ring_diff,
int& fan_size,
const ProjDataInfo& proj_data_info);
93make_det_pair_data(DetPairData& det_pair_data,
const ProjDataInfo& proj_data_info,
const int segment_num,
const int ax_pos_num);
95void make_det_pair_data(DetPairData& det_pair_data,
const ProjData& proj_data,
const int segment_num,
const int ax_pos_num);
97void set_det_pair_data(ProjData& proj_data,
const DetPairData& det_pair_data,
const int segment_num,
const int ax_pos_num);
98void apply_block_norm(DetPairData& det_pair_data,
const BlockData& geo_data,
const bool apply =
true);
99void apply_geo_norm(DetPairData& det_pair_data,
const GeoData& geo_data,
const bool apply =
true);
100void apply_efficiencies(DetPairData& det_pair_data,
const Array<1, float>& efficiencies,
const bool apply =
true);
102void make_fan_sum_data(Array<1, float>& data_fan_sums,
const DetPairData& det_pair_data);
104void make_geo_data(GeoData& geo_data,
const DetPairData& det_pair_data);
106void make_block_data(BlockData& block_data,
const DetPairData& det_pair_data);
108void iterate_efficiencies(Array<1, float>& efficiencies,
const Array<1, float>& data_fan_sums,
const DetPairData& model);
110void iterate_geo_norm(GeoData& norm_geo_data,
const GeoData& measured_geo_data,
const DetPairData& model);
112void iterate_block_norm(BlockData& norm_block_data,
const BlockData& measured_block_data,
const DetPairData& model);
115typedef Array<2, float> DetectorEfficiencies;
117class GeoData3D :
public Array<4, float>
121 GeoData3D(
const int num_axial_crystals_per_block,
122 const int half_num_transaxial_crystals_per_block,
124 const int num_detectors_per_ring);
125 ~GeoData3D()
override;
126 GeoData3D& operator=(
const GeoData3D&);
128 float& operator()(
const int ra,
const int a,
const int rb,
const int b);
130 float operator()(
const int ra,
const int a,
const int rb,
const int b)
const;
132 void fill(
const float d);
134 int get_min_ra()
const;
135 int get_max_ra()
const;
136 int get_min_a()
const;
137 int get_max_a()
const;
138 int get_min_b(
const int a)
const;
139 int get_max_b(
const int a)
const;
140 int get_min_rb(
const int ra)
const;
141 int get_max_rb(
const int ra)
const;
142 bool is_in_data(
const int ra,
const int a,
const int rb,
const int b)
const;
144 float sum(
const int ra,
const int a)
const;
145 float find_max()
const;
146 float find_min()
const;
147 int get_num_axial_crystals_per_block()
const;
148 int get_half_num_transaxial_crystals_per_block()
const;
151 friend std::ostream& operator<<(std::ostream&,
const GeoData3D&);
152 friend std::istream& operator>>(std::istream&, GeoData3D&);
153 typedef Array<4, float> base_type;
154 int num_axial_crystals_per_block;
155 int half_num_transaxial_crystals_per_block;
157 int num_detectors_per_ring;
160class FanProjData :
private Array<4, float>
164 FanProjData(
const int num_rings,
const int num_detectors_per_ring,
const int max_ring_diff,
const int fan_size);
165 ~FanProjData()
override;
166 FanProjData& operator=(
const FanProjData&);
168 float& operator()(
const int ra,
const int a,
const int rb,
const int b);
170 float operator()(
const int ra,
const int a,
const int rb,
const int b)
const;
172 void fill(
const float d);
174 int get_min_ra()
const;
175 int get_max_ra()
const;
176 int get_min_a()
const;
177 int get_max_a()
const;
178 int get_min_b(
const int a)
const;
179 int get_max_b(
const int a)
const;
180 int get_min_rb(
const int ra)
const;
181 int get_max_rb(
const int ra)
const;
182 bool is_in_data(
const int ra,
const int a,
const int rb,
const int b)
const;
184 float sum(
const int ra,
const int a)
const;
185 float find_max()
const;
186 float find_min()
const;
187 int get_max_delta()
const;
188 int get_num_detectors_per_ring()
const;
189 int get_num_rings()
const;
192 friend std::ostream& operator<<(std::ostream&,
const FanProjData&);
193 friend std::istream& operator>>(std::istream&, FanProjData&);
194 typedef Array<4, float> base_type;
198 int num_detectors_per_ring;
203typedef FanProjData BlockData3D;
205void make_fan_data_remove_gaps(FanProjData& fan_data,
const ProjData& proj_data);
207void set_fan_data_add_gaps(ProjData& proj_data,
const FanProjData& fan_data,
const float gap_value = 0.F);
209void apply_block_norm(FanProjData& fan_data,
const BlockData3D& block_data,
const bool apply =
true);
211void apply_geo_norm(FanProjData& fan_data,
const GeoData3D& geo_data,
const bool apply =
true);
213void apply_efficiencies(FanProjData& fan_data,
const DetectorEfficiencies& efficiencies,
const bool apply =
true);
215void make_fan_sum_data(Array<2, float>& data_fan_sums,
const FanProjData& fan_data);
217void make_fan_sum_data(Array<2, float>& data_fan_sums,
const ProjData& proj_data);
219void make_fan_sum_data(Array<2, float>& data_fan_sums,
220 const DetectorEfficiencies& efficiencies,
221 const int max_ring_diff,
222 const int half_fan_size);
224void make_geo_data(GeoData3D& geo_data,
const FanProjData& fan_data);
226void make_block_data(BlockData3D& block_data,
const FanProjData& fan_data);
228void iterate_efficiencies(DetectorEfficiencies& efficiencies,
const Array<2, float>& data_fan_sums,
const FanProjData& model);
231void iterate_efficiencies(DetectorEfficiencies& efficiencies,
232 const Array<2, float>& data_fan_sums,
233 const int max_ring_diff,
234 const int half_fan_size);
236void iterate_geo_norm(GeoData3D& geo_data,
const GeoData3D& measured_geo_data,
const FanProjData& model);
238void iterate_block_norm(BlockData3D& norm_block_data,
const BlockData3D& measured_block_data,
const FanProjData& model);
241KL(
const double a,
const double b,
const double threshold_a = 0)
245 double res = a <= threshold_a ? b : (a * (log(a) - log(b)) + b - a);
248 warning(
"KL nan at a=%g b=%g, threshold %g\n", a, b, threshold_a);
250 warning(
"KL large at a=%g b=%g, threshold %g\n", a, b, threshold_a);
252 assert(res >= -1.e-4);
256template <
int num_dimensions,
typename elemT>
258KL(
const Array<num_dimensions, elemT>& a,
const Array<num_dimensions, elemT>& b,
const double threshold_a = 0)
261 typename Array<num_dimensions, elemT>::const_full_iterator iter_a = a.begin_all();
262 typename Array<num_dimensions, elemT>::const_full_iterator iter_b = b.begin_all();
263 while (iter_a != a.end_all())
265 sum +=
static_cast<double>(KL(*iter_a++, *iter_b++, threshold_a));
267 return static_cast<double>(
sum);
270double KL(
const DetPairData& d1,
const DetPairData& d2,
const double threshold);
272double 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::ProjData.
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:175
Declaration of stir::warning().
1.18.0