text: /*====================================================================*
   -  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.
   *====================================================================*/
  
  /*
   *  bilateral.c
   *
   *     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).
   */
  
  #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  *
   *--------------------------------------------------------------------------*/
  /*!
   *  pixBilateral()
   *
   *      Input:  pixs (8 bpp gray or 32 bpp rgb, no colormap)
   *              spatial_stdev  (of gaussian kernel; in pixels, > 0.5)
   *              range_stdev  (of gaussian range kernel; > 5.0; typ. 50.0)
   *              ncomps (number of intermediate sums J(k,x); in [4 ... 30])
   *              reduction  (1, 2 or 4)
   *      Return: pixd (bilateral filtered image), or null on error
   *
   *  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.
   */
  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;
  }
  
  
  /*!
   *  pixBilateralGray()
   *
   *      Input:  pixs (8 bpp gray)
   *              spatial_stdev  (of gaussian kernel; in pixels, > 0.5)
   *              range_stdev  (of gaussian range kernel; > 5.0; typ. 50.0)
   *              ncomps (number of intermediate sums J(k,x); in [4 ... 30])
   *              reduction  (1, 2 or 4)
   *      Return: pixd (8 bpp bilateral filtered image), or null on error
   *
   *  Notes:
   *      (1) See pixBilateral() for constraints on the input parameters.
   *      (2) See pixBilateral() for algorithm details.
   */
  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       *
   *----------------------------------------------------------------------*/
  /*!
   *  bilateralCreate()
   *
   *      Input:  pixs (8 bpp gray, no colormap)
   *              spatial_stdev  (of gaussian kernel; in pixels, > 0.5)
   *              range_stdev  (of gaussian range kernel; > 5.0; typ. 50.0)
   *              ncomps (number of intermediate sums J(k,x); in [4 ... 30])
   *              reduction  (1, 2 or 4)
   *      Return: bil, or null on error
   *
   *  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.
   */
  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 *)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 *)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 *)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 *)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 *)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 *)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 * sstdev + 0.5);
      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;
  }
  
  
  /*!
   *  bilateralApply()
   *
   *      Input:  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;
  }
  
  
  /*!
   *  bilateralDestroy()
   *
   *      Input:  &bil
   *      Return: void
   */
  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);
      FREE(bil->spatial);
      FREE(bil->range);
      FREE(bil->nc);
      FREE(bil->kindex);
      FREE(bil->kfract);
      for (i = 0; i < bil->ncomps; i++)
          FREE(bil->lineset[i]);
      FREE(bil->lineset);
      FREE(bil);
      *pbil = NULL;
      return;
  }
  
  
  /*----------------------------------------------------------------------*
   *    Exact implementation of grayscale or color bilateral filtering    *
   *----------------------------------------------------------------------*/
  /*!
   *  pixBilateralExact()
   *
   *      Input:  pixs (8 bpp gray or 32 bpp rgb)
   *              spatial_kel  (gaussian kernel)
   *              range_kel (<optional> 256 x 1, monotonically decreasing)
   *      Return: pixd (8 bpp bilateral filtered image)
   *
   *  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.
   */
  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;
      }
  }
  
  
  /*!
   *  pixBilateralGrayExact()
   *
   *      Input:  pixs (8 bpp gray)
   *              spatial_kel  (gaussian kernel)
   *              range_kel (<optional> 256 x 1, monotonically decreasing)
   *      Return: pixd (8 bpp bilateral filtered image)
   *
   *  Notes:
   *      (1) See pixBilateralExact().
   */
  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)
          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;
  }
  
  
  /*!
   *  pixBlockBilateralExact()
   *
   *      Input:  pixs (8 bpp gray or 32 bpp rgb)
   *              spatial_stdev (> 0.0)
   *              range_stdev (> 0.0)
   *      Return: pixd (8 bpp or 32 bpp bilateral filtered image)
   *
   *  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.
   */
  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                       *
   *----------------------------------------------------------------------*/
  /*!
   *  makeRangeKernel()
   *
   *      Input:  range_stdev (> 0)
   *      Return: kel, or null on error
   *
   *  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.
   */
  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;
  }
  