A.I.M.S algorithms


emthreshold.h
Go to the documentation of this file.
1 /* Copyright (C) 2000-2013 CEA
2  *
3  * This software and supporting documentation were developed by
4  * bioPICSEL
5  * CEA/DSV/I²BM/MIRCen/LMN, Batiment 61,
6  * 18, route du Panorama
7  * 92265 Fontenay-aux-Roses
8  * France
9  */
10 
11 /*
12  * EM Threshold operators
13  */
14 #ifndef AIMS_UTILITY_EMTHRESHOLD_H
15 #define AIMS_UTILITY_EMTHRESHOLD_H
16 
18 #include <aims/data/data_g.h>
19 #include <aims/math/math_g.h>
20 #include <aims/io/writer.h>
21 #include <aims/io/reader.h>
24 #include <aims/utility/threshold.h>
25 
26 
27 template <class T, class U>
29 {
30  public:
31  inline
33  T backgd = 0,
34  float foregd = 32767,
35  int classes = 2,
36  int levelindex = 0, /* Level expressed as index */
37  int level2index = 1, /* Level expressed as index */
38  int bins = 0,
39  int binstart = 0, /* First bin considered during thresholding process */
40  int binend = 0 /* Last bin considered during thresholding process */
41  );
42 
43  virtual ~AimsEMThreshold() {}
44 
45  inline
46  std::vector<T> processEMThresholds(const AimsData<T> &image);
47 
49  inline AimsData<T> operator () (const AimsData<T> &image);
50 
52  inline
53  AimsData<U> bin(const AimsData<T> &image);
54 
55  private :
56  int _classes;
57  int _levelindex;
58  int _level2index;
59  int _bins;
60  int _binstart;
61  int _binend;
62 };
63 
64 template <class T, class U> inline
66  T backgd, float foregd,
67  int classes,
68  int levelindex, int levelindex2,
69  int bins,
70  int binstart,
71  int binend )
72  : _classes(classes),
73  _levelindex(levelindex),
74  _level2index(levelindex2),
75  _bins(bins),
76  _binstart(binstart),
77  _binend(binend),
78  AimsThreshold<T, U>( type, 0, 0, backgd, foregd )
79 {
80  if (carto::verbose)
81  {
82  std::cout << "Type " << carto::toString(this->_type) << " ..." << std::endl << std::flush;
83  std::cout << "Backgd " << carto::toString(this->_backgd) << " ..." << std::endl << std::flush;
84  std::cout << "Foregd " << carto::toString(this->_foregd) << " ..." << std::endl << std::flush;
85  std::cout << "Classes " << carto::toString(_classes) << " ..." << std::endl << std::flush;
86  std::cout << "Level index " << carto::toString(_levelindex) << " ..." << std::endl << std::flush;
87  std::cout << "Level 2 index " << carto::toString(_level2index) << " ..." << std::endl << std::flush;
88 // std::cout << "Bins " << carto::toString(_bins) << " ..." << std::endl << std::flush;
89  }
90 }
91 // EM
92 template <class T, class U> inline
94 {
95  std::vector<T> thresholds, thresholds_mu;
96  int bins = _bins;
97  int binend;
98 
99 
100  T maxi = image.maximum(), mini = image.minimum();
101  if ( bins <= 0 )
102  bins = int(maxi - mini);
103 
104  if (!_binend)
105  binend = bins;
106  else
107  binend = _binend + 1;
108 
109  if (carto::verbose)
110  std::cout << "Bins " << carto::toString(bins) << " ..." << std::endl << std::flush;
111 
112  // Histogram
114  his.doit(image);
115 
116  // Declarations
117  AimsData<double> th_pi(_classes);
118  AimsData<double> th_mu(_classes);
119  AimsData<double> th_sigma(_classes);
120  AimsData<double> th_tau(bins, _classes);
121  AimsData<double> aux(_classes);
122  int n, i, j;
123 
124  // Initialisation of the EM algorithm
125  double scl = (double) bins / (double)( maxi - mini );
126 // std::cout << "scl : " << carto::toString(scl) << std::endl << std::flush;
127  double points_count = 0.;
128  for (n = _binstart; n < binend; ++n)
129  points_count += his.data()(n);
130 
131 // std::cout << "points_count : " << carto::toString(points_count) << std::endl << std::flush;
132  // Same number of points for each EM class
133  double class_points_count = 0.;
134  for (n = _binstart; n < binend; ++n)
135  {
136  class_points_count += his.data()(n);
137  for (i = 0; i < _classes; ++i)
138  {
139  if ((i <= (_classes * class_points_count / points_count))
140  && (i >= (_classes * class_points_count / points_count) - 1))
141  th_tau(n, i) = 1;
142  else
143  th_tau(n, i) = 0;
144  }
145  }
146 
147  // MAIN interation
148  int count = 0;
149  double sum1 = 0, sum2 = 0, sum = 0, sum_aux;
150  while ( count <= (_classes + 1) )
151  {
152  count++;
153 
154  // EM STEP E : calculate expect classification variable tau for each class
155  // the first time jump this step to initialize mu, sigma and pi.
156  if (count > 1)
157  {
158  th_tau = 0.;
159  for (n = _binstart; n < binend; ++n)
160  {
161  for (i = 0; i < _classes; ++i)
162  {
163  aux = 0.;
164  sum_aux = 0.;
165  for (j = 0; j < _classes; ++j)
166  {
167  aux(j) = th_pi(j) * (1 / sqrt(2 * M_PI * th_sigma(j)))
168  * exp(-(1 / (2 * th_sigma(j))) * (n - th_mu(j)) * (n - th_mu(j)));
169  sum_aux += aux(j);
170  }
171  th_tau(n, i) = aux(i) / sum_aux;
172  }
173  }
174  }
175 
176  // EM STEP M : calculate new threshold means and variances for classified data
177  th_mu = 0.; th_sigma = 0.; th_pi = 0.;
178 
179  for (j = 0; j < _classes; ++j)
180  {
181  // Mean and covariance
182  sum1 = 0.;
183  sum2 = 0.;
184  sum = 0.;
185  for (n = _binstart; n < binend; ++n) {
186  sum += th_tau(n, j) * his.data()(n);
187  sum1 += th_tau(n, j) * n * his.data()(n);
188  }
189  th_mu(j) = sum1 / sum;
190 
191  for (n = _binstart; n < binend; ++n)
192  sum2 += th_tau(n, j) * (n - th_mu(j)) * (n - th_mu(j)) * his.data()(n);
193  th_sigma(j) = sum2 / sum;
194 
195  // Classifying variable pi
196  sum1 = 0.;
197  sum = 0.;
198  for (n = _binstart; n < binend; ++n)
199  {
200  sum1 += th_tau(n, j) * his.data()(n);
201  sum += his.data()(n);
202  }
203  th_pi(j) = sum1 / sum;
204  }
205  }
206 
207  // Compute thresholds
208  T level;
209  for (j = 0; j < _classes; ++j)
210  {
211  level = (T)th_mu(j);
212 
213  if (level < 0)
214  level = 0;
215 
216  else if (level > (binend - 1))
217  level = binend - 1;
218 
219  thresholds_mu.push_back(level);
220  }
221 
222  std::sort(thresholds_mu.begin(), thresholds_mu.end());
223 
224  if (carto::verbose) {
225  std::cout << "Found class means : [";
226 
227  for (i = 0; i < _classes; ++i)
228  {
229  // Search for last minimum between class means
230  std::cout << carto::toString((T)(thresholds_mu[i] / scl) + mini);
231  if (i < (_classes - 1))
232  std::cout << ", ";
233  }
234  std::cout << "]" << std::endl << std::flush;
235  }
236 
237  if (carto::verbose)
238  std::cout << "Found thresholds : [";
239 
240  for (i = 0; i < (_classes - 1); ++i)
241  {
242  // Search for last minimum between class means
243  level = thresholds_mu[i];
244  for(int j = thresholds_mu[i]; j <= thresholds_mu[i + 1]; j++)
245  {
246  if (his.data()(j) < his.data()(level))
247  level = j;
248  }
249 
250  thresholds.push_back(((T)(level / scl) + mini));
251 
252  if (carto::verbose) {
253  std::cout << carto::toString(thresholds[i]);
254  if (i < (_classes - 2))
255  std::cout << ", ";
256  }
257  }
258 
259  if (carto::verbose)
260  std::cout << "]" << std::endl << std::flush;
261 
262  return thresholds;
263 }
264 
265 template <class T, class U> inline
267 {
268  // Threshold
269  std::vector<T> thresholds = processEMThresholds(image);
270  if (_levelindex < thresholds.size())
271  AimsEMThreshold<T, U>::_level = thresholds[_levelindex];
272 
273  if (_level2index < thresholds.size())
274  AimsEMThreshold<T, U>::_level2 = thresholds[_level2index];
275 
276 // std::cout << "Thresholding values ..." << std::endl << std::flush;
277  AimsData<U> thresholded = AimsThreshold<T, U>::operator ()( image );
278 
279  return thresholded;
280 }
281 
282 template <class T, class U> inline
284 {
285  // Threshold
286  std::vector<T> thresholds = processEMThresholds(image);
287  if (_levelindex < thresholds.size())
288  AimsEMThreshold<T, U>::_level = thresholds[_levelindex];
289 
290  if (_level2index < thresholds.size())
291  AimsEMThreshold<T, U>::_level2 = thresholds[_level2index];
292 
293  AimsData<U> thresholded = AimsThreshold<T, U>::bin( image );
294 
295  return thresholded;
296 }
297 
298 #endif
T minimum() const
AimsData< U > bin(const AimsData< T > &image)
Return the binary thresholded image.
Definition: emthreshold.h:283
AimsEMThreshold(threshold_t type, T backgd=0, float foregd=32767, int classes=2, int levelindex=0, int level2index=1, int bins=0, int binstart=0, int binend=0)
Definition: emthreshold.h:65
Histogram container class, with a specified number of regular bins.
AimsData< int32_t > & data()
return a reference to the data field of the histogram class.
Definition: histogram.h:66
int verbose
T maximum() const
std::vector< T > processEMThresholds(const AimsData< T > &image)
Definition: emthreshold.h:93
AimsData< T > operator()(const AimsData< T > &image)
Return the multi-level thresholded image.
Definition: emthreshold.h:266
AimsData< U > bin(const AimsData< T > &sqv)
void doit(const AimsData< T > &thing)
classical histogram computation function.
virtual ~AimsEMThreshold()
Definition: emthreshold.h:43
threshold_t
AimsData< T > operator()(const AimsData< T > &sqv)
std::string toString(const T &object)