Visual Servoing Platform version 3.6.0
Loading...
Searching...
No Matches
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
82namespace
83{
84int fastRound(float value) { return (int)(value + 0.5f); }
85
86void 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
118void 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(0, blockXCenter - blockRadius);
124 int yMin = std::max(0, blockYCenter - blockRadius);
125 int xMax = std::min((int)I.getWidth(), blockXCenter + blockRadius + 1);
126 int yMax = std::min((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
135std::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
162float 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
186float 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
194namespace vp
195{
196void 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(0, r - 1);
291 int r1 = std::min((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(0, c - 1);
309 int c1 = std::min((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(0, std::min(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)I1.getWidth(), blockRadius);
364
365 for (int y = 0; y < (int)I1.getHeight(); y++) {
366 int yMin = std::max(0, y - (int)blockRadius);
367 int yMax = std::min((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(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)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
440void 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=NULL)
Definition of the vpImage class member functions.
Definition vpImage.h:135
unsigned int getWidth() const
Definition vpImage.h:242
void resize(unsigned int h, unsigned int w)
resize the image : Image initialization
Definition vpImage.h:795
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