diff mupdf-source/thirdparty/leptonica/src/rank.c @ 2:b50eed0cc0ef upstream

ADD: MuPDF v1.26.7: the MuPDF source as downloaded by a default build of PyMuPDF 1.26.4. The directory name has changed: no version number in the expanded directory now.
author Franz Glasner <fzglas.hg@dom66.de>
date Mon, 15 Sep 2025 11:43:07 +0200
parents
children
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/mupdf-source/thirdparty/leptonica/src/rank.c	Mon Sep 15 11:43:07 2025 +0200
@@ -0,0 +1,534 @@
+/*====================================================================*
+ -  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
+ * <pre>
+ *
+ *      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.
+ * </pre>
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config_auto.h>
+#endif  /* HAVE_CONFIG_H */
+
+#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
+ *
+ * <pre>
+ * 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.
+ * </pre>
+ */
+PIX  *
+pixRankFilter(PIX       *pixs,
+              l_int32    wf,
+              l_int32    hf,
+              l_float32  rank)
+{
+l_int32  d;
+
+    if (!pixs)
+        return (PIX *)ERROR_PTR("pixs not defined", __func__, NULL);
+    if (pixGetColormap(pixs) != NULL)
+        return (PIX *)ERROR_PTR("pixs has colormap", __func__, NULL);
+    d = pixGetDepth(pixs);
+    if (d != 8 && d != 32)
+        return (PIX *)ERROR_PTR("pixs not 8 or 32 bpp", __func__, NULL);
+    if (wf < 1 || hf < 1)
+        return (PIX *)ERROR_PTR("wf < 1 || hf < 1", __func__, NULL);
+    if (rank < 0.0 || rank > 1.0)
+        return (PIX *)ERROR_PTR("rank must be in [0.0, 1.0]", __func__, 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
+ *
+ * <pre>
+ * 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.
+ * </pre>
+ */
+PIX  *
+pixRankFilterRGB(PIX       *pixs,
+                 l_int32    wf,
+                 l_int32    hf,
+                 l_float32  rank)
+{
+PIX  *pixr, *pixg, *pixb, *pixrf, *pixgf, *pixbf, *pixd;
+
+    if (!pixs)
+        return (PIX *)ERROR_PTR("pixs not defined", __func__, NULL);
+    if (pixGetDepth(pixs) != 32)
+        return (PIX *)ERROR_PTR("pixs not 32 bpp", __func__, NULL);
+    if (wf < 1 || hf < 1)
+        return (PIX *)ERROR_PTR("wf < 1 || hf < 1", __func__, NULL);
+    if (rank < 0.0 || rank > 1.0)
+        return (PIX *)ERROR_PTR("rank must be in [0.0, 1.0]", __func__, 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
+ *
+ * <pre>
+ * 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.
+ * </pre>
+ */
+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;
+
+    if (!pixs)
+        return (PIX *)ERROR_PTR("pixs not defined", __func__, NULL);
+    if (pixGetColormap(pixs) != NULL)
+        return (PIX *)ERROR_PTR("pixs has colormap", __func__, NULL);
+    pixGetDimensions(pixs, &w, &h, &d);
+    if (d != 8)
+        return (PIX *)ERROR_PTR("pixs not 8 bpp", __func__, NULL);
+    if (wf < 1 || hf < 1)
+        return (PIX *)ERROR_PTR("wf < 1 || hf < 1", __func__, NULL);
+    if (rank < 0.0 || rank > 1.0)
+        return (PIX *)ERROR_PTR("rank must be in [0.0, 1.0]", __func__, 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.0001f;
+    if (rank == 1.0) rank = 0.9999f;
+
+        /* 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", __func__, 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", __func__);
+                    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", __func__);
+                    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)
+{
+    if (!pixs)
+        return (PIX *)ERROR_PTR("pixs not defined", __func__, 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
+ *
+ * <pre>
+ * 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.
+ * </pre>
+ */
+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;
+
+    if (!pixs)
+        return (PIX *)ERROR_PTR("pixs not defined", __func__, NULL);
+    if (pixGetColormap(pixs) != NULL)
+        return (PIX *)ERROR_PTR("pixs has colormap", __func__, NULL);
+    d = pixGetDepth(pixs);
+    if (d != 8 && d != 32)
+        return (PIX *)ERROR_PTR("pixs not 8 or 32 bpp", __func__, NULL);
+    if (wf < 1 || hf < 1)
+        return (PIX *)ERROR_PTR("wf < 1 || hf < 1", __func__, NULL);
+    if (rank < 0.0 || rank > 1.0)
+        return (PIX *)ERROR_PTR("rank must be in [0.0, 1.0]", __func__, 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", __func__);
+        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;
+}