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>
26 #include <mia/core/splinekernel.hh>
27 #include <mia/core/histogram.hh>
28 
30 
45 {
46 public:
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();
125 protected:
127  void fill_histograms(const std::vector<double>& values,
128  double rmin, double rmax,
129  double mmin, double mmax);
130 private:
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 
174 template <typename MovIterator, typename RefIterator>
175 BOOST_CONCEPT_REQUIRES( ((::boost::ForwardIterator<MovIterator>))
176  ((::boost::ForwardIterator<RefIterator>)),
177  (void)
178  )
179 CSplineParzenMI::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 
248 template <typename MovIterator, typename RefIterator, typename MaskIterator>
249 void 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 
326 template <typename Iterator>
327 std::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 
340 template <typename DataIterator, typename MaskIterator>
341 std::pair<double, double>
342 CSplineParzenMI::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
cvinfo
vstream & cvinfo()
informal output that may be of interest to understand problems with a program and are of higher prior...
Definition: msgstream.hh:262
NS_MIA_BEGIN
#define NS_MIA_BEGIN
conveniance define to start the mia namespace
Definition: defines.hh:33
get_gradient
EXPORT_2D C2DFVectorfield get_gradient(const C2DImage &image)
histogram.hh
THistogram
a simple histogram that uses an instance of THistogramFeeder as input converter
Definition: histogram.hh:134
NS_MIA_END
#define NS_MIA_END
conveniance define to end the mia namespace
Definition: defines.hh:36
CSplineParzenMI
Implementation of mutual information based on B-splines.
Definition: splineparzenmi.hh:44
cvdebug
CDebugSink & cvdebug()
Definition: msgstream.hh:226
PSplineKernel
std::shared_ptr< CSplineKernel > PSplineKernel
Definition: splinekernel.hh:291
CSplineParzenMI::fill
void fill(MovIterator mov_begin, MovIterator mov_end, RefIterator ref_begin, RefIterator ref_end)
Definition: splineparzenmi.hh:179
EXPORT_CORE
#define EXPORT_CORE
Macro to manage Visual C++ style dllimport/dllexport.
Definition: defines.hh:101
THistogramFeeder
A class to normalize and quantizize input data to a given histogram range with its given number of bi...
Definition: histogram.hh:49
splinekernel.hh