| /*====================================================================* |
| - 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. |
| *====================================================================*/ |
| |
| /* |
| * rank.c |
| * |
| * Rank filter (gray and rgb) |
| * PIX *pixRankFilter() |
| * PIX *pixRankFilterRGB() |
| * PIX *pixRankFilterGray() |
| * |
| * Median filter |
| * PIX *pixMedianFilter() |
| * |
| * 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! |
| */ |
| |
| #include <stdio.h> |
| #include <stdlib.h> |
| #include "allheaders.h" |
| |
| |
| /*----------------------------------------------------------------------* |
| * Rank order filter * |
| *----------------------------------------------------------------------*/ |
| /*! |
| * pixRankFilter() |
| * |
| * Input: pixs (8 or 32 bpp; no colormap) |
| * wf, hf (width and height of filter; each is >= 1) |
| * 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); |
| } |
| |
| |
| /*! |
| * pixRankFilterRGB() |
| * |
| * Input: pixs (32 bpp) |
| * wf, hf (width and height of filter; each is >= 1) |
| * 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; |
| } |
| |
| |
| /*! |
| * pixRankFilterGray() |
| * |
| * Input: pixs (8 bpp; no colormap) |
| * wf, hf (width and height of filter; each is >= 1) |
| * 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 *)CALLOC(256, sizeof(l_int32)); |
| histo16 = (l_int32 *)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; |
| } |
| } |
| 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; |
| } |
| } |
| 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); |
| FREE(histo); |
| FREE(histo16); |
| return pixd; |
| } |
| |
| |
| /*----------------------------------------------------------------------* |
| * Median filter * |
| *----------------------------------------------------------------------*/ |
| /*! |
| * pixMedianFilter() |
| * |
| * Input: pixs (8 or 32 bpp; no colormap) |
| * 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); |
| } |
| |
| |
