[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]

histogram.hxx
1 /************************************************************************/
2 /* */
3 /* Copyright 2011-2012 by Ullrich Koethe */
4 /* */
5 /* This file is part of the VIGRA computer vision library. */
6 /* The VIGRA Website is */
7 /* http://hci.iwr.uni-heidelberg.de/vigra/ */
8 /* Please direct questions, bug reports, and contributions to */
9 /* ullrich.koethe@iwr.uni-heidelberg.de or */
10 /* vigra@informatik.uni-hamburg.de */
11 /* */
12 /* Permission is hereby granted, free of charge, to any person */
13 /* obtaining a copy of this software and associated documentation */
14 /* files (the "Software"), to deal in the Software without */
15 /* restriction, including without limitation the rights to use, */
16 /* copy, modify, merge, publish, distribute, sublicense, and/or */
17 /* sell copies of the Software, and to permit persons to whom the */
18 /* Software is furnished to do so, subject to the following */
19 /* conditions: */
20 /* */
21 /* The above copyright notice and this permission notice shall be */
22 /* included in all copies or substantial portions of the */
23 /* Software. */
24 /* */
25 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */
26 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */
27 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */
28 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */
29 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */
30 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */
31 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */
32 /* OTHER DEALINGS IN THE SOFTWARE. */
33 /* */
34 /************************************************************************/
35 
36 #ifndef VIGRA_HISTOGRAM_HXX
37 #define VIGRA_HISTOGRAM_HXX
38 
39 #include "config.hxx"
40 #include "array_vector.hxx"
41 #include <algorithm>
42 
43 namespace vigra {
44 
45 /** \brief Set histogram options.
46 
47  HistogramOptions objects are used to pass histogram options to other objects. This \ref acc_hist_options "example" shows how it is is used to pass histogram options to an accumulator chain.
48 */
50 {
51  public:
52 
53  /** \brief Lower bound for linear range mapping from values to indices. */
54  double minimum;
55 
56  /** \brief Upper bound for linear range mapping from values to indices. */
57  double maximum;
58 
59  /** \brief Total number of bins in the histogram. */
60  int binCount;
61 
62  /** \brief If true, range mapping bounds are defined by minimum and maximum of the data. */
64 
65  /** Initialize members with default values:
66 
67  - minimum, maximum = 0.0
68  - binCount = 64
69  - local_auto_init = false
70  */
72  : minimum(0.0), maximum(0.0),
73  binCount(64),
74  local_auto_init(false)
75  {}
76 
77  /** Set minimum = mi and maximum = ma. Requirement: mi < ma.
78  */
79  HistogramOptions & setMinMax(double mi, double ma)
80  {
81  vigra_precondition(mi < ma,
82  "HistogramOptions::setMinMax(): min < max required.");
83  minimum = mi;
84  maximum = ma;
85  return *this;
86  }
87 
88  /** Set binCount = c. Requirement: c > 0.
89  */
91  {
92  vigra_precondition(c > 0,
93  "HistogramOptions::setBinCount(): binCount > 0 required.");
94  binCount = c;
95  return *this;
96  }
97 
98  /** Set local_auto_init = true. Requirement: setMinMax() must not have been called before. */
100  {
101  vigra_precondition(!validMinMax(),
102  "HistogramOptions::regionAutoInit(): you must not call setMinMax() when auto initialization is desired.");
103  local_auto_init = true;
104  return *this;
105  }
106 
107  /** Set local_auto_init = false. Requirement: setMinMax() must not have been called before. */
109  {
110  vigra_precondition(!validMinMax(),
111  "HistogramOptions::globalAutoInit(): you must not call setMinMax() when auto initialization is desired.");
112  local_auto_init = false;
113  return *this;
114  }
115 
116  /** Return minimum < maximum.
117  */
118  bool validMinMax() const
119  {
120  return minimum < maximum;
121  }
122 };
123 
124 template <class DataType, class BinType>
125 class HistogramView
126 {
127  BinType * bins_;
128  int size_, stride_;
129  DataType offset_;
130  double scale_, scaleInverse_;
131 
132  public:
133  HistogramView(DataType const & min, DataType const & max, int binCount,
134  BinType * bins = 0, int stride = 1)
135  : bins_(bins),
136  size_(binCount),
137  stride_(stride),
138  offset_(min),
139  scale_(double(binCount) / (max - min)),
140  scaleInverse_(1.0 / scale_)
141  {}
142 
143  HistogramView & setData(BinType * bins , int stride = 1)
144  {
145  bins_ = bins;
146  stride_ = stride;
147  return *this;
148  }
149 
150  HistogramView & reset()
151  {
152  if(hasData())
153  for(int k=0; k<size_; ++k)
154  *(bins_ +k*stride_) = BinType();
155  return *this;
156  }
157 
158  void getBinCenters(ArrayVector<DataType> * centers) const
159  {
160  double invScale = 1.0 / scale_;
161  for(int k=0; k < size_; ++k)
162  {
163  (*centers)[k] = mapItemInverse(k + 0.5) ;
164  }
165  }
166 
167  int size() const
168  {
169  return size_;
170  }
171 
172  bool hasData() const
173  {
174  return bins_ != 0;
175  }
176 
177  BinType const & operator[](int k) const
178  {
179  return *(bins_ + k*stride_);
180  }
181 
182  double mapItem(DataType const & d) const
183  {
184  return scale_ * (d - offset_);
185  }
186 
187  DataType mapItemInverse(double d) const
188  {
189  return DataType(d * scaleInverse_ + offset_);
190  }
191 
192  void add(DataType const & d, BinType weight = NumericTraits<BinType>::one())
193  {
194  get(int(mapItem(d))) += weight;
195  }
196 
197  protected:
198 
199  BinType & get(int index)
200  {
201  if(index < 0)
202  index = 0;
203  if(index >= size_)
204  index = size_ - 1;
205  return *(bins_ + index*stride_);
206  }
207 };
208 
209 template <class T>
210 class TrapezoidKernel
211 {
212  public:
213  typedef T value_type;
214 
215  T operator[](double f) const
216  {
217  if(f < -0.5)
218  return 0.5*(f + 1.5);
219  if(f > 0.5)
220  return 0.5*(1.5 - f);
221  return 0.5;
222  }
223 
224  double radius() const
225  {
226  return 1.5;
227  }
228 
229  T findMaximum(double l, double c, double r) const
230  {
231  double curv = -2.0*c + r + l;
232  if(curv == 0.0)
233  return T(-0.5);
234  double extr = 0.5*(l-r) / curv;
235  if(curv < 0.0)
236  {
237  return extr < -0.5
238  ? T(-0.5)
239  : extr > 0.5
240  ? T(0.5)
241  : T(extr);
242  }
243  else
244  {
245  return extr < 0.0
246  ? T(0.5)
247  : T(-0.5);
248  }
249  }
250 
251  bool findMode(double l, double c, double r, double * m) const
252  {
253  double curv = -2.0*c + r + l;
254  if(curv >= 0.0)
255  return false;
256  *m = 0.5*(l-r) / curv;
257  if(*m < -0.5 || *m > 0.5)
258  return false;
259  return true;
260  }
261 };
262 
263 template <class DataType, class KernelType>
264 class KernelHistogramView
265 : public HistogramView<DataType, typename KernelType::value_type>
266 {
267  KernelType kernel_;
268  int radius_;
269 
270  public:
271 
272  typedef typename KernelType::value_type BinType;
273  typedef HistogramView<DataType, BinType> BaseType;
274 
275  KernelHistogramView(DataType const & min, DataType const & max, int binCount,
276  BinType * bins = 0, int stride = 1)
277  : BaseType(min, max, binCount, bins, stride),
278  radius_(kernel_.radius()-0.5) // FIXME: this needs generalization
279  {}
280 
281  void add(DataType const & d, BinType weight = NumericTraits<BinType>::one())
282  {
283  double mapped = this->mapItem(d);
284  double f = mapped - std::floor(mapped) - kernel_.radius();
285  int center = int(mapped);
286 
287  for(int k=center+radius_; k>=center-radius_; --k, f += 1.0)
288  {
289  this->get(k) += weight*kernel_[f];
290  }
291  }
292 
293  DataType findMode() const
294  {
295  double mmax = 0, vmax = 0, m;
296 
297  for(int k=0; k<this->size(); ++k)
298  {
299  double l = k > 0
300  ? (*this)[k-1]
301  : 0.0;
302  double c = (*this)[k];
303  double r = k < this->size() - 1
304  ? (*this)[k+1]
305  : 0.0;
306  if(kernel_.findMode(l, c, r, &m))
307  {
308  double v = l*kernel_[m+1.0] + c*kernel_[m] + r*kernel_[m-1.0];
309  if(vmax < v)
310  {
311  mmax = m + k + 0.5;
312  vmax = v;
313  }
314  }
315  }
316  return this->mapItemInverse(mmax);
317  }
318 
319  template <class Array>
320  void findModes(Array * modes)
321  {
322  double m;
323  for(int k=0; k<this->size(); ++k)
324  {
325  double l = k > 0
326  ? (*this)[k-1]
327  : 0.0;
328  double c = (*this)[k];
329  double r = k < this->size() - 1
330  ? (*this)[k+1]
331  : 0.0;
332  if(kernel_.findMode(l, c, r, &m))
333  {
334  double v = l*kernel_[m+1.0] + c*kernel_[m] + r*kernel_[m-1.0];
335  modes->push_back(std::make_pair(this->mapItemInverse(m + k + 0.5), v));
336  }
337  }
338  }
339 };
340 
341 template <class DataType, class BinType>
342 class Histogram
343 : public HistogramView<DataType, BinType>
344 {
345  public:
346  typedef HistogramView<DataType, BinType> BaseType;
347  ArrayVector<BinType> data_;
348 
349  public:
350  Histogram(DataType const & min, DataType const & max, int binCount,
351  BinType * bins = 0, int stride = 1)
352  : BaseType(min, max, binCount),
353  data_(binCount)
354  {
355  this->setData(&data_[0]);
356  }
357 
358  Histogram const & reset()
359  {
360  this->setData(&data_[0]);
361  BaseType::reset();
362  return *this;
363  }
364  };
365 
366 } // namespace vigra
367 
368 #endif // VIGRA_HISTOGRAM_HXX

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
vigra 1.9.0 (Sat Oct 5 2013)