Loading...
Searching...
No Matches
Histogram.h
1/*
2
3MIT License
4
5Copyright (c) 2017 FMI Open Development / Markus Peura, first.last@fmi.fi
6
7Permission is hereby granted, free of charge, to any person obtaining a copy
8of this software and associated documentation files (the "Software"), to deal
9in the Software without restriction, including without limitation the rights
10to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
11copies of the Software, and to permit persons to whom the Software is
12furnished to do so, subject to the following conditions:
13
14The above copyright notice and this permission notice shall be included in all
15copies or substantial portions of the Software.
16
17THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
18IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
19FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
20AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
21LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
22OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
23SOFTWARE.
24
25*/
26/*
27Part of Rack development has been done in the BALTRAD projects part-financed
28by the European Union (European Regional Development Fund and European
29Neighbourhood Partnership Instrument, Baltic Sea Region Programme 2007-2013)
30*/
31#ifndef HISTOGRAM_H_
32#define HISTOGRAM_H_
33
34//
35#include <drain/Log.h>
36#include <drain/TypeUtils.h>
37#include <typeinfo>
38#include <cmath>
39
40#include <vector>
41#include <iostream>
42#include <stdexcept>
43#include <map>
44
45#include "ValueScaling.h"
46
47
48
49namespace drain
50{
51
52// // // using namespace std;
53
54//std::cerr;
55
57
61class Histogram : protected std::vector<unsigned long> {
62public:
63
64 typedef unsigned long count_t;
65 typedef std::vector<count_t> vect_t;
66
67 Histogram(size_t size=256);
68
69 Histogram(const Histogram & histogram);
70
71 virtual ~Histogram(){};
72
73 ValueScaling scaling;
74
76 void setSize(size_t s);
77
78 static
79 std::size_t recommendSizeByType(const std::type_info & type, std::size_t defaultValue = 256);
80
81 inline
82 int getSize() const { return size(); };
83
84 inline
85 int autoSize(const std::type_info & type){
86 if (empty()){
87 resize(recommendSizeByType(type, 256));
88 }
89 return size();
90 };
91
93 inline
94 void clearBins(){
95 std::fill(begin(), end(), 0);
96 }
97
99
105 template <class T>
106 void compute(const T & src, const std::type_info & type = typeid(double), const UniTuple<double,2> & scaling = {1.0, 0.0});
107
109 //void clearBins();
110
111
113 void setSampleCount(long int n){
114 //sampleCount = n;
115 sampleCount = n;
116 sampleCountMedian = static_cast<size_t>(weight * static_cast<double>(n));
117 }
118
119 size_t getSampleCount() const {
120 return sampleCount;
121 }
122
123
125 void setRange(double dataMin, double dataMax);
126
128 inline
129 void setRange(const Range<double> & range){
130 setRange(range.min, range.max);
131 }
132
133 inline
134 void deriveScaling(const ValueScaling & s, const Type & type){
135 setRange(s.fwd(0.0), s.fwd(1 << (8*drain::Type::call<drain::sizeGetter>(type))));
136 }
137
138 inline
139 void setScale(const ValueScaling & scaling){
140 drain::Logger mout(__FILE__, __FUNCTION__);
141 mout.discouraged("use setRange instead (scaling=", scaling, ")");
142 mout.advice("range perhaps: [", this->scaling.getPhysicalRange());
143 this->scaling.assignSequence(scaling);
144 //int s = 1 << (8*drain::Type::call<drain::sizeGetter>(type));
145 }
146
148 /*
149 inline
150 void setScale(const ValueScaling & scaling){
151 //scaling.setRange(0.0, bins-1.0, dataMin, dataMax);
152 };
153 */
154
155 inline
156 int getInMin() const { return 0; };
157
159 inline
160 int getUpperBoundIn() const { return size(); };
161
162
163 inline
164 int getOutMin() const { return scaling.fwd(0); };
165
167 inline
168 int getUpperBoundOut() const { return scaling.fwd(size()); };
169
170
171
173
176 inline
177 void setMedianPosition(double pos){
178 weight = pos;
179 sampleCountMedian = static_cast<size_t>(weight * static_cast<double>(sampleCount));
180 };
181
182 /*
183 template <class T>
184 inline
185 bool withinLimits(T i) const {
186 return ((i >= inMin) && (i <= inMax)); // check if ok inMax is too much
187 }
188 */
189
190 template <class T>
191 inline
192 void increment(T i){
193 /*
194 size_t i2 = scaling.inv(i);
195 if (i2 >= this->size()){
196 std::cerr << "fail: " << i << "\t => \t" << i2 << '\t';
197 return;
198 }*/
199 ++(*this)[scaling.inv(i)];
200 ++sampleCount; // TODO: slow?
201 }
202
203 template <class T>
204 inline
205 void increment(T i, int count){
206 /*
207 size_t i2 = scaling.inv(i);
208 if (i2 >= this->size()){
209 std::cerr << "fail: " << i << "\t => \t" << i2 << '\t';
210 return;
211 }
212 */
213 (*this)[scaling.inv(i)] += count;
214 //(*this)[((i-inMin)*bins)/inSpan] += count;
215 sampleCount += count;
216 }
217
218 template <class T>
219 inline
220 void incrementRaw(T i){
221 /*
222 if (i >= this->size()){
223 std::cerr << "failRaw: " << i << '\t';
224 return;
225 }*/
226 ++(*this)[i];
227 ++sampleCount; // TODO: slow?
228 }
229
230 template <class T>
231 inline
232 void incrementRaw(T i, int count){
233 if (i > this->size()){
234 std::cerr << "failRaw: " << i << '\t';
235 return;
236 }
237 (*this)[i] += count;
238 sampleCount += count;
239 }
240
241 template <class T>
242 inline
243 void decrement(T i){
244 --(*this)[scaling.inv(i)]; // [((i-inMin)*bins)/inSpan];
245 --sampleCount;
246 }
247
248 template <class T>
249 inline
250 void decrement(T i, int count){
251 (*this)[scaling.inv(i)] -= count; //[((i-inMin)*bins)/inSpan] -= count;
252 sampleCount -= count;
253 }
254
255 template <class T>
256 inline
257 void decrementRaw(T i){
258 --(*this)[i];
259 --sampleCount;
260 }
261
262 template <class T>
263 inline
264 void decrementRaw(T i, int count){
265 (*this)[i] -= count;
266 sampleCount -= count;
267 }
268
269
270 template <class T>
271 inline
272 T scaleOut(size_type i) const {
273 //return static_cast<T>(outMin + (i*outSpan)/bins);
274 return static_cast<T>(scaling.fwd(i));
275 }
276
277
278
279
281
283 // @param p applies to weighted median; for standard median, p = 0.5.
284 template <class T>
285 inline
286 T getMax() const {
287 for (size_type i = size()-1; i > 0; --i)
288 if ((*this)[i] > 0)
289 return scaleOut<T>(i);
290 // cerr << "Warning: Histogram empty.\n";
291 return getUpperBoundOut(); //static_cast<T>(outMax);
292 }
293
294
295 template <class T>
296 inline
297 T getMin() const {
298 for (size_type i = 0; i < size(); ++i)
299 if ((*this)[i] > 0)
300 return scaleOut<T>(i);
301 //cerr << "Warning: Histogram empty.\n";
302 return getUpperBoundOut(); //static_cast<T>(outMax);
303 }
304
306
309 template <class T>
310 //inline
311 T getSum() const {
312 double sum = 0;
313 //std::cout << __FUNCTION__ << ':' << sum << std::endl;
314 for (size_type i = 0; i < size(); i++){
315 sum += ((*this)[i] * scaleOut<T>(i)); // count * nominal (should be scaled to middle)
316 //if ((*this)[i])
317 // std::cout << i << '\t' << (*this)[i] << '\t' << scaleOut<T>(i) << '\t' << sum << '\n';
318 }
319 return sum;
320 }
321
322
324 template <class T>
325 inline
326 T getMean() const {
327 //return scaleOut(getSum()/sampleCount);
328 if (sampleCount > 0)
329 return getSum<T>()/sampleCount;
330 else
331 return 0;
332 }
333
334
335 template <class T>
336 inline
337 T getMedian() const {
338 size_type sum = 0;
339 //const
340 //size_type limit = static_cast<size_t>(weight*sampleCount);
341
342 for (size_type i = 0; i < size(); ++i){
343 sum += (*this)[i];
344 if (sum >= sampleCountMedian){ // limit){
345 //if (sum >= limit){
346 return scaleOut<T>(i);
347 }
348 }
349 //return (245*sum)/sampleCount;
350
351 //return ::rand(); //static_cast<T>(outMax);
352 return getUpperBoundOut(); //static_cast<T>(outMax);
353 }
354
355
356 template <class T>
357 inline
358 T getWeightedMedian(float p) const {
359 size_type _sum;
360 if ((p < 0.0) || (p>1.0)){
361 throw std::runtime_error("Histogram<T>::getMedian: median point <0 or >1.0 .");
362 }
363 const size_t limit = sampleCountMedian; // static_cast<size_t>(p * sampleCount);
364 _sum = 0;
365 for (size_type i = 0; i < size(); i++){
366 _sum += (*this)[i];
367 if (_sum >= limit){
368 return static_cast<T>(i);
369 }
370 }
371 return static_cast<T>(size());
372 }
373
374
376 // NOTE: could be pipelined.
377 template <class T>
378 inline
379 T getVariance() const {
380 int n;
381 long sum;
382 long sum2;
383 T f;
384 T sumT;
385
386 sum = 0;
387 sum2 = 0;
388 for (size_type i = 0; i < size(); i++){
389 f = scaleOut<T>(i);
390 n = (*this)[i];
391 sum += n * f;
392 sum2 += n * (f*f);
393 }
394 sumT = static_cast<T>(sum)/sampleCount;
395 return static_cast<T>(sum2)/sampleCount - sumT*sumT;
396 }
397
398
399 template <class T>
400 inline
401 T getStdDeviation() const {
402 return static_cast<T>(sqrt(getVariance<double>()));
403 }
404
405
406
407
408
409
410 inline
411 const vect_t getVector() const {
412 return *this;
413 };
414
415 std::ostream & toStream(std::ostream & ostr) const;
416
417 void dump(std::ostream & ostr = std::cout);
418
419 std::string delimiter;
420
421 inline
422 double getValue(){
423 return (this->*statisticPtr)();
424 }
425
427
431 inline
432 double getValue(char c){
433 return (this->*getStatisticPtr(c))();
434 }
435
436
437 inline
438 void setValueFunc(char c){
439 statisticPtr = getStatisticPtr(c);
440 }
441
442
443protected:
444
445 void initialize();
446
447 typedef double (Histogram::*stat_ptr_t)() const;
448
449 double (Histogram::*statisticPtr)() const;
450
451 stat_ptr_t getStatisticPtr(char c);
452
453
454 //double (Histogram::*)getPtr const;
455
456private:
458 // size_t bins;
459
461 //size_type sampleCount;
462
464 size_type sampleCount = 0;
465
466 // Half of sampleCount.
467 size_type sampleCountMedian = 0;
468
469 // weight for weighted median;
470 float weight = 0.5;
471
473 static inline
474 bool isSmallInt(const std::type_info & type){
475 return (type == typeid(unsigned char)) || (type == typeid(unsigned short int));
476 }
477
479
482 static inline
483 short getBitShift(unsigned int value){
484 short i = 0;
485 while ((value = (value>>1)) > 0){
486 ++i;
487 }
488 return i;
489 }
490
491
492};
493
494
495template <class T>
496void Histogram::compute(const T & dst, const std::type_info & type, const UniTuple<double,2> & scaling){
497
498 drain::Logger mout(__FILE__, __FUNCTION__);
499
500 const int size = autoSize(type);
501
502 if (size <= 0){
503 mout.error(size, " bins, something went wrong");
504 return;
505 }
506
507 if (size > 0x10000){
508 mout.error(size, " bins exceeds Histogram size limit: ", 0x10000);
509 return;
510 }
511
512
513 initialize();
514
515 ValueScaling inputScaling(scaling);
516
517 const bool SAME_SCALING = (inputScaling == this->scaling);
518
519 mout.attention("scalings: ", this->scaling, ", ", inputScaling, " same? ", SAME_SCALING);
520 if (inputScaling == this->scaling){
521
522 }
523
524 if (isSmallInt(type)){
525 const unsigned short histBits = getBitShift(size);
526 if (size == (1<<histBits)){ //
527 const signed short dataBits = (8*drain::Type::call<drain::sizeGetter>(type));
528 const short int bitShift = dataBits - histBits;
529
530 mout.note("Small int (", dataBits, "b) data, histogram[", size,"] ", histBits, " b, computing with bit shift ", bitShift);
531 // setRange(inputScaling.fwd(0), inputScaling.fwd(size -1));
532 mout.info("Physical range of the histogram: ", this->scaling.getPhysicalRange());
533 // mout.warn("Histogram size ", size, "=2^N, using bitShift=", bitShift, ", computing raw values.");
534 for (typename T::const_iterator it = dst.begin(); it != dst.end(); ++it){
535 this->incrementRaw(static_cast<unsigned short int>(*it) >> bitShift);
536 }
537 }
538 else {
539 mout.advice("Consider histogram size of 2^N, e.g.", (1<<histBits));
540 }
541 return;
542 }
543
544 //mout.warn("Data range: ", range, ", mapped (physical) Range: ", rangeOut);
545 //mout.warn("Size ", s, " != 2^N (slow mode)");
546 //mout.advice("Consider histogram size of 2^N, e.g.", (1<<histBits));
547 mout.error("Skipping...");
548
549
550}
551
552inline
553std::ostream & operator<<(std::ostream &ostr,const Histogram &h){
554 return h.toStream(ostr);
555}
556
557}
558
559//std::ostream & operator<<(std::ostream &ostr,const drain::Histogram &h);
560
561//template <class T>
562
563
564#endif /*HISTOGRAM_H_*/
565
566// Drain
Class for computing a histogram and some statistics: average, min, max, mean, std....
Definition Histogram.h:61
void clearBins()
Does not change the size of the histogram.
Definition Histogram.h:94
double getValue(char c)
Return requested statistic.
Definition Histogram.h:432
void setRange(double dataMin, double dataMax)
Set range of original (physical) values to be mapped on the limited number of bins....
Definition Histogram.cpp:76
void setMedianPosition(double pos)
Set location of the median, if not in the middle (50%).
Definition Histogram.h:177
int getUpperBoundIn() const
Returns the upperLimit (exclusive)
Definition Histogram.h:160
int getInMin() const
Set range of original (physical) values to be mapped on the limited number of bins....
Definition Histogram.h:156
T getSum() const
Sum of the samples.
Definition Histogram.h:311
T getMean() const
Unscaled mean.
Definition Histogram.h:326
void setSampleCount(long int n)
Does not change the size of the histogram.
Definition Histogram.h:113
int getUpperBoundOut() const
Returns the upperLimit (exclusive)
Definition Histogram.h:168
void compute(const T &src, const std::type_info &type=typeid(double), const UniTuple< double, 2 > &scaling={1.0, 0.0})
Collect distribution.
Definition Histogram.h:496
T getVariance() const
Computes variance of the values inside the window.
Definition Histogram.h:379
T getMax() const
Statistics.
Definition Histogram.h:286
void setRange(const Range< double > &range)
Set range of original (physical) values to be mapped on the limited number of bins....
Definition Histogram.h:129
void setSize(size_t s)
Sets the number of bins; the resolution of the histogram.
Definition Histogram.cpp:62
LogSourc e is the means for a function or any program segment to "connect" to a Log.
Definition Log.h:312
Logger & note(const TT &... args)
For top-level information.
Definition Log.h:489
Logger & attention(const TT &... args)
Possible error, but execution can continue. Special type of Logger::warn().
Definition Log.h:476
Logger & error(const TT &... args)
Echoes.
Definition Log.h:416
Definition Range.h:52
Utilities related to std::type_info.
Definition Type.h:51
Tuple of N elements of type T.
Definition UniTuple.h:65
Linear scaling and physical range for image intensities.
Definition ValueScaling.h:64
double fwd(double x) const
Forward scaling: given encoded value x, returns corresponding value (possibly physically meaningful).
Definition ValueScaling.h:295
Definition DataSelector.cpp:1277