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-2013 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 public:
40  typedef typename TCost::Data Data;
41  typedef typename TCost::Force Force;
42 
43  TSSDCost();
44  TSSDCost(bool normalize);
45 private:
46  virtual double do_value(const Data& a, const Data& b) const;
47  virtual double do_evaluate_force(const Data& a, const Data& b, Force& force) const;
48  bool m_normalize;
49 };
50 
51 
52 struct FEvalSSD : public mia::TFilter<double> {
53  FEvalSSD(bool normalize):m_normalize(normalize){}
54 
55  template <typename T, typename R>
56  struct SQD {
57  double operator ()(T a, R b) const {
58  double d = (double)a - (double)b;
59  return d * d;
60  }
61  };
62 
63  template <typename T, typename R>
64  FEvalSSD::result_type operator () (const T& a, const R& b) const {
65  double scale = m_normalize ? 0.5 / a.size() : 0.5;
66  return scale * inner_product(a.begin(), a.end(), b.begin(), 0.0, ::std::plus<double>(),
67  SQD<typename T::value_type , typename R::value_type >());
68  }
69  bool m_normalize;
70 };
71 
72 
73 template <typename TCost>
74 TSSDCost<TCost>::TSSDCost():
75  m_normalize(true)
76 {
77  this->add(::mia::property_gradient);
78 }
79 
80 template <typename TCost>
81 TSSDCost<TCost>::TSSDCost(bool normalize):
82  m_normalize(normalize)
83 {
84  this->add(::mia::property_gradient);
85 }
86 
87 template <typename TCost>
88 double TSSDCost<TCost>::do_value(const Data& a, const Data& b) const
89 {
90  FEvalSSD essd(m_normalize);
91  return filter(essd, a, b);
92 }
93 
94 template <typename Force>
95 struct FEvalForce: public mia::TFilter<float> {
96  FEvalForce(Force& force, bool normalize):
97  m_force(force),
98  m_normalize(normalize)
99  {
100  }
101  template <typename T, typename R>
102  float operator ()( const T& a, const R& b) const {
103  Force gradient = get_gradient(a);
104  float cost = 0.0;
105  auto ai = a.begin();
106  auto bi = b.begin();
107  auto fi = m_force.begin();
108  auto gi = gradient.begin();
109 
110  float scale = m_normalize ? 1.0 / a.size() : 1.0;
111 
112  for (size_t i = 0; i < a.size(); ++i, ++ai, ++bi, ++fi, ++gi) {
113  float delta = float(*ai) - float(*bi);
114  *fi = *gi * delta * scale;
115  cost += delta * delta * scale;
116  }
117  return 0.5 * cost;
118  }
119 private:
120  Force& m_force;
121  float m_scale;
122  bool m_normalize;
123 };
124 
125 
129 template <typename TCost>
130 double TSSDCost<TCost>::do_evaluate_force(const Data& a, const Data& b, Force& force) const
131 {
132  assert(a.get_size() == b.get_size());
133  assert(a.get_size() == force.get_size());
134  FEvalForce<Force> ef(force, m_normalize);
135  return filter(ef, a, b);
136 }
137 
138 
144 template <typename CP, typename C>
145 class TSSDCostPlugin: public CP {
146 public:
147  TSSDCostPlugin();
148  C *do_create()const;
149 private:
150  bool m_normalize;
151 };
152 
153 
157 template <typename CP, typename C>
158 TSSDCostPlugin<CP,C>::TSSDCostPlugin():
159  CP("ssd"),
160  m_normalize(false)
161 {
162  TRACE("TSSDCostPlugin<CP,C>::TSSDCostPlugin()");
163  this->add_property(::mia::property_gradient);
164  this->add_parameter("norm", new mia::CBoolParameter(m_normalize, false,
165  "Set whether the metric should be normalized by the number of image pixels"));
166 
167 }
168 
172 template <typename CP, typename C>
173 C *TSSDCostPlugin<CP,C>::do_create() const
174 {
175  return new TSSDCost<C>(m_normalize);
176 }
177 
179 NS_END