bilateral.c

Go to the documentation of this file.

/*====================================================================*
 -  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.
 *====================================================================*/

#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  *
 *--------------------------------------------------------------------------*/
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;
}


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       *
 *----------------------------------------------------------------------*/
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;
}


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;
}


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    *
 *----------------------------------------------------------------------*/
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;
    }
}


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;
}


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                       *
 *----------------------------------------------------------------------*/
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;
}