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
48template <typename T>
50{
51public:
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;
81private:
82 T m_min;
83 T m_max;
84 size_t m_bins;
85 double m_step;
86 double m_inv_step;
87};
88
99template <>
100class THistogramFeeder<unsigned char>
101{
102public:
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
133template <typename Feeder>
135{
136public:
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;
229private:
230 Feeder m_feeder;
231 std::vector<size_t> m_histogram;
232 size_t m_n;
233};
234
235// inline inplementation
236template <typename T>
237THistogramFeeder<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
246template <typename T>
248{
249 return m_bins;
250}
251
252template <typename T>
253inline 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
266template <typename T>
268{
269 return k * m_step + m_min;
270}
271
272inline THistogramFeeder<unsigned char>::THistogramFeeder(unsigned char /*min*/, unsigned char /*max*/, size_t /*bins*/)
273{
274}
275
276inline size_t THistogramFeeder<unsigned char>::size() const
277{
278 return 256;
279}
280
281inline
282size_t THistogramFeeder<unsigned char>::index(unsigned char x) const
283{
284 return x;
285}
286
287inline
288unsigned char THistogramFeeder<unsigned char>::value(size_t k) const
289{
290 return k;
291}
292
293template <typename Feeder>
295 m_feeder(f),
296 m_histogram(f.size()),
297 m_n(0)
298{
299}
300
301template <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
318template <typename Feeder>
320{
321 return m_histogram.size();
322}
323
324template <typename Feeder>
325void THistogram<Feeder>::push(typename Feeder::value_type x)
326{
327 ++m_n;
328 ++m_histogram[m_feeder.index(x)];
329}
330
331template <typename Feeder>
332template <typename Iterator>
333void THistogram<Feeder>::push_range(Iterator begin, Iterator end)
334{
335 while (begin != end)
336 push(*begin++);
337}
338
339
340
341template <typename Feeder>
342void 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
348template <typename Feeder>
350{
351 return m_histogram.begin();
352}
353
354template <typename Feeder>
356{
357 return m_histogram.end();
358}
359
360template <typename Feeder>
361size_t THistogram<Feeder>::operator [] (size_t idx) const
362{
363 assert(idx < m_histogram.size());
364 return m_histogram[idx];
365}
366
367template <typename Feeder>
368typename 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
380template <typename Feeder>
381typename 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
397template <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
406template <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
422template <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
445template <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
466template <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
484template <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
unsigned char value_type
typedef for generic programming
Definition histogram.hh:104
THistogramFeeder(unsigned char min, unsigned char max, size_t bins)
Construct the feeder, the parameters are ignored.
size_t index(unsigned char x) const
unsigned char value(size_t k) const
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
THistogramFeeder(T min, T max, size_t bins)
Definition histogram.hh:237
size_t index(T x) const
Definition histogram.hh:253
T value(size_t k) const
Definition histogram.hh:267
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