My Project
lsd.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 struct valpos {
33  double val;
34  int pos;
35 };
36 
37 
38 template <typename TCost>
39 class TLSDImageCost: public TCost
40 {
41 public:
42  typedef typename TCost::Data Data;
43  typedef typename TCost::Force Force;
44 private:
45 
46  class CRefPrepare : public mia::TFilter<void>
47  {
48  public:
49  CRefPrepare(std::vector<double>& QtQinv, std::vector<int>& Q_mappping):
50  m_QtQinv(QtQinv),
51  m_Q_mappping(Q_mappping)
52  {
53  }
54  template <typename DataTempl>
55  void operator()(const DataTempl& ref);
56 
57  std::vector<double>& m_QtQinv;
58  std::vector<int>& m_Q_mappping;
59  };
60 
61  class RunCost : mia::TFilter<double>
62  {
63  public:
64  typedef TFilter<double>::result_type result_type;
65  RunCost( const std::vector<double>& QtQinv, const std::vector<int>& Q_mappping);
66 
67  template <typename DataTempl>
68  double operator()(const DataTempl& ref)const;
69 
70  template <typename DataTempl>
71  double operator()(const DataTempl& ref, Force& force)const;
72  private:
73  const std::vector<double>& m_QtQinv;
74  const std::vector<int>& m_Q_mappping;
75  };
76 
77 
78 
79  virtual double do_value(const Data& a, const Data& b) const;
80 
81  virtual double do_evaluate_force(const Data& a, const Data& /*b*/, Force& force) const;
82 
83  virtual void post_set_reference(const Data& ref);
84 
85  std::vector<double> m_QtQinv;
86  std::vector<int> m_Q_mappping;
87 };
88 
89 
90 inline bool operator < (const valpos& a, const valpos& b)
91 {
92  return a.val < b.val;
93 }
94 
95 template <typename CP, typename C>
96 class TLSDImageCostPlugin: public CP
97 {
98 public:
99  TLSDImageCostPlugin();
100  virtual TLSDImageCost<C> *do_create()const;
101  const std::string do_get_descr()const;
102 };
103 
104 template <typename TCost>
105 template <typename DataTempl>
106 void TLSDImageCost<TCost>::CRefPrepare::operator()(const DataTempl& ref)
107 {
108  size_t npixels = ref.get_size().x * ref.get_size().y;
109  std::vector<valpos> buffer(npixels);
110  m_Q_mappping.resize(npixels);
111  m_QtQinv.resize(npixels);
112  int idx = 0;
113 
114  for ( auto x = ref.begin(); x != ref.end(); ++x, ++idx) {
115  double value = *x;
116  valpos vp = {value, idx};
117  buffer[idx] = vp;
118  }
119 
120  std::sort(buffer.begin(), buffer.end());
121  idx = 0;
122  auto x = buffer.begin();
123  double oldv = x->val;
124  ++x;
125  ++m_QtQinv[idx];
126  // build histogram and prepare target mapping
127  m_Q_mappping[x->pos] = idx;
128 
129  while ( x != buffer.end() ) {
130  if (x->val != oldv) {
131  oldv = x->val;
132  ++idx;
133  }
134 
135  ++m_QtQinv[idx];
136  m_Q_mappping[x->pos] = idx;
137  ++x;
138  }
139 
140  ++idx;
141  m_QtQinv.resize(idx);
142  std::transform(m_QtQinv.begin(), m_QtQinv.end(), m_QtQinv.begin(),
143  [](double x) {
144  return 1.0 / x;
145  });
146 }
147 
148 template <typename TCost>
149 void TLSDImageCost<TCost>::post_set_reference(const Data& ref)
150 {
151  CRefPrepare rp(m_QtQinv, m_Q_mappping);
152  mia::accumulate(rp, ref);
153 }
154 
155 
156 template <typename TCost>
157 double TLSDImageCost<TCost>::do_value(const Data& a, const Data& /*b*/) const
158 {
159  RunCost rf(m_QtQinv, m_Q_mappping);
160  return mia::filter(rf, a);
161 }
162 
163 template <typename TCost>
164 double TLSDImageCost<TCost>::do_evaluate_force(const Data& a, const Data& /*b*/, Force& force) const
165 {
166  RunCost rf(m_QtQinv, m_Q_mappping);
167  return mia::filter_and_output(rf, a, force);
168 }
169 
170 template <typename TCost>
171 TLSDImageCost<TCost>::RunCost::RunCost(const std::vector<double>& QtQinv, const std::vector<int>& Q_mappping):
172  m_QtQinv(QtQinv),
173  m_Q_mappping(Q_mappping)
174 {
175 }
176 
177 template <typename TCost>
178 template <typename DataTempl>
179 double TLSDImageCost<TCost>::RunCost::operator()(const DataTempl& a)const
180 {
181  double val1 = 0.0;
182  double val2 = 0.0;
183  std::vector<double> sums(m_QtQinv.size(), 0.0);
184  auto idx = m_Q_mappping.begin();
185 
186  for (auto ia = a.begin(); ia != a.end(); ++ia, ++idx) {
187  val1 += *ia * *ia;
188  sums[*idx] += *ia;
189  }
190 
191  for (size_t i = 0; i < sums.size(); ++i)
192  val2 += sums[i] * sums[i] * m_QtQinv[i];
193 
194  return 0.5 * (val1 - val2);
195 }
196 
197 template <typename TCost>
198 template <typename DataTempl>
199 double TLSDImageCost<TCost>::RunCost::operator()(const DataTempl& a, Force& force)const
200 {
201  double value = 0.0;
202  std::vector<double> sums(m_QtQinv.size(), 0.0);
203  auto idx = m_Q_mappping.begin();
204 
205  for (auto ia = a.begin(); ia != a.end(); ++ia, ++idx) {
206  value += *ia * *ia;
207  sums[*idx] += *ia;
208  }
209 
210  auto s = sums.begin();
211 
212  for (auto q = m_QtQinv.begin(); q != m_QtQinv.end(); ++q, ++s) {
213  value -= *q * *s * *s;
214  *s *= *q;
215  }
216 
217  auto gradient = get_gradient(a);
218  auto iforce = force.begin();
219  auto igrad = gradient.begin();
220  idx = m_Q_mappping.begin();
221 
222  for (auto ia = a.begin(); ia != a.end(); ++ia, ++iforce, ++igrad, ++idx)
223  *iforce = *igrad * (*ia - sums[*idx]);
224 
225  return 0.5 * value;
226 }
227 
228 template <typename CP, typename C>
229 TLSDImageCostPlugin<CP, C>::TLSDImageCostPlugin():
230  CP("lsd")
231 {
232 }
233 
234 template <typename CP, typename C>
235 TLSDImageCost<C> *TLSDImageCostPlugin<CP, C>::do_create()const
236 {
237  return new TLSDImageCost<C>();
238 }
239 
240 template <typename CP, typename C>
241 const std::string TLSDImageCostPlugin<CP, C>::do_get_descr()const
242 {
243  return "Least-Squares Distance measure";
244 }
245 
247 NS_END
get_gradient
EXPORT_2D C2DFVectorfield get_gradient(const C2DImage &image)
parameter.hh
TFilter::result_type
R result_type
defines the return type of the filter function
Definition: core/filter.hh:72
TCost::Data
T Data
typedef for generic programming: The data type used by the cost function
Definition: core/cost.hh:68
filter
static F::result_type filter(const F &f, const B &b)
Definition: core/filter.hh:258
msgstream.hh
filter_and_output
static F::result_type filter_and_output(const F &f, const B &a, O &b)
Definition: core/filter.hh:647
NS_BEGIN
#define NS_BEGIN(NS)
conveniance define to start a namespace
Definition: defines.hh:42
operator<
bool operator<(const T2DVector< T > &a, const T2DVector< S > &b)
Definition: 2d/vector.hh:495
NS_END
#define NS_END
conveniance define to end a namespace
Definition: defines.hh:45
accumulate
static F::result_type accumulate(F &f, const B &data)
Definition: core/filter.hh:371
TCost
The generic cost function interface.
Definition: core/cost.hh:65
property_flags.hh
TCost::Force
V Force
typedef for generic programming: The gradient forca type create by the cost function
Definition: core/cost.hh:71
filter.hh