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