21 #ifndef mia_core_splineparzenmi_hh
22 #define mia_core_splineparzenmi_hh
24 #include <boost/concept/requires.hpp>
25 #include <boost/concept_check.hpp>
69 template <
typename MovIterator,
typename RefIterator>
70 BOOST_CONCEPT_REQUIRES( ((::boost::ForwardIterator<MovIterator>))
71 ((::boost::ForwardIterator<RefIterator>)),
74 fill(MovIterator mov_begin, MovIterator mov_end,
75 RefIterator ref_begin, RefIterator ref_end);
91 double get_gradient(
double moving,
double reference)
const;
100 double get_gradient_slow(
double moving,
double reference)
const;
107 void fill_histograms(
const std::vector<double>& values,
108 double rmin,
double rmax,
109 double mmin,
double mmax);
112 double scale_moving(
double x)
const;
113 double scale_reference(
double x)
const;
115 void evaluate_histograms();
116 void evaluate_log_cache();
121 size_t m_ref_real_bins;
130 size_t m_mov_real_bins;
136 std::vector<double> m_joined_histogram;
137 std::vector<double> m_ref_histogram;
138 std::vector<double> m_mov_histogram;
140 std::vector<std::vector<double> > m_pdfLogCache;
141 double m_cut_histogram;
143 template <
typename Iterator>
144 std::pair<double,double> get_reduced_range(Iterator begin, Iterator end)
const;
147 template <
typename MovIterator,
typename RefIterator>
148 BOOST_CONCEPT_REQUIRES( ((::boost::ForwardIterator<MovIterator>))
149 ((::boost::ForwardIterator<RefIterator>)),
153 RefIterator ref_begin, RefIterator ref_end)
155 std::fill(m_joined_histogram.begin(), m_joined_histogram.end(), 0.0);
157 if (m_mov_max < m_mov_min) {
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";
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");
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";
181 std::vector<double> mweights(m_mov_kernel->size());
182 std::vector<double> rweights(m_ref_kernel->size());
185 while (ref_begin != ref_end && mov_begin != mov_end) {
187 const double mov = scale_moving(*mov_begin);
188 const double ref = scale_reference(*ref_begin);
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;
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;});
206 cvdebug() <<
"CSplineParzenMI::fill: counted " << N <<
" pixels\n";
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;});
212 evaluate_histograms();
213 evaluate_log_cache();
216 template <
typename Iterator>
217 std::pair<double,double> CSplineParzenMI::get_reduced_range(Iterator begin, Iterator end)
const
219 auto range = std::minmax_element(begin, end);
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);