My Project
ssd.hh
Go to the documentation of this file.
1/* -*- mia-c++ -*-
2 *
3 * This file is part of MIA - a toolbox for medical image analysis
4 * Copyright (c) Leipzig, Madrid 1999-2017 Gert Wollny
5 *
6 * MIA is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with MIA; if not, see <http://www.gnu.org/licenses/>.
18 *
19 */
20
21#include <mia/core/filter.hh>
22#include <mia/core/msgstream.hh>
23#include <mia/core/parameter.hh>
25
26#include <numeric>
27#include <limits>
28
29NS_BEGIN(NS)
30
31
32
33
37template <typename TCost>
38class TSSDCost: public TCost
39{
40public:
41 typedef typename TCost::Data Data;
42 typedef typename TCost::Force Force;
43
44 TSSDCost();
45 TSSDCost(bool normalize, float automask_thresh);
46private:
47 virtual double do_value(const Data& a, const Data& b) const;
48 virtual double do_evaluate_force(const Data& a, const Data& b, Force& force) const;
49 bool m_normalize;
50 float m_automask_thresh;
51};
52
53
54struct FEvalSSD : public mia::TFilter<double> {
55 FEvalSSD(bool normalize, float automask_thresh):
56 m_normalize(normalize), m_automask_thresh(automask_thresh) {}
57
58 template <typename T, typename R>
59 struct SQD {
60 double operator ()(T a, R b) const
61 {
62 double d = (double)a - (double)b;
63 return d * d;
64 }
65 };
66
67 template <typename T, typename R>
68 FEvalSSD::result_type operator () (const T& a, const R& b) const
69 {
70 double scale = m_normalize ? 0.5 / a.size() : 0.5;
71
72 if (m_automask_thresh == 0.0f)
73 return scale * inner_product(a.begin(), a.end(), b.begin(), 0.0, ::std::plus<double>(),
74 SQD<typename T::value_type, typename R::value_type >());
75 else {
76 auto ia = a.begin();
77 auto ib = b.begin();
78 auto ie = a.end();
79 long n = 0;
80 double sum = 0.0;
81
82 while (ia != ie) {
83 if (*ia > m_automask_thresh) {
84 double d = (double) * ia - (double) * ib;
85 sum += d * d;
86 ++n;
87 }
88
89 ++ia;
90 ++ib;
91 }
92
93 // high penalty if the mask don't overlap at all
94 return n > 0 ? 0.5 * sum / n : std::numeric_limits<float>::max();
95 }
96 }
97 bool m_normalize;
98 float m_automask_thresh;
99};
100
101
102template <typename TCost>
103TSSDCost<TCost>::TSSDCost():
104 m_normalize(true),
105 m_automask_thresh(0.0)
106{
107 this->add(::mia::property_gradient);
108}
109
110template <typename TCost>
111TSSDCost<TCost>::TSSDCost(bool normalize, float automask_thresh):
112 m_normalize(normalize),
113 m_automask_thresh(automask_thresh)
114{
115 this->add(::mia::property_gradient);
116}
117
118template <typename TCost>
119double TSSDCost<TCost>::do_value(const Data& a, const Data& b) const
120{
121 FEvalSSD essd(m_normalize, m_automask_thresh);
122 return filter(essd, a, b);
123}
124
125template <typename Force>
126struct FEvalForce: public mia::TFilter<float> {
127 FEvalForce(Force& force, bool normalize, float automask_thresh):
128 m_force(force),
129 m_normalize(normalize),
130 m_automask_thresh(automask_thresh)
131 {
132 }
133 template <typename T, typename R>
134 float operator ()( const T& a, const R& b) const
135 {
136 Force gradient = get_gradient(a);
137 float cost = 0.0;
138 auto ai = a.begin();
139 auto bi = b.begin();
140 auto fi = m_force.begin();
141 auto gi = gradient.begin();
142
143 if (m_automask_thresh == 0.0f) {
144 float scale = m_normalize ? 1.0 / a.size() : 1.0;
145
146 for (size_t i = 0; i < a.size(); ++i, ++ai, ++bi, ++fi, ++gi) {
147 float delta = float(*ai) - float(*bi);
148 *fi = *gi * delta * scale;
149 cost += delta * delta * scale;
150 }
151
152 return 0.5 * cost;
153 } else {
154 long n = 0;
155
156 for (size_t i = 0; i < a.size(); ++i, ++ai, ++bi, ++fi, ++gi) {
157 if (*ai > m_automask_thresh) {
158 float delta = float(*ai) - float(*bi);
159 *fi = *gi * delta;
160 cost += delta * delta;
161 ++n;
162 }
163 }
164
165 if ( n > 0) {
166 float scale = 1.0f / n;
167 transform(m_force.begin(), m_force.end(), m_force.begin(),
168 [scale](const typename Force::value_type & x) {
169 return scale * x;
170 });
171 return 0.5 * cost * scale;
172 } else {
173 return std::numeric_limits<float>::max();
174 }
175 }
176 }
177private:
178 Force& m_force;
179 bool m_normalize;
180 float m_automask_thresh;
181};
182
183
187template <typename TCost>
188double TSSDCost<TCost>::do_evaluate_force(const Data& a, const Data& b, Force& force) const
189{
190 assert(a.get_size() == b.get_size());
191 assert(a.get_size() == force.get_size());
192 FEvalForce<Force> ef(force, m_normalize, m_automask_thresh);
193 return filter(ef, a, b);
194}
195
196
202template <typename CP, typename C>
203class TSSDCostPlugin: public CP
204{
205public:
206 TSSDCostPlugin();
207 C *do_create()const;
208private:
209 bool m_normalize;
210 float m_automask_thresh;
211};
212
213
217template <typename CP, typename C>
218TSSDCostPlugin<CP, C>::TSSDCostPlugin():
219 CP("ssd"),
220 m_normalize(false),
221 m_automask_thresh(0)
222{
223 TRACE("TSSDCostPlugin<CP,C>::TSSDCostPlugin()");
224 this->add_property(::mia::property_gradient);
225 this->add_parameter("norm", new mia::CBoolParameter(m_normalize, false,
226 "Set whether the metric should be normalized by the number of image pixels")
227 );
228 this->add_parameter("autothresh", mia::make_ci_param(m_automask_thresh, 0.0f, 1000.0f, false,
229 "Use automatic masking of the moving image by only takeing "
230 "intensity values into accound that are larger than the given threshold"));
231}
232
236template <typename CP, typename C>
237C *TSSDCostPlugin<CP, C>::do_create() const
238{
239 return new TSSDCost<C>(m_normalize, m_automask_thresh);
240}
241
243NS_END
EXPORT_2D C2DFVectorfield get_gradient(const C2DImage &image)
The generic cost function interface.
Definition: core/cost.hh:65
V Force
typedef for generic programming: The gradient forca type create by the cost function
Definition: core/cost.hh:71
T Data
typedef for generic programming: The data type used by the cost function
Definition: core/cost.hh:68
#define NS_BEGIN(NS)
conveniance define to start a namespace
Definition: defines.hh:42
#define NS_END
conveniance define to end a namespace
Definition: defines.hh:45
static F::result_type filter(const F &f, const B &b)
Definition: core/filter.hh:258
#define TRACE(DOMAIN)
Definition: msgstream.hh:208
std::shared_ptr< Image > normalize(const Image &image)
a normalizer for image intensities
Definition: normalize.hh:135
CParameter * make_ci_param(T &value, S1 lower_bound, S2 upper_bound, bool required, const char *descr)
Definition: parameter.hh:329
CTParameter< bool > CBoolParameter
boolean parameter
Definition: parameter.hh:561
EXPORT_CORE const char * property_gradient
constant defining the gradient property