ssd.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)
template <typename TCost>
class TSSDCost: public TCost
{
public:
typedef typename TCost::Data Data;
typedef typename TCost::Force Force;
TSSDCost();
TSSDCost(bool normalize, float automask_thresh);
private:
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;
bool m_normalize;
float m_automask_thresh;
};
struct FEvalSSD : public mia::TFilter<double> {
FEvalSSD(bool normalize, float automask_thresh):
m_normalize(normalize), m_automask_thresh(automask_thresh) {}
template <typename T, typename R>
struct SQD {
double operator ()(T a, R b) const
{
double d = (double)a - (double)b;
return d * d;
}
};
template <typename T, typename R>
FEvalSSD::result_type operator () (const T& a, const R& b) const
{
double scale = m_normalize ? 0.5 / a.size() : 0.5;
if (m_automask_thresh == 0.0f)
return scale * inner_product(a.begin(), a.end(), b.begin(), 0.0, ::std::plus<double>(),
SQD<typename T::value_type, typename R::value_type >());
else {
auto ia = a.begin();
auto ib = b.begin();
auto ie = a.end();
long n = 0;
double sum = 0.0;
while (ia != ie) {
if (*ia > m_automask_thresh) {
double d = (double) * ia - (double) * ib;
sum += d * d;
++n;
}
++ia;
++ib;
}
// high penalty if the mask don't overlap at all
return n > 0 ? 0.5 * sum / n : std::numeric_limits<float>::max();
}
}
bool m_normalize;
float m_automask_thresh;
};
template <typename TCost>
TSSDCost<TCost>::TSSDCost():
m_normalize(true),
m_automask_thresh(0.0)
{
this->add(::mia::property_gradient);
}
template <typename TCost>
TSSDCost<TCost>::TSSDCost(bool normalize, float automask_thresh):
m_normalize(normalize),
m_automask_thresh(automask_thresh)
{
this->add(::mia::property_gradient);
}
template <typename TCost>
double TSSDCost<TCost>::do_value(const Data& a, const Data& b) const
{
FEvalSSD essd(m_normalize, m_automask_thresh);
return filter(essd, a, b);
}
template <typename Force>
struct FEvalForce: public mia::TFilter<float> {
FEvalForce(Force& force, bool normalize, float automask_thresh):
m_force(force),
m_normalize(normalize),
m_automask_thresh(automask_thresh)
{
}
template <typename T, typename R>
float operator ()( const T& a, const R& b) const
{
Force gradient = get_gradient(a);
float cost = 0.0;
auto ai = a.begin();
auto bi = b.begin();
auto fi = m_force.begin();
auto gi = gradient.begin();
if (m_automask_thresh == 0.0f) {
float scale = m_normalize ? 1.0 / a.size() : 1.0;
for (size_t i = 0; i < a.size(); ++i, ++ai, ++bi, ++fi, ++gi) {
float delta = float(*ai) - float(*bi);
*fi = *gi * delta * scale;
cost += delta * delta * scale;
}
return 0.5 * cost;
} else {
long n = 0;
for (size_t i = 0; i < a.size(); ++i, ++ai, ++bi, ++fi, ++gi) {
if (*ai > m_automask_thresh) {
float delta = float(*ai) - float(*bi);
*fi = *gi * delta;
cost += delta * delta;
++n;
}
}
if ( n > 0) {
float scale = 1.0f / n;
transform(m_force.begin(), m_force.end(), m_force.begin(),
[scale](const typename Force::value_type & x) {
return scale * x;
});
return 0.5 * cost * scale;
} else {
return std::numeric_limits<float>::max();
}
}
}
private:
Force& m_force;
bool m_normalize;
float m_automask_thresh;
};
template <typename TCost>
double TSSDCost<TCost>::do_evaluate_force(const Data& a, const Data& b, Force& force) const
{
assert(a.get_size() == b.get_size());
assert(a.get_size() == force.get_size());
FEvalForce<Force> ef(force, m_normalize, m_automask_thresh);
return filter(ef, a, b);
}
template <typename CP, typename C>
class TSSDCostPlugin: public CP
{
public:
TSSDCostPlugin();
C *do_create()const;
private:
bool m_normalize;
float m_automask_thresh;
};
template <typename CP, typename C>
TSSDCostPlugin<CP, C>::TSSDCostPlugin():
CP("ssd"),
m_normalize(false),
m_automask_thresh(0)
{
TRACE("TSSDCostPlugin<CP,C>::TSSDCostPlugin()");
this->add_property(::mia::property_gradient);
this->add_parameter("norm", new mia::CBoolParameter(m_normalize, false,
"Set whether the metric should be normalized by the number of image pixels")
);
this->add_parameter("autothresh", mia::make_ci_param(m_automask_thresh, 0.0f, 1000.0f, false,
"Use automatic masking of the moving image by only takeing "
"intensity values into accound that are larger than the given threshold"));
}
template <typename CP, typename C>
C *TSSDCostPlugin<CP, C>::do_create() const
{
return new TSSDCost<C>(m_normalize, m_automask_thresh);
}
NS_END
get_gradient
EXPORT_2D C2DFVectorfield get_gradient(const C2DImage &image)
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
filter
static F::result_type filter(const F &f, const B &b)
Definition core/filter.hh:258
msgstream.hh
TRACE
#define TRACE(DOMAIN)
Definition msgstream.hh:208
normalize
std::shared_ptr< Image > normalize(const Image &image)
a normalizer for image intensities
Definition normalize.hh:135
parameter.hh
property_flags.hh