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-2013 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 public:
52  typedef T value_type;
53 
61  THistogramFeeder(T min, T max, size_t bins);
62 
64  size_t size() const;
65 
72  size_t index(T x) const;
73 
79  T value(size_t k) const;
80 private:
81  T m_min;
82  T m_max;
83  size_t m_bins;
84  double m_step;
85  double m_inv_step;
86 };
87 
98 template <>
99 class THistogramFeeder<unsigned char> {
100 public:
102  typedef unsigned char value_type;
103 
105  THistogramFeeder(unsigned char min, unsigned char max, size_t bins);
106 
108  size_t size() const;
109 
111  size_t index(unsigned char x) const;
112 
114  unsigned char value(size_t k) const;
115 };
116 
119 
131 template <typename Feeder>
132 class THistogram {
133 public:
134 
136  typedef std::vector<size_t>::const_iterator const_iterator;
137 
139  typedef std::pair<typename Feeder::value_type, size_t> value_type;
140 
141  typedef std::pair<typename Feeder::value_type, typename Feeder::value_type> range_type;
142 
146  THistogram(const Feeder& f);
147 
157  THistogram(const THistogram<Feeder>& org, double perc);
158 
163  void push(typename Feeder::value_type x);
164 
170  void push(typename Feeder::value_type x, size_t count);
171 
178  template <typename Iterator>
179  void push_range(Iterator begin, Iterator end);
180 
182  size_t size() const;
183 
185  const_iterator begin() const;
186 
188  const_iterator end() const;
189 
191  size_t operator [] (size_t idx) const;
192 
198  const value_type at(size_t idx) const;
199 
200 
202  typename Feeder::value_type median() const;
203 
205  typename Feeder::value_type MAD() const;
206 
208  double average() const;
209 
211  double deviation() const;
212 
219  range_type get_reduced_range(double remove) const;
220 private:
221  Feeder m_feeder;
222  std::vector<size_t> m_histogram;
223  size_t m_n;
224 };
225 
226 // inline inplementation
227 template <typename T>
228 THistogramFeeder<T>::THistogramFeeder(T min, T max, size_t bins):
229  m_min(min),
230  m_max(max),
231  m_bins(bins),
232  m_step(( double(max) - double(min) ) / double(bins - 1)),
233  m_inv_step(double(bins - 1) / (double(max) - double(min)))
234 {
235 }
236 
237 template <typename T>
239 {
240  return m_bins;
241 }
242 
243 template <typename T>
244 inline size_t THistogramFeeder<T>::index(T x) const
245 {
246  double val = floor(m_inv_step * (x - m_min) + 0.5);
247  if (val < 0)
248  return 0;
249  if (val < m_bins)
250  return val;
251  return m_bins - 1;
252 }
253 
254 template <typename T>
255 T THistogramFeeder<T>::value(size_t k) const
256 {
257  return k * m_step + m_min;
258 }
259 
260 inline THistogramFeeder<unsigned char>::THistogramFeeder(unsigned char /*min*/, unsigned char /*max*/, size_t /*bins*/)
261 {
262 }
263 
265 {
266  return 256;
267 }
268 
269 inline
270 size_t THistogramFeeder<unsigned char>::index(unsigned char x) const
271 {
272  return x;
273 }
274 
275 inline
276 unsigned char THistogramFeeder<unsigned char>::value(size_t k) const
277 {
278  return k;
279 }
280 
281 template <typename Feeder>
283  m_feeder(f),
284  m_histogram(f.size()),
285  m_n(0)
286 {
287 }
288 
289 template <typename Feeder>
291  m_feeder(org.m_feeder),
292  m_histogram(m_feeder.size()),
293  m_n(0)
294 {
295  size_t n = (size_t)(org.m_n * (1.0 - perc));
296 
297  size_t i = 0;
298  while (n > m_n && i < m_histogram.size()) {
299  m_n += org.m_histogram[i];
300  m_histogram[i] = org.m_histogram[i];
301  ++i;
302  }
303 }
304 
305 
306 template <typename Feeder>
308 {
309  return m_histogram.size();
310 }
311 
312 template <typename Feeder>
313 void THistogram<Feeder>::push(typename Feeder::value_type x)
314 {
315  ++m_n;
316  ++m_histogram[m_feeder.index(x)];
317 }
318 
319 template <typename Feeder>
320 template <typename Iterator>
321 void THistogram<Feeder>::push_range(Iterator begin, Iterator end)
322 {
323  while (begin != end)
324  push(*begin++);
325 }
326 
327 
328 
329 template <typename Feeder>
330 void THistogram<Feeder>::push(typename Feeder::value_type x, size_t count)
331 {
332  m_n += count;
333  m_histogram[m_feeder.index(x)] += count;
334 }
335 
336 template <typename Feeder>
338 {
339  return m_histogram.begin();
340 }
341 
342 template <typename Feeder>
344 {
345  return m_histogram.end();
346 }
347 
348 template <typename Feeder>
349 size_t THistogram<Feeder>::operator [] (size_t idx) const
350 {
351  assert(idx < m_histogram.size());
352  return m_histogram[idx];
353 }
354 
355 template <typename Feeder>
356 typename Feeder::value_type THistogram<Feeder>::median() const
357 {
358  float n_2 = m_n / 2.0f;
359  float sum = 0;
360  size_t k = 0;
361  while ( sum < n_2 )
362  sum += m_histogram[k++];
363 
364  return m_feeder.value(k-1);
365 }
366 
367 template <typename Feeder>
368 typename Feeder::value_type THistogram<Feeder>::MAD() const
369 {
370  typedef typename Feeder::value_type T;
371  T m = median();
372 
373  THistogram<Feeder> help(m_feeder);
374 
375  ;
376  for (size_t k = 0; k < size(); ++k) {
377  T v = m_feeder.value(k);
378  help.push(v > m ? v - m : m -v, m_histogram[k]);
379  }
380  return help.median();
381 }
382 
383 
384 template <typename Feeder>
386 {
387  if (idx < m_histogram.size())
388  return value_type(m_feeder.value(idx), m_histogram[idx]);
389  else
390  return value_type(m_feeder.value(idx), 0);
391 }
392 
393 template <typename Feeder>
395 {
396  if (m_n < 1)
397  return 0.0;
398  double sum = 0.0;
399  for (size_t i = 0; i < size(); ++i) {
400  const typename THistogram<Feeder>::value_type value = at(i);
401  sum += value.first * value.second;
402  }
403  return sum / m_n;
404 }
405 
406 template <typename Feeder>
408 {
409  if (m_n < 2)
410  return 0.0;
411  double sum = 0.0;
412  double sum2 = 0.0;
413  for (size_t i = 0; i < size(); ++i) {
414  const typename THistogram<Feeder>::value_type value = at(i);
415  sum += value.first * value.second;
416  sum2 += value.first * value.first * value.second;
417  }
418  return sqrt((sum2 - sum * sum / m_n) / (m_n - 1));
419 }
420 
421 template <typename Feeder>
424 {
425  assert(remove >= 0.0 && remove < 49.0);
426  long remove_count = static_cast<long>(remove * m_n / 100.0);
427 
428  range_type result(m_feeder.value(0), m_feeder.value(m_histogram.size() - 1));
429 
430  if (remove_count > 0) {
431  long low_end = -1;
432  long counted_pixels_low = 0;
433 
434  while (counted_pixels_low < remove_count && low_end < (long)m_histogram.size())
435  counted_pixels_low += m_histogram[++low_end];
436 
437  result.first = m_feeder.value(low_end);
438 
439  long high_end = m_histogram.size();
440  long counted_pixels_high = 0;
441  while (counted_pixels_high <= remove_count && high_end > 0)
442  counted_pixels_high += m_histogram[--high_end];
443  cvdebug() << " int range = " << low_end << ", " << high_end << " removing " << remove_count << " pixels at each end\n";
444  result.second = m_feeder.value(high_end);
445  }
446  return result;
447 }
448 
450 
451 #endif
452 
453