[C++] "dead reckoning" SDF generator

deadreckoning.cpp

// "dead reckoning" transform function which generates signed distance field
// from binary image.
//
// "The Dead Reckoning algorithm" is described in Fig.3 of the following paper:
//   The "dead reckoning" signed distance transform by George J. Grevera
//   http://perso.ensta-paristech.fr/~manzaner/Download/IAD/Grevera_04.pdf
//
// License:
//   Copyright (C) 2021, Takayuki Matsuoka.
//   SPDX-License-Identifier: CC0-1.0

#define DEAD_RECKONING_TEST // Define this symbol if you need test driver.
#include <functional>
#include <vector>
#include <algorithm>

//  width               Width of the binary image in pixels.
//  height              Height of the binary image in pixels.
//  getBinaryPixel      A function which returns pixel of binary image.  true:inside, false:outside.
//  setSignedDistance   A function which sets signed distance to the signed distance image.
//                      - Distance has sign, in pixels.
//                      - Positive means inside of binary image.  Negative means outside.
//                      - Distance is not normalized.  Therefore, -6.5f means "6.5 pixel, outside"
//                        +1.2f means "1.2 pixels, inside".  You can find example of "normalization"
//                        of signed distance in proc().
void deadReckoning(
      int width
    , int height
    , const std::function<bool(int x, int y)>& getBinaryPixel
    , const std::function<void(int x, int y, float signedDistance)>& setSignedDistance
) {
    const auto contains = [&](int x, int y) -> bool {
        return x >= 0 && x < width && y >= 0 && y < height;
    };

    // I - a 2D binary image
    const bool outside = false;
    const bool inside = true;
    const auto I = [&](int x, int y) -> bool { return getBinaryPixel(x, y) ? inside : outside; };

    // d - a 2D grey image representing the distance image
    std::vector<float> ds(width * height);
    const float inf = static_cast<float>(256 * 256 * 256);
    const auto getd = [&](int x, int y) -> float { return ds[y*width + x]; };
    const auto setd = [&](int x, int y, float v) { ds[y*width + x] = v; };
    const auto distance = [](int x, int y) -> float {
        return hypotf(static_cast<float>(x), static_cast<float>(y));
    };

    // p - for each pixel, the corresponding border point
    struct Point { int x, y; };
    std::vector<Point> ps(width * height);
    const Point outOfBounds = { -1, -1 };
    const auto getp = [&](int x, int y) -> Point { return ps[y*width + x]; };
    const auto setp = [&](int x, int y, const Point& p) { ps[y*width + x] = p; };

    // initialize d
    // initialize immediate interior & exterior elements
    for(int y = 0; y < height; ++y) {
        for(int x = 0; x < width; ++x) {
            const bool c = I(x, y);
            const bool w = contains(x-1, y) ? I(x-1, y) : outside;
            const bool e = contains(x+1, y) ? I(x+1, y) : outside;
            const bool n = contains(x, y-1) ? I(x, y-1) : outside;
            const bool s = contains(x, y+1) ? I(x, y+1) : outside;
            if((w != c) || (e != c) || (n != c) || (s != c)) {
                setd(x, y, 0);
                setp(x, y, { x, y });
            } else {
                setd(x, y, inf);
                setp(x, y, outOfBounds);
            }
        }
    }

    // Perform minimum distance choice for single pixel, single direction.
    enum class Dir {
        NW, N, NE,      // NW=(x-1, y-1), N=(x, y-1), NE=(x+1, y-1)
        W,     E,       // W =(x-1, y),               E =(x+1, y)
        SW, S, SE       // SW=(x-1, y+1), S=(x, y+1), SE=(x+1, y+1)
    };
    const auto f = [&](int x, int y, Dir dir) {
        // d1 - distance between two adjacent pixels in either the x or y direction
        const float d1 = 1.0f;

        // d2 - distance between two diagonally adjacent pixels (sqrt(2))
        const float d2 = d1 * 1.4142135623730950488f;

        int dx, dy;
        float od;
        switch(dir) {
        default:
        case Dir::NW: dx = -1; dy = -1; od = d2; break; // first pass
        case Dir::N:  dx =  0; dy = -1; od = d1; break; // first pass
        case Dir::NE: dx =  1; dy = -1; od = d2; break; // first pass
        case Dir::W:  dx = -1; dy =  0; od = d1; break; // first pass
        case Dir::E:  dx =  1; dy =  0; od = d1; break; // final pass
        case Dir::SW: dx = -1; dy =  1; od = d2; break; // final pass
        case Dir::S:  dx =  0; dy =  1; od = d1; break; // final pass
        case Dir::SE: dx =  1; dy =  1; od = d2; break; // final pass
        }
        const bool b = contains(x+dx, y+dy);
        const float cd = b ? getd(x+dx, y+dy) : inf;
        if(cd + od < getd(x, y)) {
            const Point p = b ? getp(x+dx, y+dy) : outOfBounds;
            setp(x, y, p);
            const int xx = x - p.x;
            const int yy = y - p.y;
            const float nd = distance(xx, yy);
            setd(x, y, nd);
        }
    };

    // perform the first pass
    for(int y = 0; y < height; ++y) {           // top to bottom
        for(int x = 0; x < width; ++x) {        // left to right
            static const Dir dirs[] = { Dir::NW, Dir::N, Dir::NE, Dir::W };
            for(const Dir dir : dirs) {
                f(x, y, dir);
            }
        }
    }

    // perform the final pass
    for(int y = height-1; y >= 0; --y) {        // bottom to top
        for(int x = width-1; x >= 0; --x) {     // right to left
            static const Dir dirs[] = { Dir::E, Dir::SW, Dir::S, Dir::SE };
            for(const Dir dir : dirs) {
                f(x, y, dir);
            }
        }
    }

    // indicate inside & outside
    for(int y = 0; y < height; ++y) {
        for(int x = 0; x < width; ++x) {
            float d = getd(x, y);
            if(I(x, y) == outside) {
                d = -d;                         // Negative distance means outside.
            }
            setSignedDistance(x, y, d);
        }
    }
}


#if defined(DEAD_RECKONING_TEST)
#include <string>
#include <lodepng.h>        // https://lodev.org/lodepng/

static bool proc(int spread, bool useAlpha, const std::string& srcFilename, const std::string& outFilename) {
    printf("Processing (%s), spread=%d, useAlpha=%s\n", srcFilename.c_str(), spread, useAlpha ? "true" : "false");
    if(spread <= 0) {
        printf("ERROR: spread must be greater than 0\n");
        return false;
    }

    std::vector<unsigned char> srcRgbaImage;
    unsigned int width, height;
    {
        const unsigned int decoderError = lodepng::decode(srcRgbaImage, width, height, srcFilename);
        if(decoderError) {
            printf("ERROR: lodepng decoder error : %s\n", lodepng_error_text(decoderError));
            return false;
        }
    }

    std::vector<unsigned char> outRgbaImage(width * height * 4);
    {
        // Convert an RGBA pixel to binary.
        const auto getBinaryPixel = [&](int x, int y) -> bool {
            const int o = (y * width + x) * 4;
            const uint8_t r = srcRgbaImage[o + 0];
            const uint8_t g = srcRgbaImage[o + 1];
            const uint8_t b = srcRgbaImage[o + 2];
            const uint8_t a = srcRgbaImage[o + 3];

            if(useAlpha) {
                return (a >= 128);
            } else {
                return ((r + g + b) >= 128);
            }
        };

        // Normalize signed distance to [0,255] and write it as a RGBA pixel.
        const auto setSignedDistance = [&](int x, int y, float signedDistance) {
            const float d = ((signedDistance / spread) + 1.0f) * 128.0f;
            const auto c = static_cast<uint8_t>(std::clamp(d, 0.0f, 255.0f));
            const int o = (y * width + x) * 4;
            outRgbaImage[o + 0] = c;    // R
            outRgbaImage[o + 1] = c;    // G
            outRgbaImage[o + 2] = c;    // B
            outRgbaImage[o + 3] = c;    // A
        };

        deadReckoning(width, height, getBinaryPixel, setSignedDistance);
    }

    {
        const unsigned int encoderError = lodepng::encode(outFilename, outRgbaImage, width, height);
        if(encoderError) {
            printf("ERROR: lodepng encoder error : %s\n", lodepng_error_text(encoderError));
            return false;
        }
    }
    printf("Done (%s)\n", outFilename.c_str());
    return true;
}


int main(int argc, const char* argv[]) {
    int spread = 8;
    bool help = false;
    bool useAlpha = false;
    std::string inpFilename;
    std::string outFilename;
    for(int i = 1; i < argc; ++i) {
        std::string arg = argv[i];
        std::string nextArg = (i+1 >= argc) ? "" : argv[i+1];

        if(arg == "--help") {
            help = true;
            continue;
        }
        if(arg == "--use-alpha") {
            useAlpha = true;
            continue;
        }
        if(arg == "--spread") {
            spread = std::stoi(nextArg);
            ++i;
            continue;
        }
        if(arg == "--input") {
            inpFilename = nextArg;
            ++i;
            continue;
        }
        if(arg == "--output") {
            outFilename = nextArg;
            ++i;
            continue;
        }
    }
    if(inpFilename.empty()) {
        printf("--input is empty\n");
    }
    if(help || inpFilename.empty()) {
        printf("Usage: %s --input=<INP.png> --output=<OUT.png> --spread=<SPREAD>\n", argv[0]);
        printf("  --input <INP.png>  : Input image filename\n");
        printf("  --output <OUT.png> : Output signed distance image filename\n");
        printf("  --spread <SPREAD>  : Spread of signed distance in pixels\n");
        printf("  --use-alpha        : Use alpha channel to make binary image\n");
        printf("\n");
        printf("Example:\n");
        printf("  %s --input my-file.png --output my-sdf.png --spread 8\n", argv[0]);
        return EXIT_FAILURE;
    }
    if(outFilename.empty()) {
        outFilename = inpFilename + ".spread=" + std::to_string(spread) + ".out.png";
    }
    const bool result = proc(spread, useAlpha, inpFilename, outFilename);
    return result ? EXIT_SUCCESS : EXIT_FAILURE;
}
#endif // DEAD_RECKONING_TEST