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 
29 NS_BEGIN(NS)
30 
31 
37 template <typename TCost>
38 class TSSDCost: public TCost
39 {
40 public:
41  typedef typename TCost::Data Data;
42  typedef typename TCost::Force Force;
43 
44  TSSDCost();
45  TSSDCost(bool normalize, float automask_thresh);
46 private:
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 
54 struct 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 
102 template <typename TCost>
103 TSSDCost<TCost>::TSSDCost():
104  m_normalize(true),
105  m_automask_thresh(0.0)
106 {
107  this->add(::mia::property_gradient);
108 }
109 
110 template <typename TCost>
111 TSSDCost<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 
118 template <typename TCost>
119 double 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 
125 template <typename Force>
126 struct 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  }
177 private:
178  Force& m_force;
179  bool m_normalize;
180  float m_automask_thresh;
181 };
182 
183 
187 template <typename TCost>
188 double 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 
202 template <typename CP, typename C>
203 class TSSDCostPlugin: public CP
204 {
205 public:
206  TSSDCostPlugin();
207  C *do_create()const;
208 private:
209  bool m_normalize;
210  float m_automask_thresh;
211 };
212 
213 
217 template <typename CP, typename C>
218 TSSDCostPlugin<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 
236 template <typename CP, typename C>
237 C *TSSDCostPlugin<CP, C>::do_create() const
238 {
239  return new TSSDCost<C>(m_normalize, m_automask_thresh);
240 }
241 
243 NS_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