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
29NS_BEGIN(NS)
30
31
32struct valpos {
33 double val;
34 int pos;
35};
36
37
38template <typename TCost>
39class TLSDImageCost: public TCost
40{
41public:
42 typedef typename TCost::Data Data;
43 typedef typename TCost::Force Force;
44private:
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
90inline bool operator < (const valpos& a, const valpos& b)
91{
92 return a.val < b.val;
93}
94
95template <typename CP, typename C>
96class TLSDImageCostPlugin: public CP
97{
98public:
99 TLSDImageCostPlugin();
100 virtual TLSDImageCost<C> *do_create()const;
101 const std::string do_get_descr()const;
102};
103
104template <typename TCost>
105template <typename DataTempl>
106void 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
148template <typename TCost>
149void 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
156template <typename TCost>
157double 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
163template <typename TCost>
164double 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
170template <typename TCost>
171TLSDImageCost<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
177template <typename TCost>
178template <typename DataTempl>
179double 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
197template <typename TCost>
198template <typename DataTempl>
199double 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
228template <typename CP, typename C>
229TLSDImageCostPlugin<CP, C>::TLSDImageCostPlugin():
230 CP("lsd")
231{
232}
233
234template <typename CP, typename C>
235TLSDImageCost<C> *TLSDImageCostPlugin<CP, C>::do_create()const
236{
237 return new TLSDImageCost<C>();
238}
239
240template <typename CP, typename C>
241const std::string TLSDImageCostPlugin<CP, C>::do_get_descr()const
242{
243 return "Least-Squares Distance measure";
244}
245
247NS_END
EXPORT_2D C2DFVectorfield get_gradient(const C2DImage &image)
bool operator<(const T2DVector< T > &a, const T2DVector< S > &b)
Definition 2d/vector.hh:495
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
R result_type
defines the return type of the filter function