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-2013 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 public:
55  CSplineParzenMI(size_t rbins, PSplineKernel rkernel,
56  size_t mbins, PSplineKernel mkernel, double cut_histogram);
57 
58 
69  template <typename MovIterator, typename RefIterator>
70  BOOST_CONCEPT_REQUIRES( ((::boost::ForwardIterator<MovIterator>))
71  ((::boost::ForwardIterator<RefIterator>)),
72  (void)
73  )
74  fill(MovIterator mov_begin, MovIterator mov_end,
75  RefIterator ref_begin, RefIterator ref_end);
76 
77 
82  double value() const;
83 
91  double get_gradient(double moving, double reference) const;
92 
100  double get_gradient_slow(double moving, double reference) const;
104  void reset();
105 protected:
107  void fill_histograms(const std::vector<double>& values,
108  double rmin, double rmax,
109  double mmin, double mmax);
110 private:
111 
112  double scale_moving(double x) const;
113  double scale_reference(double x) const;
114 
115  void evaluate_histograms();
116  void evaluate_log_cache();
117 
118  size_t m_ref_bins;
119  PSplineKernel m_ref_kernel;
120  size_t m_ref_border;
121  size_t m_ref_real_bins;
122  double m_ref_max;
123  double m_ref_min;
124  double m_ref_scale;
125 
126  size_t m_mov_bins;
127 
128  PSplineKernel m_mov_kernel;
129  size_t m_mov_border;
130  size_t m_mov_real_bins;
131  double m_mov_max;
132  double m_mov_min;
133  double m_mov_scale;
134 
135 
136  std::vector<double> m_joined_histogram;
137  std::vector<double> m_ref_histogram;
138  std::vector<double> m_mov_histogram;
139 
140  std::vector<std::vector<double> > m_pdfLogCache;
141  double m_cut_histogram;
142 
143  template <typename Iterator>
144  std::pair<double,double> get_reduced_range(Iterator begin, Iterator end)const;
145 };
146 
147 template <typename MovIterator, typename RefIterator>
148 BOOST_CONCEPT_REQUIRES( ((::boost::ForwardIterator<MovIterator>))
149  ((::boost::ForwardIterator<RefIterator>)),
150  (void)
151  )
152  CSplineParzenMI::fill(MovIterator mov_begin, MovIterator mov_end,
153  RefIterator ref_begin, RefIterator ref_end)
154 {
155  std::fill(m_joined_histogram.begin(), m_joined_histogram.end(), 0.0);
156 
157  if (m_mov_max < m_mov_min) {
158  // (re)evaluate the ranges
159  auto mov_range = get_reduced_range(mov_begin, mov_end);
160  if (mov_range.second == mov_range.first)
161  throw std::invalid_argument("relevant moving image intensity range is zero");
162  m_mov_min = mov_range.first;
163  m_mov_max = mov_range.second;
164  m_mov_scale = (m_mov_bins - 1) / (m_mov_max - m_mov_min);
165  cvdebug() << "Mov Range = [" << m_mov_min << ", " << m_mov_max << "]\n";
166  }
167 
168 
169  if (m_ref_max < m_ref_min) {
170  auto ref_range = get_reduced_range(ref_begin, ref_end);
171  if (ref_range.second == ref_range.first)
172  throw std::invalid_argument("relevant reference image intensity range is zero");
173 
174  m_ref_min = ref_range.first;
175  m_ref_max = ref_range.second;
176  m_ref_scale = (m_ref_bins - 1) / (m_ref_max - m_ref_min);
177  cvdebug() << "Ref Range = [" << m_ref_min << ", " << m_ref_max << "]\n";
178  }
179 
180 
181  std::vector<double> mweights(m_mov_kernel->size());
182  std::vector<double> rweights(m_ref_kernel->size());
183 
184  size_t N = 0;
185  while (ref_begin != ref_end && mov_begin != mov_end) {
186 
187  const double mov = scale_moving(*mov_begin);
188  const double ref = scale_reference(*ref_begin);
189 
190  const int mov_start = m_mov_kernel->get_start_idx_and_value_weights(mov, mweights) + m_mov_border;
191  const int ref_start = m_ref_kernel->get_start_idx_and_value_weights(ref, rweights) + m_ref_border;
192 
193  for (size_t r = 0; r < m_ref_kernel->size(); ++r) {
194  auto inbeg = m_joined_histogram.begin() +
195  m_mov_real_bins * (ref_start + r) + mov_start;
196  double rw = rweights[r];
197  std::transform(mweights.begin(), mweights.end(), inbeg, inbeg,
198  [rw](double mw, double jhvalue){ return mw * rw + jhvalue;});
199  }
200 
201  ++N;
202  ++mov_begin;
203  ++ref_begin;
204  }
205 
206  cvdebug() << "CSplineParzenMI::fill: counted " << N << " pixels\n";
207  // normalize joined histogram
208  const double nscale = 1.0/N;
209  transform(m_joined_histogram.begin(), m_joined_histogram.end(), m_joined_histogram.begin(),
210  [&nscale](double jhvalue){return jhvalue * nscale;});
211 
212  evaluate_histograms();
213  evaluate_log_cache();
214 }
215 
216 template <typename Iterator>
217 std::pair<double,double> CSplineParzenMI::get_reduced_range(Iterator begin, Iterator end)const
218 {
219  auto range = std::minmax_element(begin, end);
221  THistogram<Feeder> h(Feeder(*range.first, *range.second, 4096));
222  h.push_range(begin, end);
223  auto reduced_range = h.get_reduced_range(m_cut_histogram);
224  cvinfo() << "CSplineParzenMI: reduce range by "<< m_cut_histogram
225  <<"% from [" << *range.first << ", " << *range.second
226  << "] to [" << reduced_range.first << ", " << reduced_range.second << "]\n";
227  return std::pair<double,double>(reduced_range.first, reduced_range.second);
228 
229 }
230 
232 #endif