My Project
mi.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>
26
27#include <numeric>
28#include <limits>
29
30NS_BEGIN(NS)
31
32
33
34template <typename T>
35class TMIImageCost: public T
36{
37public:
38 typedef typename T::Data Data;
39 typedef typename T::Force Force;
40
41 TMIImageCost(size_t fbins, mia::PSplineKernel fkernel, size_t rbins, mia::PSplineKernel rkernel, double cut);
42private:
43 virtual double do_value(const Data& a, const Data& b) const;
44 virtual double do_evaluate_force(const Data& a, const Data& b, Force& force) const;
45 virtual void post_set_reference(const Data& ref);
46 mutable mia::CSplineParzenMI m_parzen_mi;
47
48};
49
50
51struct FEvalMI : public mia::TFilter<double> {
52 FEvalMI( mia::CSplineParzenMI& parzen_mi):
53 m_parzen_mi(parzen_mi)
54 {}
55
56
57 template <typename T, typename R>
58 FEvalMI::result_type operator () (const T& a, const R& b) const
59 {
60 m_parzen_mi.fill(a.begin(), a.end(), b.begin(), b.end());
61 return m_parzen_mi.value();
62 }
63 mia::CSplineParzenMI& m_parzen_mi;
64};
65
66
67template <typename T>
68TMIImageCost<T>::TMIImageCost(size_t rbins, mia::PSplineKernel rkernel, size_t mbins,
69 mia::PSplineKernel mkernel, double cut):
70 m_parzen_mi(rbins, rkernel, mbins, mkernel, cut)
71
72{
73 this->add(::mia::property_gradient);
74}
75
76template <typename T>
77double TMIImageCost<T>::do_value(const Data& a, const Data& b) const
78{
79 FEvalMI essd(m_parzen_mi);
80 return filter(essd, a, b);
81}
82
83template <typename Force>
84struct FEvalForce: public mia::TFilter<float> {
85 FEvalForce(Force& force, mia::CSplineParzenMI& parzen_mi):
86 m_force(force),
87 m_parzen_mi(parzen_mi)
88 {
89 }
90 template <typename T, typename R>
91 float operator ()( const T& a, const R& b) const
92 {
93 Force gradient = get_gradient(a);
94 m_parzen_mi.fill(a.begin(), a.end(), b.begin(), b.end());
95 typename T::const_iterator ai = a.begin();
96 typename R::const_iterator bi = b.begin();
97
98 for (size_t i = 0; i < a.size(); ++i, ++ai, ++bi) {
99 float delta = -m_parzen_mi.get_gradient_slow(*ai, *bi);
100 m_force[i] = gradient[i] * delta;
101 }
102
103 return m_parzen_mi.value();
104 }
105private:
106 Force& m_force;
107 mia::CSplineParzenMI& m_parzen_mi;
108};
109
110
114template <typename T>
115double TMIImageCost<T>::do_evaluate_force(const Data& a, const Data& b, Force& force) const
116{
117 assert(a.get_size() == b.get_size());
118 assert(a.get_size() == force.get_size());
119 FEvalForce<Force> ef(force, m_parzen_mi);
120 return filter(ef, a, b);
121}
122
123template <typename T>
124void TMIImageCost<T>::post_set_reference(const Data& MIA_PARAM_UNUSED(ref))
125{
126 m_parzen_mi.reset();
127}
128
134template <typename CP, typename C>
135class TMIImageCostPlugin: public CP
136{
137public:
138 TMIImageCostPlugin();
139 C *do_create()const;
140private:
141 const std::string do_get_descr() const;
142 unsigned int m_rbins;
143 unsigned int m_mbins;
144 mia::PSplineKernel m_mkernel;
145 mia::PSplineKernel m_rkernel;
146 float m_histogram_cut;
147};
148
149
153template <typename CP, typename C>
154TMIImageCostPlugin<CP, C>::TMIImageCostPlugin():
155 CP("mi"),
156 m_rbins(64),
157 m_mbins(64),
158 m_histogram_cut(0.0)
159{
160 TRACE("TMIImageCostPlugin<CP,C>::TMIImageCostPlugin()");
161 this->add_property(::mia::property_gradient);
162 this->add_parameter("rbins", mia::make_ci_param(m_rbins, 1u, 256u, false,
163 "Number of histogram bins used for the reference image"));
164 this->add_parameter("mbins", mia::make_ci_param(m_mbins, 1u, 256u, false,
165 "Number of histogram bins used for the moving image"));
166 this->add_parameter("rkernel", mia::make_param(m_rkernel, "bspline:d=0", false,
167 "Spline kernel for reference image parzen hinstogram"));
168 this->add_parameter("mkernel", mia::make_param(m_mkernel, "bspline:d=3", false,
169 "Spline kernel for moving image parzen hinstogram"));
170 this->add_parameter("cut", mia::make_ci_param(m_histogram_cut, 0.0f, 40.0f,
171 false, "Percentage of pixels to cut at high and low "
172 "intensities to remove outliers"));
173}
174
178template <typename CP, typename C>
179C *TMIImageCostPlugin<CP, C>::do_create() const
180{
181 return new TMIImageCost<C>(m_rbins, m_rkernel, m_mbins, m_mkernel, m_histogram_cut);
182}
183
184template <typename CP, typename C>
185const std::string TMIImageCostPlugin<CP, C>::do_get_descr() const
186{
187 return "Spline parzen based mutual information.";
188}
189
191NS_END
EXPORT_2D C2DFVectorfield get_gradient(const C2DImage &image)
#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)
#define TRACE(DOMAIN)
Definition msgstream.hh:208