Visual Servoing Platform version 3.5.0
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
87namespace
88{
89int fastRound(float value) { return (int)(value + 0.5f); }
90
91void 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
123void 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
140std::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
167float 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
191float 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
228void 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
493void 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}
static void split(const vpImage< vpRGBa > &src, vpImage< unsigned char > *pR, vpImage< unsigned char > *pG, vpImage< unsigned char > *pB, vpImage< unsigned char > *pa=NULL)
unsigned int getWidth() const
Definition: vpImage.h:246
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
unsigned int getHeight() const
Definition: vpImage.h:188
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