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 assert(mov_begin != mov_end);
158 assert(ref_begin != ref_end);
160 if (m_mov_max < m_mov_min) {
162 auto mov_range = get_reduced_range(mov_begin, mov_end);
163 if (mov_range.second == mov_range.first)
164 throw std::invalid_argument(
"relevant moving image intensity range is zero");
165 m_mov_min = mov_range.first;
166 m_mov_max = mov_range.second;
167 m_mov_scale = (m_mov_bins - 1) / (m_mov_max - m_mov_min);
168 cvdebug() <<
"Mov Range = [" << m_mov_min <<
", " << m_mov_max <<
"]\n";
172 if (m_ref_max < m_ref_min) {
173 auto ref_range = get_reduced_range(ref_begin, ref_end);
174 if (ref_range.second == ref_range.first)
175 throw std::invalid_argument(
"relevant reference image intensity range is zero");
177 m_ref_min = ref_range.first;
178 m_ref_max = ref_range.second;
179 m_ref_scale = (m_ref_bins - 1) / (m_ref_max - m_ref_min);
180 cvdebug() <<
"Ref Range = [" << m_ref_min <<
", " << m_ref_max <<
"]\n";
184 std::vector<double> mweights(m_mov_kernel->size());
185 std::vector<double> rweights(m_ref_kernel->size());
188 while (ref_begin != ref_end && mov_begin != mov_end) {
190 const double mov = scale_moving(*mov_begin);
191 const double ref = scale_reference(*ref_begin);
193 const int mov_start = m_mov_kernel->get_start_idx_and_value_weights(mov, mweights) + m_mov_border;
194 const int ref_start = m_ref_kernel->get_start_idx_and_value_weights(ref, rweights) + m_ref_border;
196 for (
size_t r = 0; r < m_ref_kernel->size(); ++r) {
197 auto inbeg = m_joined_histogram.begin() +
198 m_mov_real_bins * (ref_start + r) + mov_start;
199 double rw = rweights[r];
200 std::transform(mweights.begin(), mweights.end(), inbeg, inbeg,
201 [rw](
double mw,
double jhvalue){
return mw * rw + jhvalue;});
209 cvdebug() <<
"CSplineParzenMI::fill: counted " << N <<
" pixels\n";
211 const double nscale = 1.0/N;
212 transform(m_joined_histogram.begin(), m_joined_histogram.end(), m_joined_histogram.begin(),
213 [&nscale](
double jhvalue){
return jhvalue * nscale;});
215 evaluate_histograms();
216 evaluate_log_cache();
219 template <
typename Iterator>
220 std::pair<double,double> CSplineParzenMI::get_reduced_range(Iterator begin, Iterator end)
const
222 auto range = std::minmax_element(begin, end);
225 h.push_range(begin, end);
226 auto reduced_range = h.get_reduced_range(m_cut_histogram);
227 cvinfo() <<
"CSplineParzenMI: reduce range by "<< m_cut_histogram
228 <<
"% from [" << *range.first <<
", " << *range.second
229 <<
"] to [" << reduced_range.first <<
", " << reduced_range.second <<
"]\n";
230 return std::pair<double,double>(reduced_range.first, reduced_range.second);
vstream & cvinfo()
informal output that may be of interest to understand problems with a program and are of higher prior...
#define NS_MIA_BEGIN
conveniance define to start the mia namespace
Implementation of mutual information based on B-splines.
a simple histogram that uses an instance of THistogramFeeder as input converter
std::shared_ptr< CSplineKernel > PSplineKernel
EXPORT_2D C2DFVectorfield get_gradient(const C2DImage &image)
A class to normalize and quantizize input data to a given histogram range with its given number of bi...
vstream & cvdebug()
Short for debug output in non-debug build output send to this will be ignored.
#define EXPORT_CORE
Macro to manage Visual C++ style dllimport/dllexport.
#define NS_MIA_END
conveniance define to end the mia namespace