Fast algorithm for finding the minimum distances to a set

MJW · December 29, 2015

Some time ago, Red ochre and I discussed ideas for make his Edge Shader plugin run faster. While looking for a paper on a different subject, I stumbled upon a fascinating (and surprising) paper by Pedro F. Felzenszwalb and Daniel P. Huttenlocher on finding the minimum Euclidean distances to a set of points in O(n) time. I converted the C++ code provided by one of the authors into C# and used it in a new version on the Edge Shader.

Here is the plugin: EdgeShaderMJW.zip

Here is the plugin code.

Spoiler


        #region User Entered Code
        // Name: EdgeShaderMJW
        // Author: Red ochre (John Robbins) and MJW
        // Submenu: Advanced
        // URL: http://www.getpaint.net/redirect/plugins.html
        // Title:EdgeShader                       Red ochre Dec 2015
        
        #region UICode
        int Amount1 = 255; // [0,255] Opacity threshold
        int Amount2 = 50; // [0,100] Maximum distance
        double Amount3 = 50; // [0,100] Effect mix %
        bool Amount4 = false; // [0,1] Swap Colors
        #endregion
             
        int [,] distanceArray = null;

        int threshA;
        int maxdistance;
        double effectmix;
        bool swapcols;
             
        // Call from OnSetRenderInfo
        void InitializeRendering(Surface dst, Surface src)
        {         
            threshA = Amount1; // Opacity threshold
            maxdistance = Amount2; // Maximum distance
            effectmix = Amount3; // Effect mix %
            swapcols = Amount4;

            if (distanceArray == null)
            {
                distanceArray = new int[src.Width, src.Height];
            }
         
            // Fill in the starting distance values: 0.0 for transparent (inside), MaxValue for opaque (outside).
            for (int y = 0; y < src.Height; y++)
                for (int x = 0; x < src.Width; x++)
                    distanceArray[x, y] = (src[x, y].A < threshA) ? 0 : DistanceTransform.MaxValue;

            // Fill in the minimum squared distances.
            DistanceTransform.TransformColumns(distanceArray);
        }
        
        void Render(Surface dst, Surface src, Rectangle rect)
        {
            // For CodeLab processing.
            if (distanceArray == null)
            {
                dst.CopySurface(src, rect.Location, rect);// copy surface quicker than looping through
                return;
            }
            
            ColorBgra PC = (ColorBgra)EnvironmentParameters.PrimaryColor;
            ColorBgra SC = (ColorBgra)EnvironmentParameters.SecondaryColor;
            if (swapcols) { SC = (ColorBgra)EnvironmentParameters.PrimaryColor; PC = (ColorBgra)EnvironmentParameters.SecondaryColor; }
            int H = rect.Bottom - rect.Top;
            int W = rect.Right - rect.Left;
            ColorBgra cp;
            int B, G, R, A;
            double shadrat = 0, ishadrat;
            double emix = effectmix / 100; double iemix = 1 - emix;
            double distScale = 1.0 / (double)maxdistance;

            DistanceTransform.TransformRows(distanceArray, rect.Top, rect.Bottom);
            for (int y = rect.Top; y < rect.Bottom; y++)
            {
                if (IsCancelRequested) return;
                for (int x = rect.Left; x < rect.Right; x++)
                {        
                    cp = src[x, y];
                    if (cp.A < threshA)
                    {
                        dst[x, y] = cp;
                    }
                    else
                    {
                        shadrat = distScale * Math.Sqrt((double)distanceArray[x, y]);
                        ishadrat = 1 - shadrat;
        
                        double Bnew = (shadrat * SC. + (ishadrat * PC.;
                        double Gnew = (shadrat * SC.G) + (ishadrat * PC.G);
                        double Rnew = (shadrat * SC.R) + (ishadrat * PC.R);
        
                        B = (int)((iemix * cp. + (emix * Bnew));
                        G = (int)((iemix * cp.G) + (emix * Gnew));
                        R = (int)((iemix * cp.R) + (emix * Rnew));
                        A = cp.A;
        
                        cp = ColorBgra.FromBgra(Int32Util.ClampToByte(,
                                                Int32Util.ClampToByte(G),
                                                Int32Util.ClampToByte(R), Int32Util.ClampToByte(A));
                        dst[x, y] = cp;
                    }
                }// end of x loop
            }//end of y loop
        }//end render

        #endregion

Here is the distance transform code:

Spoiler


namespace DistanceTransformation
{
    // DistanceTransform by MJW.
    // Based on an algorithm by A. Meijster, et al., and later by
    // Pedro F. Felzenszwalb and Daniel P. Huttenlocher.
    // This version combines elements from both papers along with some modifications.
    //
    // Suppose we have an image consisting of two disjoint subsets: "set" pixels and "not set"
    // pixels. For each not-set pixel, find the nearest set pixel, or the squared Euclidean
    // distance to the nearest set pixel.
    //
    // This problem can be solved using an O(n) algorithm which relies on the following facts:
    //
    // First, if (X, Y) is a pixel, and (Xs, Ys) is the nearest set pixel, there is no
    // pixel, (Xs, Ys'), such that |Ys' - Y| < |Ys - Y|; in other words, for the pixels within
    // a given row, it is only necessary to check the set of pixels which, for some X, minimizes
    // the Y distance. This set of pixels -- one for each X -- will be referred to as the
    // "candidate" pixels. (For a given X, there may be two pixels which are the same mininmal Y
    // distance. Since these two pixels are the same distance from every pixel in the row, either
    // can be used.)
    //
    // Second, the set of pixels in a row that are nearest to a given candidate pixel is
    // contiguous. This set will be referred to as the candidate's "span."
    //
    // Third, on any given row, the spans for two differnt candidate pixels occur in the same
    // X order as the pixels; that is, if S the span for P and S' is the span for P', and P is
    // left of P', then S is left of S'. (A pixel may be the same distance from two candidates;
    // in that case, it can either be the rightmost pixel of the left candidate's span, or the
    // leftmost pixel of the right candidate's span.)
    //
    // The first fact follows from observing that if (Xs, Ys) is the nearest pixel to (X, Y),
    // there can be no pixel (X, Ys') where |Ys' - Y| < |Ys - Y|, since such a pixel would be
    // nearer to (X, Y).
    //
    // The second and third facts follow from observing that if (Xs, Ys) and (Xs', Ys') are two
    // candidates, Xs < Xs', then there is some (real-valued) point on the Y row which is
    // equidistance from the two points. All the pixels on the row to the left are nearer to
    // (Xs, Ys) and all the pixels to the right are nearer to (Xs', Ys'). The spans for each
    // candidate must lie on the respective side of the equidistance point.
    class DistanceTransform
    {
        public static void Transform(int[,] image)
        {
            TransformColumns(image);
            TransformRows(image);
        }

        public const int MaxValue = 32767;
        const int MaxValueSquared = MaxValue * MaxValue;

        public static void TransformColumns(int[,] image)
        {
            TransformColumns(image, 0, image.GetLength(0));
        }

        public static void TransformColumns(int[,] image, int left, int right)
        {
            int height = image.GetLength(1);

            // The array elements are initially 0 if set and nonzero if not set.
            // Replace each element with the distance to the nearest element in the same column.
            for (int x = left; x < right; x++)
            {
                int lastSet = -MaxValue;
                for (int y = 0; y < height; y++)
                {
                    if (image[x, y] == 0)
                    {
                        lastSet = y;

                        // Check previous pixels to see if this Y is nearer.
                        int by = y;
                        while ((--by >= 0) && (image[x, by] > lastSet - by))
                            image[x, by] = lastSet - by;
                    }
                    else
                    {
                        image[x, y] = y - lastSet;
                    }
                }
            }
        }

        public static void TransformRows(int[,] image)
        {
            TransformRows(image, 0, image.GetLength(1));
        }

        struct SpanDistSq
        {
            public int NearestPixelX;
            public int SpanX;
            public int Distance2;

            public SpanDistSq(int nearestPixelX, int spanX, int distance2)
            {
                NearestPixelX = nearestPixelX;
                SpanX = spanX;
                Distance2 = distance2;
            }
        }

        public static void TransformRows(int[,] image, int top, int bottom)
        {
            int width = image.GetLength(0);
            SpanDistSq[] span = new SpanDistSq[width + 1];

            // Each entry in the image contains the distance to the nearest set pixel in the same
            // column. These are the candidate pixels.
            for (int y = top; y < bottom; y++)
            {
                // Taking each candidate pixel one at a time from left to right, determine its span. Each new
                // candidate pixel will be nearest to the pixels in the row from some point onward. Since the
                // spans are in the same X order as the candidate pixels, the previous spans are checked from
                // right to left. If the new candidate is closer than a previous candidate for that pixel's
                // entire span, the previous candidate's span is eliminated. Once a candidate is found that's
                // closer at the beginning of its span than the new candidate, the point at which the new
                // candidate becomes closer is appended to the span list as the end of the previous span and
                // the beginning of the new span.
                int spanIndex = 0;
                span[0] = new SpanDistSq(0, 0, Sq(image[0, y]));
                for (int x = 1; x < width; x++)
                {
                    // Calculate the distance squared from (0, Y) for each of the candidate pixels.
                    // If (X, Y) is a pixel in the Y row, and (Xs, Ys) is a set pixel:
                    // (X - Xs)^2 + (Y - Ys)^2
                    // X^2 - 2 * Xs * X + Xs^2 + (Y - Ys)^2
                    // [Xs^2 + (Y - Ys)^2] + X^2 - 2 * Xs * X
                    //
                    // If (Xs, Ys) and (Xs', Ys'), Xs < Xs', are two candidates for nearest pixels,
                    // (X, Y) is nearer to (Xs', Ys') then to (Xs, Ys) if,
                    // (X - Xs)^2 + (Y - Ys)^2 > (X - Xs')^2 + (Y - Ys')^2
                    // [Xs^2 + (Y - Ys)^2] + X^2 - 2 * Xs * X > [Xs'^2 + (Y - Ys')^2] + X^2 - 2 * Xs' * X
                    // [Xs'^2 + (Y - Ys')^2] - [Xs^2 + (Y - Ys)^2] < X * [2 * (Xs' - Xs)]
                    // X > {[Xs'^2 + (Y - Ys')^2] - [Xs^2 + (Y - Ys)^2]} / [2 * (Xs' - Xs)]
                    int dist2 = Sq(x) + Sq(image[x, y]);

                    // If the distance to beyond the maximum image distance, don't bother checking further,
                    // since it can't produce a valid span. (This test isn't necessary to the algorithm.)
                    if (dist2 >= MaxValueSquared)
                        continue;

                    // Loop until the previous span can't be eliminated or there are no previous spans.
                    while (true)
                    {
                        // Compare the distance of the new candidate at the beginning of the previous span
                        // to the distance of the span's associated candidate.
                        SpanDistSq prevSpan = span[spanIndex];
                        int prevPixelX = prevSpan.NearestPixelX;
                        int dist2Diff = dist2 - prevSpan.Distance2;
                        int twiceXDiff = 2 * (x - prevPixelX);
                        if (dist2Diff > twiceXDiff * prevSpan.SpanX)
                        {
                            // The span for the new candidate starts after the previous span.
                            // Add a new span to the list.
                            int spanX = (dist2Diff + twiceXDiff - 1) / twiceXDiff; // Ceiling integer division.
                            if (spanX < width)
                                span[++spanIndex] = new SpanDistSq(x, spanX, dist2);
                            break;
                        }
                        else if (--spanIndex < 0)
                        {
                            // No more previous spans, so start over.
                            spanIndex = 0;
                            span[0] = new SpanDistSq(x, 0, dist2);
                            break;
                        }
                    }
                }

                // The spans have been determined. Move through the list of spans, filling in the squared distances
                // for each pixel. The spans contain the X address of the nearest pixel and the distance squared of
                // its associated candidate from (0, Y).
                // (X - Xs)^2 + (Y - Ys)^2 =
                // X^2 - 2 * Xs * X + Xs^2 + (Y - Ys)^2 =
                // [Xs^2 + (Y - Ys)^2] + X^2 - 2 * Xs * X =
                // [Xs^2 + (Y - Ys)^2] + X * (X - 2 * Xs)
                span[spanIndex + 1].SpanX = width; // Avoid extra test for end-of-list.
                int endX = 0;
                for (int i = 0; i <= spanIndex; i++)
                {
                    int startX = endX;
                    endX = span[i + 1].SpanX;
                    int twicePixelX = 2 * span[i].NearestPixelX;
                    int pixelXDist2 = span[i].Distance2;
                    for (int x = startX; x < endX; x++)
                        image[x, y] = pixelXDist2 + x * (x - twicePixelX);
                }
            }
        }

        //---------------------------------------------------------------------------------------------
        // The following routines determine the nearest pixel rather than the distance to the nearest
        // pixel. Each entry in the 2D array is a point (using shorts to reduce memory use).
        // Entries within the set should be initialized to their own coordinates (since they are their
        // own closest points). Non-set entries can be initalized to anything else.
        public struct PointS
        {
            public short X, Y;
            public PointS(short x, short y)
            {
                X = x; Y = y;
            }

            public PointS(int x, int y)
            {
                X = (short)x; Y = (short)y;
            }
        }

        public static PointS MaxPointS = new PointS(32767, 32767);

        public static bool Transform(PointS[,] image)
        {
            TransformColumns(image);
            return TransformRows(image);
        }


        public static void TransformColumns(PointS[,] image)
        {
            TransformColumns(image, 0, image.GetLength(0));
        }

        public static void TransformColumns(PointS[,] image, int left, int right)
        {
            int height = image.GetLength(1);

            // The array pixels are initially set to their own coordinates if set, and other
            // coordinates if not set.
            // Replace each element with the distance to the nearest element in the same column
            // that's within the set.
            for (int x = left; x < right; x++)
            {
                // Set the elements to the value to the nearest set-member in the same column.
                short nearestY = -MaxValue;
                for (int y = 0; y < height; y++)
                {
                    if ((image[x, y].X == x) && (image[x, y].Y == y))
                    {
                        nearestY = (short)y;

                        // Check previous pixels to see if this Y is nearer.
                        int by = y;
                        while ((--by >= 0) && (nearestY - by < by - image[x, by].Y))
                            image[x, by].Y = nearestY;
                    }
                    else
                    {
                        image[x, y] = new PointS(x, nearestY);
                    }
                }
            }
        }

        struct SpanXY
        {
            public PointS NearestPixel;    // Nearest pixel associated with span.
            public int SpanX;               // Beginning of span.

            public SpanXY(PointS nearestPixel, int spanX)
            {
                NearestPixel = nearestPixel;
                SpanX = spanX;
            }
        }

        public static bool TransformRows(PointS[,] image)
        {
            return TransformRows(image, 0, image.GetLength(1));
        }

        public static bool TransformRows(PointS[,] image, int top, int bottom)
        {
            int width = image.GetLength(0);
            SpanXY[] span = new SpanXY[width + 1];

            // Each entry in the image contains the nearest set pixel in the same column. These are the
            // candidate pixels.
            span[0] = new SpanXY(image[0, top], 0);
            for (int y = top; y < bottom; y++)
            {
                // Taking each candidate pixel one at a time from left to right, determine its span. Each new
                // candidate pixel will be nearest to the pixels in the row from some point onward. Since the
                // spans are in the same X order as the candidate pixels, the previous spans are checked from
                // right to left. If the new candidate is closer than a previous candidate for that pixel's
                // entire span, the previous candidate's span is eliminated. Once a candidate is found that's
                // closer at the beginning of its span than the new candidate, the point at which the new
                // candidate becomes closer is appended to the span list as the end of the previous span and
                // the beginning of the new span.
                int spanIndex = 0;
                span[0] = new SpanXY(image[0, y], 0);
                for (int x = 1; x < width; x++)
                {
                    int dist2 = Sq(x) + Sq(image[x, y].Y - y);

                    // If the distance to beyond the maximum image distance, don't bother checking further,
                    // since it can't produce a valid span. (This test isn't necessary to the algorithm.)
                    if (dist2 >= MaxValueSquared)
                        continue;

                    // Loop until the previous span can't be eliminated or there are no previous spans.
                    while (true)
                    {
                        // Compare the distance of the new candidate at the beginning of the previous span
                        // to the distance of its associated candidate.
                        SpanXY prevSpan = span[spanIndex];
                        int prevPixelX = prevSpan.NearestPixel.X;
                        int dist2Diff = dist2 - Sq(prevPixelX) - Sq(y - prevSpan.NearestPixel.Y);
                        int twiceXDiff = 2 * (x - prevPixelX);
                        if (dist2Diff > twiceXDiff * prevSpan.SpanX)
                        {
                            // The span for the new candidate starts after the previous span.
                            // Add a new span to the list.
                            int spanX = (dist2Diff + twiceXDiff - 1) / twiceXDiff; // Ceiling integer division.
                            if (spanX < width)
                                span[++spanIndex] = new SpanXY(image[x, y], spanX);
                            break;
                        }
                        else if (--spanIndex < 0)
                        {
                            // No more previous spans, so start over.
                            spanIndex = 0;
                            span[0] = new SpanXY(image[x, y], 0);
                            break;
                        }
                    }
                }

                // The spans have been determined. Move through the list of spans, filling in nearest pixel for each pixel.
                span[spanIndex + 1].SpanX = width; // Avoid extra test for end-of-list.
                int endX = 0;
                for (int i = 0; i <= spanIndex; i++)
                {
                    int startX = endX;
                    endX = span[i + 1].SpanX;
                    PointS nearestPixel = span[i].NearestPixel;
                    for (int x = startX; x < endX; x++)
                        image[x, y] = nearestPixel;
                }
            }
           
            // Return a flag which is false if there were no set pixels in the original image.
            // If there weren't, the nearest-pixel coordinates will be invalid and outside the image range.
            return (image[0, top].Y >= 0);
        }

        static int Sq(int x)
        {
            return x * x;
        }
    }
}

EDIT: I found that the same idea was published earlier by A. Meijster, et al. in a paper called "A General Algorithm for Computing Distance Transforms in Linear Time."

EDIT 2: Because the algorithm processes the columns and then the rows independently, the process could be parallelized to some extent in PDN by performing the column processing in OnSetRenderInfo and the row processing in Render. It shouldn't be too difficult, but for right now, I'll leave that as an exercise for the reader.

EDIT 3: I improved (I believe) the code by combining elements from Meijster's version with the original. I made the routines static. I also processed the rows in the Render code for parallelization.

EDIT4: Put the distance-transform class in a separate file. Rewrote much of the code, using integer arithmetic. Added comments to explain how it works. Added routines to find the nearest points' coordinates rather than the distances. These routines are probably slightly slower than the distance version, but more general, since the distance can be easily calculated from the nearest pixel's coordinates. Added routines to transform a range of columns to aid in parallelization of the column as well as the row transformations (I don't know how practical that would be within PDN, but it could be used in other situations.)

EDIT 5: Replaced the plugin zip file with the current version. It does the same thing as the previous version; it just uses the new routines.

BoltBait · December 29, 2015

Example images?

MJW · December 29, 2015

There are several on the linked-to thread. The result is the same, the method is improved. Here is another example:

Eli · December 29, 2015

Thanks MJW, it works great.

MJW · December 29, 2015

That's quite flashy looking, Eli! I especially like the magenta glow.

(Credit for the Edge Shader goes to Red ochre, of course. I hope by appending the MJW to the end of the name on the posted version, it doesn't appear I'm trying steal credit for it. I just didn't want it to overwrite other versions.)

Red ochre · December 29, 2015

(Credit for the Edge Shader goes to Red ochre, of course. I hope by appending the MJW to the end of the name on the posted version, it doesn't appear I'm trying steal credit for it. I just didn't want it to overwrite other versions.)

Absoloutely not a problem MJW! - your help with this one and the sub pixel sampling method (smoothing) for clipdisplace are very much appreciated. :beer:

Seerose · December 29, 2015

MJW and Red ochre! Thank you so much.

Rick Brewster · December 29, 2015

const float Infinity = float.MaxValue;
const float NegativeInfinity = float.MinValue;

You could also use float.PositiveInfinity and float.NegativeInfinity.

MJW · December 30, 2015

I rewrote some of code using ideas from Meijster's paper. His method of transforming the columns is less general, but more efficient. I'm not sure the generality was of much use. I'd like to do everything "in place" without the input array. Meijster does that, but I couldn't get his version to work, and I don't yet understand the algorithm well enough to know which elements I can safely write without overwriting something that will be used later. Meijster's version uses ints instead of floats (which might be desirable). That was, I think, part of the problem with making it work as is.

I also parallelized it somewhat. The rows are now processed in the Render calls.

One thing that's interesting about the algorithm is that I'm pretty sure it could be modified to record for each pixel the coordinates of nearest pixel in the set, That might be useful for something or other.

Edited December 30, 2015 by MJW

MJW · January 7, 2016

I've updated the distance transform code, using integer, instead of floating point, arithmetic. I made various other changes, as described in the EDIT 4 comment.

MJW · January 7, 2016

For some reason, I can no longer edit the initial comment. Selecting Edit brings up an empty edit box, and none of the buttons -- even Cancel -- work. Maybe the forum has decided five edits is plenty. In any case:

EDIT 6: Renamed Point16 to PointS and SpanDist2 to SpanDistSq. I just decided I didn't like the original names.

Edited January 7, 2016 by MJW

Rick Brewster · January 8, 2016

If you want to follow the naming conventions I use in Paint.NET itself, you could call that PointInt16.

Int16 because it's comprised of a pair of System.Int16's. "short" is the C# name, "Int16" is the .NET and thus language-neutral name.

MJW · January 8, 2016

That sounds like a good name for it. I changed my version to use "PointInt16," but I can't seem to update my original comment to change it there.

MJW · January 20, 2016

After a little experimentation, I realized something I should have already known, which makes parallelizing the row transform more difficult, and perhaps not worth the effort. I'd assumed each row would only be processed once in each render pass. I didn't consider the case of a non-convex selection, which could break rows into multiple segments. That would result in the rows being transformed multiple times. One possible solution is to only do the parallelized row transform if the selection is a rectangle, which is the common case.

Edited January 20, 2016 by MJW

MJW · January 21, 2016

I've revised the distance transform code. I can't edit my initial comment, so I'll post it here.

Spoiler

Hidden Content:


namespace DistanceTransformation
{
    // DistanceTransform by MJW.
    // Based on an algorithm by A. Meijster, et al., and later by
    // Pedro F. Felzenszwalb and Daniel P. Huttenlocher.
    // This version combines elements from both papers along with some modifications.
    // 
    // Suppose we have an image consisting of two disjoint subsets: "set" pixels and "not set"
    // pixels. For each not-set pixel, find the nearest set pixel, or the squared Euclidean
    // distance to the nearest set pixel.
    //
    // This problem can be solved using an O(n) algorithm which relies on the following facts:
    //
    // First, if (X, Y) is a pixel, and (Xs, Ys) is the nearest set pixel, there is no
    // pixel, (Xs, Ys'), such that |Ys' - Y| < |Ys - Y|; in other words, for the pixels within
    // a given row, it is only necessary to check the set of pixels which, for some X, minimizes
    // the Y distance. This set of pixels -- one for each X -- will be referred to as the
    // "candidate" pixels. (For a given X, there may be two pixels which are the same mininmal Y
    // distance. Since these two pixels are the same distance from every pixel in the row, either
    // can be used.)
    //
    // Second, the set of pixels in a row that are nearest to a given candidate pixel is
    // contiguous. This set will be referred to as the candidate's "span."
    //
    // Third, on any given row, the spans for two differnt candidate pixels occur in the same
    // X order as the pixels; that is, if S the span for P and S' is the span for P', and P is
    // left of P', then S is left of S'. (A pixel may be the same distance from two candidates;
    // in that case, it can either be the rightmost pixel of the left candidate's span, or the
    // leftmost pixel of the right candidate's span.)
    //
    // The first fact follows from observing that if (Xs, Ys) is the nearest pixel to (X, Y),
    // there can be no pixel (X, Ys') where |Ys' - Y| < |Ys - Y|, since such a pixel would be
    // nearer to (X, Y).
    //
    // The second and third facts follow from observing that if (Xs, Ys) and (Xs', Ys') are two
    // candidates, Xs < Xs', then there is some (real-valued) point on the Y row which is
    // equidistance from the two points. All the pixels on the row to the left are nearer to
    // (Xs, Ys) and all the pixels to the right are nearer to (Xs', Ys'). The spans for each
    // candidate must lie on the respective side of the equidistance point. 
    class DistanceTransform : Table2D<int>
    {
        public DistanceTransform() : base()
        {
        }

        public DistanceTransform(int width, int height) : base(width, height)
        {
        }

        public DistanceTransform(int left, int top, int width, int height) :
            base(left, top, width, height)
        {
        }

        public const int MaxValue = 32767;
        const int MaxValueSquared = MaxValue * MaxValue;

        public void Transform()
        {
            TransformColumns();
            TransformRows();
        }

        public void Transform(IsInSetCallback isInSet)
        {
            Include(isInSet);
            Transform();
        }

        public void TransformColumns()
        {
            TransformColumns(Left, Right);
        }

        public void TransformColumns(IsInSetCallback isInSet)
        {
            Include(isInSet);
            TransformColumns();
        }

        public DistanceTransform(int left, int top, int width, int height, IsInSetCallback isInSet) :
            base(left, top, width, height)
        {
            Transform(isInSet);
        }

        public DistanceTransform(int width, int height, IsInSetCallback isInSet) :
            base(width, height)
        {
            Transform(isInSet);
        }

        public void TransformColumns(int left, int right)
        {
            // The array elements are initially 0 if set and nonzero if not set.
            // Replace each element with the distance to the nearest element in the same column.
            for (int x = left; x < right; x++)
            {
                int nearestY = -MaxValue;
                for (int y = Top; y < Bottom; y++)
                {
                    if (this[x, y] == 0)
                    {
                        nearestY = y;

                        // Check previous pixels to see if this Y is nearer.
                        int by = y;
                        while ((--by >= Top) && (this[x, by] > nearestY - by))
                            this[x, by] = nearestY - by;
                    }
                    else
                    {
                        this[x, y] = y - nearestY;
                    }
                }
            }
        }

        struct SpanDistSq
        {
            public int NearestPixelX;
            public int SpanX;
            public int Distance2;

            public SpanDistSq(int nearestPixelX, int spanX, int distance2)
            {
                NearestPixelX = nearestPixelX;
                SpanX = spanX;
                Distance2 = distance2;
            }
        }

        public void TransformRows()
        {
            TransformRows(Top, Bottom);
        }

        public void TransformRows(int top, int bottom)
        {
            SpanDistSq[] span = new SpanDistSq[Width + 1];

            // Each entry in the image contains the distance to the nearest set pixel in the same
            // column. These are the candidate pixels.
            for (int y = top; y < bottom; y++)
            {
                // Taking each candidate pixel one at a time from left to right, determine its span. Each new
                // candidate pixel will be nearest to the pixels in the row from some point onward. Since the
                // spans are in the same X order as the candidate pixels, the previous spans are checked from
                // right to left. If the new candidate is closer than a previous candidate for that pixel's
                // entire span, the previous candidate's span is eliminated. Once a candidate is found that's
                // closer at the beginning of its span than the new candidate, the point at which the new
                // candidate becomes closer is appended to the span list as the end of the previous span and
                // the beginning of the new span.
                int spanIndex = 0;
                span[0] = new SpanDistSq(Left, Left, Sq(Left) + Sq(this[Left, y]));
                for (int x = Left + 1; x < Right; x++)
                {
                    // Calculate the distance squared from (0, Y) for each of the candidate pixels.
                    // If (X, Y) is a pixel in the Y row, and (Xs, Ys) is a set pixel:
                    // (X - Xs)^2 + (Y - Ys)^2
                    // X^2 - 2 * Xs * X + Xs^2 + (Y - Ys)^2
                    // [Xs^2 + (Y - Ys)^2] + X^2 - 2 * Xs * X
                    //
                    // If (Xs, Ys) and (Xs', Ys'), Xs < Xs', are two candidates for nearest pixels,
                    // (X, Y) is nearer to (Xs', Ys') then to (Xs, Ys) if,
                    // (X - Xs)^2 + (Y - Ys)^2 > (X - Xs')^2 + (Y - Ys')^2
                    // [Xs^2 + (Y - Ys)^2] + X^2 - 2 * Xs * X > [Xs'^2 + (Y - Ys')^2] + X^2 - 2 * Xs' * X
                    // [Xs'^2 + (Y - Ys')^2] - [Xs^2 + (Y - Ys)^2] < X * [2 * (Xs' - Xs)]
                    // X > {[Xs'^2 + (Y - Ys')^2] - [Xs^2 + (Y - Ys)^2]} / [2 * (Xs' - Xs)]
                    int dist2 = Sq(x) + Sq(this[x, y]);

                    // If the distance to beyond the maximum image distance, don't bother checking further,
                    // since it can't produce a valid span. (This test isn't necessary to the algorithm.)
                    if (dist2 >= MaxValueSquared)
                        continue;

                    // Loop until the previous span can't be eliminated or there are no previous spans.
                    while (true)
                    {
                        // Compare the distance of the new candidate at the beginning of the previous span
                        // to the distance of the span's associated candidate.
                        SpanDistSq prevSpan = span[spanIndex];
                        int prevPixelX = prevSpan.NearestPixelX;
                        int dist2Diff = dist2 - prevSpan.Distance2;
                        int twiceXDiff = 2 * (x - prevPixelX);
                        if (dist2Diff > twiceXDiff * prevSpan.SpanX)
                        {
                            // The span for the new candidate starts after the previous span.
                            // Add a new span to the list.
                            int spanX = (dist2Diff + twiceXDiff - 1) / twiceXDiff; // Ceiling integer division.
                            if (spanX < Right)
                                span[++spanIndex] = new SpanDistSq(x, spanX, dist2);
                            break;
                        }
                        else if (--spanIndex < 0)
                        {
                            // No more previous spans, so start over.
                            spanIndex = 0;
                            span[0] = new SpanDistSq(x, Left, dist2);
                            break;
                        }
                    }
                }

                // The spans have been determined. Move through the list of spans, filling in the squared distances
                // for each pixel. The spans contain the X address of the nearest pixel and the distance squared of
                // its associated candidate from (0, Y).
                // (X - Xs)^2 + (Y - Ys)^2 =
                // X^2 - 2 * Xs * X + Xs^2 + (Y - Ys)^2 =
                // [Xs^2 + (Y - Ys)^2] + X^2 - 2 * Xs * X =
                // [Xs^2 + (Y - Ys)^2] + X * (X - 2 * Xs)
                span[spanIndex + 1].SpanX = Right; // Avoid extra test for end-of-list.
                int endX = Left;
                for (int i = 0; i <= spanIndex; i++)
                {
                    int startX = endX;
                    endX = span[i + 1].SpanX;
                    int twicePixelX = 2 * span[i].NearestPixelX;
                    int pixelXDist2 = span[i].Distance2;
                    for (int x = startX; x < endX; x++)
                        this[x, y] = pixelXDist2 + x * (x - twicePixelX);
                }
            }
        }

        static int Sq(int x)
        {
            return x * x;
        }

        public delegate bool IsInSetCallback(int x, int y);

        public void Include(IsInSetCallback isInSet)
        {
            Include(Left, Top, Width, Height, isInSet);
        }

        public void Include(int x, int y) { this[x, y] = 0; }
        public void Exclude(int x, int y) { this[x, y] = MaxValue; }
        public void Include(int x, int y, bool inSet) { this[x, y] = inSet ? 0 : MaxValue; }
        public bool IsInSet(int x, int y) { return this[x, y] == 0; }

        public int DistanceSquared(int x, int y)
        {
            return this[x, y];
        }

        public double Distance(int x, int y)
        {
            return System.Math.Sqrt((double)DistanceSquared(x, y));
        }

        public void ExcludeAll()
        {
            // Set to "ouside of set" value.
            Clear(MaxValue);
        }

        public void Include(int left, int top, int width, int height, IsInSetCallback isInSet)
        {
            int right = left + width, bottom = top + height;
            for (int y = top; y < bottom; y++)
                for (int x = left; x < right; x++)
                    this[x, y] = isInSet(x, y) ? 0 : MaxValue;
        }
    }

    //----------------------------------------------------------------------------------------
    // This class determines the nearest pixel rather than the distance to the nearest pixel.
    // The simplest use is to create a new instance specifying the size, then call Transform
    // with a callback delegate that determines which pixels belong in the set to be measured
    // from. The class can then be indexed like a two-dimensional to get the coordinates of
    // the nearest pixels.
    //
    // Note: This class and the distance class share many methods which could be implemented
    // in an intermediate class which they both could be derived from. At least for now, I've
    // chosen not to do that. I'm not sure the distance class is that useful. It's less
    // flexible, and I'm not sure it's much faster. If I decide to eliminate it, the
    // intermediate class would just add complexity.
    class NearestPixelTransform : Table2D<PointInt16>
    {
        public const int MaxValue = 32767;
        const int MaxValueSquared = MaxValue * MaxValue;

        public NearestPixelTransform() : base()
        {
        }

        public NearestPixelTransform(int width, int height) : base(width, height)
        {
        }

        public NearestPixelTransform(int left, int top, int width, int height) :
            base(left, top, width, height)
        {
        }

        public NearestPixelTransform(int left, int top, int width, int height, IsInSetCallback isInSet) :
            base(left, top, width, height)
        {
            Transform(isInSet);
        }

        public NearestPixelTransform(int width, int height, IsInSetCallback isInSet) :
            base(width, height)
        {
            Transform(isInSet);
        }

        public bool Transform()
        {
            TransformColumns();
            return TransformRows();
        }

        public bool Transform(IsInSetCallback isInSet)
        {
            Include(isInSet);
            return Transform();
        }

        public void TransformColumns()
        {
            TransformColumns(Left, Right);
        }

        public void TransformColumns(IsInSetCallback isInSet)
        {
            Include(isInSet);
            TransformColumns();
        }

        public void TransformColumns(int left, int right)
        {
            // The array pixels are initially set to their own coordinates if set, and other
            // coordinates if not set.
            // Replace each element with the distance to the nearest element in the same column
            // that's within the set.
            for (int x = left; x < right; x++)
            {
                // Set the elements to the value to the nearest set-member in the same column.
                short nearestY = -MaxValue;
                for (int y = Top; y < Bottom; y++)
                {
                    if ((this[x, y].X == x) && (this[x, y].Y == y))
                    {
                        nearestY = (short)y;

                        // Check previous pixels to see if this Y is nearer.
                        int by = y;
                        while ((--by >= Top) && (nearestY - by < by - this[x, by].Y))
                            this[x, by] = new PointInt16(x, nearestY);
                    }
                    else
                    {
                        this[x, y] = new PointInt16(x, nearestY);
                    }
                }
            }
        }

        struct SpanXY
        {
            public PointInt16 NearestPixel;    // Nearest pixel associated with span.
            public int SpanX;                  // Beginning of span.

            public SpanXY(PointInt16 nearestPixel, int spanX)
            {
                NearestPixel = nearestPixel;
                SpanX = spanX;
            }
        }

        public bool TransformRows()
        {
            return TransformRows(Top, Bottom);
        }

        public bool TransformRows(int top, int bottom)
        {
            SpanXY[] span = new SpanXY[Width + 1];

            // Each entry in the image contains the nearest set pixel in the same column. These are the
            // candidate pixels.
            for (int y = top; y < bottom; y++)
            {
                // Taking each candidate pixel one at a time from left to right, determine its span. Each new
                // candidate pixel will be nearest to the pixels in the row from some point onward. Since the
                // spans are in the same X order as the candidate pixels, the previous spans are checked from
                // right to left. If the new candidate is closer than a previous candidate for that pixel's
                // entire span, the previous candidate's span is eliminated. Once a candidate is found that's
                // closer at the beginning of its span than the new candidate, the point at which the new
                // candidate becomes closer is appended to the span list as the end of the previous span and
                // the beginning of the new span.
                int spanIndex = 0;
                span[0] = new SpanXY(this[Left, y], Left);
                for (int x = Left + 1; x < Right; x++)
                {
                    int dist2 = Sq(x) + Sq(this[x, y].Y - y);

                    // If the distance to beyond the maximum image distance, don't bother checking further,
                    // since it can't produce a valid span. (This test isn't necessary to the algorithm.)
                    if (dist2 >= MaxValueSquared)
                        continue;

                    // Loop until the previous span can't be eliminated or there are no previous spans.
                    while (true)
                    {
                        // Compare the distance of the new candidate at the beginning of the previous span
                        // to the distance of its associated candidate.
                        SpanXY prevSpan = span[spanIndex];
                        int prevPixelX = prevSpan.NearestPixel.X;
                        int dist2Diff = dist2 - Sq(prevPixelX) - Sq(y - prevSpan.NearestPixel.Y);
                        int twiceXDiff = 2 * (x - prevPixelX);
                        if (dist2Diff > twiceXDiff * prevSpan.SpanX)
                        {
                            // The span for the new candidate starts after the previous span.
                            // Add a new span to the list.
                            int spanX = (dist2Diff + twiceXDiff - 1) / twiceXDiff; // Ceiling integer division.
                            if (spanX < Right)
                                span[++spanIndex] = new SpanXY(this[x, y], spanX);
                            break;
                        }
                        else if (--spanIndex < 0)
                        {
                            // No more previous spans, so start over.
                            spanIndex = 0;
                            span[0] = new SpanXY(this[x, y], Left);
                            break;
                        }
                    }
                }

                // The spans have been determined. Move through the list of spans, filling in nearest pixel for each pixel.
                span[spanIndex + 1].SpanX = Right; // Avoid extra test for end-of-list.
                int endX = Left;
                for (int i = 0; i <= spanIndex; i++)
                {
                    int startX = endX;
                    endX = span[i + 1].SpanX;
                    PointInt16 nearestPixel = span[i].NearestPixel;
                    for (int x = startX; x < endX; x++)
                        this[x, y] = nearestPixel;
                }
            }

            // Return a flag which is false if there were no set pixels in the original image.
            // If there weren't, the nearest-pixel coordinates will be invalid and outside the image range.
            return (this[Left, top].Y >= 0);
        }

        static int Sq(int x)
        {
            return x * x;
        }

        public delegate bool IsInSetCallback(int x, int y);

        public void Include(IsInSetCallback isIncluded)
        {
            Include(Left, Top, Width, Height, isIncluded);
        }

        public void Include(int x, int y) { this[x, y] = new PointInt16(x, y); }
        public void Exclude(int x, int y) { this[x, y] = PointInt16.MaxValue; }
        public void Include(int x, int y, bool inSet)
        {
            this[x, y] = inSet ? new PointInt16(x, y) : PointInt16.MaxValue;
        }

        public bool IsInSet(int x, int y, bool inSet)
        {
            return this[x, y] == new PointInt16(x, y);
        }

        public int DistanceSquared(int x, int y)
        {
            return this[x, y].DistanceSquared(x, y);
        }

        public double Distance(int x, int y)
        {
            return System.Math.Sqrt((double)DistanceSquared(x, y));
        }

        public void Exclude()
        {
            // Set to "ouside of set" value.
            Clear(PointInt16.MaxValue);
        }

        public void Include(int left, int top, int width, int height, IsInSetCallback isIncluded)
        {
            int right = left + width, bottom = top + height;
            for (int y = top; y < bottom; y++)
                for (int x = left; x < right; x++)
                    this[x, y] = isIncluded(x, y) ? new PointInt16(x, y) : PointInt16.MaxValue;
        }
    }

    public struct PointInt16
    {
        public short x, y;
        public PointInt16(short x, short y)
        {
            this.x = x; this.y = y;
        }

        public PointInt16(int x, int y)
        {
            this.x = (short)x; this.y = (short)y;
        }

        public short X
        {
            get { return x; }
            set { x = value; }
        }

        public short Y
        {
            get { return y; }
            set { y = value; }
        }

        public override bool Equals(object obj)
        {
            return (obj is PointInt16) && (this == (PointInt16)obj);
        }
        public override int GetHashCode()
        {
            return (0xffff & x) | (y << 16);
        }
        public static bool operator ==(PointInt16 p0, PointInt16 p1)
        {
            return (p0.x == p1.x) && (p0.y == p1.y);
        }
        public static bool operator !=(PointInt16 p0, PointInt16 p1)
        {
            return !(p0 == p1);
        }
        public static implicit operator System.Drawing.Point(PointInt16 p)
        {
            return new System.Drawing.Point(p.X, p.Y);
        }
        public static explicit operator PointInt16(System.Drawing.Point p)
        {
            return new PointInt16(p.X, p.Y);
        }

        public int DistanceSquared(int x, int y)
        {
            return Sq(this.x - x) + Sq(this.y  - y);
        }
        public int DistanceSquared(PointInt16 p)
        {
            return DistanceSquared(p.X, p.Y);
        }
        public int DistanceSquared(System.Drawing.Point p)
        {
            return DistanceSquared(p.X, p.Y);
        }
        int Sq(int x)
        {
            return x * x;
        }

        public static readonly PointInt16 MaxValue = new PointInt16(32767, 32767);
    }

    // Addressed like an array, but the starting points can be nonzero.
    // The indices aren't range checked.
    public class Table2D<ElementType>
    {
        int left, top, width, height;
        int right, bottom;
        int size, arraySize;
        int offset;

        public ElementType[] Array;

        public Table2D()
        {
        }

        public Table2D(int width, int height)
        {
            Resize(0, 0, width, height);
        }

        public Table2D(int left, int top, int width, int height)
        {
            Resize(left, top, width, height);
        }

        public void Resize(int width, int height, bool reallocate = false)
        {
            Resize(0, 0, width, height, reallocate);
        }

        public void Resize(int left, int top, int width, int height, bool reallocate = false)
        {
            this.left = left; this.top = top; this.width = width; this.height = height;
            right = left + width; bottom = top + height;

            size = width * height;
            offset = left + width * top;
            if (reallocate || (Array == null) || (arraySize < size))
                CreateArray(size);
        }

        void CreateArray(int arraySize)
        {
            this.arraySize = arraySize;
            Array = new ElementType[arraySize];
        }

        public ElementType this[int x, int y]
        {
            get { return Array[x + y * width - offset]; }
            set { Array[x + y * width - offset] = value; }
        }

        public void Clear(ElementType clearValue)
        {
            for (int i = 0; i < size; i++)
                Array[i] = clearValue;
        }

        public int ToIndex(int x, int y) { return x + y * width; }
        public int ToArrayIndex(int x, int y) { return x + y * width - offset; }

        public int Left { get { return left;  } }
        public int Top { get { return top; } }
        public int Width { get { return width; } }
        public int Height { get { return height; } }
        public int Right { get { return right; } }
        public int Bottom { get { return bottom; } }
        public int Size { get { return size; } }
        public int ArraySize { get { return arraySize; } }
    }
}

I've made many changes with the objective of making it easier to use in PDN (at the expense of making it somewhat more complex and obscure).

I changed back to using non-static methods, and replaced the arrays required by the old version with indexing that makes the classes themselves addressable like arrays. I'll explain the NearestPixelTransform version, which works almost identically to the DistanceTransform version.

In order to use it, an instance must be created in the usual way. To work with selections, the rectangle it covers can have a nonzero origin.

NearestPixelTransform NearestPixels = new NearestPixelTransform(left, top, width, height);

The next step is to tell the class which pixels are inside the set from which the nearest pixels are selected. This can be done by calling the Include method with a delegate that returns true for inside pixels and false for outside pixels. For example (using a lambda expression for the delegate):

NearestPixels.Include((x, y) => (src[x, y].A >= OpacityThreshold));

Then the transform must be performed:

NearestPixels.Transform();

The class can then be indexed (with indexes within the originally specified rectangle). It returns a special PointInt16 value, which is just a point with short components. It can be implicitly converted to a System.Drawing.Point.

dst[x, y] = src[NearestPixels[x, y]];

There a variety of other methods and such for flexibility or convenience. For example, the transform can be computed with:

NearestPixelTransform NearestPixels = new NearestPixelTransform(src.Width, src.Height,
    delegate(int x, int y) { return src[x, y].A >= OpacityThreshold; });

lynxster4 · July 16, 2019

Hi @MJW!

Is it possible to do this? I use EdgeShader so much, I've often wished I could offset the middle color a bit.

This is just a simple example, but it would help soooo much if it's possible.

Thanks!

MJW · July 16, 2019

Unfortunately, I don't think that's possible. I'll give it some thought, just in case I've overlooked something, but I don't believe the algorithm allows that.

Red ochre · July 16, 2019

Possibly add an offset angle chooser and a '% of Max distance offset slider. Use the offset angle added to the angle to the edge pixel (sine ... or cosine) to multiply the % distance (so always smaller than max dist) and subtract/add that from the max distance divisor for the shading?
... Just thinking out loud ... quite possible that approach wouldn't work with the fast algorithm! Hope that made some sense?
Good luck MJW.

Was thinking about multiplying all distances by 255 then adding the object edge pixel alpha before dividing by 255, to get smoother results too?

Just ideas to try .. if MJW hasn't explored these already.

MJW · July 16, 2019

4 hours ago, Red ochre said:

Was thinking about multiplying all distances by 255 then adding the object edge pixel alpha before dividing by 255, to get smoother results too?

That's an excellent suggestion, and one I've been thinking about for a while. Whether or not it will work will need to be seen. I have reasons it might, and reasons it might not.

I plan to work on an updated version of the Edge Shader as soon as I finish the plugin I'm currently working on (which is nearly complete, and which took a lot longer than I anticipated -- more time than it's worth, I'd say.)

The problem with the angle offset is that the distance measure must have certain characteristics in order for the fast algorithm to work, and I don't believe a non-symmetric measure that depends on the angle to the pixel qualifies. On the other hand, I can't say for sure without thinking through the details.

EDIT: Let me add a few details about the transparency suggestion.

The reason it might not work is that the transparency of an edge doesn't reveal which direction the edge is shifted. That is, you can't tell from the transparency whether the "actual" edge that was antialiased to produce the non-opaqueness is horizontal, and close to the X side, vertical, and close to the Y side, or angled. This was a major problem I encountered when I tried to use the idea in the Object Rounder plugin. However, for the nearest pixel, the distance from a pixel to an edge is generally minimized when the edge is perpendicular to the segment from the pixel to the edge, so it seems like we can generally assume the offset to the actual edge is in that direction. That may help make the opacity a valid measure of the extra distance.

Also, height errors when lighted are much more obvious than color errors, so what didn't work for the Object Rounder might work fine for the Edge Shader. (However, I hope to produce a separate height version of the Edge Shader, so I hope the method works well for that.)

lynxster4 · July 17, 2019

@Red ochre I'm so pleased to see you on here again! You've been sorely missed.

Hope all is well with you and you continue to 'pop' in!

@MJW Hopefully some of Red's suggestions will pan out and I appreciate the time you put into these plugins.

lynxster4 · October 28, 2020

@MJW Any hopes of jumping back on this? This is one of my 'go-to' plugins. Thanks!

MJW · October 28, 2020

5 hours ago, lynxster4 said:

@MJW Any hopes of jumping back on this? This is one of my 'go-to' plugins. Thanks!

Sorry about the long, long delay. I will get back on it.

I really wanted to add the type of smoothing Red Ochre referred to here:

On 7/16/2019 at 3:47 PM, Red ochre said:

Was thinking about multiplying all distances by 255 then adding the object edge pixel alpha before dividing by 255, to get smoother results too?

I tried a number of approaches, but so far haven't found anything satisfactory. As I unfortunately sometimes do in that situation, I sort of put it aside in frustration, and worked on other things.

If I can't come up with a solution to that problem soon, I'll release a version without it, then try to solve the smoothing problem, so I can release an update with it.

lynxster4 · October 29, 2020

On 7/16/2019 at 4:51 PM, lynxster4 said:

Is it possible to do this? I use EdgeShader so much, I've often wished I could offset the middle color a bit.

This is just a simple example, but it would help soooo much if it's possible.

Thanks!

Then, this might be impossible to implement?

MJW · October 29, 2020

2 hours ago, lynxster4 said:

Then, this might be impossible to implement?

I believe offsetting is possible. Whether it will work exactly as you wish is hard to say.

The part I was having difficulty with is using the edge transparency of antialiased objects to adjust the distance in a way that reduces jagginess. It could be that when I take a fresh look at it, a solution will be clear. That sometimes happens.

MJW · December 31, 2020

On 10/28/2020 at 11:50 AM, MJW said:

Sorry about the long, long delay. I will get back on it.

Just want to apologize for the further delay and assure you I haven't forgotten about this. I have not, for whatever reason, been too industrious lately, but hope I can release a version of this plugin soon.

Sign In

Fast algorithm for finding the minimum distances to a set

Recommended Posts

MJW

BoltBait

MJW

Eli

MJW

Red ochre

Seerose

Rick Brewster

MJW

MJW

MJW

Rick Brewster

MJW

MJW

MJW

lynxster4

MJW

Red ochre

MJW

lynxster4

lynxster4

MJW

lynxster4

MJW

MJW

Join the conversation

Browse

Activity