Visual Servoing Platform  version 3.5.0 under development (2022-02-15)
vpCLAHE.cpp
1 /****************************************************************************
2  *
3  * ViSP, open source Visual Servoing Platform software.
4  * Copyright (C) 2005 - 2019 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(float value) { return (int)(value + 0.5f); }
90 
91 void clipHistogram(const std::vector<int> &hist, std::vector<int> &clippedHist, 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(int blockRadius, int bins, int blockXCenter, 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, 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(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(int v, const std::vector<int> &hist, std::vector<int> &clippedHist, int limit)
192 {
193  clipHistogram(hist, clippedHist, limit);
194 
195  return transferValue(v, clippedHist);
196 }
197 }
198 
228 void vp::clahe(const vpImage<unsigned char> &I1, vpImage<unsigned char> &I2, int blockRadius, int bins,
229  float slope, bool fast)
230 {
231  if (blockRadius < 0) {
232  std::cerr << "Error: blockRadius < 0!" << std::endl;
233  return;
234  }
235 
236  if (bins < 0 || bins > 256) {
237  std::cerr << "Error: (bins < 0 || bins > 256)!" << std::endl;
238  return;
239  }
240 
241  if ((unsigned int)(2 * blockRadius + 1) > I1.getWidth() || (unsigned int)(2 * blockRadius + 1) > I1.getHeight()) {
242  std::cerr << "Error: (unsigned int) (2*blockRadius+1) > I1.getWidth() || "
243  "(unsigned int) (2*blockRadius+1) > I1.getHeight()!"
244  << std::endl;
245  return;
246  }
247 
248  I2.resize(I1.getHeight(), I1.getWidth());
249 
250  if (fast) {
251  int blockSize = 2 * blockRadius + 1;
252  int limit = (int)(slope * blockSize * blockSize / bins + 0.5);
253 
254  /* div */
255  int nc = I1.getWidth() / blockSize;
256  int nr = I1.getHeight() / blockSize;
257 
258  /* % */
259  int cm = I1.getWidth() - nc * blockSize;
260  std::vector<int> cs;
261 
262  switch (cm) {
263  case 0:
264  cs.resize(nc);
265  for (int i = 0; i < nc; ++i) {
266  cs[i] = i * blockSize + blockRadius + 1;
267  }
268  break;
269 
270  case 1:
271  cs.resize(nc + 1);
272  for (int i = 0; i < nc; ++i) {
273  cs[i] = i * blockSize + blockRadius + 1;
274  }
275  cs[nc] = I1.getWidth() - blockRadius - 1;
276  break;
277 
278  default:
279  cs.resize(nc + 2);
280  cs[0] = blockRadius + 1;
281  for (int i = 0; i < nc; ++i) {
282  cs[i + 1] = i * blockSize + blockRadius + 1 + cm / 2;
283  }
284  cs[nc + 1] = I1.getWidth() - blockRadius - 1;
285  }
286 
287  int rm = I1.getHeight() - nr * blockSize;
288  std::vector<int> rs;
289 
290  switch (rm) {
291  case 0:
292  rs.resize((size_t)nr);
293  for (int i = 0; i < nr; ++i) {
294  rs[i] = i * blockSize + blockRadius + 1;
295  }
296  break;
297 
298  case 1:
299  rs.resize((size_t)(nr + 1));
300  for (int i = 0; i < nr; ++i) {
301  rs[i] = i * blockSize + blockRadius + 1;
302  }
303  rs[nr] = I1.getHeight() - blockRadius - 1;
304  break;
305 
306  default:
307  rs.resize((size_t)(nr + 2));
308  rs[0] = blockRadius + 1;
309  for (int i = 0; i < nr; ++i) {
310  rs[i + 1] = i * blockSize + blockRadius + 1 + rm / 2;
311  }
312  rs[nr + 1] = I1.getHeight() - blockRadius - 1;
313  }
314 
315  std::vector<int> hist((size_t)(bins + 1));
316  std::vector<int> cdfs((size_t)(bins + 1));
317  std::vector<float> tl;
318  std::vector<float> tr;
319  std::vector<float> br;
320  std::vector<float> bl;
321 
322  for (int r = 0; r <= (int)rs.size(); ++r) {
323  int r0 = std::max(0, r - 1);
324  int r1 = std::min((int)rs.size() - 1, r);
325  int dr = rs[r1] - rs[r0];
326 
327  createHistogram(blockRadius, bins, cs[0], rs[r0], I1, hist);
328  tr = createTransfer(hist, limit, cdfs);
329  if (r0 == r1) {
330  br = tr;
331  } else {
332  createHistogram(blockRadius, bins, cs[0], rs[r1], I1, hist);
333  br = createTransfer(hist, limit, cdfs);
334  }
335 
336  int yMin = (r == 0 ? 0 : rs[r0]);
337  int yMax = (r < (int)rs.size() ? rs[r1] : I1.getHeight());
338 
339  for (int c = 0; c <= (int)cs.size(); ++c) {
340  int c0 = std::max(0, c - 1);
341  int c1 = std::min((int)cs.size() - 1, c);
342  int dc = cs[c1] - cs[c0];
343 
344  tl = tr;
345  bl = br;
346 
347  if (c0 != c1) {
348  createHistogram(blockRadius, bins, cs[c1], rs[r0], I1, hist);
349  tr = createTransfer(hist, limit, cdfs);
350  if (r0 == r1) {
351  br = tr;
352  } else {
353  createHistogram(blockRadius, bins, cs[c1], rs[r1], I1, hist);
354  br = createTransfer(hist, limit, cdfs);
355  }
356  }
357 
358  int xMin = (c == 0 ? 0 : cs[c0]);
359  int xMax = (c < (int)cs.size() ? cs[c1] : I1.getWidth());
360  for (int y = yMin; y < yMax; ++y) {
361  float wy = (float)(rs[r1] - y) / dr;
362 
363  for (int x = xMin; x < xMax; ++x) {
364  float wx = (float)(cs[c1] - x) / dc;
365  int v = fastRound(I1[y][x] / 255.0f * bins);
366  float t00 = tl[v];
367  float t01 = tr[v];
368  float t10 = bl[v];
369  float t11 = br[v];
370  float t0 = 0.0f, t1 = 0.0f;
371 
372  if (c0 == c1) {
373  t0 = t00;
374  t1 = t10;
375  } else {
376  t0 = wx * t00 + (1.0f - wx) * t01;
377  t1 = wx * t10 + (1.0f - wx) * t11;
378  }
379 
380  float t = (r0 == r1) ? t0 : wy * t0 + (1.0f - wy) * t1;
381  I2[y][x] = std::max(0, std::min(255, fastRound(t * 255.0f)));
382  }
383  }
384  }
385  }
386  } else {
387  std::vector<int> hist(bins + 1), prev_hist(bins + 1);
388  std::vector<int> clippedHist(bins + 1);
389 
390  bool first = true;
391  int xMin0 = 0;
392  int xMax0 = std::min((int)I1.getWidth(), blockRadius);
393 
394  for (int y = 0; y < (int)I1.getHeight(); y++) {
395  int yMin = std::max(0, y - (int)blockRadius);
396  int yMax = std::min((int)I1.getHeight(), y + blockRadius + 1);
397  int h = yMax - yMin;
398 
399 #if 0
400  std::fill(hist.begin(), hist.end(), 0);
401  // Compute histogram for the current block
402  for (int yi = yMin; yi < yMax; yi++) {
403  for (int xi = xMin0; xi < xMax0; xi++) {
404  ++hist[fastRound(I1[yi][xi] / 255.0f * bins)];
405  }
406  }
407 #else
408  if (first) {
409  first = false;
410  // Compute histogram for the block at (0,0)
411  for (int yi = yMin; yi < yMax; yi++) {
412  for (int xi = xMin0; xi < xMax0; xi++) {
413  ++hist[fastRound(I1[yi][xi] / 255.0f * bins)];
414  }
415  }
416  } else {
417  hist = prev_hist;
418 
419  if (yMin > 0) {
420  int yMin1 = yMin - 1;
421  // Sliding histogram, remove top
422  for (int xi = xMin0; xi < xMax0; xi++) {
423  --hist[fastRound(I1[yMin1][xi] / 255.0f * bins)];
424  }
425  }
426 
427  if (y + blockRadius < (int)I1.getHeight()) {
428  int yMax1 = yMax - 1;
429  // Sliding histogram, add bottom
430  for (int xi = xMin0; xi < xMax0; xi++) {
431  ++hist[fastRound(I1[yMax1][xi] / 255.0f * bins)];
432  }
433  }
434  }
435  prev_hist = hist;
436 #endif
437 
438  for (int x = 0; x < (int)I1.getWidth(); x++) {
439  int xMin = std::max(0, x - (int)blockRadius);
440  int xMax = x + blockRadius + 1;
441 
442  if (xMin > 0) {
443  int xMin1 = xMin - 1;
444  // Sliding histogram, remove left
445  for (int yi = yMin; yi < yMax; yi++) {
446  --hist[fastRound(I1[yi][xMin1] / 255.0f * bins)];
447  }
448  }
449 
450  if (xMax <= (int)I1.getWidth()) {
451  int xMax1 = xMax - 1;
452  // Sliding histogram, add right
453  for (int yi = yMin; yi < yMax; yi++) {
454  ++hist[fastRound(I1[yi][xMax1] / 255.0f * bins)];
455  }
456  }
457 
458  int v = fastRound(I1[y][x] / 255.0f * bins);
459  int w = std::min((int)I1.getWidth(), xMax) - xMin;
460  int n = h * w;
461  int limit = (int)(slope * n / bins + 0.5f);
462  I2[y][x] = fastRound(transferValue(v, hist, clippedHist, limit) * 255.0f);
463  }
464  }
465  }
466 }
467 
493 void vp::clahe(const vpImage<vpRGBa> &I1, vpImage<vpRGBa> &I2, int blockRadius, int bins, float slope,
494  bool fast)
495 {
496  // Split
501 
502  vpImageConvert::split(I1, &pR, &pG, &pB, &pa);
503 
504  // Apply CLAHE independently on RGB channels
505  vpImage<unsigned char> resR, resG, resB;
506  clahe(pR, resR, blockRadius, bins, slope, fast);
507  clahe(pG, resG, blockRadius, bins, slope, fast);
508  clahe(pB, resB, blockRadius, bins, slope, fast);
509 
510  I2.resize(I1.getHeight(), I1.getWidth());
511  unsigned int size = I2.getWidth() * I2.getHeight();
512  unsigned char *ptrStart = (unsigned char *)I2.bitmap;
513  unsigned char *ptrEnd = ptrStart + size * 4;
514  unsigned char *ptrCurrent = ptrStart;
515 
516  unsigned int cpt = 0;
517  while (ptrCurrent != ptrEnd) {
518  *ptrCurrent = resR.bitmap[cpt];
519  ++ptrCurrent;
520 
521  *ptrCurrent = resG.bitmap[cpt];
522  ++ptrCurrent;
523 
524  *ptrCurrent = resB.bitmap[cpt];
525  ++ptrCurrent;
526 
527  *ptrCurrent = pa.bitmap[cpt];
528  ++ptrCurrent;
529 
530  cpt++;
531  }
532 }
void resize(unsigned int h, unsigned int w)
resize the image : Image initialization
Definition: vpImage.h:800
Type * bitmap
points toward the bitmap
Definition: vpImage.h:143
static void split(const vpImage< vpRGBa > &src, vpImage< unsigned char > *pR, vpImage< unsigned char > *pG, vpImage< unsigned char > *pB, vpImage< unsigned char > *pa=NULL)
VISP_EXPORT void clahe(const vpImage< unsigned char > &I1, vpImage< unsigned char > &I2, int blockRadius=150, int bins=256, float slope=3.0f, bool fast=true)
Definition: vpCLAHE.cpp:228
unsigned int getHeight() const
Definition: vpImage.h:188
unsigned int getWidth() const
Definition: vpImage.h:246