splineparzenmi.hh
Go to the documentation of this file.
/* -*- mia-c++ -*-
*
* This file is part of MIA - a toolbox for medical image analysis
* Copyright (c) Leipzig, Madrid 1999-2017 Gert Wollny
*
* MIA is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with MIA; if not, see <http://www.gnu.org/licenses/>.
*
*/
#ifndef mia_core_splineparzenmi_hh
#define mia_core_splineparzenmi_hh
#include <boost/concept/requires.hpp>
#include <boost/concept_check.hpp>
#include <mia/core/splinekernel.hh>
#include <mia/core/histogram.hh>
NS_MIA_BEGIN
class EXPORT_CORE CSplineParzenMI
{
public:
CSplineParzenMI(size_t rbins, PSplineKernel rkernel,
size_t mbins, PSplineKernel mkernel, double cut_histogram);
template <typename MovIterator, typename RefIterator>
BOOST_CONCEPT_REQUIRES( ((::boost::ForwardIterator<MovIterator>))
((::boost::ForwardIterator<RefIterator>)),
(void)
)
fill(MovIterator mov_begin, MovIterator mov_end,
RefIterator ref_begin, RefIterator ref_end);
template <typename MovIterator, typename RefIterator, typename MaskIterator>
void fill(MovIterator mov_begin, MovIterator mov_end,
RefIterator ref_begin, RefIterator ref_end,
MaskIterator mask_begin, MaskIterator mask_end);
double value() const;
double get_gradient(double moving, double reference) const;
double get_gradient_slow(double moving, double reference) const;
void reset();
protected:
void fill_histograms(const std::vector<double>& values,
double rmin, double rmax,
double mmin, double mmax);
private:
double scale_moving(double x) const;
double scale_reference(double x) const;
void evaluate_histograms();
void evaluate_log_cache();
size_t m_ref_bins;
PSplineKernel m_ref_kernel;
size_t m_ref_border;
size_t m_ref_real_bins;
double m_ref_max;
double m_ref_min;
double m_ref_scale;
size_t m_mov_bins;
PSplineKernel m_mov_kernel;
size_t m_mov_border;
size_t m_mov_real_bins;
double m_mov_max;
double m_mov_min;
double m_mov_scale;
std::vector<double> m_joined_histogram;
std::vector<double> m_ref_histogram;
std::vector<double> m_mov_histogram;
std::vector<std::vector<double>> m_pdfLogCache;
double m_cut_histogram;
template <typename Iterator>
std::pair<double, double> get_reduced_range(Iterator begin, Iterator end)const;
template <typename DataIterator, typename MaskIterator>
std::pair<double, double>
get_reduced_range_masked(DataIterator dbegin,
DataIterator dend,
MaskIterator mbegin)const;
};
template <typename MovIterator, typename RefIterator>
BOOST_CONCEPT_REQUIRES( ((::boost::ForwardIterator<MovIterator>))
((::boost::ForwardIterator<RefIterator>)),
(void)
)
CSplineParzenMI::fill(MovIterator mov_begin, MovIterator mov_end,
RefIterator ref_begin, RefIterator ref_end)
{
std::fill(m_joined_histogram.begin(), m_joined_histogram.end(), 0.0);
assert(mov_begin != mov_end);
assert(ref_begin != ref_end);
if (m_mov_max < m_mov_min) {
// (re)evaluate the ranges
auto mov_range = get_reduced_range(mov_begin, mov_end);
if (mov_range.second == mov_range.first)
throw std::invalid_argument("relevant moving image intensity range is zero");
m_mov_min = mov_range.first;
m_mov_max = mov_range.second;
m_mov_scale = (m_mov_bins - 1) / (m_mov_max - m_mov_min);
cvdebug() << "Mov Range = [" << m_mov_min << ", " << m_mov_max << "]\n";
}
if (m_ref_max < m_ref_min) {
auto ref_range = get_reduced_range(ref_begin, ref_end);
if (ref_range.second == ref_range.first)
throw std::invalid_argument("relevant reference image intensity range is zero");
m_ref_min = ref_range.first;
m_ref_max = ref_range.second;
m_ref_scale = (m_ref_bins - 1) / (m_ref_max - m_ref_min);
cvdebug() << "Ref Range = [" << m_ref_min << ", " << m_ref_max << "]\n";
}
std::vector<double> mweights(m_mov_kernel->size());
std::vector<double> rweights(m_ref_kernel->size());
size_t N = 0;
while (ref_begin != ref_end && mov_begin != mov_end) {
const double mov = scale_moving(*mov_begin);
const double ref = scale_reference(*ref_begin);
const int mov_start = m_mov_kernel->get_start_idx_and_value_weights(mov, mweights) + m_mov_border;
const int ref_start = m_ref_kernel->get_start_idx_and_value_weights(ref, rweights) + m_ref_border;
for (size_t r = 0; r < m_ref_kernel->size(); ++r) {
auto inbeg = m_joined_histogram.begin() +
m_mov_real_bins * (ref_start + r) + mov_start;
double rw = rweights[r];
std::transform(mweights.begin(), mweights.end(), inbeg, inbeg,
[rw](double mw, double jhvalue) {
return mw * rw + jhvalue;
});
}
++N;
++mov_begin;
++ref_begin;
}
cvdebug() << "CSplineParzenMI::fill: counted " << N << " pixels\n";
// normalize joined histogram
const double nscale = 1.0 / N;
transform(m_joined_histogram.begin(), m_joined_histogram.end(), m_joined_histogram.begin(),
[&nscale](double jhvalue) {
return jhvalue * nscale;
});
evaluate_histograms();
evaluate_log_cache();
}
template <typename MovIterator, typename RefIterator, typename MaskIterator>
void CSplineParzenMI::fill(MovIterator mov_begin, MovIterator mov_end,
RefIterator ref_begin, RefIterator ref_end,
MaskIterator mask_begin, MaskIterator mask_end)
{
std::fill(m_joined_histogram.begin(), m_joined_histogram.end(), 0.0);
assert(mov_begin != mov_end);
assert(ref_begin != ref_end);
// no mask
if (mask_begin == mask_end)
fill(mov_begin, mov_end, ref_begin, ref_end);
if (m_mov_max < m_mov_min) {
// (re)evaluate the ranges
auto mov_range = get_reduced_range_masked(mov_begin, mov_end, mask_begin);
if (mov_range.second == mov_range.first)
throw std::invalid_argument("relevant moving image intensity range is zero");
m_mov_min = mov_range.first;
m_mov_max = mov_range.second;
m_mov_scale = (m_mov_bins - 1) / (m_mov_max - m_mov_min);
cvdebug() << "Mov Range = [" << m_mov_min << ", " << m_mov_max << "]\n";
}
if (m_ref_max < m_ref_min) {
auto ref_range = get_reduced_range_masked(ref_begin, ref_end, mask_begin);
if (ref_range.second == ref_range.first)
throw std::invalid_argument("relevant reference image intensity range is zero");
m_ref_min = ref_range.first;
m_ref_max = ref_range.second;
m_ref_scale = (m_ref_bins - 1) / (m_ref_max - m_ref_min);
cvdebug() << "Ref Range = [" << m_ref_min << ", " << m_ref_max << "]\n";
}
std::vector<double> mweights(m_mov_kernel->size());
std::vector<double> rweights(m_ref_kernel->size());
size_t N = 0;
while (ref_begin != ref_end && mov_begin != mov_end) {
if (*mask_begin) {
const double mov = scale_moving(*mov_begin);
const double ref = scale_reference(*ref_begin);
const int mov_start = m_mov_kernel->get_start_idx_and_value_weights(mov, mweights) + m_mov_border;
const int ref_start = m_ref_kernel->get_start_idx_and_value_weights(ref, rweights) + m_ref_border;
for (size_t r = 0; r < m_ref_kernel->size(); ++r) {
auto inbeg = m_joined_histogram.begin() +
m_mov_real_bins * (ref_start + r) + mov_start;
double rw = rweights[r];
std::transform(mweights.begin(), mweights.end(), inbeg, inbeg,
[rw](double mw, double jhvalue) {
return mw * rw + jhvalue;
});
}
++N;
}
++mask_begin;
++mov_begin;
++ref_begin;
}
cvdebug() << "CSplineParzenMI::fill: counted " << N << " pixels\n";
// normalize joined histogram
const double nscale = 1.0 / N;
transform(m_joined_histogram.begin(), m_joined_histogram.end(), m_joined_histogram.begin(),
[&nscale](double jhvalue) {
return jhvalue * nscale;
});
evaluate_histograms();
evaluate_log_cache();
}
template <typename Iterator>
std::pair<double, double> CSplineParzenMI::get_reduced_range(Iterator begin, Iterator end)const
{
auto range = std::minmax_element(begin, end);
typedef THistogramFeeder<typename Iterator::value_type> Feeder;
THistogram<Feeder> h(Feeder(*range.first, *range.second, 4096));
h.push_range(begin, end);
auto reduced_range = h.get_reduced_range(m_cut_histogram);
cvinfo() << "CSplineParzenMI: reduce range by " << m_cut_histogram
<< "% from [" << *range.first << ", " << *range.second
<< "] to [" << reduced_range.first << ", " << reduced_range.second << "]\n";
return std::make_pair(reduced_range.first, reduced_range.second);
}
template <typename DataIterator, typename MaskIterator>
std::pair<double, double>
CSplineParzenMI::get_reduced_range_masked(DataIterator dbegin,
DataIterator dend,
MaskIterator mbegin)const
{
auto ib = dbegin;
auto im = mbegin;
while (! *im && ib != dend) {
++im;
++ib;
}
if (ib == dend)
throw std::runtime_error("CSplineParzenMI: empty mask");
double range_max = *ib;
double range_min = *ib;
++ib;
++im;
while (ib != dend) {
if (*im) {
if (*ib < range_min)
range_min = *ib;
if (*ib > range_max)
range_max = *ib;
}
++ib;
++im;
}
typedef THistogramFeeder<typename DataIterator::value_type> Feeder;
THistogram<Feeder> h(Feeder(range_min, range_max, 4096));
ib = dbegin;
im = mbegin;
while (ib != dend) {
if (*im)
h.push(*ib);
++ib;
++im;
}
auto reduced_range = h.get_reduced_range(m_cut_histogram);
cvinfo() << "CSplineParzenMI: reduce range by " << m_cut_histogram
<< "% from [" << range_min << ", " << range_max
<< "] to [" << reduced_range.first << ", " << reduced_range.second << "]\n";
return std::make_pair(reduced_range.first, reduced_range.second);
}
NS_MIA_END
#endif
CSplineParzenMI::fill
void fill(MovIterator mov_begin, MovIterator mov_end, RefIterator ref_begin, RefIterator ref_end)
Definition splineparzenmi.hh:179
CSplineParzenMI::get_gradient_slow
double get_gradient_slow(double moving, double reference) const
CSplineParzenMI::value
double value() const
CSplineParzenMI::fill_histograms
void fill_histograms(const std::vector< double > &values, double rmin, double rmax, double mmin, double mmax)
CSplineParzenMI::reset
void reset()
CSplineParzenMI::get_gradient
double get_gradient(double moving, double reference) const
CSplineParzenMI::CSplineParzenMI
CSplineParzenMI(size_t rbins, PSplineKernel rkernel, size_t mbins, PSplineKernel mkernel, double cut_histogram)
THistogramFeeder
A class to normalize and quantizize input data to a given histogram range with its given number of bi...
Definition histogram.hh:50
THistogram
a simple histogram that uses an instance of THistogramFeeder as input converter
Definition histogram.hh:135
NS_MIA_BEGIN
#define NS_MIA_BEGIN
conveniance define to start the mia namespace
Definition defines.hh:33
EXPORT_CORE
#define EXPORT_CORE
Macro to manage Visual C++ style dllimport/dllexport.
Definition defines.hh:101
NS_MIA_END
#define NS_MIA_END
conveniance define to end the mia namespace
Definition defines.hh:36
PSplineKernel
std::shared_ptr< CSplineKernel > PSplineKernel
Definition splinekernel.hh:291
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
histogram.hh
cvdebug
CDebugSink & cvdebug()
Definition msgstream.hh:226
splinekernel.hh