mirror of http://192.168.1.51:8099/lmh188/twain3.0
811 lines
30 KiB
C
811 lines
30 KiB
C
/*====================================================================*
|
|
- Copyright (C) 2001 Leptonica. All rights reserved.
|
|
-
|
|
- Redistribution and use in source and binary forms, with or without
|
|
- modification, are permitted provided that the following conditions
|
|
- are met:
|
|
- 1. Redistributions of source code must retain the above copyright
|
|
- notice, this list of conditions and the following disclaimer.
|
|
- 2. Redistributions in binary form must reproduce the above
|
|
- copyright notice, this list of conditions and the following
|
|
- disclaimer in the documentation and/or other materials
|
|
- provided with the distribution.
|
|
-
|
|
- THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
|
|
- ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
|
|
- LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
|
|
- A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL ANY
|
|
- CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
|
|
- EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
|
|
- PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
|
|
- PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY
|
|
- OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
|
|
- NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
|
|
- SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
|
|
*====================================================================*/
|
|
|
|
/*!
|
|
* \file bilateral.c
|
|
* <pre>
|
|
*
|
|
* Top level approximate separable grayscale or color bilateral filtering
|
|
* PIX *pixBilateral()
|
|
* PIX *pixBilateralGray()
|
|
*
|
|
* Implementation of approximate separable bilateral filter
|
|
* static L_BILATERAL *bilateralCreate()
|
|
* static void *bilateralDestroy()
|
|
* static PIX *bilateralApply()
|
|
*
|
|
* Slow, exact implementation of grayscale or color bilateral filtering
|
|
* PIX *pixBilateralExact()
|
|
* PIX *pixBilateralGrayExact()
|
|
* PIX *pixBlockBilateralExact()
|
|
*
|
|
* Kernel helper function
|
|
* L_KERNEL *makeRangeKernel()
|
|
*
|
|
* This includes both a slow, exact implementation of the bilateral
|
|
* filter algorithm (given by Sylvain Paris and Frédo Durand),
|
|
* and a fast, approximate and separable implementation (following
|
|
* Yang, Tan and Ahuja). See bilateral.h for algorithmic details.
|
|
*
|
|
* The bilateral filter has the nice property of applying a gaussian
|
|
* filter to smooth parts of the image that don't vary too quickly,
|
|
* while at the same time preserving edges. The filter is nonlinear
|
|
* and cannot be decomposed into two separable filters; however,
|
|
* there exists an approximate method that is separable. To further
|
|
* speed up the separable implementation, you can generate the
|
|
* intermediate data at reduced resolution.
|
|
*
|
|
* The full kernel is composed of two parts: a spatial gaussian filter
|
|
* and a nonlinear "range" filter that depends on the intensity difference
|
|
* between the reference pixel at the spatial kernel origin and any other
|
|
* pixel within the kernel support.
|
|
*
|
|
* In our implementations, the range filter is a parameterized,
|
|
* one-sided, 256-element, monotonically decreasing gaussian function
|
|
* of the absolute value of the difference between pixel values; namely,
|
|
* abs(I2 - I1). In general, any decreasing function can be used,
|
|
* and more generally, any two-dimensional kernel can be used if
|
|
* you wish to relax the 'abs' condition. (In that case, the range
|
|
* filter can be 256 x 256).
|
|
* </pre>
|
|
*/
|
|
|
|
#include <math.h>
|
|
#include "allheaders.h"
|
|
#include "bilateral.h"
|
|
|
|
static L_BILATERAL *bilateralCreate(PIX *pixs, l_float32 spatial_stdev,
|
|
l_float32 range_stdev, l_int32 ncomps,
|
|
l_int32 reduction);
|
|
static PIX *bilateralApply(L_BILATERAL *bil);
|
|
static void bilateralDestroy(L_BILATERAL **pbil);
|
|
|
|
|
|
#ifndef NO_CONSOLE_IO
|
|
#define DEBUG_BILATERAL 0
|
|
#endif /* ~NO_CONSOLE_IO */
|
|
|
|
|
|
/*--------------------------------------------------------------------------*
|
|
* Top level approximate separable grayscale or color bilateral filtering *
|
|
*--------------------------------------------------------------------------*/
|
|
/*!
|
|
* \brief pixBilateral()
|
|
*
|
|
* \param[in] pixs 8 bpp gray or 32 bpp rgb, no colormap
|
|
* \param[in] spatial_stdev of gaussian kernel; in pixels, > 0.5
|
|
* \param[in] range_stdev of gaussian range kernel; > 5.0; typ. 50.0
|
|
* \param[in] ncomps number of intermediate sums J(k,x);
|
|
* in [4 ... 30]
|
|
* \param[in] reduction 1, 2 or 4
|
|
* \return pixd bilateral filtered image, or NULL on error
|
|
*
|
|
* <pre>
|
|
* Notes:
|
|
* (1) This performs a relatively fast, separable bilateral
|
|
* filtering operation. The time is proportional to ncomps
|
|
* and varies inversely approximately as the cube of the
|
|
* reduction factor. See bilateral.h for algorithm details.
|
|
* (2) We impose minimum values for range_stdev and ncomps to
|
|
* avoid nasty artifacts when either are too small. We also
|
|
* impose a constraint on their product:
|
|
* ncomps * range_stdev >= 100.
|
|
* So for values of range_stdev >= 25, ncomps can be as small as 4.
|
|
* Here is a qualitative, intuitive explanation for this constraint.
|
|
* Call the difference in k values between the J(k) == 'delta', where
|
|
* 'delta' ~ 200 / ncomps
|
|
* Then this constraint is roughly equivalent to the condition:
|
|
* 'delta' < 2 * range_stdev
|
|
* Note that at an intensity difference of (2 * range_stdev), the
|
|
* range part of the kernel reduces the effect by the factor 0.14.
|
|
* This constraint requires that we have a sufficient number of
|
|
* PCBs (i.e, a small enough 'delta'), so that for any value of
|
|
* image intensity I, there exists a k (and a PCB, J(k), such that
|
|
* |I - k| < range_stdev
|
|
* Any fewer PCBs and we don't have enough to support this condition.
|
|
* (3) The upper limit of 30 on ncomps is imposed because the
|
|
* gain in accuracy is not worth the extra computation.
|
|
* (4) The size of the gaussian kernel is twice the spatial_stdev
|
|
* on each side of the origin. The minimum value of
|
|
* spatial_stdev, 0.5, is required to have a finite sized
|
|
* spatial kernel. In practice, a much larger value is used.
|
|
* (5) Computation of the intermediate images goes inversely
|
|
* as the cube of the reduction factor. If you can use a
|
|
* reduction of 2 or 4, it is well-advised.
|
|
* (6) The range kernel is defined over the absolute value of pixel
|
|
* grayscale differences, and hence must have size 256 x 1.
|
|
* Values in the array represent the multiplying weight
|
|
* depending on the absolute gray value difference between
|
|
* the source pixel and the neighboring pixel, and should
|
|
* be monotonically decreasing.
|
|
* (7) Interesting observation. Run this on prog/fish24.jpg, with
|
|
* range_stdev = 60, ncomps = 6, and spatial_dev = {10, 30, 50}.
|
|
* As spatial_dev gets larger, we get the counter-intuitive
|
|
* result that the body of the red fish becomes less blurry.
|
|
* </pre>
|
|
*/
|
|
PIX *
|
|
pixBilateral(PIX *pixs,
|
|
l_float32 spatial_stdev,
|
|
l_float32 range_stdev,
|
|
l_int32 ncomps,
|
|
l_int32 reduction)
|
|
{
|
|
l_int32 d;
|
|
l_float32 sstdev; /* scaled spatial stdev */
|
|
PIX *pixt, *pixr, *pixg, *pixb, *pixd;
|
|
|
|
PROCNAME("pixBilateral");
|
|
|
|
if (!pixs || pixGetColormap(pixs))
|
|
return (PIX *)ERROR_PTR("pixs not defined or cmapped", procName, NULL);
|
|
d = pixGetDepth(pixs);
|
|
if (d != 8 && d != 32)
|
|
return (PIX *)ERROR_PTR("pixs not 8 or 32 bpp", procName, NULL);
|
|
if (reduction != 1 && reduction != 2 && reduction != 4)
|
|
return (PIX *)ERROR_PTR("reduction invalid", procName, NULL);
|
|
sstdev = spatial_stdev / (l_float32)reduction; /* reduced spat. stdev */
|
|
if (sstdev < 0.5)
|
|
return (PIX *)ERROR_PTR("sstdev < 0.5", procName, NULL);
|
|
if (range_stdev <= 5.0)
|
|
return (PIX *)ERROR_PTR("range_stdev <= 5.0", procName, NULL);
|
|
if (ncomps < 4 || ncomps > 30)
|
|
return (PIX *)ERROR_PTR("ncomps not in [4 ... 30]", procName, NULL);
|
|
if (ncomps * range_stdev < 100.0)
|
|
return (PIX *)ERROR_PTR("ncomps * range_stdev < 100.0", procName, NULL);
|
|
|
|
if (d == 8)
|
|
return pixBilateralGray(pixs, spatial_stdev, range_stdev,
|
|
ncomps, reduction);
|
|
|
|
pixt = pixGetRGBComponent(pixs, COLOR_RED);
|
|
pixr = pixBilateralGray(pixt, spatial_stdev, range_stdev, ncomps,
|
|
reduction);
|
|
pixDestroy(&pixt);
|
|
pixt = pixGetRGBComponent(pixs, COLOR_GREEN);
|
|
pixg = pixBilateralGray(pixt, spatial_stdev, range_stdev, ncomps,
|
|
reduction);
|
|
pixDestroy(&pixt);
|
|
pixt = pixGetRGBComponent(pixs, COLOR_BLUE);
|
|
pixb = pixBilateralGray(pixt, spatial_stdev, range_stdev, ncomps,
|
|
reduction);
|
|
pixDestroy(&pixt);
|
|
pixd = pixCreateRGBImage(pixr, pixg, pixb);
|
|
pixDestroy(&pixr);
|
|
pixDestroy(&pixg);
|
|
pixDestroy(&pixb);
|
|
return pixd;
|
|
}
|
|
|
|
|
|
/*!
|
|
* \brief pixBilateralGray()
|
|
*
|
|
* \param[in] pixs 8 bpp gray
|
|
* \param[in] spatial_stdev of gaussian kernel; in pixels, > 0.5
|
|
* \param[in] range_stdev of gaussian range kernel; > 5.0; typ. 50.0
|
|
* \param[in] ncomps number of intermediate sums J(k,x);
|
|
* in [4 ... 30]
|
|
* \param[in] reduction 1, 2 or 4
|
|
* \return pixd 8 bpp bilateral filtered image, or NULL on error
|
|
*
|
|
* <pre>
|
|
* Notes:
|
|
* (1) See pixBilateral() for constraints on the input parameters.
|
|
* (2) See pixBilateral() for algorithm details.
|
|
* </pre>
|
|
*/
|
|
PIX *
|
|
pixBilateralGray(PIX *pixs,
|
|
l_float32 spatial_stdev,
|
|
l_float32 range_stdev,
|
|
l_int32 ncomps,
|
|
l_int32 reduction)
|
|
{
|
|
l_float32 sstdev; /* scaled spatial stdev */
|
|
PIX *pixd;
|
|
L_BILATERAL *bil;
|
|
|
|
PROCNAME("pixBilateralGray");
|
|
|
|
if (!pixs || pixGetColormap(pixs))
|
|
return (PIX *)ERROR_PTR("pixs not defined or cmapped", procName, NULL);
|
|
if (pixGetDepth(pixs) != 8)
|
|
return (PIX *)ERROR_PTR("pixs not 8 bpp gray", procName, NULL);
|
|
if (reduction != 1 && reduction != 2 && reduction != 4)
|
|
return (PIX *)ERROR_PTR("reduction invalid", procName, NULL);
|
|
sstdev = spatial_stdev / (l_float32)reduction; /* reduced spat. stdev */
|
|
if (sstdev < 0.5)
|
|
return (PIX *)ERROR_PTR("sstdev < 0.5", procName, NULL);
|
|
if (range_stdev <= 5.0)
|
|
return (PIX *)ERROR_PTR("range_stdev <= 5.0", procName, NULL);
|
|
if (ncomps < 4 || ncomps > 30)
|
|
return (PIX *)ERROR_PTR("ncomps not in [4 ... 30]", procName, NULL);
|
|
if (ncomps * range_stdev < 100.0)
|
|
return (PIX *)ERROR_PTR("ncomps * range_stdev < 100.0", procName, NULL);
|
|
|
|
bil = bilateralCreate(pixs, spatial_stdev, range_stdev, ncomps, reduction);
|
|
if (!bil) return (PIX *)ERROR_PTR("bil not made", procName, NULL);
|
|
pixd = bilateralApply(bil);
|
|
bilateralDestroy(&bil);
|
|
return pixd;
|
|
}
|
|
|
|
|
|
/*----------------------------------------------------------------------*
|
|
* Implementation of approximate separable bilateral filter *
|
|
*----------------------------------------------------------------------*/
|
|
/*!
|
|
* \brief bilateralCreate()
|
|
*
|
|
* \param[in] pixs 8 bpp gray, no colormap
|
|
* \param[in] spatial_stdev of gaussian kernel; in pixels, > 0.5
|
|
* \param[in] range_stdev of gaussian range kernel; > 5.0; typ. 50.0
|
|
* \param[in] ncomps number of intermediate sums J(k,x);
|
|
* in [4 ... 30]
|
|
* \param[in] reduction 1, 2 or 4
|
|
* \return bil, or NULL on error
|
|
*
|
|
* <pre>
|
|
* Notes:
|
|
* (1) This initializes a bilateral filtering operation, generating all
|
|
* the data required. It takes most of the time in the bilateral
|
|
* filtering operation.
|
|
* (2) See bilateral.h for details of the algorithm.
|
|
* (3) See pixBilateral() for constraints on input parameters, which
|
|
* are not checked here.
|
|
* </pre>
|
|
*/
|
|
static L_BILATERAL *
|
|
bilateralCreate(PIX *pixs,
|
|
l_float32 spatial_stdev,
|
|
l_float32 range_stdev,
|
|
l_int32 ncomps,
|
|
l_int32 reduction)
|
|
{
|
|
l_int32 w, ws, wd, h, hs, hd, i, j, k, index;
|
|
l_int32 border, minval, maxval, spatial_size;
|
|
l_int32 halfwidth, wpls, wplt, wpld, kval, nval, dval;
|
|
l_float32 sstdev, fval1, fval2, denom, sum, norm, kern;
|
|
l_int32 *nc, *kindex;
|
|
l_float32 *kfract, *range, *spatial;
|
|
l_uint32 *datas, *datat, *datad, *lines, *linet, *lined;
|
|
L_BILATERAL *bil;
|
|
PIX *pixt, *pixt2, *pixsc, *pixd;
|
|
PIXA *pixac;
|
|
|
|
PROCNAME("bilateralCreate");
|
|
|
|
sstdev = spatial_stdev / (l_float32)reduction; /* reduced spat. stdev */
|
|
if ((bil = (L_BILATERAL *)LEPT_CALLOC(1, sizeof(L_BILATERAL))) == NULL)
|
|
return (L_BILATERAL *)ERROR_PTR("bil not made", procName, NULL);
|
|
bil->spatial_stdev = sstdev;
|
|
bil->range_stdev = range_stdev;
|
|
bil->reduction = reduction;
|
|
bil->ncomps = ncomps;
|
|
|
|
if (reduction == 1) {
|
|
pixt = pixClone(pixs);
|
|
} else if (reduction == 2) {
|
|
pixt = pixScaleAreaMap2(pixs);
|
|
} else { /* reduction == 4) */
|
|
pixt2 = pixScaleAreaMap2(pixs);
|
|
pixt = pixScaleAreaMap2(pixt2);
|
|
pixDestroy(&pixt2);
|
|
}
|
|
|
|
pixGetExtremeValue(pixt, 1, L_SELECT_MIN, NULL, NULL, NULL, &minval);
|
|
pixGetExtremeValue(pixt, 1, L_SELECT_MAX, NULL, NULL, NULL, &maxval);
|
|
bil->minval = minval;
|
|
bil->maxval = maxval;
|
|
|
|
border = (l_int32)(2 * sstdev + 1);
|
|
pixsc = pixAddMirroredBorder(pixt, border, border, border, border);
|
|
bil->pixsc = pixsc;
|
|
pixDestroy(&pixt);
|
|
bil->pixs = pixClone(pixs);
|
|
|
|
|
|
/* -------------------------------------------------------------------- *
|
|
* Generate arrays for interpolation of J(k,x):
|
|
* (1.0 - kfract[.]) * J(kindex[.], x) + kfract[.] * J(kindex[.] + 1, x),
|
|
* where I(x) is the index into kfract[] and kindex[],
|
|
* and x is an index into the 2D image array.
|
|
* -------------------------------------------------------------------- */
|
|
/* nc is the set of k values to be used in J(k,x) */
|
|
nc = (l_int32 *)LEPT_CALLOC(ncomps, sizeof(l_int32));
|
|
for (i = 0; i < ncomps; i++)
|
|
nc[i] = minval + i * (maxval - minval) / (ncomps - 1);
|
|
bil->nc = nc;
|
|
|
|
/* kindex maps from intensity I(x) to the lower k index for J(k,x) */
|
|
kindex = (l_int32 *)LEPT_CALLOC(256, sizeof(l_int32));
|
|
for (i = minval, k = 0; i <= maxval && k < ncomps - 1; k++) {
|
|
fval2 = nc[k + 1];
|
|
while (i < fval2) {
|
|
kindex[i] = k;
|
|
i++;
|
|
}
|
|
}
|
|
kindex[maxval] = ncomps - 2;
|
|
bil->kindex = kindex;
|
|
|
|
/* kfract maps from intensity I(x) to the fraction of J(k+1,x) used */
|
|
kfract = (l_float32 *)LEPT_CALLOC(256, sizeof(l_float32)); /* from lower */
|
|
for (i = minval, k = 0; i <= maxval && k < ncomps - 1; k++) {
|
|
fval1 = nc[k];
|
|
fval2 = nc[k + 1];
|
|
while (i < fval2) {
|
|
kfract[i] = (l_float32)(i - fval1) / (l_float32)(fval2 - fval1);
|
|
i++;
|
|
}
|
|
}
|
|
kfract[maxval] = 1.0;
|
|
bil->kfract = kfract;
|
|
|
|
#if DEBUG_BILATERAL
|
|
for (i = minval; i <= maxval; i++)
|
|
fprintf(stderr, "kindex[%d] = %d; kfract[%d] = %5.3f\n",
|
|
i, kindex[i], i, kfract[i]);
|
|
for (i = 0; i < ncomps; i++)
|
|
fprintf(stderr, "nc[%d] = %d\n", i, nc[i]);
|
|
#endif /* DEBUG_BILATERAL */
|
|
|
|
|
|
/* -------------------------------------------------------------------- *
|
|
* Generate 1-D kernel arrays (spatial and range) *
|
|
* -------------------------------------------------------------------- */
|
|
spatial_size = 2 * sstdev + 1;
|
|
spatial = (l_float32 *)LEPT_CALLOC(spatial_size, sizeof(l_float32));
|
|
denom = 2. * sstdev * sstdev;
|
|
for (i = 0; i < spatial_size; i++)
|
|
spatial[i] = expf(-(l_float32)(i * i) / denom);
|
|
bil->spatial = spatial;
|
|
|
|
range = (l_float32 *)LEPT_CALLOC(256, sizeof(l_float32));
|
|
denom = 2. * range_stdev * range_stdev;
|
|
for (i = 0; i < 256; i++)
|
|
range[i] = expf(-(l_float32)(i * i) / denom);
|
|
bil->range = range;
|
|
|
|
|
|
/* -------------------------------------------------------------------- *
|
|
* Generate principal bilateral component images *
|
|
* -------------------------------------------------------------------- */
|
|
pixac = pixaCreate(ncomps);
|
|
pixGetDimensions(pixsc, &ws, &hs, NULL);
|
|
datas = pixGetData(pixsc);
|
|
wpls = pixGetWpl(pixsc);
|
|
pixGetDimensions(pixs, &w, &h, NULL);
|
|
wd = (w + reduction - 1) / reduction;
|
|
hd = (h + reduction - 1) / reduction;
|
|
halfwidth = (l_int32)(2.0 * sstdev);
|
|
for (index = 0; index < ncomps; index++) {
|
|
pixt = pixCopy(NULL, pixsc);
|
|
datat = pixGetData(pixt);
|
|
wplt = pixGetWpl(pixt);
|
|
kval = nc[index];
|
|
/* Separable convolutions: horizontal first */
|
|
for (i = 0; i < hd; i++) {
|
|
lines = datas + (border + i) * wpls;
|
|
linet = datat + (border + i) * wplt;
|
|
for (j = 0; j < wd; j++) {
|
|
sum = 0.0;
|
|
norm = 0.0;
|
|
for (k = -halfwidth; k <= halfwidth; k++) {
|
|
nval = GET_DATA_BYTE(lines, border + j + k);
|
|
kern = spatial[L_ABS(k)] * range[L_ABS(kval - nval)];
|
|
sum += kern * nval;
|
|
norm += kern;
|
|
}
|
|
dval = (l_int32)((sum / norm) + 0.5);
|
|
SET_DATA_BYTE(linet, border + j, dval);
|
|
}
|
|
}
|
|
/* Vertical convolution */
|
|
pixd = pixCreate(wd, hd, 8);
|
|
datad = pixGetData(pixd);
|
|
wpld = pixGetWpl(pixd);
|
|
for (i = 0; i < hd; i++) {
|
|
linet = datat + (border + i) * wplt;
|
|
lined = datad + i * wpld;
|
|
for (j = 0; j < wd; j++) {
|
|
sum = 0.0;
|
|
norm = 0.0;
|
|
for (k = -halfwidth; k <= halfwidth; k++) {
|
|
nval = GET_DATA_BYTE(linet + k * wplt, border + j);
|
|
kern = spatial[L_ABS(k)] * range[L_ABS(kval - nval)];
|
|
sum += kern * nval;
|
|
norm += kern;
|
|
}
|
|
dval = (l_int32)((sum / norm) + 0.5);
|
|
SET_DATA_BYTE(lined, j, dval);
|
|
}
|
|
}
|
|
pixDestroy(&pixt);
|
|
pixaAddPix(pixac, pixd, L_INSERT);
|
|
}
|
|
bil->pixac = pixac;
|
|
bil->lineset = (l_uint32 ***)pixaGetLinePtrs(pixac, NULL);
|
|
|
|
return bil;
|
|
}
|
|
|
|
|
|
/*!
|
|
* \brief bilateralApply()
|
|
*
|
|
* \param[in] bil
|
|
* \return pixd
|
|
*/
|
|
static PIX *
|
|
bilateralApply(L_BILATERAL *bil)
|
|
{
|
|
l_int32 i, j, k, ired, jred, w, h, wpls, wpld, ncomps, reduction;
|
|
l_int32 vals, vald, lowval, hival;
|
|
l_int32 *kindex;
|
|
l_float32 fract;
|
|
l_float32 *kfract;
|
|
l_uint32 *lines, *lined, *datas, *datad;
|
|
l_uint32 ***lineset = NULL; /* for set of PBC */
|
|
PIX *pixs, *pixd;
|
|
PIXA *pixac;
|
|
|
|
PROCNAME("bilateralApply");
|
|
|
|
if (!bil)
|
|
return (PIX *)ERROR_PTR("bil not defined", procName, NULL);
|
|
pixs = bil->pixs;
|
|
ncomps = bil->ncomps;
|
|
kindex = bil->kindex;
|
|
kfract = bil->kfract;
|
|
reduction = bil->reduction;
|
|
pixac = bil->pixac;
|
|
lineset = bil->lineset;
|
|
if (pixaGetCount(pixac) != ncomps)
|
|
return (PIX *)ERROR_PTR("PBC images do not exist", procName, NULL);
|
|
|
|
if ((pixd = pixCreateTemplate(pixs)) == NULL)
|
|
return (PIX *)ERROR_PTR("pixd not made", procName, NULL);
|
|
datas = pixGetData(pixs);
|
|
wpls = pixGetWpl(pixs);
|
|
datad = pixGetData(pixd);
|
|
wpld = pixGetWpl(pixd);
|
|
pixGetDimensions(pixs, &w, &h, NULL);
|
|
for (i = 0; i < h; i++) {
|
|
lines = datas + i * wpls;
|
|
lined = datad + i * wpld;
|
|
ired = i / reduction;
|
|
for (j = 0; j < w; j++) {
|
|
jred = j / reduction;
|
|
vals = GET_DATA_BYTE(lines, j);
|
|
k = kindex[vals];
|
|
lowval = GET_DATA_BYTE(lineset[k][ired], jred);
|
|
hival = GET_DATA_BYTE(lineset[k + 1][ired], jred);
|
|
fract = kfract[vals];
|
|
vald = (l_int32)((1.0 - fract) * lowval + fract * hival + 0.5);
|
|
SET_DATA_BYTE(lined, j, vald);
|
|
}
|
|
}
|
|
|
|
return pixd;
|
|
}
|
|
|
|
|
|
/*!
|
|
* \brief bilateralDestroy()
|
|
*
|
|
* \param[in,out] pbil will be set to null before returning
|
|
*/
|
|
static void
|
|
bilateralDestroy(L_BILATERAL **pbil)
|
|
{
|
|
l_int32 i;
|
|
L_BILATERAL *bil;
|
|
|
|
PROCNAME("bilateralDestroy");
|
|
|
|
if (pbil == NULL) {
|
|
L_WARNING("ptr address is null!\n", procName);
|
|
return;
|
|
}
|
|
|
|
if ((bil = *pbil) == NULL)
|
|
return;
|
|
|
|
pixDestroy(&bil->pixs);
|
|
pixDestroy(&bil->pixsc);
|
|
pixaDestroy(&bil->pixac);
|
|
LEPT_FREE(bil->spatial);
|
|
LEPT_FREE(bil->range);
|
|
LEPT_FREE(bil->nc);
|
|
LEPT_FREE(bil->kindex);
|
|
LEPT_FREE(bil->kfract);
|
|
for (i = 0; i < bil->ncomps; i++)
|
|
LEPT_FREE(bil->lineset[i]);
|
|
LEPT_FREE(bil->lineset);
|
|
LEPT_FREE(bil);
|
|
*pbil = NULL;
|
|
return;
|
|
}
|
|
|
|
|
|
/*----------------------------------------------------------------------*
|
|
* Exact implementation of grayscale or color bilateral filtering *
|
|
*----------------------------------------------------------------------*/
|
|
/*!
|
|
* \brief pixBilateralExact()
|
|
*
|
|
* \param[in] pixs 8 bpp gray or 32 bpp rgb
|
|
* \param[in] spatial_kel gaussian kernel
|
|
* \param[in] range_kel [optional] 256 x 1, monotonically decreasing
|
|
* \return pixd 8 bpp bilateral filtered image
|
|
*
|
|
* <pre>
|
|
* Notes:
|
|
* (1) The spatial_kel is a conventional smoothing kernel, typically a
|
|
* 2-d Gaussian kernel or other block kernel. It can be either
|
|
* normalized or not, but must be everywhere positive.
|
|
* (2) The range_kel is defined over the absolute value of pixel
|
|
* grayscale differences, and hence must have size 256 x 1.
|
|
* Values in the array represent the multiplying weight for each
|
|
* gray value difference between the target pixel and center of the
|
|
* kernel, and should be monotonically decreasing.
|
|
* (3) If range_kel == NULL, a constant weight is applied regardless
|
|
* of the range value difference. This degenerates to a regular
|
|
* pixConvolve() with a normalized kernel.
|
|
* </pre>
|
|
*/
|
|
PIX *
|
|
pixBilateralExact(PIX *pixs,
|
|
L_KERNEL *spatial_kel,
|
|
L_KERNEL *range_kel)
|
|
{
|
|
l_int32 d;
|
|
PIX *pixt, *pixr, *pixg, *pixb, *pixd;
|
|
|
|
PROCNAME("pixBilateralExact");
|
|
if (!pixs)
|
|
return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
|
|
if (pixGetColormap(pixs) != NULL)
|
|
return (PIX *)ERROR_PTR("pixs is cmapped", procName, NULL);
|
|
d = pixGetDepth(pixs);
|
|
if (d != 8 && d != 32)
|
|
return (PIX *)ERROR_PTR("pixs not 8 or 32 bpp", procName, NULL);
|
|
if (!spatial_kel)
|
|
return (PIX *)ERROR_PTR("spatial_ke not defined", procName, NULL);
|
|
|
|
if (d == 8) {
|
|
return pixBilateralGrayExact(pixs, spatial_kel, range_kel);
|
|
} else { /* d == 32 */
|
|
pixt = pixGetRGBComponent(pixs, COLOR_RED);
|
|
pixr = pixBilateralGrayExact(pixt, spatial_kel, range_kel);
|
|
pixDestroy(&pixt);
|
|
pixt = pixGetRGBComponent(pixs, COLOR_GREEN);
|
|
pixg = pixBilateralGrayExact(pixt, spatial_kel, range_kel);
|
|
pixDestroy(&pixt);
|
|
pixt = pixGetRGBComponent(pixs, COLOR_BLUE);
|
|
pixb = pixBilateralGrayExact(pixt, spatial_kel, range_kel);
|
|
pixDestroy(&pixt);
|
|
pixd = pixCreateRGBImage(pixr, pixg, pixb);
|
|
|
|
pixDestroy(&pixr);
|
|
pixDestroy(&pixg);
|
|
pixDestroy(&pixb);
|
|
return pixd;
|
|
}
|
|
}
|
|
|
|
|
|
/*!
|
|
* \brief pixBilateralGrayExact()
|
|
*
|
|
* \param[in] pixs 8 bpp gray
|
|
* \param[in] spatial_kel gaussian kernel
|
|
* \param[in] range_kel [optional] 256 x 1, monotonically decreasing
|
|
* \return pixd 8 bpp bilateral filtered image
|
|
*
|
|
* <pre>
|
|
* Notes:
|
|
* (1) See pixBilateralExact().
|
|
* </pre>
|
|
*/
|
|
PIX *
|
|
pixBilateralGrayExact(PIX *pixs,
|
|
L_KERNEL *spatial_kel,
|
|
L_KERNEL *range_kel)
|
|
{
|
|
l_int32 i, j, id, jd, k, m, w, h, d, sx, sy, cx, cy, wplt, wpld;
|
|
l_int32 val, center_val;
|
|
l_uint32 *datat, *datad, *linet, *lined;
|
|
l_float32 sum, weight_sum, weight;
|
|
L_KERNEL *keli;
|
|
PIX *pixt, *pixd;
|
|
|
|
PROCNAME("pixBilateralGrayExact");
|
|
|
|
if (!pixs)
|
|
return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
|
|
if (pixGetDepth(pixs) != 8)
|
|
return (PIX *)ERROR_PTR("pixs must be gray", procName, NULL);
|
|
pixGetDimensions(pixs, &w, &h, &d);
|
|
if (!spatial_kel)
|
|
return (PIX *)ERROR_PTR("spatial kel not defined", procName, NULL);
|
|
|
|
if (!range_kel)
|
|
return pixConvolve(pixs, spatial_kel, 8, 1);
|
|
if (range_kel->sx != 256 || range_kel->sy != 1)
|
|
return (PIX *)ERROR_PTR("range kel not {256 x 1", procName, NULL);
|
|
|
|
keli = kernelInvert(spatial_kel);
|
|
kernelGetParameters(keli, &sy, &sx, &cy, &cx);
|
|
if ((pixt = pixAddMirroredBorder(pixs, cx, sx - cx, cy, sy - cy)) == NULL) {
|
|
kernelDestroy(&keli);
|
|
return (PIX *)ERROR_PTR("pixt not made", procName, NULL);
|
|
}
|
|
|
|
pixd = pixCreate(w, h, 8);
|
|
datat = pixGetData(pixt);
|
|
datad = pixGetData(pixd);
|
|
wplt = pixGetWpl(pixt);
|
|
wpld = pixGetWpl(pixd);
|
|
for (i = 0, id = 0; id < h; i++, id++) {
|
|
lined = datad + id * wpld;
|
|
for (j = 0, jd = 0; jd < w; j++, jd++) {
|
|
center_val = GET_DATA_BYTE(datat + (i + cy) * wplt, j + cx);
|
|
weight_sum = 0.0;
|
|
sum = 0.0;
|
|
for (k = 0; k < sy; k++) {
|
|
linet = datat + (i + k) * wplt;
|
|
for (m = 0; m < sx; m++) {
|
|
val = GET_DATA_BYTE(linet, j + m);
|
|
weight = keli->data[k][m] *
|
|
range_kel->data[0][L_ABS(center_val - val)];
|
|
weight_sum += weight;
|
|
sum += val * weight;
|
|
}
|
|
}
|
|
SET_DATA_BYTE(lined, jd, (l_int32)(sum / weight_sum + 0.5));
|
|
}
|
|
}
|
|
|
|
kernelDestroy(&keli);
|
|
pixDestroy(&pixt);
|
|
return pixd;
|
|
}
|
|
|
|
|
|
/*!
|
|
* \brief pixBlockBilateralExact()
|
|
*
|
|
* \param[in] pixs 8 bpp gray or 32 bpp rgb
|
|
* \param[in] spatial_stdev must be > 0.0
|
|
* \param[in] range_stdev must be > 0.0
|
|
* \return pixd 8 bpp or 32 bpp bilateral filtered image
|
|
*
|
|
* <pre>
|
|
* Notes:
|
|
* (1) See pixBilateralExact(). This provides an interface using
|
|
* the standard deviations of the spatial and range filters.
|
|
* (2) The convolution window halfwidth is 2 * spatial_stdev,
|
|
* and the square filter size is 4 * spatial_stdev + 1.
|
|
* The kernel captures 95% of total energy. This is compensated
|
|
* by normalization.
|
|
* (3) The range_stdev is analogous to spatial_halfwidth in the
|
|
* grayscale domain [0...255], and determines how much damping of the
|
|
* smoothing operation is applied across edges. The larger this
|
|
* value is, the smaller the damping. The smaller the value, the
|
|
* more edge details are preserved. These approximations are useful
|
|
* for deciding the appropriate cutoff.
|
|
* kernel[1 * stdev] ~= 0.6 * kernel[0]
|
|
* kernel[2 * stdev] ~= 0.14 * kernel[0]
|
|
* kernel[3 * stdev] ~= 0.01 * kernel[0]
|
|
* If range_stdev is infinite there is no damping, and this
|
|
* becomes a conventional gaussian smoothing.
|
|
* This value does not affect the run time.
|
|
* (4) If range_stdev is negative or zero, the range kernel is
|
|
* ignored and this degenerates to a straight gaussian convolution.
|
|
* (5) This is very slow for large spatial filters. The time
|
|
* on a 3GHz pentium is roughly
|
|
* T = 1.2 * 10^-8 * (A * sh^2) sec
|
|
* where A = # of pixels, sh = spatial halfwidth of filter.
|
|
* </pre>
|
|
*/
|
|
PIX*
|
|
pixBlockBilateralExact(PIX *pixs,
|
|
l_float32 spatial_stdev,
|
|
l_float32 range_stdev)
|
|
{
|
|
l_int32 d, halfwidth;
|
|
L_KERNEL *spatial_kel, *range_kel;
|
|
PIX *pixd;
|
|
|
|
PROCNAME("pixBlockBilateralExact");
|
|
|
|
if (!pixs)
|
|
return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
|
|
d = pixGetDepth(pixs);
|
|
if (d != 8 && d != 32)
|
|
return (PIX *)ERROR_PTR("pixs not 8 or 32 bpp", procName, NULL);
|
|
if (pixGetColormap(pixs) != NULL)
|
|
return (PIX *)ERROR_PTR("pixs is cmapped", procName, NULL);
|
|
if (spatial_stdev <= 0.0)
|
|
return (PIX *)ERROR_PTR("invalid spatial stdev", procName, NULL);
|
|
if (range_stdev <= 0.0)
|
|
return (PIX *)ERROR_PTR("invalid range stdev", procName, NULL);
|
|
|
|
halfwidth = 2 * spatial_stdev;
|
|
spatial_kel = makeGaussianKernel(halfwidth, halfwidth, spatial_stdev, 1.0);
|
|
range_kel = makeRangeKernel(range_stdev);
|
|
pixd = pixBilateralExact(pixs, spatial_kel, range_kel);
|
|
kernelDestroy(&spatial_kel);
|
|
kernelDestroy(&range_kel);
|
|
return pixd;
|
|
}
|
|
|
|
|
|
/*----------------------------------------------------------------------*
|
|
* Kernel helper function *
|
|
*----------------------------------------------------------------------*/
|
|
/*!
|
|
* \brief makeRangeKernel()
|
|
*
|
|
* \param[in] range_stdev must be > 0.0
|
|
* \return kel, or NULL on error
|
|
*
|
|
* <pre>
|
|
* Notes:
|
|
* (1) Creates a one-sided Gaussian kernel with the given
|
|
* standard deviation. At grayscale difference of one stdev,
|
|
* the kernel falls to 0.6, and to 0.01 at three stdev.
|
|
* (2) A typical input number might be 20. Then pixels whose
|
|
* value differs by 60 from the center pixel have their
|
|
* weight in the convolution reduced by a factor of about 0.01.
|
|
* </pre>
|
|
*/
|
|
L_KERNEL *
|
|
makeRangeKernel(l_float32 range_stdev)
|
|
{
|
|
l_int32 x;
|
|
l_float32 val, denom;
|
|
L_KERNEL *kel;
|
|
|
|
PROCNAME("makeRangeKernel");
|
|
|
|
if (range_stdev <= 0.0)
|
|
return (L_KERNEL *)ERROR_PTR("invalid stdev <= 0", procName, NULL);
|
|
|
|
denom = 2. * range_stdev * range_stdev;
|
|
if ((kel = kernelCreate(1, 256)) == NULL)
|
|
return (L_KERNEL *)ERROR_PTR("kel not made", procName, NULL);
|
|
kernelSetOrigin(kel, 0, 0);
|
|
for (x = 0; x < 256; x++) {
|
|
val = expf(-(l_float32)(x * x) / denom);
|
|
kernelSetElement(kel, 0, x, val);
|
|
}
|
|
return kel;
|
|
}
|