My Project
histogram.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 #ifndef mia_core_histogram_hh
22 #define mia_core_histogram_hh
23 
24 #include <mia/core/defines.hh>
25 
26 #include <cmath>
27 #include <cassert>
28 #include <vector>
29 #include <mia/core/defines.hh>
30 #include <mia/core/msgstream.hh>
31 #include <boost/type_traits.hpp>
32 
36 
48 template <typename T>
50 {
51 public:
53  typedef T value_type;
54 
62  THistogramFeeder(T min, T max, size_t bins);
63 
65  size_t size() const;
66 
73  size_t index(T x) const;
74 
80  T value(size_t k) const;
81 private:
82  T m_min;
83  T m_max;
84  size_t m_bins;
85  double m_step;
86  double m_inv_step;
87 };
88 
99 template <>
100 class THistogramFeeder<unsigned char>
101 {
102 public:
104  typedef unsigned char value_type;
105 
107  THistogramFeeder(unsigned char min, unsigned char max, size_t bins);
108 
110  size_t size() const;
111 
113  size_t index(unsigned char x) const;
114 
116  unsigned char value(size_t k) const;
117 };
118 
121 
133 template <typename Feeder>
135 {
136 public:
137 
139  typedef std::vector<size_t>::const_iterator const_iterator;
140 
142  typedef std::pair<typename Feeder::value_type, size_t> value_type;
143 
144  typedef std::pair<typename Feeder::value_type, typename Feeder::value_type> range_type;
145 
149  THistogram(const Feeder& f);
150 
160  THistogram(const THistogram<Feeder>& org, double perc);
161 
166  void push(typename Feeder::value_type x);
167 
173  void push(typename Feeder::value_type x, size_t count);
174 
181  template <typename Iterator>
182  void push_range(Iterator begin, Iterator end);
183 
185  size_t size() const;
186 
188  const_iterator begin() const;
189 
191  const_iterator end() const;
192 
194  size_t operator [] (size_t idx) const;
195 
201  const value_type at(size_t idx) const;
202 
203 
205  typename Feeder::value_type median() const;
206 
208  typename Feeder::value_type MAD() const;
209 
211  double average() const;
212 
214  double deviation() const;
215 
217  double excess_kurtosis() const;
218 
220  double skewness() const;
221 
228  range_type get_reduced_range(double remove) const;
229 private:
230  Feeder m_feeder;
231  std::vector<size_t> m_histogram;
232  size_t m_n;
233 };
234 
235 // inline inplementation
236 template <typename T>
237 THistogramFeeder<T>::THistogramFeeder(T min, T max, size_t bins):
238  m_min(min),
239  m_max(max),
240  m_bins(bins),
241  m_step(( double(max) - double(min) ) / double(bins - 1)),
242  m_inv_step(double(bins - 1) / (double(max) - double(min)))
243 {
244 }
245 
246 template <typename T>
248 {
249  return m_bins;
250 }
251 
252 template <typename T>
253 inline size_t THistogramFeeder<T>::index(T x) const
254 {
255  double val = floor(m_inv_step * (x - m_min) + 0.5);
256 
257  if (val < 0)
258  return 0;
259 
260  if (val < m_bins)
261  return val;
262 
263  return m_bins - 1;
264 }
265 
266 template <typename T>
267 T THistogramFeeder<T>::value(size_t k) const
268 {
269  return k * m_step + m_min;
270 }
271 
272 inline THistogramFeeder<unsigned char>::THistogramFeeder(unsigned char /*min*/, unsigned char /*max*/, size_t /*bins*/)
273 {
274 }
275 
277 {
278  return 256;
279 }
280 
281 inline
282 size_t THistogramFeeder<unsigned char>::index(unsigned char x) const
283 {
284  return x;
285 }
286 
287 inline
288 unsigned char THistogramFeeder<unsigned char>::value(size_t k) const
289 {
290  return k;
291 }
292 
293 template <typename Feeder>
295  m_feeder(f),
296  m_histogram(f.size()),
297  m_n(0)
298 {
299 }
300 
301 template <typename Feeder>
303  m_feeder(org.m_feeder),
304  m_histogram(m_feeder.size()),
305  m_n(0)
306 {
307  size_t n = (size_t)(org.m_n * (1.0 - perc));
308  size_t i = 0;
309 
310  while (n > m_n && i < m_histogram.size()) {
311  m_n += org.m_histogram[i];
312  m_histogram[i] = org.m_histogram[i];
313  ++i;
314  }
315 }
316 
317 
318 template <typename Feeder>
320 {
321  return m_histogram.size();
322 }
323 
324 template <typename Feeder>
325 void THistogram<Feeder>::push(typename Feeder::value_type x)
326 {
327  ++m_n;
328  ++m_histogram[m_feeder.index(x)];
329 }
330 
331 template <typename Feeder>
332 template <typename Iterator>
333 void THistogram<Feeder>::push_range(Iterator begin, Iterator end)
334 {
335  while (begin != end)
336  push(*begin++);
337 }
338 
339 
340 
341 template <typename Feeder>
342 void THistogram<Feeder>::push(typename Feeder::value_type x, size_t count)
343 {
344  m_n += count;
345  m_histogram[m_feeder.index(x)] += count;
346 }
347 
348 template <typename Feeder>
350 {
351  return m_histogram.begin();
352 }
353 
354 template <typename Feeder>
356 {
357  return m_histogram.end();
358 }
359 
360 template <typename Feeder>
361 size_t THistogram<Feeder>::operator [] (size_t idx) const
362 {
363  assert(idx < m_histogram.size());
364  return m_histogram[idx];
365 }
366 
367 template <typename Feeder>
368 typename Feeder::value_type THistogram<Feeder>::median() const
369 {
370  float n_2 = m_n / 2.0f;
371  float sum = 0;
372  size_t k = 0;
373 
374  while ( sum < n_2 )
375  sum += m_histogram[k++];
376 
377  return m_feeder.value(k > 0 ? k - 1 : k);
378 }
379 
380 template <typename Feeder>
381 typename Feeder::value_type THistogram<Feeder>::MAD() const
382 {
383  typedef typename Feeder::value_type T;
384  T m = median();
385  THistogram<Feeder> help(m_feeder);
386  ;
387 
388  for (size_t k = 0; k < size(); ++k) {
389  T v = m_feeder.value(k);
390  help.push(v > m ? v - m : m - v, m_histogram[k]);
391  }
392 
393  return help.median();
394 }
395 
396 
397 template <typename Feeder>
399 {
400  if (idx < m_histogram.size())
401  return value_type(m_feeder.value(idx), m_histogram[idx]);
402  else
403  return value_type(m_feeder.value(idx), 0);
404 }
405 
406 template <typename Feeder>
408 {
409  if (m_n < 1)
410  return 0.0;
411 
412  double sum = 0.0;
413 
414  for (size_t i = 0; i < size(); ++i) {
415  const typename THistogram<Feeder>::value_type value = at(i);
416  sum += value.first * value.second;
417  }
418 
419  return sum / m_n;
420 }
421 
422 template <typename Feeder>
424 {
425  double mu = average();
426  double sum1 = 0.0;
427  double sum2 = 0.0;
428  double sum3 = 0.0;
429 
430  for (size_t i = 0; i < size(); ++i) {
431  const auto value = at(i);
432  sum1 += value.first * value.second;
433  sum2 += value.first * value.first * value.second;
434  double h = (value.first - mu);
435  h *= h;
436  h *= h;
437  sum3 += h * value.second;
438  }
439 
440  double sigma2 = (sum2 - sum1 * sum1 / m_n) / m_n;
441  double mu4 = sum3 / m_n;
442  return mu4 / (sigma2 * sigma2) - 3;
443 }
444 
445 template <typename Feeder>
447 {
448  double mu = average();
449  double sum1 = 0.0;
450  double sum2 = 0.0;
451  double sum3 = 0.0;
452 
453  for (size_t i = 0; i < size(); ++i) {
454  const auto value = at(i);
455  sum1 += value.first * value.second;
456  sum2 += value.first * value.first * value.second;
457  double h = (value.first - mu);
458  sum3 += h * h * h * value.second;
459  }
460 
461  double sigma2 = (sum2 - sum1 * sum1 / m_n) / m_n;
462  double mu4 = sum3 / m_n;
463  return mu4 / (sigma2 * sqrt(sigma2));
464 }
465 
466 template <typename Feeder>
468 {
469  if (m_n < 2)
470  return 0.0;
471 
472  double sum = 0.0;
473  double sum2 = 0.0;
474 
475  for (size_t i = 0; i < size(); ++i) {
476  const typename THistogram<Feeder>::value_type value = at(i);
477  sum += value.first * value.second;
478  sum2 += value.first * value.first * value.second;
479  }
480 
481  return sqrt((sum2 - sum * sum / m_n) / (m_n - 1));
482 }
483 
484 template <typename Feeder>
487 {
488  assert(remove >= 0.0 && remove < 49.0);
489  long remove_count = static_cast<long>(remove * m_n / 100.0);
490  range_type result(m_feeder.value(0), m_feeder.value(m_histogram.size() - 1));
491 
492  if (remove_count > 0) {
493  long low_end = -1;
494  long counted_pixels_low = 0;
495 
496  while (counted_pixels_low < remove_count && low_end < (long)m_histogram.size())
497  counted_pixels_low += m_histogram[++low_end];
498 
499  result.first = m_feeder.value(low_end);
500  long high_end = m_histogram.size();
501  long counted_pixels_high = 0;
502 
503  while (counted_pixels_high <= remove_count && high_end > 0)
504  counted_pixels_high += m_histogram[--high_end];
505 
506  cvdebug() << " int range = " << low_end << ", " << high_end << " removing " << remove_count << " pixels at each end\n";
507  result.second = m_feeder.value(high_end);
508  }
509 
510  return result;
511 }
512 
514 
515 #endif
516 
517 
specialization of the THistogramFeeder for unsigned byte input data
Definition: histogram.hh:101
unsigned char value_type
typedef for generic programming
Definition: histogram.hh:104
A class to normalize and quantizize input data to a given histogram range with its given number of bi...
Definition: histogram.hh:50
T value_type
typedef for generic programming
Definition: histogram.hh:53
size_t size() const
Definition: histogram.hh:247
T value(size_t k) const
Definition: histogram.hh:267
THistogramFeeder(T min, T max, size_t bins)
Definition: histogram.hh:237
size_t index(T x) const
Definition: histogram.hh:253
a simple histogram that uses an instance of THistogramFeeder as input converter
Definition: histogram.hh:135
Feeder::value_type median() const
Definition: histogram.hh:368
std::vector< size_t >::const_iterator const_iterator
STL iterator.
Definition: histogram.hh:139
void push_range(Iterator begin, Iterator end)
Definition: histogram.hh:333
range_type get_reduced_range(double remove) const
Definition: histogram.hh:486
const_iterator end() const
Definition: histogram.hh:355
THistogram(const Feeder &f)
Definition: histogram.hh:294
std::pair< typename Feeder::value_type, size_t > value_type
A type for the value-index pair.
Definition: histogram.hh:142
double deviation() const
Definition: histogram.hh:467
double average() const
Definition: histogram.hh:407
void push(typename Feeder::value_type x)
Definition: histogram.hh:325
const_iterator begin() const
Definition: histogram.hh:349
double skewness() const
Definition: histogram.hh:446
double excess_kurtosis() const
Definition: histogram.hh:423
Feeder::value_type MAD() const
Definition: histogram.hh:381
const value_type at(size_t idx) const
Definition: histogram.hh:398
size_t operator[](size_t idx) const
Definition: histogram.hh:361
size_t size() const
Definition: histogram.hh:319
std::pair< typename Feeder::value_type, typename Feeder::value_type > range_type
Definition: histogram.hh:144
#define NS_MIA_BEGIN
conveniance define to start the mia namespace
Definition: defines.hh:33
#define NS_MIA_END
conveniance define to end the mia namespace
Definition: defines.hh:36
THistogramFeeder< unsigned char > CUBHistogramFeeder
typedef for the unsigned byte histogram feeder specialization
Definition: histogram.hh:120
CDebugSink & cvdebug()
Definition: msgstream.hh:226