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 
30 NS_BEGIN(NS)
31 
32 
34 template <typename T>
35 class TMIImageCost: public T
36 {
37 public:
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);
42 private:
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 
51 struct 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 
67 template <typename T>
68 TMIImageCost<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 
76 template <typename T>
77 double 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 
83 template <typename Force>
84 struct 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  }
105 private:
106  Force& m_force;
107  mia::CSplineParzenMI& m_parzen_mi;
108 };
109 
110 
114 template <typename T>
115 double 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 
123 template <typename T>
124 void TMIImageCost<T>::post_set_reference(const Data& MIA_PARAM_UNUSED(ref))
125 {
126  m_parzen_mi.reset();
127 }
128 
134 template <typename CP, typename C>
135 class TMIImageCostPlugin: public CP
136 {
137 public:
138  TMIImageCostPlugin();
139  C *do_create()const;
140 private:
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 
153 template <typename CP, typename C>
154 TMIImageCostPlugin<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 
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 
191 NS_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)
Definition: core/filter.hh:258
std::shared_ptr< CSplineKernel > PSplineKernel
#define TRACE(DOMAIN)
Definition: msgstream.hh:208
CParameter * make_param(T &value, bool required, const char *descr)
Definition: parameter.hh:280
CParameter * make_ci_param(T &value, S1 lower_bound, S2 upper_bound, bool required, const char *descr)
Definition: parameter.hh:329
EXPORT_CORE const char * property_gradient
constant defining the gradient property