/*====================================================================* - 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 rank.c *
* * Rank filter (gray and rgb) * PIX *pixRankFilter() * PIX *pixRankFilterRGB() * PIX *pixRankFilterGray() * * Median filter * PIX *pixMedianFilter() * * Rank filter (accelerated with downscaling) * PIX *pixRankFilterWithScaling() * * What is a brick rank filter? * * A brick rank order filter evaluates, for every pixel in the image, * a rectangular set of n = wf x hf pixels in its neighborhood (where the * pixel in question is at the "center" of the rectangle and is * included in the evaluation). It determines the value of the * neighboring pixel that is the r-th smallest in the set, * where r is some integer between 1 and n. The input rank parameter * is a fraction between 0.0 and 1.0, where 0.0 represents the * smallest value (r = 1) and 1.0 represents the largest value (r = n). * A median filter is a rank filter where rank = 0.5. * * It is important to note that grayscale erosion is equivalent * to rank = 0.0, and grayscale dilation is equivalent to rank = 1.0. * These are much easier to calculate than the general rank value, * thanks to the van Herk/Gil-Werman algorithm: * http://www.leptonica.com/grayscale-morphology.html * so you should use pixErodeGray() and pixDilateGray() for * rank 0.0 and 1.0, rsp. See notes below in the function header. * * How is a rank filter implemented efficiently on an image? * * Sorting will not work. * * * The best sort algorithms are O(n*logn), where n is the number * of values to be sorted (the area of the filter). For large * filters this is an impractically large number. * * * Selection of the rank value is O(n). (To understand why it's not * O(n*logn), see Numerical Recipes in C, 2nd edition, 1992, p. 355ff). * This also still far too much computation for large filters. * * * Suppose we get clever. We really only need to do an incremental * selection or sorting, because, for example, moving the filter * down by one pixel causes one filter width of pixels to be added * and another to be removed. Can we do this incrementally in * an efficient way? Unfortunately, no. The sorted values will be * in an array. Even if the filter width is 1, we can expect to * have to move O(n) pixels, because insertion and deletion can happen * anywhere in the array. By comparison, heapsort is excellent for * incremental sorting, where the cost for insertion or deletion * is O(logn), because the array itself doesn't need to * be sorted into strictly increasing order. However, heapsort * only gives the max (or min) value, not the general rank value. * * This leaves histograms. * * * Represented as an array. The problem with an array of 256 * bins is that, in general, a significant fraction of the * entire histogram must be summed to find the rank value bin. * Suppose the filter size is 5x5. You spend most of your time * adding zeroes. Ouch! * * * Represented as a linked list. This would overcome the * summing-over-empty-bin problem, but you lose random access * for insertions and deletions. No way. * * * Two histogram solution. Maintain two histograms with * bin sizes of 1 and 16. Proceed from coarse to fine. * First locate the coarse bin for the given rank, of which * there are only 16. Then, in the 256 entry (fine) histogram, * you need look at a maximum of 16 bins. For each output * pixel, the average number of bins summed over, both in the * coarse and fine histograms, is thus 16. * * If someone has a better method, please let me know! * * The rank filtering operation is relatively expensive, compared to most * of the other imaging operations. The speed is only weakly dependent * on the size of the rank filter. On standard hardware, it runs at * about 10 Mpix/sec for a 50 x 50 filter, and 25 Mpix/sec for * a 5 x 5 filter. For applications where the rank filter can be * performed on a downscaled image, significant speedup can be * achieved because the time goes as the square of the scaling factor. * We provide an interface that handles the details, and only * requires the amount of downscaling to be input. **/ #include "allheaders.h" /*----------------------------------------------------------------------* * Rank order filter * *----------------------------------------------------------------------*/ /*! * \brief pixRankFilter() * * \param[in] pixs 8 or 32 bpp; no colormap * \param[in] wf, hf width and height of filter; each is >= 1 * \param[in] rank in [0.0 ... 1.0] * \return pixd of rank values, or NULL on error * *
* Notes: * (1) This defines, for each pixel in pixs, a neighborhood of * pixels given by a rectangle "centered" on the pixel. * This set of wf*hf pixels has a distribution of values. * For each component, if the values are sorted in increasing * order, we choose the component such that rank*(wf*hf-1) * pixels have a lower or equal value and * (1-rank)*(wf*hf-1) pixels have an equal or greater value. * (2) See notes in pixRankFilterGray() for further details. **/ PIX * pixRankFilter(PIX *pixs, l_int32 wf, l_int32 hf, l_float32 rank) { l_int32 d; PROCNAME("pixRankFilter"); if (!pixs) return (PIX *)ERROR_PTR("pixs not defined", procName, NULL); if (pixGetColormap(pixs) != NULL) return (PIX *)ERROR_PTR("pixs has colormap", procName, NULL); d = pixGetDepth(pixs); if (d != 8 && d != 32) return (PIX *)ERROR_PTR("pixs not 8 or 32 bpp", procName, NULL); if (wf < 1 || hf < 1) return (PIX *)ERROR_PTR("wf < 1 || hf < 1", procName, NULL); if (rank < 0.0 || rank > 1.0) return (PIX *)ERROR_PTR("rank must be in [0.0, 1.0]", procName, NULL); if (wf == 1 && hf == 1) /* no-op */ return pixCopy(NULL, pixs); if (d == 8) return pixRankFilterGray(pixs, wf, hf, rank); else /* d == 32 */ return pixRankFilterRGB(pixs, wf, hf, rank); } /*! * \brief pixRankFilterRGB() * * \param[in] pixs 32 bpp * \param[in] wf, hf width and height of filter; each is >= 1 * \param[in] rank in [0.0 ... 1.0] * \return pixd of rank values, or NULL on error * *
* Notes: * (1) This defines, for each pixel in pixs, a neighborhood of * pixels given by a rectangle "centered" on the pixel. * This set of wf*hf pixels has a distribution of values. * For each component, if the values are sorted in increasing * order, we choose the component such that rank*(wf*hf-1) * pixels have a lower or equal value and * (1-rank)*(wf*hf-1) pixels have an equal or greater value. * (2) Apply gray rank filtering to each component independently. * (3) See notes in pixRankFilterGray() for further details. **/ PIX * pixRankFilterRGB(PIX *pixs, l_int32 wf, l_int32 hf, l_float32 rank) { PIX *pixr, *pixg, *pixb, *pixrf, *pixgf, *pixbf, *pixd; PROCNAME("pixRankFilterRGB"); if (!pixs) return (PIX *)ERROR_PTR("pixs not defined", procName, NULL); if (pixGetDepth(pixs) != 32) return (PIX *)ERROR_PTR("pixs not 32 bpp", procName, NULL); if (wf < 1 || hf < 1) return (PIX *)ERROR_PTR("wf < 1 || hf < 1", procName, NULL); if (rank < 0.0 || rank > 1.0) return (PIX *)ERROR_PTR("rank must be in [0.0, 1.0]", procName, NULL); if (wf == 1 && hf == 1) /* no-op */ return pixCopy(NULL, pixs); pixr = pixGetRGBComponent(pixs, COLOR_RED); pixg = pixGetRGBComponent(pixs, COLOR_GREEN); pixb = pixGetRGBComponent(pixs, COLOR_BLUE); pixrf = pixRankFilterGray(pixr, wf, hf, rank); pixgf = pixRankFilterGray(pixg, wf, hf, rank); pixbf = pixRankFilterGray(pixb, wf, hf, rank); pixd = pixCreateRGBImage(pixrf, pixgf, pixbf); pixDestroy(&pixr); pixDestroy(&pixg); pixDestroy(&pixb); pixDestroy(&pixrf); pixDestroy(&pixgf); pixDestroy(&pixbf); return pixd; } /*! * \brief pixRankFilterGray() * * \param[in] pixs 8 bpp; no colormap * \param[in] wf, hf width and height of filter; each is >= 1 * \param[in] rank in [0.0 ... 1.0] * \return pixd of rank values, or NULL on error * *
* Notes: * (1) This defines, for each pixel in pixs, a neighborhood of * pixels given by a rectangle "centered" on the pixel. * This set of wf*hf pixels has a distribution of values, * and if they are sorted in increasing order, we choose * the pixel such that rank*(wf*hf-1) pixels have a lower * or equal value and (1-rank)*(wf*hf-1) pixels have an equal * or greater value. * (2) By this definition, the rank = 0.0 pixel has the lowest * value, and the rank = 1.0 pixel has the highest value. * (3) We add mirrored boundary pixels to avoid boundary effects, * and put the filter center at (0, 0). * (4) This dispatches to grayscale erosion or dilation if the * filter dimensions are odd and the rank is 0.0 or 1.0, rsp. * (5) Returns a copy if both wf and hf are 1. * (6) Uses row-major or column-major incremental updates to the * histograms depending on whether hf > wf or hv <= wf, rsp. **/ PIX * pixRankFilterGray(PIX *pixs, l_int32 wf, l_int32 hf, l_float32 rank) { l_int32 w, h, d, i, j, k, m, n, rankloc, wplt, wpld, val, sum; l_int32 *histo, *histo16; l_uint32 *datat, *linet, *datad, *lined; PIX *pixt, *pixd; PROCNAME("pixRankFilterGray"); if (!pixs) return (PIX *)ERROR_PTR("pixs not defined", procName, NULL); if (pixGetColormap(pixs) != NULL) return (PIX *)ERROR_PTR("pixs has colormap", procName, NULL); pixGetDimensions(pixs, &w, &h, &d); if (d != 8) return (PIX *)ERROR_PTR("pixs not 8 bpp", procName, NULL); if (wf < 1 || hf < 1) return (PIX *)ERROR_PTR("wf < 1 || hf < 1", procName, NULL); if (rank < 0.0 || rank > 1.0) return (PIX *)ERROR_PTR("rank must be in [0.0, 1.0]", procName, NULL); if (wf == 1 && hf == 1) /* no-op */ return pixCopy(NULL, pixs); /* For rank = 0.0, this is a grayscale erosion, and for rank = 1.0, * a dilation. Grayscale morphology operations are implemented * for filters of odd dimension, so we dispatch to grayscale * morphology if both wf and hf are odd. Otherwise, we * slightly adjust the rank (to get the correct behavior) and * use the slower rank filter here. */ if (wf % 2 && hf % 2) { if (rank == 0.0) return pixErodeGray(pixs, wf, hf); else if (rank == 1.0) return pixDilateGray(pixs, wf, hf); } if (rank == 0.0) rank = 0.0001; if (rank == 1.0) rank = 0.9999; /* Add wf/2 to each side, and hf/2 to top and bottom of the * image, mirroring for accuracy and to avoid special-casing * the boundary. */ if ((pixt = pixAddMirroredBorder(pixs, wf / 2, wf / 2, hf / 2, hf / 2)) == NULL) return (PIX *)ERROR_PTR("pixt not made", procName, NULL); /* Set up the two histogram arrays. */ histo = (l_int32 *)LEPT_CALLOC(256, sizeof(l_int32)); histo16 = (l_int32 *)LEPT_CALLOC(16, sizeof(l_int32)); rankloc = (l_int32)(rank * wf * hf); /* Place the filter center at (0, 0). This is just a * convenient location, because it allows us to perform * the rank filter over x:(0 ... w - 1) and y:(0 ... h - 1). */ pixd = pixCreateTemplate(pixs); datat = pixGetData(pixt); wplt = pixGetWpl(pixt); datad = pixGetData(pixd); wpld = pixGetWpl(pixd); /* If hf > wf, it's more efficient to use row-major scanning. * Otherwise, traverse the image in use column-major order. */ if (hf > wf) { for (j = 0; j < w; j++) { /* row-major */ /* Start each column with clean histogram arrays. */ for (n = 0; n < 256; n++) histo[n] = 0; for (n = 0; n < 16; n++) histo16[n] = 0; for (i = 0; i < h; i++) { /* fast scan on columns */ /* Update the histos for the new location */ lined = datad + i * wpld; if (i == 0) { /* do full histo */ for (k = 0; k < hf; k++) { linet = datat + (i + k) * wplt; for (m = 0; m < wf; m++) { val = GET_DATA_BYTE(linet, j + m); histo[val]++; histo16[val >> 4]++; } } } else { /* incremental update */ linet = datat + (i - 1) * wplt; for (m = 0; m < wf; m++) { /* remove top line */ val = GET_DATA_BYTE(linet, j + m); histo[val]--; histo16[val >> 4]--; } linet = datat + (i + hf - 1) * wplt; for (m = 0; m < wf; m++) { /* add bottom line */ val = GET_DATA_BYTE(linet, j + m); histo[val]++; histo16[val >> 4]++; } } /* Find the rank value */ sum = 0; for (n = 0; n < 16; n++) { /* search over coarse histo */ sum += histo16[n]; if (sum > rankloc) { sum -= histo16[n]; break; } } if (n == 16) { /* avoid accessing out of bounds */ L_WARNING("n = 16; reducing\n", procName); n = 15; sum -= histo16[n]; } k = 16 * n; /* starting value in fine histo */ for (m = 0; m < 16; m++) { sum += histo[k]; if (sum > rankloc) { SET_DATA_BYTE(lined, j, k); break; } k++; } } } } else { /* wf >= hf */ for (i = 0; i < h; i++) { /* column-major */ /* Start each row with clean histogram arrays. */ for (n = 0; n < 256; n++) histo[n] = 0; for (n = 0; n < 16; n++) histo16[n] = 0; lined = datad + i * wpld; for (j = 0; j < w; j++) { /* fast scan on rows */ /* Update the histos for the new location */ if (j == 0) { /* do full histo */ for (k = 0; k < hf; k++) { linet = datat + (i + k) * wplt; for (m = 0; m < wf; m++) { val = GET_DATA_BYTE(linet, j + m); histo[val]++; histo16[val >> 4]++; } } } else { /* incremental update at left and right sides */ for (k = 0; k < hf; k++) { linet = datat + (i + k) * wplt; val = GET_DATA_BYTE(linet, j - 1); histo[val]--; histo16[val >> 4]--; val = GET_DATA_BYTE(linet, j + wf - 1); histo[val]++; histo16[val >> 4]++; } } /* Find the rank value */ sum = 0; for (n = 0; n < 16; n++) { /* search over coarse histo */ sum += histo16[n]; if (sum > rankloc) { sum -= histo16[n]; break; } } if (n == 16) { /* avoid accessing out of bounds */ L_WARNING("n = 16; reducing\n", procName); n = 15; sum -= histo16[n]; } k = 16 * n; /* starting value in fine histo */ for (m = 0; m < 16; m++) { sum += histo[k]; if (sum > rankloc) { SET_DATA_BYTE(lined, j, k); break; } k++; } } } } pixDestroy(&pixt); LEPT_FREE(histo); LEPT_FREE(histo16); return pixd; } /*----------------------------------------------------------------------* * Median filter * *----------------------------------------------------------------------*/ /*! * \brief pixMedianFilter() * * \param[in] pixs 8 or 32 bpp; no colormap * \param[in] wf, hf width and height of filter; each is >= 1 * \return pixd of median values, or NULL on error */ PIX * pixMedianFilter(PIX *pixs, l_int32 wf, l_int32 hf) { PROCNAME("pixMedianFilter"); if (!pixs) return (PIX *)ERROR_PTR("pixs not defined", procName, NULL); return pixRankFilter(pixs, wf, hf, 0.5); } /*----------------------------------------------------------------------* * Rank filter (accelerated with downscaling) * *----------------------------------------------------------------------*/ /*! * \brief pixRankFilterWithScaling() * * \param[in] pixs 8 or 32 bpp; no colormap * \param[in] wf, hf width and height of filter; each is >= 1 * \param[in] rank in [0.0 ... 1.0] * \param[in] scalefactor scale factor; must be >= 0.2 and <= 0.7 * \return pixd of rank values, or NULL on error * *
* Notes: * (1) This is a convenience function that downscales, does * the rank filtering, and upscales. Because the down- * and up-scaling functions are very fast compared to * rank filtering, the time it takes is reduced from that * for the simple rank filtering operation by approximately * the square of the scaling factor. **/ PIX * pixRankFilterWithScaling(PIX *pixs, l_int32 wf, l_int32 hf, l_float32 rank, l_float32 scalefactor) { l_int32 w, h, d, wfs, hfs; PIX *pix1, *pix2, *pixd; PROCNAME("pixRankFilterWithScaling"); if (!pixs) return (PIX *)ERROR_PTR("pixs not defined", procName, NULL); if (pixGetColormap(pixs) != NULL) return (PIX *)ERROR_PTR("pixs has colormap", procName, NULL); d = pixGetDepth(pixs); if (d != 8 && d != 32) return (PIX *)ERROR_PTR("pixs not 8 or 32 bpp", procName, NULL); if (wf < 1 || hf < 1) return (PIX *)ERROR_PTR("wf < 1 || hf < 1", procName, NULL); if (rank < 0.0 || rank > 1.0) return (PIX *)ERROR_PTR("rank must be in [0.0, 1.0]", procName, NULL); if (wf == 1 && hf == 1) /* no-op */ return pixCopy(NULL, pixs); if (scalefactor < 0.2 || scalefactor > 0.7) { L_ERROR("invalid scale factor; no scaling used\n", procName); return pixRankFilter(pixs, wf, hf, rank); } pix1 = pixScaleAreaMap(pixs, scalefactor, scalefactor); wfs = L_MAX(1, (l_int32)(scalefactor * wf + 0.5)); hfs = L_MAX(1, (l_int32)(scalefactor * hf + 0.5)); pix2 = pixRankFilter(pix1, wfs, hfs, rank); pixGetDimensions(pixs, &w, &h, NULL); pixd = pixScaleToSize(pix2, w, h); pixDestroy(&pix1); pixDestroy(&pix2); return pixd; }