SMIL  0.9.1
DImageHistogram.hpp
1 /*
2  * Copyright (c) 2011-2016, Matthieu FAESSEL and ARMINES
3  * All rights reserved.
4  *
5  * Redistribution and use in source and binary forms, with or without
6  * modification, are permitted provided that the following conditions are met:
7  *
8  * * Redistributions of source code must retain the above copyright
9  * notice, this list of conditions and the following disclaimer.
10  * * Redistributions in binary form must reproduce the above copyright
11  * notice, this list of conditions and the following disclaimer in the
12  * documentation and/or other materials provided with the distribution.
13  * * Neither the name of Matthieu FAESSEL, or ARMINES nor the
14  * names of its contributors may be used to endorse or promote products
15  * derived from this software without specific prior written permission.
16  *
17  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ``AS IS'' AND ANY
18  * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
19  * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
20  * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS AND CONTRIBUTORS BE LIABLE FOR ANY
21  * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
22  * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
23  * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
24  * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
25  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
26  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
27  */
28 
29 
30 #ifndef _D_IMAGE_HISTOGRAM_HPP
31 #define _D_IMAGE_HISTOGRAM_HPP
32 
33 #include <climits>
34 
35 #include "DLineHistogram.hpp"
36 #include "DImageArith.hpp"
37 
38 
39 namespace smil
40 {
41 
48 #ifndef SWIG
49  template <class T>
50  ENABLE_IF( !IS_FLOAT(T), RES_T )
51  histogram(const Image<T> &imIn, size_t *h)
52  {
53  for (size_t i=0;i<ImDtTypes<T>::cardinal();i++)
54  h[i] = 0;
55 
56  typename Image<T>::lineType pixels = imIn.getPixels();
57  for (size_t i=0;i<imIn.getPixelCount();i++)
58  h[size_t(pixels[i]-ImDtTypes<T>::min())]++;
59 
60  return RES_OK;
61  }
62 
63  template <class T>
64  ENABLE_IF( IS_FLOAT(T), RES_T )
65  histogram(const Image<T> &/*imIn*/, size_t * /*h*/)
66  {
67  return RES_ERR_NOT_IMPLEMENTED;
68  }
69 
70  template <class T>
71  RES_T histogram(const Image<T> &imIn, const Image<T> &imMask, size_t *h)
72  {
73  ASSERT(haveSameSize(&imIn, &imMask, NULL));
74 
75  for (size_t i=0;i<ImDtTypes<T>::cardinal();i++)
76  h[i] = 0;
77 
78  typename Image<T>::lineType pixels = imIn.getPixels();
79  typename Image<T>::lineType maskPix = imMask.getPixels();
80 
81  for (size_t i=0;i<imIn.getPixelCount();i++)
82  if (maskPix[i]!=0)
83  h[size_t(pixels[i]-ImDtTypes<T>::min())]++;
84 
85  return RES_OK;
86  }
87 #endif // SWIG
88 
92  template <class T>
93  std::map<T, UINT> histogram(const Image<T> &imIn, bool fullRange=false)
94  {
95  vector<T> rVals = rangeVal(imIn);
96  size_t card = rVals[1]-rVals[0]+1;
97 
98  size_t *buf = new size_t[card];
99  for (size_t i=0;i<card;i++)
100  buf[i] = 0;
101 
102  typename Image<T>::lineType pixels = imIn.getPixels();
103  for (size_t i=0;i<imIn.getPixelCount();i++)
104  buf[size_t(pixels[i]-rVals[0])]++;
105 
106  map<T, UINT> h;
107 
108  if (fullRange)
109  for (T i=ImDtTypes<T>::min();i<rVals[0];i++)
110  h.insert(pair<T,UINT>(i, 0));
111 
112  for (size_t i=0;i<card;i++)
113  h.insert(pair<T,UINT>(i+rVals[0], buf[i]));
114 
115  if (fullRange)
116  for (T i=rVals[1];i<=ImDtTypes<T>::max() && i!=ImDtTypes<T>::min();i++)
117  h.insert(pair<T,UINT>(i, 0));
118 
119  delete[] buf;
120 
121  return h;
122  }
123 
124 
130  template <class T>
131  std::map<T, UINT> histogram(const Image<T> &imIn, const Image<T> &imMask, bool fullRange=false)
132  {
133  map<T, UINT> h;
134 
135  ASSERT(haveSameSize(&imIn, &imMask, NULL), h);
136 
137  vector<T> rVals = rangeVal(imIn);
138  size_t card = rVals[1]-rVals[0]+1;
139 
140  size_t *buf = new size_t[card];
141  for (size_t i=0;i<card;i++)
142  buf[i] = 0;
143 
144  typename Image<T>::lineType pixels = imIn.getPixels();
145  typename Image<T>::lineType maskPixels = imMask.getPixels();
146  for (size_t i=0;i<imIn.getPixelCount();i++)
147  if (maskPixels[i]!=0)
148  buf[size_t(pixels[i]-rVals[0])]++;
149 
150  if (fullRange)
151  for (T i=ImDtTypes<T>::min();i<rVals[0];i++)
152  h.insert(pair<T,UINT>(i, 0));
153 
154  for (size_t i=0;i<card;i++)
155  h.insert(pair<T,UINT>(i+rVals[0], buf[i]));
156 
157  if (fullRange)
158  for (T i=rVals[1];i<=ImDtTypes<T>::max() && i!=ImDtTypes<T>::min();i++)
159  h.insert(pair<T,UINT>(i, 0));
160 
161  delete[] buf;
162 
163  return h;
164  }
165 
169  template <class T, class T_out>
170  RES_T threshold(const Image<T> &imIn, T minVal, T maxVal, T_out trueVal, T_out falseVal, Image<T_out> &imOut)
171  {
172  if (minVal>maxVal) // Loop threshold
173  return threshold(imIn, maxVal, minVal, falseVal, trueVal, imOut);
174 
175  ASSERT_ALLOCATED(&imIn, &imOut);
176  ASSERT_SAME_SIZE(&imIn, &imOut);
177 
179 
180  iFunc.lineFunction.minVal = minVal;
181  iFunc.lineFunction.maxVal = maxVal;
182  iFunc.lineFunction.trueVal = trueVal;
183  iFunc.lineFunction.falseVal = falseVal;
184 
185  return iFunc(imIn, imOut);
186  }
187 
188  template <class T, class T_out>
189  RES_T threshold(const Image<T> &imIn, T minVal, T maxVal, Image<T_out> &imOut)
190  {
191  return threshold<T>(imIn, minVal, maxVal, ImDtTypes<T_out>::max(), ImDtTypes<T_out>::min(), imOut);
192  }
193 
194  template <class T, class T_out>
195  RES_T threshold(const Image<T> &imIn, T minVal, Image<T_out> &imOut)
196  {
197  return threshold<T>(imIn, minVal, ImDtTypes<T>::max(), ImDtTypes<T_out>::max(), ImDtTypes<T_out>::min(), imOut);
198  }
199 
200  template <class T, class T_out>
201  RES_T threshold(const Image<T> &imIn, Image<T_out> &imOut)
202  {
203  T tVal = otsuThreshold(imIn, imOut);
204  if (tVal==ImDtTypes<T>::min())
205  return RES_ERR;
206  else
207  return RES_OK;
208  }
209 
213  template <class T1, class T2>
214  RES_T stretchHist(const Image<T1> &imIn, T1 inMinVal, T1 inMaxVal, Image<T2> &imOut, T2 outMinVal=numeric_limits<T2>::min(), T2 outMaxVal=numeric_limits<T2>::max())
215  {
216  ASSERT_ALLOCATED(&imIn);
217  ASSERT_SAME_SIZE(&imIn, &imOut);
218 
220  iFunc.lineFunction.coeff = double (outMaxVal-outMinVal) / double (inMaxVal-inMinVal);
221  iFunc.lineFunction.inOrig = inMinVal;
222  iFunc.lineFunction.outOrig = outMinVal;
223 
224  return iFunc(imIn, imOut);
225  }
226 
227  template <class T1, class T2>
228  RES_T stretchHist(const Image<T1> &imIn, Image<T2> &imOut, T2 outMinVal, T2 outMaxVal)
229  {
230  ASSERT_ALLOCATED(&imIn);
231  ASSERT_SAME_SIZE(&imIn, &imOut);
232 
234  vector<T1> rangeV = rangeVal(imIn);
235  iFunc.lineFunction.coeff = double (outMaxVal-outMinVal) / double (rangeV[1]-rangeV[0]);
236  iFunc.lineFunction.inOrig = rangeV[0];
237  iFunc.lineFunction.outOrig = outMinVal;
238 
239  return iFunc(imIn, imOut);
240  }
241 
242  template <class T1, class T2>
243  RES_T stretchHist(const Image<T1> &imIn, Image<T2> &imOut)
244  {
245  return stretchHist<T1,T2>(imIn, imOut, numeric_limits<T2>::min(), numeric_limits<T2>::max());
246  }
247 
248 
255  template <class T>
256  vector<T> histogramRange(const Image<T> &imIn, double ignorePercent, bool cumulative=true)
257  {
258  vector<T> rVect;
259 
260  ASSERT(imIn.isAllocated(), rVect);
261 
262  size_t *h = new size_t[ImDtTypes<T>::cardinal()];
263  histogram(imIn, h);
264 
265  double imVol;
266  if (cumulative)
267  imVol = imIn.getPixelCount();
268  else
269  imVol = *std::max_element(h, h+ImDtTypes<T>::cardinal());
270 
271  double satVol;
272  double curVol;
273  vector<T> rangeV = rangeVal(imIn);
274  T threshValLeft = rangeV[0];
275  T threshValRight = rangeV[1];
276 
277  satVol = imVol * ignorePercent / 100.;
278 
279  // left
280  curVol=0;
281  for (size_t i=rangeV[0]; i<rangeV[1]; i++)
282  {
283  if (cumulative)
284  curVol += double(h[i]);
285  else
286  curVol = double(h[i]);
287 
288  if (curVol>satVol)
289  break;
290  threshValLeft = i;
291  }
292 
293  // Right
294  curVol=0;
295  for (size_t i=rangeV[1]; i>size_t(rangeV[0]); i--)
296  {
297  if (cumulative)
298  curVol += double(h[i]);
299  else
300  curVol = double(h[i]);
301 
302  if (curVol>satVol)
303  break;
304  threshValRight = i;
305  }
306 
307  delete[] h;
308 
309  rVect.push_back(threshValLeft);
310  rVect.push_back(threshValRight);
311 
312  return rVect;
313  }
314 
318  template <class T>
319  RES_T enhanceContrast(const Image<T> &imIn, Image<T> &imOut, double sat=0.25)
320  {
321  vector<T> rangeV = histogramRange(imIn, sat, true);
322 
323  ASSERT(rangeV.size()==2);
324 
325  stretchHist(imIn, rangeV[0], rangeV[1], imOut);
326  imOut.modified();
327 
328  return RES_OK;
329  }
330 
331 
332  template <class T>
333  bool IncrementThresholds(vector<double> &thresholdIndexes, map<T, UINT> &hist, UINT threshLevels, double totalFrequency, double &globalMean, vector<double> &classMean, vector<double> &classFrequency)
334  {
335  unsigned long numberOfHistogramBins = hist.size();
336  unsigned long numberOfClasses = classMean.size();
337 
338  typedef double MeanType;
339  typedef double FrequencyType;
340 
341  MeanType meanOld;
342  FrequencyType freqOld;
343 
344  unsigned int k;
345  int j;
346 
347  // from the upper threshold down
348  for(j=static_cast<int>(threshLevels-1); j>=0; j--)
349  {
350  // if this threshold can be incremented (i.e. we're not at the end of the histogram)
351  if (thresholdIndexes[j] < numberOfHistogramBins - 2 - (threshLevels-1 - j) )
352  {
353  // increment it and update mean and frequency of the class bounded by the threshold
354  thresholdIndexes[j] += 1;
355 
356  meanOld = classMean[j];
357  freqOld = classFrequency[j];
358 
359  classFrequency[j] += hist[UINT(thresholdIndexes[j])];
360 
361  if (classFrequency[j]>0)
362  {
363  classMean[j] = (meanOld * static_cast<MeanType>(freqOld)
364  + static_cast<MeanType>(thresholdIndexes[j])
365  * static_cast<MeanType>(hist[UINT(thresholdIndexes[j])]))
366  / static_cast<MeanType>(classFrequency[j]);
367  }
368  else
369  {
370  classMean[j] = 0;
371  }
372 
373  // set higher thresholds adjacent to their previous ones, and update mean and frequency of the respective classes
374  for (k=j+1; k<threshLevels; k++)
375  {
376  thresholdIndexes[k] = thresholdIndexes[k-1] + 1;
377  classFrequency[k] = hist[UINT(thresholdIndexes[k])];
378  if (classFrequency[k]>0)
379  {
380  classMean[k] = static_cast<MeanType>(thresholdIndexes[k]);
381  }
382  else
383  {
384  classMean[k] = 0;
385  }
386  }
387 
388  // update mean and frequency of the highest class
389  classFrequency[numberOfClasses-1] = totalFrequency;
390  classMean[numberOfClasses-1] = globalMean * totalFrequency;
391 
392  for(k=0; k<numberOfClasses-1; k++)
393  {
394  classFrequency[numberOfClasses-1] -= classFrequency[k];
395  classMean[numberOfClasses-1] -= classMean[k] * static_cast<MeanType>(classFrequency[k]);
396  }
397 
398  if (classFrequency[numberOfClasses-1]>0)
399  {
400  classMean[numberOfClasses-1] /= static_cast<MeanType>(classFrequency[numberOfClasses-1]);
401  }
402  else
403  {
404  classMean[numberOfClasses-1] = 0;
405  }
406 
407  // exit the for loop if a threshold has been incremented
408  break;
409  }
410  else // if this threshold can't be incremented
411  {
412  // if it's the lowest threshold
413  if (j==0)
414  {
415  // we couldn't increment because we're done
416  return false;
417  }
418  }
419  }
420  // we incremented
421  return true;
422  }
423 
427  template <class T>
428  vector<T> otsuThresholdValues(map<T, UINT> &hist, UINT threshLevels=1)
429  {
430 
431  typedef double MeanType;
432  typedef vector<MeanType> MeanVectorType;
433 
434 
435  double totalFrequency = 0;
436  MeanType globalMean = 0;
437 
438  for (typename map<T, UINT>::iterator it=hist.begin();it!=hist.end();it++)
439  {
440  globalMean += double((*it).first) * double((*it).second);
441  totalFrequency += (*it).second;
442  }
443 
444  globalMean /= totalFrequency;
445 
446 
447  unsigned long numberOfClasses = threshLevels + 1;
448  MeanVectorType thresholdIndexes(threshLevels, 0);
449 
450  for(unsigned long j=0; j<threshLevels; j++)
451  {
452  thresholdIndexes[j] = j;
453  }
454 
455  MeanVectorType maxVarThresholdIndexes = thresholdIndexes;
456  double freqSum = 0;
457  MeanVectorType classFrequency(numberOfClasses, 0);
458 
459  for (unsigned long j=0; j<numberOfClasses-1; j++)
460  {
461  classFrequency[j] = hist[T(thresholdIndexes[j])];
462  freqSum += classFrequency[j];
463  }
464 
465  classFrequency[numberOfClasses-1] = totalFrequency - freqSum;
466 
467  double meanSum = 0;
468  MeanVectorType classMean(numberOfClasses, 0);
469 
470  for (unsigned long j=0; j < numberOfClasses-1; j++)
471  {
472  if (classFrequency[j]>0)
473  {
474  classMean[j] = j;
475  }
476  else
477  {
478  classMean[j] = 0;
479  }
480  meanSum += classMean[j] * classFrequency[j];
481  }
482 
483  if (classFrequency[numberOfClasses-1]>0)
484  {
485  classMean[numberOfClasses-1] = double(globalMean * totalFrequency - meanSum) / double(classFrequency[numberOfClasses-1]);
486  }
487  else
488  {
489  classMean[numberOfClasses-1] = 0;
490  }
491 
492  double maxVarBetween = 0;
493  for (unsigned long j=0; j<numberOfClasses; j++)
494  {
495  maxVarBetween += classFrequency[j] * (globalMean - classMean[j]) * (globalMean - classMean[j]);
496  }
497 
498  // explore all possible threshold configurations and choose the one that yields maximum between-class variance
499  while (IncrementThresholds(thresholdIndexes, hist, threshLevels, totalFrequency, globalMean, classMean, classFrequency))
500  {
501  double varBetween = 0;
502  for (unsigned long j=0; j<numberOfClasses; j++)
503  {
504  varBetween += classFrequency[j] * (globalMean - classMean[j]) * (globalMean - classMean[j]);
505  }
506 
507  if (varBetween > maxVarBetween)
508  {
509  maxVarBetween = varBetween;
510  maxVarThresholdIndexes = thresholdIndexes;
511  }
512  }
513 
514  vector<T> threshVals;
515 
516  for (unsigned long j=0; j<threshLevels; j++)
517  {
518  threshVals.push_back(T(maxVarThresholdIndexes[j])); //= histogram->GetBinMax(0,maxVarThresholdIndexes[j]);
519  }
520  threshVals.push_back(maxVarBetween);
521 
522  return threshVals;
523  }
524 
525  template <class T>
526  vector<T> otsuThresholdValues(const Image<T> &im, UINT threshLevels=1)
527  {
528  map<T, UINT> hist = histogram(im);
529  return otsuThresholdValues(hist, threshLevels);
530  }
531 
532  template <class T>
533  vector<T> otsuThresholdValues(const Image<T> &im, const Image<T> &imMask, UINT threshLevels=1)
534  {
535  map<T, UINT> hist = histogram(im, imMask);
536  return otsuThresholdValues(hist, threshLevels);
537  }
538 
539 
545  template <class T, class T_out>
546  vector<T> otsuThreshold(const Image<T> &imIn, Image<T_out> &imOut, UINT nbrThresholds)
547  {
548  if (!areAllocated(&imIn, &imOut, NULL))
549  return vector<T>();
550 
551  vector<T> tVals = otsuThresholdValues<T>(imIn, nbrThresholds);
552  map<T, T_out> lut;
553  T i = ImDtTypes<T>::min();
554  T_out lbl = 0;
555  for (typename vector<T>::iterator it=tVals.begin();it!=tVals.end();it++,lbl++)
556  {
557  while(i<(*it))
558  {
559  lut[i] = lbl;
560  i++;
561  }
562  }
563  while(i<=ImDtTypes<T>::max() && i>ImDtTypes<T>::min())
564  {
565  lut[i] = lbl;
566  i++;
567  }
568  applyLookup(imIn, lut, imOut);
569 
570  return tVals;
571 
572  }
573 
574  template <class T, class T_out>
575  T otsuThreshold(const Image<T> &imIn, Image<T_out> &imOut)
576  {
577  if (!areAllocated(&imIn, &imOut, NULL))
578  return ImDtTypes<T>::min();
579 
580  vector<T> tVals = otsuThresholdValues<T>(imIn, 1);
581  threshold(imIn, tVals[0], imOut);
582  return tVals[0];
583  }
584 
585 
586  template <class T, class T_out>
587  vector<T> otsuThreshold(const Image<T> &imIn, const Image<T> &imMask, Image<T_out> &imOut, UINT nbrThresholds=1)
588  {
589  if (!areAllocated(&imIn, &imOut, NULL))
590  return vector<T>();
591 
592  vector<T> tVals = otsuThresholdValues<T>(imIn, imMask, nbrThresholds);
593  threshold(imIn, tVals[0], imOut);
594 
595  return tVals;
596 
597  }
598 
601 } // namespace smil
602 
603 
604 #endif // _D_IMAGE_HISTOGRAM_HPP
605 
std::map< T, UINT > histogram(const Image< T > &imIn, const Image< T > &imMask, bool fullRange=false)
Image histogram with a mask image.
Definition: DImageHistogram.hpp:131
RES_T threshold(const Image< T > &imIn, T minVal, T maxVal, T_out trueVal, T_out falseVal, Image< T_out > &imOut)
Image threshold.
Definition: DImageHistogram.hpp:170
Definition: DColorConvert.h:38
Definition: DBaseImageOperations.hpp:72
lineType getPixels() const
Get the pixels as a 1D array.
Definition: DImage.hpp:110
bool isAllocated() const
Check if the image is allocated.
Definition: DBaseImage.h:163
RES_T enhanceContrast(const Image< T > &imIn, Image< T > &imOut, double sat=0.25)
Enhance contrast.
Definition: DImageHistogram.hpp:319
vector< T > rangeVal(const Image< T > &imIn, bool onlyNonZero=false)
Min and Max values of an image.
Definition: DMeasures.hpp:313
vector< T > otsuThresholdValues(map< T, UINT > &hist, UINT threshLevels=1)
Return the different threshold values and the value of the resulting variance between classes...
Definition: DImageHistogram.hpp:428
virtual void modified()
Trigger modified event (allows to force display update)
Definition: DImage.hxx:223
vector< T > histogramRange(const Image< T > &imIn, double ignorePercent, bool cumulative=true)
Min and Max values of an histogram ignoring left/right low values (lower than a given height/cumulati...
Definition: DImageHistogram.hpp:256
size_t getPixelCount() const
Get the number of pixels.
Definition: DBaseImage.h:150
RES_T stretchHist(const Image< T1 > &imIn, T1 inMinVal, T1 inMaxVal, Image< T2 > &imOut, T2 outMinVal=numeric_limits< T2 >::min(), T2 outMaxVal=numeric_limits< T2 >::max())
Stretch histogram.
Definition: DImageHistogram.hpp:214
Definition: DTypes.hpp:78
Main Image class.
Definition: DQVtkViewer.hpp:44
RES_T applyLookup(const Image< T1 > &imIn, const mapT &_map, Image< T2 > &imOut, T2 defaultValue=T2(0))
Apply a lookup map.
Definition: DImageArith.hpp:1092
vector< T > otsuThreshold(const Image< T > &imIn, Image< T_out > &imOut, UINT nbrThresholds)
Otsu Threshold.
Definition: DImageHistogram.hpp:546
T maxVal(const Image< T > &imIn, bool onlyNonZero=false)
Max value of an image.
Definition: DMeasures.hpp:260
T minVal(const Image< T > &imIn, bool onlyNonZero=false)
Min value of an image.
Definition: DMeasures.hpp:198