My Project
splineparzenmi.hh
Go to the documentation of this file.
1/* -*- mia-c++ -*-
2 *
3 * This file is part of MIA - a toolbox for medical image analysis
4 * Copyright (c) Leipzig, Madrid 1999-2017 Gert Wollny
5 *
6 * MIA is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with MIA; if not, see <http://www.gnu.org/licenses/>.
18 *
19 */
20
21#ifndef mia_core_splineparzenmi_hh
22#define mia_core_splineparzenmi_hh
23
24#include <boost/concept/requires.hpp>
25#include <boost/concept_check.hpp>
27#include <mia/core/histogram.hh>
28
30
45{
46public:
56 CSplineParzenMI(size_t rbins, PSplineKernel rkernel,
57 size_t mbins, PSplineKernel mkernel, double cut_histogram);
58
59
70 template <typename MovIterator, typename RefIterator>
71 BOOST_CONCEPT_REQUIRES( ((::boost::ForwardIterator<MovIterator>))
72 ((::boost::ForwardIterator<RefIterator>)),
73 (void)
74 )
75 fill(MovIterator mov_begin, MovIterator mov_end,
76 RefIterator ref_begin, RefIterator ref_end);
77
78
92 template <typename MovIterator, typename RefIterator, typename MaskIterator>
93 void fill(MovIterator mov_begin, MovIterator mov_end,
94 RefIterator ref_begin, RefIterator ref_end,
95 MaskIterator mask_begin, MaskIterator mask_end);
96
97
102 double value() const;
103
111 double get_gradient(double moving, double reference) const;
112
120 double get_gradient_slow(double moving, double reference) const;
124 void reset();
125protected:
127 void fill_histograms(const std::vector<double>& values,
128 double rmin, double rmax,
129 double mmin, double mmax);
130private:
131
132 double scale_moving(double x) const;
133 double scale_reference(double x) const;
134
135 void evaluate_histograms();
136 void evaluate_log_cache();
137
138 size_t m_ref_bins;
139 PSplineKernel m_ref_kernel;
140 size_t m_ref_border;
141 size_t m_ref_real_bins;
142 double m_ref_max;
143 double m_ref_min;
144 double m_ref_scale;
145
146 size_t m_mov_bins;
147
148 PSplineKernel m_mov_kernel;
149 size_t m_mov_border;
150 size_t m_mov_real_bins;
151 double m_mov_max;
152 double m_mov_min;
153 double m_mov_scale;
154
155
156 std::vector<double> m_joined_histogram;
157 std::vector<double> m_ref_histogram;
158 std::vector<double> m_mov_histogram;
159
160 std::vector<std::vector<double>> m_pdfLogCache;
161 double m_cut_histogram;
162
163 template <typename Iterator>
164 std::pair<double, double> get_reduced_range(Iterator begin, Iterator end)const;
165
166 template <typename DataIterator, typename MaskIterator>
167 std::pair<double, double>
168 get_reduced_range_masked(DataIterator dbegin,
169 DataIterator dend,
170 MaskIterator mbegin)const;
171
172};
173
174template <typename MovIterator, typename RefIterator>
175BOOST_CONCEPT_REQUIRES( ((::boost::ForwardIterator<MovIterator>))
176 ((::boost::ForwardIterator<RefIterator>)),
177 (void)
178 )
179CSplineParzenMI::fill(MovIterator mov_begin, MovIterator mov_end,
180 RefIterator ref_begin, RefIterator ref_end)
181{
182 std::fill(m_joined_histogram.begin(), m_joined_histogram.end(), 0.0);
183 assert(mov_begin != mov_end);
184 assert(ref_begin != ref_end);
185
186 if (m_mov_max < m_mov_min) {
187 // (re)evaluate the ranges
188 auto mov_range = get_reduced_range(mov_begin, mov_end);
189
190 if (mov_range.second == mov_range.first)
191 throw std::invalid_argument("relevant moving image intensity range is zero");
192
193 m_mov_min = mov_range.first;
194 m_mov_max = mov_range.second;
195 m_mov_scale = (m_mov_bins - 1) / (m_mov_max - m_mov_min);
196 cvdebug() << "Mov Range = [" << m_mov_min << ", " << m_mov_max << "]\n";
197 }
198
199 if (m_ref_max < m_ref_min) {
200 auto ref_range = get_reduced_range(ref_begin, ref_end);
201
202 if (ref_range.second == ref_range.first)
203 throw std::invalid_argument("relevant reference image intensity range is zero");
204
205 m_ref_min = ref_range.first;
206 m_ref_max = ref_range.second;
207 m_ref_scale = (m_ref_bins - 1) / (m_ref_max - m_ref_min);
208 cvdebug() << "Ref Range = [" << m_ref_min << ", " << m_ref_max << "]\n";
209 }
210
211 std::vector<double> mweights(m_mov_kernel->size());
212 std::vector<double> rweights(m_ref_kernel->size());
213 size_t N = 0;
214
215 while (ref_begin != ref_end && mov_begin != mov_end) {
216 const double mov = scale_moving(*mov_begin);
217 const double ref = scale_reference(*ref_begin);
218 const int mov_start = m_mov_kernel->get_start_idx_and_value_weights(mov, mweights) + m_mov_border;
219 const int ref_start = m_ref_kernel->get_start_idx_and_value_weights(ref, rweights) + m_ref_border;
220
221 for (size_t r = 0; r < m_ref_kernel->size(); ++r) {
222 auto inbeg = m_joined_histogram.begin() +
223 m_mov_real_bins * (ref_start + r) + mov_start;
224 double rw = rweights[r];
225 std::transform(mweights.begin(), mweights.end(), inbeg, inbeg,
226 [rw](double mw, double jhvalue) {
227 return mw * rw + jhvalue;
228 });
229 }
230
231 ++N;
232 ++mov_begin;
233 ++ref_begin;
234 }
235
236 cvdebug() << "CSplineParzenMI::fill: counted " << N << " pixels\n";
237 // normalize joined histogram
238 const double nscale = 1.0 / N;
239 transform(m_joined_histogram.begin(), m_joined_histogram.end(), m_joined_histogram.begin(),
240 [&nscale](double jhvalue) {
241 return jhvalue * nscale;
242 });
243 evaluate_histograms();
244 evaluate_log_cache();
245}
246
247
248template <typename MovIterator, typename RefIterator, typename MaskIterator>
249void CSplineParzenMI::fill(MovIterator mov_begin, MovIterator mov_end,
250 RefIterator ref_begin, RefIterator ref_end,
251 MaskIterator mask_begin, MaskIterator mask_end)
252{
253 std::fill(m_joined_histogram.begin(), m_joined_histogram.end(), 0.0);
254 assert(mov_begin != mov_end);
255 assert(ref_begin != ref_end);
256
257 // no mask
258 if (mask_begin == mask_end)
259 fill(mov_begin, mov_end, ref_begin, ref_end);
260
261 if (m_mov_max < m_mov_min) {
262 // (re)evaluate the ranges
263 auto mov_range = get_reduced_range_masked(mov_begin, mov_end, mask_begin);
264
265 if (mov_range.second == mov_range.first)
266 throw std::invalid_argument("relevant moving image intensity range is zero");
267
268 m_mov_min = mov_range.first;
269 m_mov_max = mov_range.second;
270 m_mov_scale = (m_mov_bins - 1) / (m_mov_max - m_mov_min);
271 cvdebug() << "Mov Range = [" << m_mov_min << ", " << m_mov_max << "]\n";
272 }
273
274 if (m_ref_max < m_ref_min) {
275 auto ref_range = get_reduced_range_masked(ref_begin, ref_end, mask_begin);
276
277 if (ref_range.second == ref_range.first)
278 throw std::invalid_argument("relevant reference image intensity range is zero");
279
280 m_ref_min = ref_range.first;
281 m_ref_max = ref_range.second;
282 m_ref_scale = (m_ref_bins - 1) / (m_ref_max - m_ref_min);
283 cvdebug() << "Ref Range = [" << m_ref_min << ", " << m_ref_max << "]\n";
284 }
285
286 std::vector<double> mweights(m_mov_kernel->size());
287 std::vector<double> rweights(m_ref_kernel->size());
288 size_t N = 0;
289
290 while (ref_begin != ref_end && mov_begin != mov_end) {
291 if (*mask_begin) {
292 const double mov = scale_moving(*mov_begin);
293 const double ref = scale_reference(*ref_begin);
294 const int mov_start = m_mov_kernel->get_start_idx_and_value_weights(mov, mweights) + m_mov_border;
295 const int ref_start = m_ref_kernel->get_start_idx_and_value_weights(ref, rweights) + m_ref_border;
296
297 for (size_t r = 0; r < m_ref_kernel->size(); ++r) {
298 auto inbeg = m_joined_histogram.begin() +
299 m_mov_real_bins * (ref_start + r) + mov_start;
300 double rw = rweights[r];
301 std::transform(mweights.begin(), mweights.end(), inbeg, inbeg,
302 [rw](double mw, double jhvalue) {
303 return mw * rw + jhvalue;
304 });
305 }
306
307 ++N;
308 }
309
310 ++mask_begin;
311 ++mov_begin;
312 ++ref_begin;
313 }
314
315 cvdebug() << "CSplineParzenMI::fill: counted " << N << " pixels\n";
316 // normalize joined histogram
317 const double nscale = 1.0 / N;
318 transform(m_joined_histogram.begin(), m_joined_histogram.end(), m_joined_histogram.begin(),
319 [&nscale](double jhvalue) {
320 return jhvalue * nscale;
321 });
322 evaluate_histograms();
323 evaluate_log_cache();
324}
325
326template <typename Iterator>
327std::pair<double, double> CSplineParzenMI::get_reduced_range(Iterator begin, Iterator end)const
328{
329 auto range = std::minmax_element(begin, end);
331 THistogram<Feeder> h(Feeder(*range.first, *range.second, 4096));
332 h.push_range(begin, end);
333 auto reduced_range = h.get_reduced_range(m_cut_histogram);
334 cvinfo() << "CSplineParzenMI: reduce range by " << m_cut_histogram
335 << "% from [" << *range.first << ", " << *range.second
336 << "] to [" << reduced_range.first << ", " << reduced_range.second << "]\n";
337 return std::make_pair(reduced_range.first, reduced_range.second);
338}
339
340template <typename DataIterator, typename MaskIterator>
341std::pair<double, double>
342CSplineParzenMI::get_reduced_range_masked(DataIterator dbegin,
343 DataIterator dend,
344 MaskIterator mbegin)const
345{
346 auto ib = dbegin;
347 auto im = mbegin;
348
349 while (! *im && ib != dend) {
350 ++im;
351 ++ib;
352 }
353
354 if (ib == dend)
355 throw std::runtime_error("CSplineParzenMI: empty mask");
356
357 double range_max = *ib;
358 double range_min = *ib;
359 ++ib;
360 ++im;
361
362 while (ib != dend) {
363 if (*im) {
364 if (*ib < range_min)
365 range_min = *ib;
366
367 if (*ib > range_max)
368 range_max = *ib;
369 }
370
371 ++ib;
372 ++im;
373 }
374
376 THistogram<Feeder> h(Feeder(range_min, range_max, 4096));
377 ib = dbegin;
378 im = mbegin;
379
380 while (ib != dend) {
381 if (*im)
382 h.push(*ib);
383
384 ++ib;
385 ++im;
386 }
387
388 auto reduced_range = h.get_reduced_range(m_cut_histogram);
389 cvinfo() << "CSplineParzenMI: reduce range by " << m_cut_histogram
390 << "% from [" << range_min << ", " << range_max
391 << "] to [" << reduced_range.first << ", " << reduced_range.second << "]\n";
392 return std::make_pair(reduced_range.first, reduced_range.second);
393}
394
396#endif
Implementation of mutual information based on B-splines.
void fill(MovIterator mov_begin, MovIterator mov_end, RefIterator ref_begin, RefIterator ref_end)
double get_gradient_slow(double moving, double reference) const
double value() const
void fill_histograms(const std::vector< double > &values, double rmin, double rmax, double mmin, double mmax)
double get_gradient(double moving, double reference) const
CSplineParzenMI(size_t rbins, PSplineKernel rkernel, size_t mbins, PSplineKernel mkernel, double cut_histogram)
A class to normalize and quantizize input data to a given histogram range with its given number of bi...
Definition: histogram.hh:50
a simple histogram that uses an instance of THistogramFeeder as input converter
Definition: histogram.hh:135
#define NS_MIA_BEGIN
conveniance define to start the mia namespace
Definition: defines.hh:33
#define EXPORT_CORE
Macro to manage Visual C++ style dllimport/dllexport.
Definition: defines.hh:101
#define NS_MIA_END
conveniance define to end the mia namespace
Definition: defines.hh:36
std::shared_ptr< CSplineKernel > PSplineKernel
vstream & cvinfo()
informal output that may be of interest to understand problems with a program and are of higher prior...
Definition: msgstream.hh:262
CDebugSink & cvdebug()
Definition: msgstream.hh:226