My Project
cmeans.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 
22 #include <mia/core/probmap.hh>
24 #include <mia/core/factory.hh>
25 
27 
28 
29 
31 {
32 public:
33  typedef std::vector<double> DVector;
35  typedef std::vector<std::pair<double, double>> NormalizedHistogram;
36 
38  {
39  public:
40  typedef std::pair<unsigned short, DVector> value_type;
41  typedef std::vector<value_type> Map;
42 
43  SparseProbmap() = delete;
44  SparseProbmap (size_t size);
45  SparseProbmap (const std::string& filename);
46 
47  value_type& operator [](int i)
48  {
49  return m_map[i];
50  }
51 
52  const value_type& operator [](int i) const
53  {
54  return m_map[i];
55  }
56 
57  Map::const_iterator begin() const
58  {
59  return m_map.begin();
60  }
61  Map::iterator begin()
62  {
63  return m_map.begin();
64  }
65 
66  Map::const_iterator end() const
67  {
68  return m_map.end();
69  }
70  Map::iterator end()
71  {
72  return m_map.end();
73  }
74 
75  bool save(const std::string& filename) const;
76 
77 
78  DVector get_fuzzy(double x) const;
79 
80  size_t size() const
81  {
82  return m_map.size();
83  }
84  private:
85  Map m_map;
86 
87  };
88 
89 
91  {
92  public:
95  static const char *data_descr;
96  static const char *type_descr;
97  virtual DVector run(const NormalizedHistogram& nh) const = 0;
98  };
99  typedef std::shared_ptr<Initializer> PInitializer;
100 
101  CMeans(double epsilon, PInitializer class_center_initializer);
102 
103  ~CMeans();
104 
105  SparseProbmap run(const SparseHistogram& histogram, DVector& class_centers) const;
106 
107  SparseProbmap run(const SparseHistogram& histogram, DVector& class_centers, bool de_normalize_results) const;
108 
109 private:
110  PInitializer m_cci;
111  struct CMeansImpl *impl;
112 
113 };
114 
143 template <typename T, template <class > class Field>
144 void cmeans_evaluate_probabilities(const Field<T>& image, const Field<float>& gain,
145  const std::vector<double>& class_centers,
146  std::vector<Field<float>>& pv)
147 {
148  assert(image.size() == gain.size());
149  assert(class_centers.size() == pv.size());
150 #ifndef NDEBUG
151 
152  for (auto i : pv)
153  assert(image.size() == i.size());
154 
155 #endif
156  auto ii = image.begin();
157  auto ie = image.end();
158  auto ig = gain.begin();
159  typedef typename Field<float>::iterator prob_iterator;
160  std::vector<prob_iterator> ipv(pv.size());
161  transform(pv.begin(), pv.end(), ipv.begin(), [](Field<float>& p) {
162  return p.begin();
163  });
164  std::vector<double> gain_class_centers(class_centers.size());
165 
166  while (ii != ie) {
167  double x = *ii;
168 
169  for (auto iipv : ipv)
170  *iipv = 0.0;
171 
172  const double vgain = *ig;
173  transform(class_centers.begin(), class_centers.end(), gain_class_centers.begin(),
174  [vgain](double x) {
175  return vgain * x;
176  });
177 
178  if ( x < gain_class_centers[0]) {
179  *ipv[0] = 1.0;
180  } else {
181  unsigned j = 1;
182  bool value_set = false;
183 
184  while (!value_set && (j < class_centers.size()) ) {
185  // between two centers
186  if (x < gain_class_centers[j]) {
187  double p0 = x - gain_class_centers[j - 1];
188  double p1 = x - gain_class_centers[j];
189  double p02 = p0 * p0;
190  double p12 = p1 * p1;
191  double normalizer = 1.0 / (p02 + p12);
192  *ipv[j] = p02 * normalizer;
193  *ipv[j - 1] = p12 * normalizer;
194  value_set = true;
195  }
196 
197  ++j;
198  }
199 
200  if (!value_set)
201  *ipv[class_centers.size() - 1] = 1.0;
202  }
203 
204  ++ii;
205  ++ig;
206 
207  for (unsigned i = 0; i < class_centers.size(); ++i)
208  ++ipv[i];
209  }
210 }
211 
233 template <typename T, template <class> class Field>
234 double cmeans_update_class_centers(const Field<T>& image, const Field<float>& gain,
235  const std::vector<Field<float>>& pv,
236  std::vector<double>& class_centers)
237 {
238  double residuum = 0.0;
239 
240  for (size_t i = 0; i < class_centers.size(); ++i) {
241  double cc = class_centers[i];
242  double sum_prob = 0.0;
243  double sum_weight = 0.0;
244  auto ie = image.end();
245  auto ii = image.begin();
246  auto ig = gain.begin();
247  auto ip = pv[i].begin();
248 
249  while (ii != ie) {
250  if (*ip > 0.0) {
251  auto v = *ip * *ip * *ig;
252  sum_prob += v * *ig;
253  sum_weight += v * *ii;
254  }
255 
256  ++ii;
257  ++ig;
258  ++ip;
259  }
260 
261  if (sum_prob != 0.0) // move slowly in the direction of new center
262  cc = sum_weight / sum_prob;
263  else {
264  cvwarn() << "class[" << i << "] has no probable members, keeping old value:" <<
265  sum_prob << ":" << sum_weight << "\n";
266  }
267 
268  double delta = (cc - class_centers[i]) * 0.5;
269  residuum += delta * delta;
270  class_centers[i] += delta;
271  }// end update class centers
272 
273  return sqrt(residuum);
274 }
275 
276 
278 
279 // the class that has only the size as a parameter
281 {
282 public:
283  CMeansInitializerSizedPlugin(const char *name);
284 protected:
285  size_t get_size_param() const;
286 private:
287  size_t m_size;
288 
289 };
290 
292 #ifdef __GNUC__
293 #pragma GCC diagnostic push
294 #pragma GCC diagnostic ignored "-Wattributes"
295 #endif
297 extern template class EXPORT_CORE TFactory<CMeans::Initializer>;
298 #ifdef __GNUC__
299 #pragma GCC diagnostic pop
300 #endif
304 
305  template <> const char *const TPluginHandler<TFactory<CMeans::Initializer>>::m_help;
306 
307 
309 
313 
314 
315  NS_MIA_END
316 
cmeans_evaluate_probabilities
void cmeans_evaluate_probabilities(const Field< T > &image, const Field< float > &gain, const std::vector< double > &class_centers, std::vector< Field< float >> &pv)
evaluate the probabilities for a c-means classification with gain field
Definition: cmeans.hh:144
CMeans::SparseProbmap::end
Map::const_iterator end() const
Definition: cmeans.hh:66
factory.hh
CMeans::PInitializer
std::shared_ptr< Initializer > PInitializer
Definition: cmeans.hh:99
CMeansInitializerPluginHandler
THandlerSingleton< TFactoryPluginHandler< CMeansInitializerPlugin > > CMeansInitializerPluginHandler
Definition: cmeans.hh:308
NS_MIA_BEGIN
#define NS_MIA_BEGIN
conveniance define to start the mia namespace
Definition: defines.hh:33
cmeans_update_class_centers
double cmeans_update_class_centers(const Field< T > &image, const Field< float > &gain, const std::vector< Field< float >> &pv, std::vector< double > &class_centers)
Definition: cmeans.hh:234
CMeans::SparseProbmap
Definition: cmeans.hh:37
CMeansInitializerPlugin
TFactory< CMeans::Initializer > CMeansInitializerPlugin
Definition: cmeans.hh:277
CMeans::Initializer::plugin_data
Initializer plugin_data
Definition: cmeans.hh:93
CMeans::SparseProbmap::size
size_t size() const
Definition: cmeans.hh:80
NS_MIA_END
#define NS_MIA_END
conveniance define to end the mia namespace
Definition: defines.hh:36
CMeans::SparseProbmap::value_type
std::pair< unsigned short, DVector > value_type
Definition: cmeans.hh:40
THandlerSingleton
the singleton that a plug-in handler really is
Definition: handler.hh:158
CSparseHistogram::Compressed
std::vector< std::pair< int, unsigned long > > Compressed
Definition: sparse_histogram.hh:42
TPluginHandler
The basic template of all plugin handlers.
Definition: handler.hh:56
sparse_histogram.hh
TPlugin
The generic base for all plug-ins.
Definition: plugin_base.hh:171
cvwarn
vstream & cvwarn()
send warnings to this stream adapter
Definition: msgstream.hh:321
CMeans::SparseHistogram
CSparseHistogram::Compressed SparseHistogram
Definition: cmeans.hh:34
CProductBase
The base class for all plug-in created object.
Definition: product_base.hh:40
CMeans::SparseProbmap::Map
std::vector< value_type > Map
Definition: cmeans.hh:41
CMeans::Initializer::type_descr
static const char * type_descr
Definition: cmeans.hh:96
probmap.hh
CMeansInitializerSizedPlugin
Definition: cmeans.hh:280
TFactoryPluginHandler
the Base class for all plugn handlers that deal with factory plugins.
Definition: factory.hh:94
EXPORT_CORE
#define EXPORT_CORE
Macro to manage Visual C++ style dllimport/dllexport.
Definition: defines.hh:101
CMeans::Initializer::plugin_type
Initializer plugin_type
Definition: cmeans.hh:94
CMeans::SparseProbmap::end
Map::iterator end()
Definition: cmeans.hh:70
CMeans::DVector
std::vector< double > DVector
Definition: cmeans.hh:33
FACTORY_TRAIT
#define FACTORY_TRAIT(F)
Definition: factory_trait.hh:69
CMeans
Definition: cmeans.hh:30
CMeans::Initializer
Definition: cmeans.hh:90
CMeans::Initializer::data_descr
static const char * data_descr
Definition: cmeans.hh:95
CMeans::NormalizedHistogram
std::vector< std::pair< double, double > > NormalizedHistogram
Definition: cmeans.hh:35
CMeans::SparseProbmap::begin
Map::const_iterator begin() const
Definition: cmeans.hh:57
CMeans::SparseProbmap::begin
Map::iterator begin()
Definition: cmeans.hh:61
TFactory
This is tha base of all plugins that create "things", like filters, cost functions time step operator...
Definition: factory.hh:49