// David Eberly, Geometric Tools, Redmond WA 98052
// Copyright (c) 1998-2020
// Distributed under the Boost Software License, Version 1.0.
// https://www.boost.org/LICENSE_1_0.txt
// https://www.geometrictools.com/License/Boost/LICENSE_1_0.txt
// Version: 4.0.2019.08.13

#pragma once

#include <Mathematics/Image2.h>
#include <cmath>
#include <functional>
#include <limits>

// Image utilities for Image2<int> objects.  TODO: Extend this to a template
// class to allow the pixel type to be int*_t and uint*_t for * in
// {8,16,32,64}.
//
// All but the Draw* functions are operations on binary images.  Let the image
// have d0 columns and d1 rows.  The input image must have zeros on its
// boundaries x = 0, x = d0-1, y = 0 and y = d1-1.  The 0-valued pixels are
// considered to be background.  The 1-valued pixels are considered to be
// foreground.  In some of the operations, to save memory and time the input
// image is modified by the algorithms.  If you need to preserve the input
// image, make a copy of it before calling these functions.  Dilation and
// erosion functions do not have the requirement that the boundary pixels of
// the binary image inputs be zero.

namespace WwiseGTE
{
    class ImageUtility2
    {
    public:
        // Compute the 4-connected components of a binary image.  The input
        // image is modified to avoid the cost of making a copy.  On output,
        // the image values are the labels for the components.  The array
        // components[k], k >= 1, contains the indices for the k-th component.
        static void GetComponents4(Image2<int>& image,
            std::vector<std::vector<size_t>>& components)
        {
            std::array<int, 4> neighbors;
            image.GetNeighborhood(neighbors);
            GetComponents(4, &neighbors[0], image, components);
        }

        // Compute the 8-connected components of a binary image.  The input
        // image is modified to avoid the cost of making a copy.  On output,
        // the image values are the labels for the components.  The array
        // components[k], k >= 1, contains the indices for the k-th component.
        static void GetComponents8(Image2<int>& image,
            std::vector<std::vector<size_t>>& components)
        {
            std::array<int, 8> neighbors;
            image.GetNeighborhood(neighbors);
            GetComponents(8, &neighbors[0], image, components);
        }

        // Compute a dilation with a structuring element consisting of the
        // 4-connected neighbors of each pixel.  The input image is binary
        // with 0 for background and 1 for foreground.  The output image must
        // be an object different from the input image.
        static void Dilate4(Image2<int> const& input, Image2<int>& output)
        {
            std::array<std::array<int, 2>, 4> neighbors;
            input.GetNeighborhood(neighbors);
            Dilate(input, 4, &neighbors[0], output);
        }

        // Compute a dilation with a structuring element consisting of the
        // 8-connected neighbors of each pixel.  The input image is binary
        // with 0 for background and 1 for foreground.  The output image must
        // be an object different from the input image.
        static void Dilate8(Image2<int> const& input, Image2<int>& output)
        {
            std::array<std::array<int, 2>, 8> neighbors;
            input.GetNeighborhood(neighbors);
            Dilate(input, 8, &neighbors[0], output);
        }

        // Compute a dilation with a structing element consisting of neighbors
        // specified by offsets relative to the pixel.  The input image is
        // binary with 0 for background and 1 for foreground.  The output
        // image must be an object different from the input image.
        static void Dilate(Image2<int> const& input, int numNeighbors,
            std::array<int, 2> const* neighbors, Image2<int>& output)
        {
            LogAssert(&output != &input, "Input and output must be different.");

            output = input;

            // If the pixel at (x,y) is 1, then the pixels at (x+dx,y+dy) are
            // set to 1 where (dx,dy) is in the 'neighbors' array.  Boundary
            // testing is used to avoid accessing out-of-range pixels.
            int const dim0 = input.GetDimension(0);
            int const dim1 = input.GetDimension(1);
            for (int y = 0; y < dim1; ++y)
            {
                for (int x = 0; x < dim0; ++x)
                {
                    if (input(x, y) == 1)
                    {
                        for (int j = 0; j < numNeighbors; ++j)
                        {
                            int xNbr = x + neighbors[j][0];
                            int yNbr = y + neighbors[j][1];
                            if (0 <= xNbr && xNbr < dim0 && 0 <= yNbr && yNbr < dim1)
                            {
                                output(xNbr, yNbr) = 1;
                            }
                        }
                    }
                }
            }
        }

        // Compute an erosion with a structuring element consisting of the
        // 4-connected neighbors of each pixel.  The input image is binary 
        // with 0 for background and 1 for foreground.  The output image must
        // be an object different from the input image.  If zeroExterior is
        // true, the image exterior is assumed to be 0, so 1-valued boundary
        // pixels are set to 0; otherwise, boundary pixels are set to 0 only
        // when they have neighboring image pixels that are 0.
        static void Erode4(Image2<int> const& input, bool zeroExterior, Image2<int>& output)
        {
            std::array<std::array<int, 2>, 4> neighbors;
            input.GetNeighborhood(neighbors);
            Erode(input, zeroExterior, 4, &neighbors[0], output);
        }

        // Compute an erosion with a structuring element consisting of the
        // 8-connected neighbors of each pixel.  The input image is binary
        // with 0 for background and 1 for foreground.  The output image must
        // be an object different from the input image.  If zeroExterior is
        // true, the image exterior is assumed to be 0, so 1-valued boundary
        // pixels are set to 0; otherwise, boundary pixels are set to 0 only
        // when they have neighboring image pixels that are 0.
        static void Erode8(Image2<int> const& input, bool zeroExterior, Image2<int>& output)
        {
            std::array<std::array<int, 2>, 8> neighbors;
            input.GetNeighborhood(neighbors);
            Erode(input, zeroExterior, 8, &neighbors[0], output);
        }

        // Compute an erosion with a structuring element consisting of
        // neighbors specified by offsets relative to the pixel.  The input
        // image is binary with 0 for background and 1 for foreground.  The
        // output image must be an object different from the input image.  If
        // zeroExterior is true, the image exterior is assumed to be 0, so
        // 1-valued boundary pixels are set to 0; otherwise, boundary pixels
        // are set to 0 only when they have neighboring image pixels that
        // are 0.
        static void Erode(Image2<int> const& input, bool zeroExterior,
            int numNeighbors, std::array<int, 2> const* neighbors, Image2<int>& output)
        {
            LogAssert(&output != &input, "Input and output must be different.");

            output = input;

            // If the pixel at (x,y) is 1, it is changed to 0 when at least
            // one neighbor (x+dx,y+dy) is 0, where (dx,dy) is in the
            // 'neighbors' array.
            int const dim0 = input.GetDimension(0);
            int const dim1 = input.GetDimension(1);
            for (int y = 0; y < dim1; ++y)
            {
                for (int x = 0; x < dim0; ++x)
                {
                    if (input(x, y) == 1)
                    {
                        for (int j = 0; j < numNeighbors; ++j)
                        {
                            int xNbr = x + neighbors[j][0];
                            int yNbr = y + neighbors[j][1];
                            if (0 <= xNbr && xNbr < dim0 && 0 <= yNbr && yNbr < dim1)
                            {
                                if (input(xNbr, yNbr) == 0)
                                {
                                    output(x, y) = 0;
                                    break;
                                }
                            }
                            else if (zeroExterior)
                            {
                                output(x, y) = 0;
                                break;
                            }
                        }
                    }
                }
            }
        }

        // Compute an opening with a structuring element consisting of the
        // 4-connected neighbors of each pixel.  The input image is binary
        // with 0 for background and 1 for foreground.  The output image must
        // be an object different from the input image.  If zeroExterior is
        // true, the image exterior is assumed to consist of 0-valued pixels;
        // otherwise, the image exterior is assumed to consist of 1-valued
        // pixels.
        static void Open4(Image2<int> const& input, bool zeroExterior, Image2<int>& output)
        {
            Image2<int> temp(input.GetDimension(0), input.GetDimension(1));
            Erode4(input, zeroExterior, temp);
            Dilate4(temp, output);
        }

        // Compute an opening with a structuring element consisting of the
        // 8-connected neighbors of each pixel.  The input image is binary
        // with 0 for background and 1 for foreground.  The output image must
        // be an object different from the input image.  If zeroExterior is
        // true, the image exterior is assumed to consist of 0-valued pixels;
        // otherwise, the image exterior is assumed to consist of 1-valued
        // pixels.
        static void Open8(Image2<int> const& input, bool zeroExterior, Image2<int>& output)
        {
            Image2<int> temp(input.GetDimension(0), input.GetDimension(1));
            Erode8(input, zeroExterior, temp);
            Dilate8(temp, output);
        }

        // Compute an opening with a structuring element consisting of
        // neighbors specified by offsets relative to the pixel.  The input
        // image is binary with 0 for background and 1 for foreground.  The
        // output image must be an object different from the input image.  If
        // zeroExterior is true, the image exterior is assumed to consist of
        // 0-valued pixels; otherwise, the image exterior is assumed to
        // consist of 1-valued pixels.
        static void Open(Image2<int> const& input, bool zeroExterior,
            int numNeighbors, std::array<int, 2> const* neighbors, Image2<int>& output)
        {
            Image2<int> temp(input.GetDimension(0), input.GetDimension(1));
            Erode(input, zeroExterior, numNeighbors, neighbors, temp);
            Dilate(temp, numNeighbors, neighbors, output);
        }

        // Compute a closing with a structuring element consisting of the
        // 4-connected neighbors of each pixel.  The input image is binary
        // with 0 for background and 1 for foreground.  The output image must
        // be an object different from the input image.  If zeroExterior is
        // true, the image exterior is assumed to consist of 0-valued pixels;
        // otherwise, the image exterior is assumed to consist of 1-valued
        // pixels.
        static void Close4(Image2<int> const& input, bool zeroExterior, Image2<int>& output)
        {
            Image2<int> temp(input.GetDimension(0), input.GetDimension(1));
            Dilate4(input, temp);
            Erode4(temp, zeroExterior, output);
        }

        // Compute a closing with a structuring element consisting of the
        // 8-connected neighbors of each pixel.  The input image is binary
        // with 0 for background and 1 for foreground.  The output image must
        // be an object different from the input image.  If zeroExterior is
        // true, the image exterior is assumed to consist of 0-valued pixels;
        // otherwise, the image exterior is assumed to consist of 1-valued
        // pixels.
        static void Close8(Image2<int> const& input, bool zeroExterior, Image2<int>& output)
        {
            Image2<int> temp(input.GetDimension(0), input.GetDimension(1));
            Dilate8(input, temp);
            Erode8(temp, zeroExterior, output);
        }

        // Compute a closing with a structuring element consisting of
        // neighbors specified by offsets relative to the pixel.  The input
        // image is binary with 0 for background and 1 for foreground.  The
        // output image must be an object different from the input image.  If
        // zeroExterior is true, the image exterior is assumed to consist of
        // 0-valued pixels; otherwise, the image exterior is assumed to
        // consist of 1-valued pixels.
        static void Close(Image2<int> const& input, bool zeroExterior,
            int numNeighbors, std::array<int, 2> const* neighbors, Image2<int>& output)
        {
            Image2<int> temp(input.GetDimension(0), input.GetDimension(1));
            Dilate(input, numNeighbors, neighbors, temp);
            Erode(temp, zeroExterior, numNeighbors, neighbors, output);
        }

        // Locate a pixel and walk around the edge of a component.  The input
        // (x,y) is where the search starts for a nonzero pixel.  If (x,y) is
        // outside the component, the walk is around the outside the
        // component.  If the component has a hole and (x,y) is inside that
        // hole, the walk is around the boundary surrounding the hole.  The
        // function returns 'true' on a success walk.  The return value is
        // 'false' when no boundary was found from the starting (x,y).
        static bool ExtractBoundary(int x, int y, Image2<int>& image, std::vector<size_t>& boundary)
        {
            // Find a first boundary pixel.
            size_t const numPixels = image.GetNumPixels();
            size_t i;
            for (i = image.GetIndex(x, y); i < numPixels; ++i)
            {
                if (image[i])
                {
                    break;
                }
            }
            if (i == numPixels)
            {
                // No boundary pixel found.
                return false;
            }

            std::array<int, 8> const dx = { -1,  0, +1, +1, +1,  0, -1, -1 };
            std::array<int, 8> const dy = { -1, -1, -1,  0, +1, +1, +1,  0 };

            // Create a new point list that contains the first boundary point.
            boundary.push_back(i);

            // The direction from background 0 to boundary pixel 1 is
            // (dx[7],dy[7]).
            std::array<int, 2> coord = image.GetCoordinates(i);
            int x0 = coord[0], y0 = coord[1];
            int cx = x0, cy = y0;
            int nx = x0 - 1, ny = y0, dir = 7;

            // Traverse the boundary in clockwise order.  Mark visited pixels
            // as 2.
            image(cx, cy) = 2;
            bool notDone = true;
            while (notDone)
            {
                int j, nbr;
                for (j = 0, nbr = dir; j < 8; ++j, nbr = (nbr + 1) % 8)
                {
                    nx = cx + dx[nbr];
                    ny = cy + dy[nbr];
                    if (image(nx, ny))  // next boundary pixel found
                    {
                        break;
                    }
                }

                if (j == 8)
                {
                    // (cx,cy) is isolated
                    notDone = false;
                    continue;
                }

                if (nx == x0 && ny == y0)
                {
                    // boundary traversal completed
                    notDone = false;
                    continue;
                }

                // (nx,ny) is next boundary point, add point to list.  Note
                // that the index for the pixel is computed for the original
                // image, not for the larger temporary image.
                boundary.push_back(image.GetIndex(nx, ny));

                // Mark visited pixels as 2.
                image(nx, ny) = 2;

                // Start search for next point.
                cx = nx;
                cy = ny;
                dir = (j + 5 + dir) % 8;
            }

            return true;
        }

        // Use a depth-first search for filling a 4-connected region.  This is
        // nonrecursive, simulated by using a heap-allocated "stack".  The
        // input (x,y) is the seed point that starts the fill.
        template <typename PixelType>
        static void FloodFill4(Image2<PixelType>& image, int x, int y,
            PixelType foreColor, PixelType backColor)
        {
            // Test for a valid seed.
            int const dim0 = image.GetDimension(0);
            int const dim1 = image.GetDimension(1);
            if (x < 0 || x >= dim0 || y < 0 || y >= dim1)
            {
                // The seed point is outside the image domain, so there is
                // nothing to fill.
                return;
            }

            // Allocate the maximum amount of space needed for the stack.
            // An empty stack has top == -1.
            size_t const numPixels = image.GetNumPixels();
            std::vector<int> xStack(numPixels), yStack(numPixels);

            // Push seed point onto stack if it has the background color.  All
            // points pushed onto stack have background color backColor.
            int top = 0;
            xStack[top] = x;
            yStack[top] = y;

            while (top >= 0)  // stack is not empty
            {
                // Read top of stack.  Do not pop since we need to return to
                // this top value later to restart the fill in a different
                // direction.
                x = xStack[top];
                y = yStack[top];

                // Fill the pixel.
                image(x, y) = foreColor;

                int xp1 = x + 1;
                if (xp1 < dim0 && image(xp1, y) == backColor)
                {
                    // Push pixel with background color.
                    ++top;
                    xStack[top] = xp1;
                    yStack[top] = y;
                    continue;
                }

                int xm1 = x - 1;
                if (0 <= xm1 && image(xm1, y) == backColor)
                {
                    // Push pixel with background color.
                    ++top;
                    xStack[top] = xm1;
                    yStack[top] = y;
                    continue;
                }

                int yp1 = y + 1;
                if (yp1 < dim1 && image(x, yp1) == backColor)
                {
                    // Push pixel with background color.
                    ++top;
                    xStack[top] = x;
                    yStack[top] = yp1;
                    continue;
                }

                int ym1 = y - 1;
                if (0 <= ym1 && image(x, ym1) == backColor)
                {
                    // Push pixel with background color.
                    ++top;
                    xStack[top] = x;
                    yStack[top] = ym1;
                    continue;
                }

                // Done in all directions, pop and return to search other
                // directions of predecessor.
                --top;
            }
        }

        // Compute the L1-distance transform of the binary image. The function
        // returns the maximum distance and a point at which the maximum
        // distance is attained.
        static void GetL1Distance(Image2<int>& image, int& maxDistance, int& xMax, int& yMax)
        {
            int const dim0 = image.GetDimension(0);
            int const dim1 = image.GetDimension(1);
            int const dim0m1 = dim0 - 1;
            int const dim1m1 = dim1 - 1;

            // Use a grass-fire approach, computing distance from boundary to
            // interior one pass at a time.
            bool changeMade = true;
            int distance;
            for (distance = 1, xMax = 0, yMax = 0; changeMade; ++distance)
            {
                changeMade = false;
                int distanceP1 = distance + 1;
                for (int y = 1; y < dim1m1; ++y)
                {
                    for (int x = 1; x < dim0m1; ++x)
                    {
                        if (image(x, y) == distance)
                        {
                            if (image(x - 1, y) >= distance
                                && image(x + 1, y) >= distance
                                && image(x, y - 1) >= distance
                                && image(x, y + 1) >= distance)
                            {
                                image(x, y) = distanceP1;
                                xMax = x;
                                yMax = y;
                                changeMade = true;
                            }
                        }
                    }
                }
            }

            maxDistance = --distance;
        }

        // Compute the L2-distance transform of the binary image.  The maximum
        // distance should not be larger than 100, so you have to ensure this
        // is the case for the input image.  The function returns the maximum
        // distance and a point at which the maximum distance is attained.
        // Comments about the algorithm are in the source file.
        static void GetL2Distance(Image2<int> const& image, float& maxDistance,
            int& xMax, int& yMax, Image2<float>& transform)
        {
            // This program calculates the Euclidean distance transform of a
            // binary input image.  The adaptive algorithm is guaranteed to
            // give exact distances for all distances < 100.  The algorithm
            // was provided John Gauch at University of Kansas.  The following
            // is a quote:
            ///
            // The basic idea is similar to a EDT described recently in PAMI
            // by Laymarie from McGill.  By keeping the dx and dy offset to
            // the nearest edge (feature) point in the image, we can search to
            // see which dx dy is closest to a given point by examining a set
            // of neighbors.  The Laymarie method (and Borgfors) look at a
            // fixed 3x3 or 5x5 neighborhood and call it a day.  What we did
            // was calculate (painfully) what neighborhoods you need to look
            // at to guarentee that the exact distance is obtained.  Thus,
            // you will see in the code, that we L2Check the current distance
            // and depending on what we have so far, we extend the search
            // region.  Since our algorithm for L2Checking the exactness of
            // each neighborhood is on the order N^4, we have only gone to
            // N=100.  In theory, you could make this large enough to get all
            // distances exact.  We have implemented the algorithm to get all
            // distances < 100 to be exact. 
            int const dim0 = image.GetDimension(0);
            int const dim1 = image.GetDimension(1);
            int const dim0m1 = dim0 - 1;
            int const dim1m1 = dim1 - 1;
            int x, y, distance;

            // Create and initialize intermediate images.
            Image2<int> xNear(dim0, dim1);
            Image2<int> yNear(dim0, dim1);
            Image2<int> dist(dim0, dim1);
            for (y = 0; y < dim1; ++y)
            {
                for (x = 0; x < dim0; ++x)
                {
                    if (image(x, y) != 0)
                    {
                        xNear(x, y) = 0;
                        yNear(x, y) = 0;
                        dist(x, y) = std::numeric_limits<int>::max();
                    }
                    else
                    {
                        xNear(x, y) = x;
                        yNear(x, y) = y;
                        dist(x, y) = 0;
                    }
                }
            }

            int const K1 = 1;
            int const K2 = 169;   // 13^2
            int const K3 = 961;   // 31^2
            int const K4 = 2401;  // 49^2
            int const K5 = 5184;  // 72^2

            // Pass in the ++ direction.
            for (y = 0; y < dim1; ++y)
            {
                for (x = 0; x < dim0; ++x)
                {
                    distance = dist(x, y);
                    if (distance > K1)
                    {
                        L2Check(x, y, -1, 0, xNear, yNear, dist);
                        L2Check(x, y, -1, -1, xNear, yNear, dist);
                        L2Check(x, y, 0, -1, xNear, yNear, dist);
                    }
                    if (distance > K2)
                    {
                        L2Check(x, y, -2, -1, xNear, yNear, dist);
                        L2Check(x, y, -1, -2, xNear, yNear, dist);
                    }
                    if (distance > K3)
                    {
                        L2Check(x, y, -3, -1, xNear, yNear, dist);
                        L2Check(x, y, -3, -2, xNear, yNear, dist);
                        L2Check(x, y, -2, -3, xNear, yNear, dist);
                        L2Check(x, y, -1, -3, xNear, yNear, dist);
                    }
                    if (distance > K4)
                    {
                        L2Check(x, y, -4, -1, xNear, yNear, dist);
                        L2Check(x, y, -4, -3, xNear, yNear, dist);
                        L2Check(x, y, -3, -4, xNear, yNear, dist);
                        L2Check(x, y, -1, -4, xNear, yNear, dist);
                    }
                    if (distance > K5)
                    {
                        L2Check(x, y, -5, -1, xNear, yNear, dist);
                        L2Check(x, y, -5, -2, xNear, yNear, dist);
                        L2Check(x, y, -5, -3, xNear, yNear, dist);
                        L2Check(x, y, -5, -4, xNear, yNear, dist);
                        L2Check(x, y, -4, -5, xNear, yNear, dist);
                        L2Check(x, y, -2, -5, xNear, yNear, dist);
                        L2Check(x, y, -3, -5, xNear, yNear, dist);
                        L2Check(x, y, -1, -5, xNear, yNear, dist);
                    }
                }
            }

            // Pass in -- direction.
            for (y = dim1m1; y >= 0; --y)
            {
                for (x = dim0m1; x >= 0; --x)
                {
                    distance = dist(x, y);
                    if (distance > K1)
                    {
                        L2Check(x, y, 1, 0, xNear, yNear, dist);
                        L2Check(x, y, 1, 1, xNear, yNear, dist);
                        L2Check(x, y, 0, 1, xNear, yNear, dist);
                    }
                    if (distance > K2)
                    {
                        L2Check(x, y, 2, 1, xNear, yNear, dist);
                        L2Check(x, y, 1, 2, xNear, yNear, dist);
                    }
                    if (distance > K3)
                    {
                        L2Check(x, y, 3, 1, xNear, yNear, dist);
                        L2Check(x, y, 3, 2, xNear, yNear, dist);
                        L2Check(x, y, 2, 3, xNear, yNear, dist);
                        L2Check(x, y, 1, 3, xNear, yNear, dist);
                    }
                    if (distance > K4)
                    {
                        L2Check(x, y, 4, 1, xNear, yNear, dist);
                        L2Check(x, y, 4, 3, xNear, yNear, dist);
                        L2Check(x, y, 3, 4, xNear, yNear, dist);
                        L2Check(x, y, 1, 4, xNear, yNear, dist);
                    }
                    if (distance > K5)
                    {
                        L2Check(x, y, 5, 1, xNear, yNear, dist);
                        L2Check(x, y, 5, 2, xNear, yNear, dist);
                        L2Check(x, y, 5, 3, xNear, yNear, dist);
                        L2Check(x, y, 5, 4, xNear, yNear, dist);
                        L2Check(x, y, 4, 5, xNear, yNear, dist);
                        L2Check(x, y, 2, 5, xNear, yNear, dist);
                        L2Check(x, y, 3, 5, xNear, yNear, dist);
                        L2Check(x, y, 1, 5, xNear, yNear, dist);
                    }
                }
            }

            // Pass in the +- direction.
            for (y = dim1m1; y >= 0; --y)
            {
                for (x = 0; x < dim0; ++x)
                {
                    distance = dist(x, y);
                    if (distance > K1)
                    {
                        L2Check(x, y, -1, 0, xNear, yNear, dist);
                        L2Check(x, y, -1, 1, xNear, yNear, dist);
                        L2Check(x, y, 0, 1, xNear, yNear, dist);
                    }
                    if (distance > K2)
                    {
                        L2Check(x, y, -2, 1, xNear, yNear, dist);
                        L2Check(x, y, -1, 2, xNear, yNear, dist);
                    }
                    if (distance > K3)
                    {
                        L2Check(x, y, -3, 1, xNear, yNear, dist);
                        L2Check(x, y, -3, 2, xNear, yNear, dist);
                        L2Check(x, y, -2, 3, xNear, yNear, dist);
                        L2Check(x, y, -1, 3, xNear, yNear, dist);
                    }
                    if (distance > K4)
                    {
                        L2Check(x, y, -4, 1, xNear, yNear, dist);
                        L2Check(x, y, -4, 3, xNear, yNear, dist);
                        L2Check(x, y, -3, 4, xNear, yNear, dist);
                        L2Check(x, y, -1, 4, xNear, yNear, dist);
                    }
                    if (distance > K5)
                    {
                        L2Check(x, y, -5, 1, xNear, yNear, dist);
                        L2Check(x, y, -5, 2, xNear, yNear, dist);
                        L2Check(x, y, -5, 3, xNear, yNear, dist);
                        L2Check(x, y, -5, 4, xNear, yNear, dist);
                        L2Check(x, y, -4, 5, xNear, yNear, dist);
                        L2Check(x, y, -2, 5, xNear, yNear, dist);
                        L2Check(x, y, -3, 5, xNear, yNear, dist);
                        L2Check(x, y, -1, 5, xNear, yNear, dist);
                    }
                }
            }

            // Pass in the -+ direction.
            for (y = 0; y < dim1; ++y)
            {
                for (x = dim0m1; x >= 0; --x)
                {
                    distance = dist(x, y);
                    if (distance > K1)
                    {
                        L2Check(x, y, 1, 0, xNear, yNear, dist);
                        L2Check(x, y, 1, -1, xNear, yNear, dist);
                        L2Check(x, y, 0, -1, xNear, yNear, dist);
                    }
                    if (distance > K2)
                    {
                        L2Check(x, y, 2, -1, xNear, yNear, dist);
                        L2Check(x, y, 1, -2, xNear, yNear, dist);
                    }
                    if (distance > K3)
                    {
                        L2Check(x, y, 3, -1, xNear, yNear, dist);
                        L2Check(x, y, 3, -2, xNear, yNear, dist);
                        L2Check(x, y, 2, -3, xNear, yNear, dist);
                        L2Check(x, y, 1, -3, xNear, yNear, dist);
                    }
                    if (distance > K4)
                    {
                        L2Check(x, y, 4, -1, xNear, yNear, dist);
                        L2Check(x, y, 4, -3, xNear, yNear, dist);
                        L2Check(x, y, 3, -4, xNear, yNear, dist);
                        L2Check(x, y, 1, -4, xNear, yNear, dist);
                    }
                    if (distance > K5)
                    {
                        L2Check(x, y, 5, -1, xNear, yNear, dist);
                        L2Check(x, y, 5, -2, xNear, yNear, dist);
                        L2Check(x, y, 5, -3, xNear, yNear, dist);
                        L2Check(x, y, 5, -4, xNear, yNear, dist);
                        L2Check(x, y, 4, -5, xNear, yNear, dist);
                        L2Check(x, y, 2, -5, xNear, yNear, dist);
                        L2Check(x, y, 3, -5, xNear, yNear, dist);
                        L2Check(x, y, 1, -5, xNear, yNear, dist);
                    }
                }
            }

            xMax = 0;
            yMax = 0;
            maxDistance = 0.0f;
            for (y = 0; y < dim1; ++y)
            {
                for (x = 0; x < dim0; ++x)
                {
                    float fdistance = std::sqrt((float)dist(x, y));
                    if (fdistance > maxDistance)
                    {
                        maxDistance = fdistance;
                        xMax = x;
                        yMax = y;
                    }
                    transform(x, y) = fdistance;
                }
            }
        }

        // Compute a skeleton of a binary image.  Boundary pixels are trimmed
        // from the object one layer at a time based on their adjacency to
        // interior pixels.  At each step the connectivity and cycles of the
        // object are preserved.  The skeleton overwrites the contents of the
        // input image.
        static void GetSkeleton(Image2<int>& image)
        {
            int const dim0 = image.GetDimension(0);
            int const dim1 = image.GetDimension(1);

            // Trim pixels, mark interior as 4.
            bool notDone = true;
            while (notDone)
            {
                if (MarkInterior(image, 4, Interior4))
                {
                    // No interior pixels, trimmed set is at most 2-pixels
                    // thick.
                    notDone = false;
                    continue;
                }

                if (ClearInteriorAdjacent(image, 4))
                {
                    // All remaining interior pixels are either articulation
                    // points or part of blobs whose boundary pixels are all
                    // articulation points.  An example of the latter case is
                    // shown below.  The background pixels are marked with '.'
                    // rather than '0' for readability.  The interior pixels
                    // are marked with '4' and the boundary pixels are marked
                    // with '1'.
                    //
                    //   .........
                    //   .....1...
                    //   ..1.1.1..
                    //   .1.141...
                    //   ..14441..
                    //   ..1441.1.
                    //   .1.11.1..
                    //   ..1..1...
                    //   .........
                    //
                    // This is a pathological problem where there are many
                    // small holes (0-pixel with north, south, west, and east
                    // neighbors all 1-pixels) that your application can try
                    // to avoid by an initial pass over the image to fill in
                    // such holes.  Of course, you do have problems with
                    // checkerboard patterns...
                    notDone = false;
                    continue;
                }
            }

            // Trim pixels, mark interior as 3.
            notDone = true;
            while (notDone)
            {
                if (MarkInterior(image, 3, Interior3))
                {
                    // No interior pixels, trimmed set is at most 2-pixels
                    // thick.
                    notDone = false;
                    continue;
                }

                if (ClearInteriorAdjacent(image, 3))
                {
                    // All remaining 3-values can be safely removed since they
                    // are not articulation points and the removal will not
                    // cause new holes.
                    for (int y = 0; y < dim1; ++y)
                    {
                        for (int x = 0; x < dim0; ++x)
                        {
                            if (image(x, y) == 3 && !IsArticulation(image, x, y))
                            {
                                image(x, y) = 0;
                            }
                        }
                    }
                    notDone = false;
                    continue;
                }
            }

            // Trim pixels, mark interior as 2.
            notDone = true;
            while (notDone)
            {
                if (MarkInterior(image, 2, Interior2))
                {
                    // No interior pixels, trimmed set is at most 1-pixel
                    // thick.  Call it a skeleton.
                    notDone = false;
                    continue;
                }

                if (ClearInteriorAdjacent(image, 2))
                {
                    // Removes 2-values that are not articulation points.
                    for (int y = 0; y < dim1; ++y)
                    {
                        for (int x = 0; x < dim0; ++x)
                        {
                            if (image(x, y) == 2 && !IsArticulation(image, x, y))
                            {
                                image(x, y) = 0;
                            }
                        }
                    }
                    notDone = false;
                    continue;
                }
            }

            // Make the skeleton a binary image.
            size_t const numPixels = image.GetNumPixels();
            for (size_t i = 0; i < numPixels; ++i)
            {
                if (image[i] != 0)
                {
                    image[i] = 1;
                }
            }
        }

        // In the remaining public member functions, the callback represents
        // the action you want applied to each pixel as it is visited.

        // Visit pixels in a (2*thick+1)x(2*thick+1) square centered at (x,y).
        static void DrawThickPixel(int x, int y, int thick,
            std::function<void(int, int)> const& callback)
        {
            for (int dy = -thick; dy <= thick; ++dy)
            {
                for (int dx = -thick; dx <= thick; ++dx)
                {
                    callback(x + dx, y + dy);
                }
            }
        }

        // Visit pixels using Bresenham's line drawing algorithm.
        static void DrawLine(int x0, int y0, int x1, int y1,
            std::function<void(int, int)> const& callback)
        {
            // Starting point of line.
            int x = x0, y = y0;

            // Direction of line.
            int dx = x1 - x0, dy = y1 - y0;

            // Increment or decrement depending on direction of line.
            int sx = (dx > 0 ? 1 : (dx < 0 ? -1 : 0));
            int sy = (dy > 0 ? 1 : (dy < 0 ? -1 : 0));

            // Decision parameters for pixel selection.
            if (dx < 0)
            {
                dx = -dx;
            }
            if (dy < 0)
            {
                dy = -dy;
            }
            int ax = 2 * dx, ay = 2 * dy;
            int decX, decY;

            // Determine largest direction component, single-step related
            // variable.
            int maxValue = dx, var = 0;
            if (dy > maxValue)
            {
                var = 1;
            }

            // Traverse Bresenham line.
            switch (var)
            {
            case 0:  // Single-step in x-direction.
                decY = ay - dx;
                for (/**/; /**/; x += sx, decY += ay)
                {
                    callback(x, y);

                    // Take Bresenham step.
                    if (x == x1)
                    {
                        break;
                    }
                    if (decY >= 0)
                    {
                        decY -= ax;
                        y += sy;
                    }
                }
                break;
            case 1:  // Single-step in y-direction.
                decX = ax - dy;
                for (/**/; /**/; y += sy, decX += ax)
                {
                    callback(x, y);

                    // Take Bresenham step.
                    if (y == y1)
                    {
                        break;
                    }
                    if (decX >= 0)
                    {
                        decX -= ay;
                        x += sx;
                    }
                }
                break;
            }
        }

        // Visit pixels using Bresenham's circle drawing algorithm.  Set
        // 'solid' to false for drawing only the circle.  Set 'solid' to true
        // to draw all pixels on and inside the circle.
        static void DrawCircle(int xCenter, int yCenter, int radius, bool solid,
            std::function<void(int, int)> const& callback)
        {
            int x, y, dec;

            if (solid)
            {
                int xValue, yMin, yMax, i;
                for (x = 0, y = radius, dec = 3 - 2 * radius; x <= y; ++x)
                {
                    xValue = xCenter + x;
                    yMin = yCenter - y;
                    yMax = yCenter + y;
                    for (i = yMin; i <= yMax; ++i)
                    {
                        callback(xValue, i);
                    }

                    xValue = xCenter - x;
                    for (i = yMin; i <= yMax; ++i)
                    {
                        callback(xValue, i);
                    }

                    xValue = xCenter + y;
                    yMin = yCenter - x;
                    yMax = yCenter + x;
                    for (i = yMin; i <= yMax; ++i)
                    {
                        callback(xValue, i);
                    }

                    xValue = xCenter - y;
                    for (i = yMin; i <= yMax; ++i)
                    {
                        callback(xValue, i);
                    }

                    if (dec >= 0)
                    {
                        dec += -4 * (y--) + 4;
                    }
                    dec += 4 * x + 6;
                }
            }
            else
            {
                for (x = 0, y = radius, dec = 3 - 2 * radius; x <= y; ++x)
                {
                    callback(xCenter + x, yCenter + y);
                    callback(xCenter + x, yCenter - y);
                    callback(xCenter - x, yCenter + y);
                    callback(xCenter - x, yCenter - y);
                    callback(xCenter + y, yCenter + x);
                    callback(xCenter + y, yCenter - x);
                    callback(xCenter - y, yCenter + x);
                    callback(xCenter - y, yCenter - x);

                    if (dec >= 0)
                    {
                        dec += -4 * (y--) + 4;
                    }
                    dec += 4 * x + 6;
                }
            }
        }

        // Visit pixels in a rectangle of the specified dimensions.  Set
        // 'solid' to false for drawing only the rectangle.  Set 'solid' to
        // true to draw all pixels on and inside the rectangle.
        static void DrawRectangle(int xMin, int yMin, int xMax, int yMax,
            bool solid, std::function<void(int, int)> const& callback)
        {
            int x, y;

            if (solid)
            {
                for (y = yMin; y <= yMax; ++y)
                {
                    for (x = xMin; x <= xMax; ++x)
                    {
                        callback(x, y);
                    }
                }
            }
            else
            {
                for (x = xMin; x <= xMax; ++x)
                {
                    callback(x, yMin);
                    callback(x, yMax);
                }
                for (y = yMin + 1; y <= yMax - 1; ++y)
                {
                    callback(xMin, y);
                    callback(xMax, y);
                }
            }
        }

        // Visit the pixels using Bresenham's algorithm for the axis-aligned
        // ellipse ((x-xc)/a)^2 + ((y-yc)/b)^2 = 1, where xCenter is xc,
        // yCenter is yc, xExtent is a, and yExtent is b.
        static void DrawEllipse(int xCenter, int yCenter, int xExtent, int yExtent,
            std::function<void(int, int)> const& callback)
        {
            int xExtSqr = xExtent * xExtent, yExtSqr = yExtent * yExtent;
            int x, y, dec;

            x = 0;
            y = yExtent;
            dec = 2 * yExtSqr + xExtSqr * (1 - 2 * yExtent);
            for (/**/; yExtSqr * x <= xExtSqr * y; ++x)
            {
                callback(xCenter + x, yCenter + y);
                callback(xCenter - x, yCenter + y);
                callback(xCenter + x, yCenter - y);
                callback(xCenter - x, yCenter - y);

                if (dec >= 0)
                {
                    dec += 4 * xExtSqr * (1 - y);
                    --y;
                }
                dec += yExtSqr * (4 * x + 6);
            }
            if (y == 0 && x < xExtent)
            {
                // The discretization caused us to reach the y-axis before the
                // x-values reached the ellipse vertices.  Draw a solid line
                // along the x-axis to those vertices.
                for (/**/; x <= xExtent; ++x)
                {
                    callback(xCenter + x, yCenter);
                    callback(xCenter - x, yCenter);
                }
                return;
            }

            x = xExtent;
            y = 0;
            dec = 2 * xExtSqr + yExtSqr * (1 - 2 * xExtent);
            for (/**/; xExtSqr * y <= yExtSqr * x; ++y)
            {
                callback(xCenter + x, yCenter + y);
                callback(xCenter - x, yCenter + y);
                callback(xCenter + x, yCenter - y);
                callback(xCenter - x, yCenter - y);

                if (dec >= 0)
                {
                    dec += 4 * yExtSqr * (1 - x);
                    --x;
                }
                dec += xExtSqr * (4 * y + 6);
            }
            if (x == 0 && y < yExtent)
            {
                // The discretization caused us to reach the x-axis before the
                // y-values reached the ellipse vertices.  Draw a solid line
                // along the y-axis to those vertices.
                for (/**/; y <= yExtent; ++y)
                {
                    callback(xCenter, yCenter + y);
                    callback(xCenter, yCenter - y);
                }
            }
        }

        // Use a depth-first search for filling a 4-connected region.  This is
        // nonrecursive, simulated by using a heap-allocated "stack".  The
        // input (x,y) is the seed point that starts the fill.  The x-value is
        // in {0..xSize-1} and the y-value is in {0..ySize-1}.
        template <typename PixelType>
        static void DrawFloodFill4(int x, int y, int xSize, int ySize,
            PixelType foreColor, PixelType backColor,
            std::function<void(int, int, PixelType)> const& setCallback,
            std::function<PixelType(int, int)> const& getCallback)
        {
            // Test for a valid seed.
            if (x < 0 || x >= xSize || y < 0 || y >= ySize)
            {
                // The seed point is outside the image domain, so nothing to
                // fill.
                return;
            }

            // Allocate the maximum amount of space needed for the stack.  An
            // empty stack has top == -1.
            int const numPixels = xSize * ySize;
            std::vector<int> xStack(numPixels), yStack(numPixels);

            // Push seed point onto stack if it has the background color.  All
            // points pushed onto stack have background color backColor.
            int top = 0;
            xStack[top] = x;
            yStack[top] = y;

            while (top >= 0)  // stack is not empty
            {
                // Read top of stack.  Do not pop since we need to return to
                // this top value later to restart the fill in a different
                // direction.
                x = xStack[top];
                y = yStack[top];

                // Fill the pixel.
                setCallback(x, y, foreColor);

                int xp1 = x + 1;
                if (xp1 < xSize && getCallback(xp1, y) == backColor)
                {
                    // Push pixel with background color.
                    ++top;
                    xStack[top] = xp1;
                    yStack[top] = y;
                    continue;
                }

                int xm1 = x - 1;
                if (0 <= xm1 && getCallback(xm1, y) == backColor)
                {
                    // Push pixel with background color.
                    ++top;
                    xStack[top] = xm1;
                    yStack[top] = y;
                    continue;
                }

                int yp1 = y + 1;
                if (yp1 < ySize && getCallback(x, yp1) == backColor)
                {
                    // Push pixel with background color.
                    ++top;
                    xStack[top] = x;
                    yStack[top] = yp1;
                    continue;
                }

                int ym1 = y - 1;
                if (0 <= ym1 && getCallback(x, ym1) == backColor)
                {
                    // Push pixel with background color.
                    ++top;
                    xStack[top] = x;
                    yStack[top] = ym1;
                    continue;
                }

                // Done in all directions, pop and return to search other
                // directions of the predecessor.
                --top;
            }
        }

    private:
        // Connected component labeling using depth-first search.
        static void GetComponents(int numNeighbors, int const* delta,
            Image2<int>& image, std::vector<std::vector<size_t>>& components)
        {
            size_t const numPixels = image.GetNumPixels();
            std::vector<int> numElements(numPixels);
            std::vector<size_t> vstack(numPixels);
            size_t i, numComponents = 0;
            int label = 2;
            for (i = 0; i < numPixels; ++i)
            {
                if (image[i] == 1)
                {
                    int top = -1;
                    vstack[++top] = i;

                    int& count = numElements[numComponents + 1];
                    count = 0;
                    while (top >= 0)
                    {
                        size_t v = vstack[top];
                        image[v] = -1;
                        int j;
                        for (j = 0; j < numNeighbors; ++j)
                        {
                            size_t adj = v + delta[j];
                            if (image[adj] == 1)
                            {
                                vstack[++top] = adj;
                                break;
                            }
                        }
                        if (j == numNeighbors)
                        {
                            image[v] = label;
                            ++count;
                            --top;
                        }
                    }

                    ++numComponents;
                    ++label;
                }
            }

            if (numComponents > 0)
            {
                components.resize(numComponents + 1);
                for (i = 1; i <= numComponents; ++i)
                {
                    components[i].resize(numElements[i]);
                    numElements[i] = 0;
                }

                for (i = 0; i < numPixels; ++i)
                {
                    int value = image[i];
                    if (value != 0)
                    {
                        // Labels started at 2 to support the depth-first
                        // search, so they need to be decremented for the
                        // correct labels.
                        image[i] = --value;
                        components[value][numElements[value]] = i;
                        ++numElements[value];
                    }
                }
            }
        }

        // Support for GetL2Distance.
        static void L2Check(int x, int y, int dx, int dy, Image2<int>& xNear,
            Image2<int>& yNear, Image2<int>& dist)
        {
            int const dim0 = dist.GetDimension(0);
            int const dim1 = dist.GetDimension(1);
            int xp = x + dx, yp = y + dy;
            if (0 <= xp && xp < dim0 && 0 <= yp && yp < dim1)
            {
                if (dist(xp, yp) < dist(x, y))
                {
                    int dx0 = xNear(xp, yp) - x;
                    int dy0 = yNear(xp, yp) - y;
                    int newDist = dx0 * dx0 + dy0 * dy0;
                    if (newDist < dist(x, y))
                    {
                        xNear(x, y) = xNear(xp, yp);
                        yNear(x, y) = yNear(xp, yp);
                        dist(x, y) = newDist;
                    }
                }
            }
        }

        // Support for GetSkeleton.
        static bool Interior2(Image2<int>& image, int x, int y)
        {
            bool b1 = (image(x, y - 1) != 0);
            bool b3 = (image(x + 1, y) != 0);
            bool b5 = (image(x, y + 1) != 0);
            bool b7 = (image(x - 1, y) != 0);
            return (b1 && b3) || (b3 && b5) || (b5 && b7) || (b7 && b1);
        }

        static bool Interior3(Image2<int>& image, int x, int y)
        {
            int numNeighbors = 0;
            if (image(x - 1, y) != 0)
            {
                ++numNeighbors;
            }
            if (image(x + 1, y) != 0)
            {
                ++numNeighbors;
            }
            if (image(x, y - 1) != 0)
            {
                ++numNeighbors;
            }
            if (image(x, y + 1) != 0)
            {
                ++numNeighbors;
            }
            return numNeighbors == 3;
        }

        static bool Interior4(Image2<int>& image, int x, int y)
        {
            return image(x - 1, y) != 0
                && image(x + 1, y) != 0
                && image(x, y - 1) != 0
                && image(x, y + 1) != 0;
        }

        static bool MarkInterior(Image2<int>& image, int value,
            bool (*function)(Image2<int>&, int, int))
        {
            int const dim0 = image.GetDimension(0);
            int const dim1 = image.GetDimension(1);
            bool noInterior = true;
            for (int y = 0; y < dim1; ++y)
            {
                for (int x = 0; x < dim0; ++x)
                {
                    if (image(x, y) > 0)
                    {
                        if (function(image, x, y))
                        {
                            image(x, y) = value;
                            noInterior = false;
                        }
                        else
                        {
                            image(x, y) = 1;
                        }
                    }
                }
            }
            return noInterior;
        }

        static bool IsArticulation(Image2<int>& image, int x, int y)
        {
            static std::array<int, 256> const articulation =
            {
                0,0,0,0,0,1,0,0,0,1,0,0,0,1,0,0,
                0,1,1,1,1,1,1,1,0,1,0,0,0,1,0,0,
                0,1,1,1,1,1,1,1,0,1,0,0,0,1,0,0,
                0,1,1,1,1,1,1,1,0,1,0,0,0,1,0,0,
                0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
                1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
                0,1,1,1,1,1,1,1,0,1,0,0,0,1,0,0,
                0,1,1,1,1,1,1,1,0,1,0,0,0,1,0,0,
                0,0,0,0,1,1,0,0,1,1,0,0,1,1,0,0,
                1,1,1,1,1,1,1,1,1,1,0,0,1,1,0,0,
                0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,
                0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,
                0,0,0,0,1,1,0,0,1,1,0,0,1,1,0,0,
                1,1,1,1,1,1,1,1,1,1,0,0,1,1,0,0,
                0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,
                0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0
            };

            // Converts 8 neighbors of pixel (x,y) to an 8-bit value,
            // bit = 1 iff pixel is set.
            int byteMask = 0;
            if (image(x - 1, y - 1) != 0)
            {
                byteMask |= 0x01;
            }
            if (image(x, y - 1) != 0)
            {
                byteMask |= 0x02;
            }
            if (image(x + 1, y - 1) != 0)
            {
                byteMask |= 0x04;
            }
            if (image(x + 1, y) != 0)
            {
                byteMask |= 0x08;
            }
            if (image(x + 1, y + 1) != 0)
            {
                byteMask |= 0x10;
            }
            if (image(x, y + 1) != 0)
            {
                byteMask |= 0x20;
            }
            if (image(x - 1, y + 1) != 0)
            {
                byteMask |= 0x40;
            }
            if (image(x - 1, y) != 0)
            {
                byteMask |= 0x80;
            }

            return articulation[byteMask] == 1;
        }

        static bool ClearInteriorAdjacent(Image2<int>& image, int value)
        {
            int const dim0 = image.GetDimension(0);
            int const dim1 = image.GetDimension(1);
            bool noRemoval = true;
            for (int y = 0; y < dim1; ++y)
            {
                for (int x = 0; x < dim0; ++x)
                {
                    if (image(x, y) == 1)
                    {
                        bool interiorAdjacent =
                            image(x - 1, y - 1) == value ||
                            image(x, y - 1) == value ||
                            image(x + 1, y - 1) == value ||
                            image(x + 1, y) == value ||
                            image(x + 1, y + 1) == value ||
                            image(x, y + 1) == value ||
                            image(x - 1, y + 1) == value ||
                            image(x - 1, y) == value;

                        if (interiorAdjacent && !IsArticulation(image, x, y))
                        {
                            image(x, y) = 0;
                            noRemoval = false;
                        }
                    }
                }
            }
            return noRemoval;
        }
    };
}