graymorphlow.c

Go to the documentation of this file.

/*====================================================================*
 -  Copyright (C) 2001 Leptonica.  All rights reserved.
 -  This software is distributed in the hope that it will be
 -  useful, but with NO WARRANTY OF ANY KIND.
 -  No author or distributor accepts responsibility to anyone for the
 -  consequences of using this software, or for whether it serves any
 -  particular purpose or works at all, unless he or she says so in
 -  writing.  Everyone is granted permission to copy, modify and
 -  redistribute this source code, for commercial or non-commercial
 -  purposes, with the following restrictions: (1) the origin of this
 -  source code must not be misrepresented; (2) modified versions must
 -  be plainly marked as such; and (3) this notice may not be removed
 -  or altered from any source or modified source distribution.
 *====================================================================*/


/*
 *  graymorphlow.c
 *
 *      Low-level grayscale morphological operations
 *
 *            void     dilateGrayLow()
 *            void     erodeGrayLow()
 *
 *      
 *      We use the van Herk/Gil-Werman (vHGW) algorithm, [van Herk,
 *      Patt. Recog. Let. 13, pp. 517-521, 1992; Gil and Werman,
 *      IEEE Trans PAMI 15(5), pp. 504-507, 1993.] 
 *      This was the first grayscale morphology
 *      algorithm to compute dilation and erosion with
 *      complexity independent of the size of the structuring
 *      element.  It is simple and elegant, and surprising that
 *      it was discovered as recently as 1992.  It works for
 *      SEs composed of horizontal and/or vertical lines.  The
 *      general case requires finding the Min or Max over an
 *      arbitrary set of pixels, and this requires a number of
 *      pixel comparisons equal to the SE "size" at each pixel
 *      in the image.  The vHGW algorithm requires not
 *      more than 3 comparisons at each point.  The algorithm has been
 *      recently refined by Gil and Kimmel ("Efficient Dilation
 *      Erosion, Opening and Closing Algorithms", in "Mathematical
 *      Morphology and its Applications to Image and Signal Processing",
 *      the proceedings of the International Symposium on Mathematical
 *      Morphology, Palo Alto, CA, June 2000, Kluwer Academic
 *      Publishers, pp. 301-310).  They bring this number down below
 *      1.5 comparisons per output pixel but at a cost of significantly
 *      increased complexity, so I don't bother with that here.
 *
 *      In brief, the method is as follows.  We evaluate the dilation
 *      in groups of "size" pixels, equal to the size of the SE.
 *      For horizontal, we start at x = "size"/2 and go
 *      (w - 2 * ("size"/2))/"size" steps.  This means that
 *      we don't evaluate the first 0.5 * "size" pixels and, worst
 *      case, the last 1.5 * "size" pixels.  Thus we embed the
 *      image in a larger image with these augmented dimensions, where
 *      the new border pixels are appropriately initialized (0 for
 *      dilation; 255 for erosion), and remove the boundary at the end.
 *      (For vertical, use h instead of w.)   Then for each group
 *      of "size" pixels, we form an array of length 2 * "size" + 1,
 *      consisting of backward and forward partial maxima (for
 *      dilation) or minima (for erosion).  This represents a
 *      jumping window computed from the source image, over which
 *      the SE will slide.  The center of the array gets the source
 *      pixel at the center of the SE.  Call this the center pixel
 *      of the window.  Array values to left of center get
 *      the maxima(minima) of the pixels from the center
 *      one and going to the left an equal distance.  Array
 *      values to the right of center get the maxima(minima) to
 *      the pixels from the center one and going to the right
 *      an equal distance.  These are computed sequentially starting
 *      from the center one.  The SE (of length "size") can slide over this 
 *      window (of length 2 * "size + 1) at "size" different places.
 *      At each place, the maxima(minima) of the values in the window
 *      that correspond to the end points of the SE give the extremal
 *      values over that interval, and these are stored at the dest
 *      pixel corresponding to the SE center.  A picture is worth
 *      at least this many words, so if this isn't clear, see the
 *      leptonica documentation on grayscale morphology.
 *      
 */

#include <stdio.h>

#include "allheaders.h"



/*-----------------------------------------------------------------*
 *              Low-level gray morphological operations            *
 *-----------------------------------------------------------------*/
/*!
 *  dilateGrayLow()
 *
 *    Input: datad, w, h, wpld (8 bpp image)
 *           datas, wpls  (8 bpp image, of same dimensions)
 *           size  (full length of SEL; restricted to odd numbers)
 *           direction  (L_HORIZ or L_VERT)
 *           buffer  (holds full line or column of src image pixels)
 *           maxarray  (array of dimension 2*size+1)
 *    Return: void
 *           
 *    Note: To eliminate border effects on the actual image, these images
 *          are prepared with an additional border of dimensions:
 *             leftpix = 0.5 * size
 *             rightpix = 1.5 * size
 *             toppix = 0.5 * size
 *             bottompix = 1.5 * size
 *          and we initialize the src border pixels to 0.
 *          This allows full processing over the actual image; at
 *          the end the border is removed.
 *
 *    Method: Algorithm by van Herk and Gil and Werman
 */
void
dilateGrayLow(l_uint32  *datad,
              l_int32    w,
              l_int32    h,
              l_int32    wpld,
              l_uint32  *datas,
              l_int32    wpls,
              l_int32    size,
              l_int32    direction,
              l_uint8   *buffer,
              l_uint8   *maxarray)
{
l_int32    i, j, k;
l_int32    hsize, nsteps, startmax, startx, starty;
l_uint8    maxval;
l_uint32  *lines, *lined;

    if (direction == L_HORIZ) {
        hsize = size / 2;
        nsteps = (w - 2 * hsize) / size;
        for (i = 0; i < h; i++) {
            lines = datas + i * wpls;
            lined = datad + i * wpld;

                /* fill buffer with pixels in byte order */
            for (j = 0; j < w; j++)
                buffer[j] = GET_DATA_BYTE(lines, j);

            for (j = 0; j < nsteps; j++)
            {
                    /* refill the minarray */
                startmax = (j + 1) * size - 1;
                maxarray[size - 1] = buffer[startmax];
                for (k = 1; k < size; k++) {
                    maxarray[size - 1 - k] =
                        L_MAX(maxarray[size - k], buffer[startmax - k]);
                    maxarray[size - 1 + k] =
                        L_MAX(maxarray[size + k - 2], buffer[startmax + k]);
                }

                    /* compute dilation values */
                startx = hsize + j * size;
                SET_DATA_BYTE(lined, startx, maxarray[0]);
                SET_DATA_BYTE(lined, startx + size - 1, maxarray[2 * size - 2]);
                for (k = 1; k < size - 1; k++) {
                    maxval = L_MAX(maxarray[k], maxarray[k + size - 1]);
                    SET_DATA_BYTE(lined, startx + k, maxval);
                }
            }
        }
    }
    else {   /* direction == L_VERT */ 
        hsize = size / 2;
        nsteps = (h - 2 * hsize) / size;
        for (j = 0; j < w; j++) {

                /* fill buffer with pixels in byte order */
            for (i = 0; i < h; i++) {
                lines = datas + i * wpls;
                buffer[i] = GET_DATA_BYTE(lines, j);
            }

            for (i = 0; i < nsteps; i++)
            {
                    /* refill the minarray */
                startmax = (i + 1) * size - 1;
                maxarray[size - 1] = buffer[startmax];
                for (k = 1; k < size; k++) {
                    maxarray[size - 1 - k] =
                        L_MAX(maxarray[size - k], buffer[startmax - k]);
                    maxarray[size - 1 + k] =
                        L_MAX(maxarray[size + k - 2], buffer[startmax + k]);
                }

                    /* compute dilation values */
                starty = hsize + i * size;
                lined = datad + starty * wpld;
                SET_DATA_BYTE(lined, j, maxarray[0]);
                SET_DATA_BYTE(lined + (size - 1) * wpld, j,
                        maxarray[2 * size - 2]);
                for (k = 1; k < size - 1; k++) {
                    maxval = L_MAX(maxarray[k], maxarray[k + size - 1]);
                    SET_DATA_BYTE(lined + wpld * k, j, maxval);
                }
            }
        }
    }
            
    return;
}


/*!
 *  erodeGrayLow()
 *
 *    Input: datad, w, h, wpld (8 bpp image)
 *           datas, wpls  (8 bpp image, of same dimensions)
 *           size  (full length of SEL; restricted to odd numbers)
 *           direction  (L_HORIZ or L_VERT)
 *           buffer  (holds full line or column of src image pixels)
 *           minarray  (array of dimension 2*size+1)
 *    Return: void
 *           
 *    Note: To eliminate border effects on the actual image, these images
 *          are prepared with an additional border of dimensions:
 *             leftpix = 0.5 * size
 *             rightpix = 1.5 * size
 *             toppix = 0.5 * size
 *             bottompix = 1.5 * size
 *          and we initialize the src border pixels to 255.
 *          This allows full processing over the actual image; at
 *          the end the border is removed.
 *
 *    Method: Algorithm by van Herk and Gil and Werman
 *
 */
void
erodeGrayLow(l_uint32  *datad,
             l_int32    w,
             l_int32    h,
             l_int32    wpld,
             l_uint32  *datas,
             l_int32    wpls,
             l_int32    size,
             l_int32    direction,
             l_uint8   *buffer,
             l_uint8   *minarray)
{
l_int32    i, j, k;
l_int32    hsize, nsteps, startmin, startx, starty;
l_uint8    minval;
l_uint32  *lines, *lined;

    if (direction == L_HORIZ) {
        hsize = size / 2;
        nsteps = (w - 2 * hsize) / size;
        for (i = 0; i < h; i++) {
            lines = datas + i * wpls;
            lined = datad + i * wpld;

                /* fill buffer with pixels in byte order */
            for (j = 0; j < w; j++)
                buffer[j] = GET_DATA_BYTE(lines, j);

            for (j = 0; j < nsteps; j++)
            {
                    /* refill the minarray */
                startmin = (j + 1) * size - 1;
                minarray[size - 1] = buffer[startmin];
                for (k = 1; k < size; k++) {
                    minarray[size - 1 - k] =
                        L_MIN(minarray[size - k], buffer[startmin - k]);
                    minarray[size - 1 + k] =
                        L_MIN(minarray[size + k - 2], buffer[startmin + k]);
                }

                    /* compute erosion values */
                startx = hsize + j * size;
                SET_DATA_BYTE(lined, startx, minarray[0]);
                SET_DATA_BYTE(lined, startx + size - 1, minarray[2 * size - 2]);
                for (k = 1; k < size - 1; k++) {
                    minval = L_MIN(minarray[k], minarray[k + size - 1]);
                    SET_DATA_BYTE(lined, startx + k, minval);
                }
            }
        }
    }
    else {   /* direction == L_VERT */ 
        hsize = size / 2;
        nsteps = (h - 2 * hsize) / size;
        for (j = 0; j < w; j++) {

                /* fill buffer with pixels in byte order */
            for (i = 0; i < h; i++) {
                lines = datas + i * wpls;
                buffer[i] = GET_DATA_BYTE(lines, j);
            }

            for (i = 0; i < nsteps; i++)
            {
                    /* refill the minarray */
                startmin = (i + 1) * size - 1;
                minarray[size - 1] = buffer[startmin];
                for (k = 1; k < size; k++) {
                    minarray[size - 1 - k] =
                        L_MIN(minarray[size - k], buffer[startmin - k]);
                    minarray[size - 1 + k] =
                        L_MIN(minarray[size + k - 2], buffer[startmin + k]);
                }

                    /* compute erosion values */
                starty = hsize + i * size;
                lined = datad + starty * wpld;
                SET_DATA_BYTE(lined, j, minarray[0]);
                SET_DATA_BYTE(lined + (size - 1) * wpld, j,
                        minarray[2 * size - 2]);
                for (k = 1; k < size - 1; k++) {
                    minval = L_MIN(minarray[k], minarray[k + size - 1]);
                    SET_DATA_BYTE(lined + wpld * k, j, minval);
                }
            }
        }
    }
            
    return;
}