lsd.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/>.
*
*/
#include <mia/core/filter.hh>
#include <mia/core/msgstream.hh>
#include <mia/core/parameter.hh>
#include <mia/core/property_flags.hh>
#include <numeric>
#include <limits>
NS_BEGIN(NS)
struct valpos {
double val;
int pos;
};
template <typename TCost>
class TLSDImageCost: public TCost
{
public:
typedef typename TCost::Data Data;
typedef typename TCost::Force Force;
private:
class CRefPrepare : public mia::TFilter<void>
{
public:
CRefPrepare(std::vector<double>& QtQinv, std::vector<int>& Q_mappping):
m_QtQinv(QtQinv),
m_Q_mappping(Q_mappping)
{
}
template <typename DataTempl>
void operator()(const DataTempl& ref);
std::vector<double>& m_QtQinv;
std::vector<int>& m_Q_mappping;
};
class RunCost : mia::TFilter<double>
{
public:
typedef TFilter<double>::result_type result_type;
RunCost( const std::vector<double>& QtQinv, const std::vector<int>& Q_mappping);
template <typename DataTempl>
double operator()(const DataTempl& ref)const;
template <typename DataTempl>
double operator()(const DataTempl& ref, Force& force)const;
private:
const std::vector<double>& m_QtQinv;
const std::vector<int>& m_Q_mappping;
};
virtual double do_value(const Data& a, const Data& b) const;
virtual double do_evaluate_force(const Data& a, const Data& /*b*/, Force& force) const;
virtual void post_set_reference(const Data& ref);
std::vector<double> m_QtQinv;
std::vector<int> m_Q_mappping;
};
inline bool operator < (const valpos& a, const valpos& b)
{
return a.val < b.val;
}
template <typename CP, typename C>
class TLSDImageCostPlugin: public CP
{
public:
TLSDImageCostPlugin();
virtual TLSDImageCost<C> *do_create()const;
const std::string do_get_descr()const;
};
template <typename TCost>
template <typename DataTempl>
void TLSDImageCost<TCost>::CRefPrepare::operator()(const DataTempl& ref)
{
size_t npixels = ref.get_size().x * ref.get_size().y;
std::vector<valpos> buffer(npixels);
m_Q_mappping.resize(npixels);
m_QtQinv.resize(npixels);
int idx = 0;
for ( auto x = ref.begin(); x != ref.end(); ++x, ++idx) {
double value = *x;
valpos vp = {value, idx};
buffer[idx] = vp;
}
std::sort(buffer.begin(), buffer.end());
idx = 0;
auto x = buffer.begin();
double oldv = x->val;
++x;
++m_QtQinv[idx];
// build histogram and prepare target mapping
m_Q_mappping[x->pos] = idx;
while ( x != buffer.end() ) {
if (x->val != oldv) {
oldv = x->val;
++idx;
}
++m_QtQinv[idx];
m_Q_mappping[x->pos] = idx;
++x;
}
++idx;
m_QtQinv.resize(idx);
std::transform(m_QtQinv.begin(), m_QtQinv.end(), m_QtQinv.begin(),
[](double x) {
return 1.0 / x;
});
}
template <typename TCost>
void TLSDImageCost<TCost>::post_set_reference(const Data& ref)
{
CRefPrepare rp(m_QtQinv, m_Q_mappping);
mia::accumulate(rp, ref);
}
template <typename TCost>
double TLSDImageCost<TCost>::do_value(const Data& a, const Data& /*b*/) const
{
RunCost rf(m_QtQinv, m_Q_mappping);
return mia::filter(rf, a);
}
template <typename TCost>
double TLSDImageCost<TCost>::do_evaluate_force(const Data& a, const Data& /*b*/, Force& force) const
{
RunCost rf(m_QtQinv, m_Q_mappping);
return mia::filter_and_output(rf, a, force);
}
template <typename TCost>
TLSDImageCost<TCost>::RunCost::RunCost(const std::vector<double>& QtQinv, const std::vector<int>& Q_mappping):
m_QtQinv(QtQinv),
m_Q_mappping(Q_mappping)
{
}
template <typename TCost>
template <typename DataTempl>
double TLSDImageCost<TCost>::RunCost::operator()(const DataTempl& a)const
{
double val1 = 0.0;
double val2 = 0.0;
std::vector<double> sums(m_QtQinv.size(), 0.0);
auto idx = m_Q_mappping.begin();
for (auto ia = a.begin(); ia != a.end(); ++ia, ++idx) {
val1 += *ia * *ia;
sums[*idx] += *ia;
}
for (size_t i = 0; i < sums.size(); ++i)
val2 += sums[i] * sums[i] * m_QtQinv[i];
return 0.5 * (val1 - val2);
}
template <typename TCost>
template <typename DataTempl>
double TLSDImageCost<TCost>::RunCost::operator()(const DataTempl& a, Force& force)const
{
double value = 0.0;
std::vector<double> sums(m_QtQinv.size(), 0.0);
auto idx = m_Q_mappping.begin();
for (auto ia = a.begin(); ia != a.end(); ++ia, ++idx) {
value += *ia * *ia;
sums[*idx] += *ia;
}
auto s = sums.begin();
for (auto q = m_QtQinv.begin(); q != m_QtQinv.end(); ++q, ++s) {
value -= *q * *s * *s;
*s *= *q;
}
auto gradient = get_gradient(a);
auto iforce = force.begin();
auto igrad = gradient.begin();
idx = m_Q_mappping.begin();
for (auto ia = a.begin(); ia != a.end(); ++ia, ++iforce, ++igrad, ++idx)
*iforce = *igrad * (*ia - sums[*idx]);
return 0.5 * value;
}
template <typename CP, typename C>
TLSDImageCostPlugin<CP, C>::TLSDImageCostPlugin():
CP("lsd")
{
}
template <typename CP, typename C>
TLSDImageCost<C> *TLSDImageCostPlugin<CP, C>::do_create()const
{
return new TLSDImageCost<C>();
}
template <typename CP, typename C>
const std::string TLSDImageCostPlugin<CP, C>::do_get_descr()const
{
return "Least-Squares Distance measure";
}
NS_END
get_gradient
EXPORT_2D C2DFVectorfield get_gradient(const C2DImage &image)
operator<
bool operator<(const T2DVector< T > &a, const T2DVector< S > &b)
Definition 2d/vector.hh:495
TCost
The generic cost function interface.
Definition core/cost.hh:65
TCost::Force
V Force
typedef for generic programming: The gradient forca type create by the cost function
Definition core/cost.hh:71
TCost::Data
T Data
typedef for generic programming: The data type used by the cost function
Definition core/cost.hh:68
filter.hh
NS_BEGIN
#define NS_BEGIN(NS)
conveniance define to start a namespace
Definition defines.hh:42
NS_END
#define NS_END
conveniance define to end a namespace
Definition defines.hh:45
msgstream.hh
parameter.hh
property_flags.hh
TFilter::result_type
R result_type
defines the return type of the filter function
Definition core/filter.hh:72