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