Mercurial > hgrepos > Python2 > PyMuPDF

comparison mupdf-source/thirdparty/leptonica/src/bilateral.c @ 2:b50eed0cc0ef upstream

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
ADD: MuPDF v1.26.7: the MuPDF source as downloaded by a default build of PyMuPDF 1.26.4. The directory name has changed: no version number in the expanded directory now.
author Franz Glasner <fzglas.hg@dom66.de>
date Mon, 15 Sep 2025 11:43:07 +0200
parents
children
comparison
equal deleted inserted replaced
1:1d09e1dec1d9 2:b50eed0cc0ef
1 /*====================================================================*
2 - Copyright (C) 2001 Leptonica. All rights reserved.
3 -
4 - Redistribution and use in source and binary forms, with or without
5 - modification, are permitted provided that the following conditions
6 - are met:
7 - 1. Redistributions of source code must retain the above copyright
8 - notice, this list of conditions and the following disclaimer.
9 - 2. Redistributions in binary form must reproduce the above
10 - copyright notice, this list of conditions and the following
11 - disclaimer in the documentation and/or other materials
12 - provided with the distribution.
13 -
14 - THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
15 - ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
16 - LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
17 - A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL ANY
18 - CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
19 - EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
20 - PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
21 - PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY
22 - OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
23 - NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
24 - SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
25 *====================================================================*/
26
27 /*!
28 * \file bilateral.c
29 * <pre>
30 *
31 * Top level approximate separable grayscale or color bilateral filtering
32 * PIX *pixBilateral()
33 * PIX *pixBilateralGray()
34 *
35 * Implementation of approximate separable bilateral filter
36 * static L_BILATERAL *bilateralCreate()
37 * static void *bilateralDestroy()
38 * static PIX *bilateralApply()
39 *
40 * Slow, exact implementation of grayscale or color bilateral filtering
41 * PIX *pixBilateralExact()
42 * PIX *pixBilateralGrayExact()
43 * PIX *pixBlockBilateralExact()
44 *
45 * Kernel helper function
46 * L_KERNEL *makeRangeKernel()
47 *
48 * This includes both a slow, exact implementation of the bilateral
49 * filter algorithm (given by Sylvain Paris and Frédo Durand),
50 * and a fast, approximate and separable implementation (following
51 * Yang, Tan and Ahuja). See bilateral.h for algorithmic details.
52 *
53 * The bilateral filter has the nice property of applying a gaussian
54 * filter to smooth parts of the image that don't vary too quickly,
55 * while at the same time preserving edges. The filter is nonlinear
56 * and cannot be decomposed into two separable filters; however,
57 * there exists an approximate method that is separable. To further
58 * speed up the separable implementation, you can generate the
59 * intermediate data at reduced resolution.
60 *
61 * The full kernel is composed of two parts: a spatial gaussian filter
62 * and a nonlinear "range" filter that depends on the intensity difference
63 * between the reference pixel at the spatial kernel origin and any other
64 * pixel within the kernel support.
65 *
66 * In our implementations, the range filter is a parameterized,
67 * one-sided, 256-element, monotonically decreasing gaussian function
68 * of the absolute value of the difference between pixel values; namely,
69 * abs(I2 - I1). In general, any decreasing function can be used,
70 * and more generally, any two-dimensional kernel can be used if
71 * you wish to relax the 'abs' condition. (In that case, the range
72 * filter can be 256 x 256).
73 * </pre>
74 */
75
76 #ifdef HAVE_CONFIG_H
77 #include <config_auto.h>
78 #endif /* HAVE_CONFIG_H */
79
80 #include <math.h>
81 #include "allheaders.h"
82 #include "bilateral.h"
83
84 static L_BILATERAL *bilateralCreate(PIX *pixs, l_float32 spatial_stdev,
85 l_float32 range_stdev, l_int32 ncomps,
86 l_int32 reduction);
87 static PIX *bilateralApply(L_BILATERAL *bil);
88 static void bilateralDestroy(L_BILATERAL **pbil);
89
90
91 #ifndef NO_CONSOLE_IO
92 #define DEBUG_BILATERAL 0
93 #endif /* ~NO_CONSOLE_IO */
94
95 /*--------------------------------------------------------------------------*
96 * Top level approximate separable grayscale or color bilateral filtering *
97 *--------------------------------------------------------------------------*/
98 /*!
99 * \brief pixBilateral()
100 *
101 * \param[in] pixs 8 bpp gray or 32 bpp rgb, no colormap
102 * \param[in] spatial_stdev of gaussian kernel; in pixels, > 0.5
103 * \param[in] range_stdev of gaussian range kernel; > 5.0; typ. 50.0
104 * \param[in] ncomps number of intermediate sums J(k,x);
105 * in [4 ... 30]
106 * \param[in] reduction 1, 2 or 4
107 * \return pixd bilateral filtered image, or NULL on error
108 *
109 * <pre>
110 * Notes:
111 * (1) This performs a relatively fast, separable bilateral
112 * filtering operation. The time is proportional to ncomps
113 * and varies inversely approximately as the cube of the
114 * reduction factor. See bilateral.h for algorithm details.
115 * (2) We impose minimum values for range_stdev and ncomps to
116 * avoid nasty artifacts when either are too small. We also
117 * impose a constraint on their product:
118 * ncomps * range_stdev >= 100.
119 * So for values of range_stdev >= 25, ncomps can be as small as 4.
120 * Here is a qualitative, intuitive explanation for this constraint.
121 * Call the difference in k values between the J(k) == 'delta', where
122 * 'delta' ~ 200 / ncomps
123 * Then this constraint is roughly equivalent to the condition:
124 * 'delta' < 2 * range_stdev
125 * Note that at an intensity difference of (2 * range_stdev), the
126 * range part of the kernel reduces the effect by the factor 0.14.
127 * This constraint requires that we have a sufficient number of
128 * PCBs (i.e, a small enough 'delta'), so that for any value of
129 * image intensity I, there exists a k (and a PCB, J(k), such that
130 * |I - k| < range_stdev
131 * Any fewer PCBs and we don't have enough to support this condition.
132 * (3) The upper limit of 30 on ncomps is imposed because the
133 * gain in accuracy is not worth the extra computation.
134 * (4) The size of the gaussian kernel is twice the spatial_stdev
135 * on each side of the origin. The minimum value of
136 * spatial_stdev, 0.5, is required to have a finite sized
137 * spatial kernel. In practice, a much larger value is used.
138 * (5) Computation of the intermediate images goes inversely
139 * as the cube of the reduction factor. If you can use a
140 * reduction of 2 or 4, it is well-advised.
141 * (6) The range kernel is defined over the absolute value of pixel
142 * grayscale differences, and hence must have size 256 x 1.
143 * Values in the array represent the multiplying weight
144 * depending on the absolute gray value difference between
145 * the source pixel and the neighboring pixel, and should
146 * be monotonically decreasing.
147 * (7) Interesting observation. Run this on prog/fish24.jpg, with
148 * range_stdev = 60, ncomps = 6, and spatial_dev = {10, 30, 50}.
149 * As spatial_dev gets larger, we get the counter-intuitive
150 * result that the body of the red fish becomes less blurry.
151 * (8) The image must be sufficiently big to get reasonable results.
152 * This requires the dimensions to be at least twice the filter size.
153 * Otherwise, return a copy of the input with warning.
154 * </pre>
155 */
156 PIX *
157 pixBilateral(PIX *pixs,
158 l_float32 spatial_stdev,
159 l_float32 range_stdev,
160 l_int32 ncomps,
161 l_int32 reduction)
162 {
163 l_int32 w, h, d, filtersize;
164 l_float32 sstdev; /* scaled spatial stdev */
165 PIX *pixt, *pixr, *pixg, *pixb, *pixd;
166
167 if (!pixs || pixGetColormap(pixs))
168 return (PIX *)ERROR_PTR("pixs not defined or cmapped", __func__, NULL);
169 pixGetDimensions(pixs, &w, &h, &d);
170 if (d != 8 && d != 32)
171 return (PIX *)ERROR_PTR("pixs not 8 or 32 bpp", __func__, NULL);
172 if (reduction != 1 && reduction != 2 && reduction != 4)
173 return (PIX *)ERROR_PTR("reduction invalid", __func__, NULL);
174 filtersize = (l_int32)(2.0 * spatial_stdev + 1.0 + 0.5);
175 if (w < 2 * filtersize || h < 2 * filtersize) {
176 L_WARNING("w = %d, h = %d; w or h < 2 * filtersize = %d; "
177 "returning copy\n", __func__, w, h, 2 * filtersize);
178 return pixCopy(NULL, pixs);
179 }
180 sstdev = spatial_stdev / (l_float32)reduction; /* reduced spat. stdev */
181 if (sstdev < 0.5)
182 return (PIX *)ERROR_PTR("sstdev < 0.5", __func__, NULL);
183 if (range_stdev <= 5.0)
184 return (PIX *)ERROR_PTR("range_stdev <= 5.0", __func__, NULL);
185 if (ncomps < 4 || ncomps > 30)
186 return (PIX *)ERROR_PTR("ncomps not in [4 ... 30]", __func__, NULL);
187 if (ncomps * range_stdev < 100.0)
188 return (PIX *)ERROR_PTR("ncomps * range_stdev < 100.0", __func__, NULL);
189
190 if (d == 8)
191 return pixBilateralGray(pixs, spatial_stdev, range_stdev,
192 ncomps, reduction);
193
194 pixt = pixGetRGBComponent(pixs, COLOR_RED);
195 pixr = pixBilateralGray(pixt, spatial_stdev, range_stdev, ncomps,
196 reduction);
197 pixDestroy(&pixt);
198 pixt = pixGetRGBComponent(pixs, COLOR_GREEN);
199 pixg = pixBilateralGray(pixt, spatial_stdev, range_stdev, ncomps,
200 reduction);
201 pixDestroy(&pixt);
202 pixt = pixGetRGBComponent(pixs, COLOR_BLUE);
203 pixb = pixBilateralGray(pixt, spatial_stdev, range_stdev, ncomps,
204 reduction);
205 pixDestroy(&pixt);
206 pixd = pixCreateRGBImage(pixr, pixg, pixb);
207 pixDestroy(&pixr);
208 pixDestroy(&pixg);
209 pixDestroy(&pixb);
210 return pixd;
211 }
212
213
214 /*!
215 * \brief pixBilateralGray()
216 *
217 * \param[in] pixs 8 bpp gray
218 * \param[in] spatial_stdev of gaussian kernel; in pixels, > 0.5
219 * \param[in] range_stdev of gaussian range kernel; > 5.0; typ. 50.0
220 * \param[in] ncomps number of intermediate sums J(k,x);
221 * in [4 ... 30]
222 * \param[in] reduction 1, 2 or 4
223 * \return pixd 8 bpp bilateral filtered image, or NULL on error
224 *
225 * <pre>
226 * Notes:
227 * (1) See pixBilateral() for constraints on the input parameters.
228 * (2) See pixBilateral() for algorithm details.
229 * </pre>
230 */
231 PIX *
232 pixBilateralGray(PIX *pixs,
233 l_float32 spatial_stdev,
234 l_float32 range_stdev,
235 l_int32 ncomps,
236 l_int32 reduction)
237 {
238 l_float32 sstdev; /* scaled spatial stdev */
239 PIX *pixd;
240 L_BILATERAL *bil;
241
242 if (!pixs || pixGetColormap(pixs))
243 return (PIX *)ERROR_PTR("pixs not defined or cmapped", __func__, NULL);
244 if (pixGetDepth(pixs) != 8)
245 return (PIX *)ERROR_PTR("pixs not 8 bpp gray", __func__, NULL);
246 if (reduction != 1 && reduction != 2 && reduction != 4)
247 return (PIX *)ERROR_PTR("reduction invalid", __func__, NULL);
248 sstdev = spatial_stdev / (l_float32)reduction; /* reduced spat. stdev */
249 if (sstdev < 0.5)
250 return (PIX *)ERROR_PTR("sstdev < 0.5", __func__, NULL);
251 if (range_stdev <= 5.0)
252 return (PIX *)ERROR_PTR("range_stdev <= 5.0", __func__, NULL);
253 if (ncomps < 4 || ncomps > 30)
254 return (PIX *)ERROR_PTR("ncomps not in [4 ... 30]", __func__, NULL);
255 if (ncomps * range_stdev < 100.0)
256 return (PIX *)ERROR_PTR("ncomps * range_stdev < 100.0", __func__, NULL);
257
258 bil = bilateralCreate(pixs, spatial_stdev, range_stdev, ncomps, reduction);
259 if (!bil) return (PIX *)ERROR_PTR("bil not made", __func__, NULL);
260 pixd = bilateralApply(bil);
261 bilateralDestroy(&bil);
262 return pixd;
263 }
264
265
266 /*----------------------------------------------------------------------*
267 * Implementation of approximate separable bilateral filter *
268 *----------------------------------------------------------------------*/
269 /*!
270 * \brief bilateralCreate()
271 *
272 * \param[in] pixs 8 bpp gray, no colormap
273 * \param[in] spatial_stdev of gaussian kernel; in pixels, > 0.5
274 * \param[in] range_stdev of gaussian range kernel; > 5.0; typ. 50.0
275 * \param[in] ncomps number of intermediate sums J(k,x);
276 * in [4 ... 30]
277 * \param[in] reduction 1, 2 or 4
278 * \return bil, or NULL on error
279 *
280 * <pre>
281 * Notes:
282 * (1) This initializes a bilateral filtering operation, generating all
283 * the data required. It takes most of the time in the bilateral
284 * filtering operation.
285 * (2) See bilateral.h for details of the algorithm.
286 * (3) See pixBilateral() for constraints on input parameters, which
287 * are not checked here.
288 * </pre>
289 */
290 static L_BILATERAL *
291 bilateralCreate(PIX *pixs,
292 l_float32 spatial_stdev,
293 l_float32 range_stdev,
294 l_int32 ncomps,
295 l_int32 reduction)
296 {
297 l_int32 w, ws, wd, h, hs, hd, i, j, k, index;
298 l_int32 border, minval, maxval, spatial_size;
299 l_int32 halfwidth, wpls, wplt, wpld, kval, nval, dval;
300 l_float32 sstdev, fval1, fval2, denom, sum, norm, kern;
301 l_int32 *nc, *kindex;
302 l_float32 *kfract, *range, *spatial;
303 l_uint32 *datas, *datat, *datad, *lines, *linet, *lined;
304 L_BILATERAL *bil;
305 PIX *pix1, *pix2, *pixt, *pixsc, *pixd;
306 PIXA *pixac;
307
308 if (reduction == 1) {
309 pix1 = pixClone(pixs);
310 } else if (reduction == 2) {
311 pix1 = pixScaleAreaMap2(pixs);
312 } else { /* reduction == 4) */
313 pix2 = pixScaleAreaMap2(pixs);
314 pix1 = pixScaleAreaMap2(pix2);
315 pixDestroy(&pix2);
316 }
317 if (!pix1)
318 return (L_BILATERAL *)ERROR_PTR("pix1 not made", __func__, NULL);
319
320 sstdev = spatial_stdev / (l_float32)reduction; /* reduced spat. stdev */
321 border = (l_int32)(2 * sstdev + 1);
322 pixsc = pixAddMirroredBorder(pix1, border, border, border, border);
323 pixGetExtremeValue(pix1, 1, L_SELECT_MIN, NULL, NULL, NULL, &minval);
324 pixGetExtremeValue(pix1, 1, L_SELECT_MAX, NULL, NULL, NULL, &maxval);
325 pixDestroy(&pix1);
326 if (!pixsc)
327 return (L_BILATERAL *)ERROR_PTR("pixsc not made", __func__, NULL);
328
329 bil = (L_BILATERAL *)LEPT_CALLOC(1, sizeof(L_BILATERAL));
330 bil->spatial_stdev = sstdev;
331 bil->range_stdev = range_stdev;
332 bil->reduction = reduction;
333 bil->ncomps = ncomps;
334 bil->minval = minval;
335 bil->maxval = maxval;
336 bil->pixsc = pixsc;
337 bil->pixs = pixClone(pixs);
338
339 /* -------------------------------------------------------------------- *
340 * Generate arrays for interpolation of J(k,x):
341 * (1.0 - kfract[.]) * J(kindex[.], x) + kfract[.] * J(kindex[.] + 1, x),
342 * where I(x) is the index into kfract[] and kindex[],
343 * and x is an index into the 2D image array.
344 * -------------------------------------------------------------------- */
345 /* nc is the set of k values to be used in J(k,x) */
346 nc = (l_int32 *)LEPT_CALLOC(ncomps, sizeof(l_int32));
347 for (i = 0; i < ncomps; i++)
348 nc[i] = minval + i * (maxval - minval) / (ncomps - 1);
349 bil->nc = nc;
350
351 /* kindex maps from intensity I(x) to the lower k index for J(k,x) */
352 kindex = (l_int32 *)LEPT_CALLOC(256, sizeof(l_int32));
353 for (i = minval, k = 0; i <= maxval && k < ncomps - 1; k++) {
354 fval2 = nc[k + 1];
355 while (i < fval2) {
356 kindex[i] = k;
357 i++;
358 }
359 }
360 kindex[maxval] = ncomps - 2;
361 bil->kindex = kindex;
362
363 /* kfract maps from intensity I(x) to the fraction of J(k+1,x) used */
364 kfract = (l_float32 *)LEPT_CALLOC(256, sizeof(l_float32)); /* from lower */
365 for (i = minval, k = 0; i <= maxval && k < ncomps - 1; k++) {
366 fval1 = nc[k];
367 fval2 = nc[k + 1];
368 while (i < fval2) {
369 kfract[i] = (l_float32)(i - fval1) / (l_float32)(fval2 - fval1);
370 i++;
371 }
372 }
373 kfract[maxval] = 1.0;
374 bil->kfract = kfract;
375
376 #if DEBUG_BILATERAL
377 for (i = minval; i <= maxval; i++)
378 lept_stderr("kindex[%d] = %d; kfract[%d] = %5.3f\n",
379 i, kindex[i], i, kfract[i]);
380 for (i = 0; i < ncomps; i++)
381 lept_stderr("nc[%d] = %d\n", i, nc[i]);
382 #endif /* DEBUG_BILATERAL */
383
384 /* -------------------------------------------------------------------- *
385 * Generate 1-D kernel arrays (spatial and range) *
386 * -------------------------------------------------------------------- */
387 spatial_size = 2 * sstdev + 1; /* same as the added border */
388 spatial = (l_float32 *)LEPT_CALLOC(spatial_size, sizeof(l_float32));
389 denom = 2. * sstdev * sstdev;
390 for (i = 0; i < spatial_size; i++)
391 spatial[i] = expf(-(l_float32)(i * i) / denom);
392 bil->spatial = spatial;
393
394 range = (l_float32 *)LEPT_CALLOC(256, sizeof(l_float32));
395 denom = 2. * range_stdev * range_stdev;
396 for (i = 0; i < 256; i++)
397 range[i] = expf(-(l_float32)(i * i) / denom);
398 bil->range = range;
399
400 /* -------------------------------------------------------------------- *
401 * Generate principal bilateral component images *
402 * -------------------------------------------------------------------- */
403 pixac = pixaCreate(ncomps);
404 pixGetDimensions(pixsc, &ws, &hs, NULL);
405 datas = pixGetData(pixsc);
406 wpls = pixGetWpl(pixsc);
407 pixGetDimensions(pixs, &w, &h, NULL);
408 wd = (w + reduction - 1) / reduction;
409 hd = (h + reduction - 1) / reduction;
410 halfwidth = (l_int32)(2.0 * sstdev);
411 for (index = 0; index < ncomps; index++) {
412 pixt = pixCopy(NULL, pixsc);
413 datat = pixGetData(pixt);
414 wplt = pixGetWpl(pixt);
415 kval = nc[index];
416 /* Separable convolutions: horizontal first */
417 for (i = 0; i < hd; i++) {
418 lines = datas + (border + i) * wpls;
419 linet = datat + (border + i) * wplt;
420 for (j = 0; j < wd; j++) {
421 sum = 0.0;
422 norm = 0.0;
423 for (k = -halfwidth; k <= halfwidth; k++) {
424 nval = GET_DATA_BYTE(lines, border + j + k);
425 kern = spatial[L_ABS(k)] * range[L_ABS(kval - nval)];
426 sum += kern * nval;
427 norm += kern;
428 }
429 if (norm > 0.0) {
430 dval = (l_int32)((sum / norm) + 0.5);
431 SET_DATA_BYTE(linet, border + j, dval);
432 }
433 }
434 }
435 /* Vertical convolution */
436 pixd = pixCreate(wd, hd, 8);
437 datad = pixGetData(pixd);
438 wpld = pixGetWpl(pixd);
439 for (i = 0; i < hd; i++) {
440 linet = datat + (border + i) * wplt;
441 lined = datad + i * wpld;
442 for (j = 0; j < wd; j++) {
443 sum = 0.0;
444 norm = 0.0;
445 for (k = -halfwidth; k <= halfwidth; k++) {
446 nval = GET_DATA_BYTE(linet + k * wplt, border + j);
447 kern = spatial[L_ABS(k)] * range[L_ABS(kval - nval)];
448 sum += kern * nval;
449 norm += kern;
450 }
451 if (norm > 0.0)
452 dval = (l_int32)((sum / norm) + 0.5);
453 else
454 dval = GET_DATA_BYTE(linet, border + j);
455 SET_DATA_BYTE(lined, j, dval);
456 }
457 }
458 pixDestroy(&pixt);
459 pixaAddPix(pixac, pixd, L_INSERT);
460 }
461 bil->pixac = pixac;
462 bil->lineset = (l_uint32 ***)pixaGetLinePtrs(pixac, NULL);
463 return bil;
464 }
465
466
467 /*!
468 * \brief bilateralApply()
469 *
470 * \param[in] bil
471 * \return pixd
472 */
473 static PIX *
474 bilateralApply(L_BILATERAL *bil)
475 {
476 l_int32 i, j, k, ired, jred, w, h, wpls, wpld, ncomps, reduction;
477 l_int32 vals, vald, lowval, hival;
478 l_int32 *kindex;
479 l_float32 fract;
480 l_float32 *kfract;
481 l_uint32 *lines, *lined, *datas, *datad;
482 l_uint32 ***lineset = NULL; /* for set of PBC */
483 PIX *pixs, *pixd;
484 PIXA *pixac;
485
486 if (!bil)
487 return (PIX *)ERROR_PTR("bil not defined", __func__, NULL);
488 pixs = bil->pixs;
489 ncomps = bil->ncomps;
490 kindex = bil->kindex;
491 kfract = bil->kfract;
492 reduction = bil->reduction;
493 pixac = bil->pixac;
494 lineset = bil->lineset;
495 if (pixaGetCount(pixac) != ncomps)
496 return (PIX *)ERROR_PTR("PBC images do not exist", __func__, NULL);
497
498 if ((pixd = pixCreateTemplate(pixs)) == NULL)
499 return (PIX *)ERROR_PTR("pixd not made", __func__, NULL);
500 datas = pixGetData(pixs);
501 wpls = pixGetWpl(pixs);
502 datad = pixGetData(pixd);
503 wpld = pixGetWpl(pixd);
504 pixGetDimensions(pixs, &w, &h, NULL);
505 for (i = 0; i < h; i++) {
506 lines = datas + i * wpls;
507 lined = datad + i * wpld;
508 ired = i / reduction;
509 for (j = 0; j < w; j++) {
510 jred = j / reduction;
511 vals = GET_DATA_BYTE(lines, j);
512 k = kindex[vals];
513 lowval = GET_DATA_BYTE(lineset[k][ired], jred);
514 hival = GET_DATA_BYTE(lineset[k + 1][ired], jred);
515 fract = kfract[vals];
516 vald = (l_int32)((1.0 - fract) * lowval + fract * hival + 0.5);
517 SET_DATA_BYTE(lined, j, vald);
518 }
519 }
520
521 return pixd;
522 }
523
524
525 /*!
526 * \brief bilateralDestroy()
527 *
528 * \param[in,out] pbil will be set to null before returning
529 */
530 static void
531 bilateralDestroy(L_BILATERAL **pbil)
532 {
533 l_int32 i;
534 L_BILATERAL *bil;
535
536 if (pbil == NULL) {
537 L_WARNING("ptr address is null!\n", __func__);
538 return;
539 }
540
541 if ((bil = *pbil) == NULL)
542 return;
543
544 pixDestroy(&bil->pixs);
545 pixDestroy(&bil->pixsc);
546 pixaDestroy(&bil->pixac);
547 LEPT_FREE(bil->spatial);
548 LEPT_FREE(bil->range);
549 LEPT_FREE(bil->nc);
550 LEPT_FREE(bil->kindex);
551 LEPT_FREE(bil->kfract);
552 for (i = 0; i < bil->ncomps; i++)
553 LEPT_FREE(bil->lineset[i]);
554 LEPT_FREE(bil->lineset);
555 LEPT_FREE(bil);
556 *pbil = NULL;
557 }
558
559
560 /*----------------------------------------------------------------------*
561 * Exact implementation of grayscale or color bilateral filtering *
562 *----------------------------------------------------------------------*/
563 /*!
564 * \brief pixBilateralExact()
565 *
566 * \param[in] pixs 8 bpp gray or 32 bpp rgb
567 * \param[in] spatial_kel gaussian kernel
568 * \param[in] range_kel [optional] 256 x 1, monotonically decreasing
569 * \return pixd 8 bpp bilateral filtered image
570 *
571 * <pre>
572 * Notes:
573 * (1) The spatial_kel is a conventional smoothing kernel, typically a
574 * 2-d Gaussian kernel or other block kernel. It can be either
575 * normalized or not, but must be everywhere positive.
576 * (2) The range_kel is defined over the absolute value of pixel
577 * grayscale differences, and hence must have size 256 x 1.
578 * Values in the array represent the multiplying weight for each
579 * gray value difference between the target pixel and center of the
580 * kernel, and should be monotonically decreasing.
581 * (3) If range_kel == NULL, a constant weight is applied regardless
582 * of the range value difference. This degenerates to a regular
583 * pixConvolve() with a normalized kernel.
584 * </pre>
585 */
586 PIX *
587 pixBilateralExact(PIX *pixs,
588 L_KERNEL *spatial_kel,
589 L_KERNEL *range_kel)
590 {
591 l_int32 d;
592 PIX *pixt, *pixr, *pixg, *pixb, *pixd;
593
594 if (!pixs)
595 return (PIX *)ERROR_PTR("pixs not defined", __func__, NULL);
596 if (pixGetColormap(pixs) != NULL)
597 return (PIX *)ERROR_PTR("pixs is cmapped", __func__, NULL);
598 d = pixGetDepth(pixs);
599 if (d != 8 && d != 32)
600 return (PIX *)ERROR_PTR("pixs not 8 or 32 bpp", __func__, NULL);
601 if (!spatial_kel)
602 return (PIX *)ERROR_PTR("spatial_ke not defined", __func__, NULL);
603
604 if (d == 8) {
605 return pixBilateralGrayExact(pixs, spatial_kel, range_kel);
606 } else { /* d == 32 */
607 pixt = pixGetRGBComponent(pixs, COLOR_RED);
608 pixr = pixBilateralGrayExact(pixt, spatial_kel, range_kel);
609 pixDestroy(&pixt);
610 pixt = pixGetRGBComponent(pixs, COLOR_GREEN);
611 pixg = pixBilateralGrayExact(pixt, spatial_kel, range_kel);
612 pixDestroy(&pixt);
613 pixt = pixGetRGBComponent(pixs, COLOR_BLUE);
614 pixb = pixBilateralGrayExact(pixt, spatial_kel, range_kel);
615 pixDestroy(&pixt);
616 pixd = pixCreateRGBImage(pixr, pixg, pixb);
617
618 pixDestroy(&pixr);
619 pixDestroy(&pixg);
620 pixDestroy(&pixb);
621 return pixd;
622 }
623 }
624
625
626 /*!
627 * \brief pixBilateralGrayExact()
628 *
629 * \param[in] pixs 8 bpp gray
630 * \param[in] spatial_kel gaussian kernel
631 * \param[in] range_kel [optional] 256 x 1, monotonically decreasing
632 * \return pixd 8 bpp bilateral filtered image
633 *
634 * <pre>
635 * Notes:
636 * (1) See pixBilateralExact().
637 * </pre>
638 */
639 PIX *
640 pixBilateralGrayExact(PIX *pixs,
641 L_KERNEL *spatial_kel,
642 L_KERNEL *range_kel)
643 {
644 l_int32 i, j, id, jd, k, m, w, h, d, sx, sy, cx, cy, wplt, wpld;
645 l_int32 val, center_val;
646 l_uint32 *datat, *datad, *linet, *lined;
647 l_float32 sum, weight_sum, weight;
648 L_KERNEL *keli;
649 PIX *pixt, *pixd;
650
651 if (!pixs)
652 return (PIX *)ERROR_PTR("pixs not defined", __func__, NULL);
653 if (pixGetDepth(pixs) != 8)
654 return (PIX *)ERROR_PTR("pixs must be gray", __func__, NULL);
655 pixGetDimensions(pixs, &w, &h, &d);
656 if (!spatial_kel)
657 return (PIX *)ERROR_PTR("spatial kel not defined", __func__, NULL);
658 kernelGetParameters(spatial_kel, &sy, &sx, NULL, NULL);
659 if (w < 2 * sx + 1 || h < 2 * sy + 1) {
660 L_WARNING("w = %d < 2 * sx + 1 = %d, or h = %d < 2 * sy + 1 = %d; "
661 "returning copy\n", __func__, w, 2 * sx + 1, h, 2 * sy + 1);
662 return pixCopy(NULL, pixs);
663 }
664 if (!range_kel)
665 return pixConvolve(pixs, spatial_kel, 8, 1);
666 if (range_kel->sx != 256 || range_kel->sy != 1)
667 return (PIX *)ERROR_PTR("range kel not {256 x 1", __func__, NULL);
668
669 keli = kernelInvert(spatial_kel);
670 kernelGetParameters(keli, &sy, &sx, &cy, &cx);
671 if ((pixt = pixAddMirroredBorder(pixs, cx, sx - cx, cy, sy - cy)) == NULL) {
672 kernelDestroy(&keli);
673 return (PIX *)ERROR_PTR("pixt not made", __func__, NULL);
674 }
675
676 pixd = pixCreate(w, h, 8);
677 datat = pixGetData(pixt);
678 datad = pixGetData(pixd);
679 wplt = pixGetWpl(pixt);
680 wpld = pixGetWpl(pixd);
681 for (i = 0, id = 0; id < h; i++, id++) {
682 lined = datad + id * wpld;
683 for (j = 0, jd = 0; jd < w; j++, jd++) {
684 center_val = GET_DATA_BYTE(datat + (i + cy) * wplt, j + cx);
685 weight_sum = 0.0;
686 sum = 0.0;
687 for (k = 0; k < sy; k++) {
688 linet = datat + (i + k) * wplt;
689 for (m = 0; m < sx; m++) {
690 val = GET_DATA_BYTE(linet, j + m);
691 weight = keli->data[k][m] *
692 range_kel->data[0][L_ABS(center_val - val)];
693 weight_sum += weight;
694 sum += val * weight;
695 }
696 }
697 SET_DATA_BYTE(lined, jd, (l_int32)(sum / weight_sum + 0.5));
698 }
699 }
700
701 kernelDestroy(&keli);
702 pixDestroy(&pixt);
703 return pixd;
704 }
705
706
707 /*!
708 * \brief pixBlockBilateralExact()
709 *
710 * \param[in] pixs 8 bpp gray or 32 bpp rgb
711 * \param[in] spatial_stdev must be > 0.0
712 * \param[in] range_stdev must be > 0.0
713 * \return pixd 8 bpp or 32 bpp bilateral filtered image
714 *
715 * <pre>
716 * Notes:
717 * (1) See pixBilateralExact(). This provides an interface using
718 * the standard deviations of the spatial and range filters.
719 * (2) The convolution window halfwidth is 2 * spatial_stdev,
720 * and the square filter size is 4 * spatial_stdev + 1.
721 * The kernel captures 95% of total energy. This is compensated
722 * by normalization.
723 * (3) The range_stdev is analogous to spatial_halfwidth in the
724 * grayscale domain [0...255], and determines how much damping of the
725 * smoothing operation is applied across edges. The larger this
726 * value is, the smaller the damping. The smaller the value, the
727 * more edge details are preserved. These approximations are useful
728 * for deciding the appropriate cutoff.
729 * kernel[1 * stdev] ~= 0.6 * kernel[0]
730 * kernel[2 * stdev] ~= 0.14 * kernel[0]
731 * kernel[3 * stdev] ~= 0.01 * kernel[0]
732 * If range_stdev is infinite there is no damping, and this
733 * becomes a conventional gaussian smoothing.
734 * This value does not affect the run time.
735 * (4) If range_stdev is negative or zero, the range kernel is
736 * ignored and this degenerates to a straight gaussian convolution.
737 * (5) This is very slow for large spatial filters. The time
738 * on a 3GHz pentium is roughly
739 * T = 1.2 * 10^-8 * (A * sh^2) sec
740 * where A = # of pixels, sh = spatial halfwidth of filter.
741 * </pre>
742 */
743 PIX*
744 pixBlockBilateralExact(PIX *pixs,
745 l_float32 spatial_stdev,
746 l_float32 range_stdev)
747 {
748 l_int32 d, halfwidth;
749 L_KERNEL *spatial_kel, *range_kel;
750 PIX *pixd;
751
752 if (!pixs)
753 return (PIX *)ERROR_PTR("pixs not defined", __func__, NULL);
754 d = pixGetDepth(pixs);
755 if (d != 8 && d != 32)
756 return (PIX *)ERROR_PTR("pixs not 8 or 32 bpp", __func__, NULL);
757 if (pixGetColormap(pixs) != NULL)
758 return (PIX *)ERROR_PTR("pixs is cmapped", __func__, NULL);
759 if (spatial_stdev <= 0.0)
760 return (PIX *)ERROR_PTR("invalid spatial stdev", __func__, NULL);
761 if (range_stdev <= 0.0)
762 return (PIX *)ERROR_PTR("invalid range stdev", __func__, NULL);
763
764 halfwidth = 2 * spatial_stdev;
765 spatial_kel = makeGaussianKernel(halfwidth, halfwidth, spatial_stdev, 1.0);
766 range_kel = makeRangeKernel(range_stdev);
767 pixd = pixBilateralExact(pixs, spatial_kel, range_kel);
768 kernelDestroy(&spatial_kel);
769 kernelDestroy(&range_kel);
770 return pixd;
771 }
772
773
774 /*----------------------------------------------------------------------*
775 * Kernel helper function *
776 *----------------------------------------------------------------------*/
777 /*!
778 * \brief makeRangeKernel()
779 *
780 * \param[in] range_stdev must be > 0.0
781 * \return kel, or NULL on error
782 *
783 * <pre>
784 * Notes:
785 * (1) Creates a one-sided Gaussian kernel with the given
786 * standard deviation. At grayscale difference of one stdev,
787 * the kernel falls to 0.6, and to 0.01 at three stdev.
788 * (2) A typical input number might be 20. Then pixels whose
789 * value differs by 60 from the center pixel have their
790 * weight in the convolution reduced by a factor of about 0.01.
791 * </pre>
792 */
793 L_KERNEL *
794 makeRangeKernel(l_float32 range_stdev)
795 {
796 l_int32 x;
797 l_float32 val, denom;
798 L_KERNEL *kel;
799
800 if (range_stdev <= 0.0)
801 return (L_KERNEL *)ERROR_PTR("invalid stdev <= 0", __func__, NULL);
802
803 denom = 2. * range_stdev * range_stdev;
804 if ((kel = kernelCreate(1, 256)) == NULL)
805 return (L_KERNEL *)ERROR_PTR("kel not made", __func__, NULL);
806 kernelSetOrigin(kel, 0, 0);
807 for (x = 0; x < 256; x++) {
808 val = expf(-(l_float32)(x * x) / denom);
809 kernelSetElement(kel, 0, x, val);
810 }
811 return kel;
812 }