Visual Servoing Platform  version 3.1.0
vpCLAHE.cpp
1 /****************************************************************************
2  *
3  * This file is part of the ViSP software.
4  * Copyright (C) 2005 - 2017 by Inria. All rights reserved.
5  *
6  * This software 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 2 of the License, or
9  * (at your option) any later version.
10  * See the file LICENSE.txt at the root directory of this source
11  * distribution for additional information about the GNU GPL.
12  *
13  * For using ViSP with software that can not be combined with the GNU
14  * GPL, please contact Inria about acquiring a ViSP Professional
15  * Edition License.
16  *
17  * See http://visp.inria.fr for more information.
18  *
19  * This software was developed at:
20  * Inria Rennes - Bretagne Atlantique
21  * Campus Universitaire de Beaulieu
22  * 35042 Rennes Cedex
23  * France
24  *
25  * If you have questions regarding the use of this file, please contact
26  * Inria at visp@inria.fr
27  *
28  * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
29  * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
30  *
31  * Description:
32  * CLAHE (Contrast Limited Adaptive Histogram Equalization) algorithm.
33  *
34  * Authors:
35  * Souriya Trinh
36  *
37  *****************************************************************************/
84 #include <visp3/core/vpImageConvert.h>
85 #include <visp3/imgproc/vpImgproc.h>
86 
87 namespace
88 {
89 int fastRound(const float value) { return (int)(value + 0.5f); }
90 
91 void clipHistogram(const std::vector<int> &hist, std::vector<int> &clippedHist, const int limit)
92 {
93  clippedHist = hist;
94  int clippedEntries = 0, clippedEntriesBefore = 0;
95  int histlength = (int)hist.size();
96 
97  do {
98  clippedEntriesBefore = clippedEntries;
99  clippedEntries = 0;
100  for (int i = 0; i < histlength; i++) {
101  int d = clippedHist[i] - limit;
102  if (d > 0) {
103  clippedEntries += d;
104  clippedHist[i] = limit;
105  }
106  }
107 
108  int d = clippedEntries / (histlength);
109  int m = clippedEntries % (histlength);
110  for (int i = 0; i < histlength; i++) {
111  clippedHist[i] += d;
112  }
113 
114  if (m != 0) {
115  int s = (histlength - 1) / m;
116  for (int i = s / 2; i < histlength; i += s) {
117  ++(clippedHist[i]);
118  }
119  }
120  } while (clippedEntries != clippedEntriesBefore);
121 }
122 
123 void createHistogram(const int blockRadius, const int bins, const int blockXCenter, const int blockYCenter,
124  const vpImage<unsigned char> &I, std::vector<int> &hist)
125 {
126  std::fill(hist.begin(), hist.end(), 0);
127 
128  int xMin = std::max(0, blockXCenter - blockRadius);
129  int yMin = std::max(0, blockYCenter - blockRadius);
130  int xMax = std::min((int)I.getWidth(), blockXCenter + blockRadius + 1);
131  int yMax = std::min((int)I.getHeight(), blockYCenter + blockRadius + 1);
132 
133  for (int y = yMin; y < yMax; ++y) {
134  for (int x = xMin; x < xMax; ++x) {
135  ++hist[fastRound(I[y][x] / 255.0f * bins)];
136  }
137  }
138 }
139 
140 std::vector<float> createTransfer(const std::vector<int> &hist, const int limit, std::vector<int> &cdfs)
141 {
142  clipHistogram(hist, cdfs, limit);
143  int hMin = (int)hist.size() - 1;
144 
145  for (int i = 0; i < hMin; ++i) {
146  if (cdfs[i] != 0) {
147  hMin = i;
148  }
149  }
150  int cdf = 0;
151  for (int i = hMin; i < (int)hist.size(); ++i) {
152  cdf += cdfs[i];
153  cdfs[i] = cdf;
154  }
155 
156  int cdfMin = cdfs[hMin];
157  int cdfMax = cdfs[hist.size() - 1];
158 
159  std::vector<float> transfer(hist.size());
160  for (int i = 0; i < (int)transfer.size(); ++i) {
161  transfer[i] = (cdfs[i] - cdfMin) / (float)(cdfMax - cdfMin);
162  }
163 
164  return transfer;
165 }
166 
167 float transferValue(const int v, std::vector<int> &clippedHist)
168 {
169  int clippedHistLength = (int)clippedHist.size();
170  int hMin = clippedHistLength - 1;
171  for (int i = 0; i < hMin; i++) {
172  if (clippedHist[i] != 0) {
173  hMin = i;
174  }
175  }
176 
177  int cdf = 0;
178  for (int i = hMin; i <= v; i++) {
179  cdf += clippedHist[i];
180  }
181 
182  int cdfMax = cdf;
183  for (int i = v + 1; i < clippedHistLength; ++i) {
184  cdfMax += clippedHist[i];
185  }
186 
187  int cdfMin = clippedHist[hMin];
188  return (cdf - cdfMin) / (float)(cdfMax - cdfMin);
189 }
190 
191 float transferValue(const int v, const std::vector<int> &hist, std::vector<int> &clippedHist, const int limit)
192 {
193  clipHistogram(hist, clippedHist, limit);
194 
195  return transferValue(v, clippedHist);
196 }
197 }
198 
224 void vp::clahe(const vpImage<unsigned char> &I1, vpImage<unsigned char> &I2, const int blockRadius, const int bins,
225  const float slope, const bool fast)
226 {
227  if (blockRadius < 0) {
228  std::cerr << "Error: blockRadius < 0!" << std::endl;
229  return;
230  }
231 
232  if (bins < 0 || bins > 256) {
233  std::cerr << "Error: (bins < 0 || bins > 256)!" << std::endl;
234  return;
235  }
236 
237  if ((unsigned int)(2 * blockRadius + 1) > I1.getWidth() || (unsigned int)(2 * blockRadius + 1) > I1.getHeight()) {
238  std::cerr << "Error: (unsigned int) (2*blockRadius+1) > I1.getWidth() || "
239  "(unsigned int) (2*blockRadius+1) > I1.getHeight()!"
240  << std::endl;
241  return;
242  }
243 
244  I2.resize(I1.getHeight(), I1.getWidth());
245 
246  if (fast) {
247  int blockSize = 2 * blockRadius + 1;
248  int limit = (int)(slope * blockSize * blockSize / bins + 0.5);
249 
250  /* div */
251  int nc = I1.getWidth() / blockSize;
252  int nr = I1.getHeight() / blockSize;
253 
254  /* % */
255  int cm = I1.getWidth() - nc * blockSize;
256  std::vector<int> cs;
257 
258  switch (cm) {
259  case 0:
260  cs.resize(nc);
261  for (int i = 0; i < nc; ++i) {
262  cs[i] = i * blockSize + blockRadius + 1;
263  }
264  break;
265 
266  case 1:
267  cs.resize(nc + 1);
268  for (int i = 0; i < nc; ++i) {
269  cs[i] = i * blockSize + blockRadius + 1;
270  }
271  cs[nc] = I1.getWidth() - blockRadius - 1;
272  break;
273 
274  default:
275  cs.resize(nc + 2);
276  cs[0] = blockRadius + 1;
277  for (int i = 0; i < nc; ++i) {
278  cs[i + 1] = i * blockSize + blockRadius + 1 + cm / 2;
279  }
280  cs[nc + 1] = I1.getWidth() - blockRadius - 1;
281  }
282 
283  int rm = I1.getHeight() - nr * blockSize;
284  std::vector<int> rs;
285 
286  switch (rm) {
287  case 0:
288  rs.resize((size_t)nr);
289  for (int i = 0; i < nr; ++i) {
290  rs[i] = i * blockSize + blockRadius + 1;
291  }
292  break;
293 
294  case 1:
295  rs.resize((size_t)(nr + 1));
296  for (int i = 0; i < nr; ++i) {
297  rs[i] = i * blockSize + blockRadius + 1;
298  }
299  rs[nr] = I1.getHeight() - blockRadius - 1;
300  break;
301 
302  default:
303  rs.resize((size_t)(nr + 2));
304  rs[0] = blockRadius + 1;
305  for (int i = 0; i < nr; ++i) {
306  rs[i + 1] = i * blockSize + blockRadius + 1 + rm / 2;
307  }
308  rs[nr + 1] = I1.getHeight() - blockRadius - 1;
309  }
310 
311  std::vector<int> hist((size_t)(bins + 1));
312  std::vector<int> cdfs((size_t)(bins + 1));
313  std::vector<float> tl;
314  std::vector<float> tr;
315  std::vector<float> br;
316  std::vector<float> bl;
317 
318  for (int r = 0; r <= (int)rs.size(); ++r) {
319  int r0 = std::max(0, r - 1);
320  int r1 = std::min((int)rs.size() - 1, r);
321  int dr = rs[r1] - rs[r0];
322 
323  createHistogram(blockRadius, bins, cs[0], rs[r0], I1, hist);
324  tr = createTransfer(hist, limit, cdfs);
325  if (r0 == r1) {
326  br = tr;
327  } else {
328  createHistogram(blockRadius, bins, cs[0], rs[r1], I1, hist);
329  br = createTransfer(hist, limit, cdfs);
330  }
331 
332  int yMin = (r == 0 ? 0 : rs[r0]);
333  int yMax = (r < (int)rs.size() ? rs[r1] : I1.getHeight());
334 
335  for (int c = 0; c <= (int)cs.size(); ++c) {
336  int c0 = std::max(0, c - 1);
337  int c1 = std::min((int)cs.size() - 1, c);
338  int dc = cs[c1] - cs[c0];
339 
340  tl = tr;
341  bl = br;
342 
343  if (c0 != c1) {
344  createHistogram(blockRadius, bins, cs[c1], rs[r0], I1, hist);
345  tr = createTransfer(hist, limit, cdfs);
346  if (r0 == r1) {
347  br = tr;
348  } else {
349  createHistogram(blockRadius, bins, cs[c1], rs[r1], I1, hist);
350  br = createTransfer(hist, limit, cdfs);
351  }
352  }
353 
354  int xMin = (c == 0 ? 0 : cs[c0]);
355  int xMax = (c < (int)cs.size() ? cs[c1] : I1.getWidth());
356  for (int y = yMin; y < yMax; ++y) {
357  float wy = (float)(rs[r1] - y) / dr;
358 
359  for (int x = xMin; x < xMax; ++x) {
360  float wx = (float)(cs[c1] - x) / dc;
361  int v = fastRound(I1[y][x] / 255.0f * bins);
362  float t00 = tl[v];
363  float t01 = tr[v];
364  float t10 = bl[v];
365  float t11 = br[v];
366  float t0 = 0.0f, t1 = 0.0f;
367 
368  if (c0 == c1) {
369  t0 = t00;
370  t1 = t10;
371  } else {
372  t0 = wx * t00 + (1.0f - wx) * t01;
373  t1 = wx * t10 + (1.0f - wx) * t11;
374  }
375 
376  float t = (r0 == r1) ? t0 : wy * t0 + (1.0f - wy) * t1;
377  I2[y][x] = std::max(0, std::min(255, fastRound(t * 255.0f)));
378  }
379  }
380  }
381  }
382  } else {
383  std::vector<int> hist(bins + 1), prev_hist(bins + 1);
384  std::vector<int> clippedHist(bins + 1);
385 
386  bool first = true;
387  int xMin0 = 0;
388  int xMax0 = std::min((int)I1.getWidth(), blockRadius);
389 
390  for (int y = 0; y < (int)I1.getHeight(); y++) {
391  int yMin = std::max(0, y - (int)blockRadius);
392  int yMax = std::min((int)I1.getHeight(), y + blockRadius + 1);
393  int h = yMax - yMin;
394 
395 #if 0
396  std::fill(hist.begin(), hist.end(), 0);
397  // Compute histogram for the current block
398  for (int yi = yMin; yi < yMax; yi++) {
399  for (int xi = xMin0; xi < xMax0; xi++) {
400  ++hist[fastRound(I1[yi][xi] / 255.0f * bins)];
401  }
402  }
403 #else
404  if (first) {
405  first = false;
406  // Compute histogram for the block at (0,0)
407  for (int yi = yMin; yi < yMax; yi++) {
408  for (int xi = xMin0; xi < xMax0; xi++) {
409  ++hist[fastRound(I1[yi][xi] / 255.0f * bins)];
410  }
411  }
412  } else {
413  hist = prev_hist;
414 
415  if (yMin > 0) {
416  int yMin1 = yMin - 1;
417  // Sliding histogram, remove top
418  for (int xi = xMin0; xi < xMax0; xi++) {
419  --hist[fastRound(I1[yMin1][xi] / 255.0f * bins)];
420  }
421  }
422 
423  if (y + blockRadius < (int)I1.getHeight()) {
424  int yMax1 = yMax - 1;
425  // Sliding histogram, add bottom
426  for (int xi = xMin0; xi < xMax0; xi++) {
427  ++hist[fastRound(I1[yMax1][xi] / 255.0f * bins)];
428  }
429  }
430  }
431  prev_hist = hist;
432 #endif
433 
434  for (int x = 0; x < (int)I1.getWidth(); x++) {
435  int xMin = std::max(0, x - (int)blockRadius);
436  int xMax = x + blockRadius + 1;
437 
438  if (xMin > 0) {
439  int xMin1 = xMin - 1;
440  // Sliding histogram, remove left
441  for (int yi = yMin; yi < yMax; yi++) {
442  --hist[fastRound(I1[yi][xMin1] / 255.0f * bins)];
443  }
444  }
445 
446  if (xMax <= (int)I1.getWidth()) {
447  int xMax1 = xMax - 1;
448  // Sliding histogram, add right
449  for (int yi = yMin; yi < yMax; yi++) {
450  ++hist[fastRound(I1[yi][xMax1] / 255.0f * bins)];
451  }
452  }
453 
454  int v = fastRound(I1[y][x] / 255.0f * bins);
455  int w = std::min((int)I1.getWidth(), xMax) - xMin;
456  int n = h * w;
457  int limit = (int)(slope * n / bins + 0.5f);
458  I2[y][x] = fastRound(transferValue(v, hist, clippedHist, limit) * 255.0f);
459  }
460  }
461  }
462 }
463 
489 void vp::clahe(const vpImage<vpRGBa> &I1, vpImage<vpRGBa> &I2, const int blockRadius, const int bins, const float slope,
490  const bool fast)
491 {
492  // Split
497 
498  vpImageConvert::split(I1, &pR, &pG, &pB, &pa);
499 
500  // Apply CLAHE independently on RGB channels
501  vpImage<unsigned char> resR, resG, resB;
502  clahe(pR, resR, blockRadius, bins, slope, fast);
503  clahe(pG, resG, blockRadius, bins, slope, fast);
504  clahe(pB, resB, blockRadius, bins, slope, fast);
505 
506  I2.resize(I1.getHeight(), I1.getWidth());
507  unsigned int size = I2.getWidth() * I2.getHeight();
508  unsigned char *ptrStart = (unsigned char *)I2.bitmap;
509  unsigned char *ptrEnd = ptrStart + size * 4;
510  unsigned char *ptrCurrent = ptrStart;
511 
512  unsigned int cpt = 0;
513  while (ptrCurrent != ptrEnd) {
514  *ptrCurrent = resR.bitmap[cpt];
515  ++ptrCurrent;
516 
517  *ptrCurrent = resG.bitmap[cpt];
518  ++ptrCurrent;
519 
520  *ptrCurrent = resB.bitmap[cpt];
521  ++ptrCurrent;
522 
523  *ptrCurrent = pa.bitmap[cpt];
524  ++ptrCurrent;
525 
526  cpt++;
527  }
528 }
Type * bitmap
points toward the bitmap
Definition: vpImage.h:133
static void split(const vpImage< vpRGBa > &src, vpImage< unsigned char > *pR, vpImage< unsigned char > *pG, vpImage< unsigned char > *pB, vpImage< unsigned char > *pa=NULL)
void resize(const unsigned int h, const unsigned int w)
resize the image : Image initialization
Definition: vpImage.h:856
VISP_EXPORT void clahe(const vpImage< unsigned char > &I1, vpImage< unsigned char > &I2, const int blockRadius=150, const int bins=256, const float slope=3.0f, const bool fast=true)
Definition: vpCLAHE.cpp:224
unsigned int getHeight() const
Definition: vpImage.h:178
unsigned int getWidth() const
Definition: vpImage.h:229