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{
32public:
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
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
109private:
110 PInitializer m_cci;
111 struct CMeansImpl *impl;
112
113};
114
143template <typename T, template <class > class Field>
144void 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
233template <typename T, template <class> class Field>
234double 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{
282public:
284protected:
285 size_t get_size_param() const;
286private:
287 size_t m_size;
288
289};
290
292#ifdef __GNUC__
293#pragma GCC diagnostic push
294#pragma GCC diagnostic ignored "-Wattributes"
295#endif
297extern 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
316
CMeansInitializerSizedPlugin(const char *name)
Initializer plugin_data
Definition: cmeans.hh:93
Initializer plugin_type
Definition: cmeans.hh:94
virtual DVector run(const NormalizedHistogram &nh) const =0
static const char * type_descr
Definition: cmeans.hh:96
static const char * data_descr
Definition: cmeans.hh:95
bool save(const std::string &filename) const
size_t size() const
Definition: cmeans.hh:80
SparseProbmap(const std::string &filename)
std::vector< value_type > Map
Definition: cmeans.hh:41
DVector get_fuzzy(double x) const
Map::iterator begin()
Definition: cmeans.hh:61
Map::iterator end()
Definition: cmeans.hh:70
Map::const_iterator end() const
Definition: cmeans.hh:66
SparseProbmap(size_t size)
Map::const_iterator begin() const
Definition: cmeans.hh:57
std::pair< unsigned short, DVector > value_type
Definition: cmeans.hh:40
Definition: cmeans.hh:31
CMeans(double epsilon, PInitializer class_center_initializer)
CSparseHistogram::Compressed SparseHistogram
Definition: cmeans.hh:34
std::vector< double > DVector
Definition: cmeans.hh:33
std::shared_ptr< Initializer > PInitializer
Definition: cmeans.hh:99
SparseProbmap run(const SparseHistogram &histogram, DVector &class_centers) const
SparseProbmap run(const SparseHistogram &histogram, DVector &class_centers, bool de_normalize_results) const
std::vector< std::pair< double, double > > NormalizedHistogram
Definition: cmeans.hh:35
The base class for all plug-in created object.
Definition: product_base.hh:41
std::vector< std::pair< int, unsigned long > > Compressed
the Base class for all plugn handlers that deal with factory plugins.
Definition: factory.hh:95
This is tha base of all plugins that create "things", like filters, cost functions time step operator...
Definition: factory.hh:51
the singleton that a plug-in handler really is
Definition: handler.hh:159
The basic template of all plugin handlers.
Definition: handler.hh:57
The generic base for all plug-ins.
Definition: plugin_base.hh:172
THandlerSingleton< TFactoryPluginHandler< CMeansInitializerPlugin > > CMeansInitializerPluginHandler
Definition: cmeans.hh:308
TFactory< CMeans::Initializer > CMeansInitializerPlugin
Definition: cmeans.hh:277
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
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
#define NS_MIA_BEGIN
conveniance define to start the mia namespace
Definition: defines.hh:33
#define EXPORT_CORE
Macro to manage Visual C++ style dllimport/dllexport.
Definition: defines.hh:101
#define NS_MIA_END
conveniance define to end the mia namespace
Definition: defines.hh:36
vstream & cvwarn()
send warnings to this stream adapter
Definition: msgstream.hh:321
#define FACTORY_TRAIT(F)