SDF.cpp

/* http://www.codersnotes.com/algorithms/signed-distance-fields */

#include <QImage>
#include <time.h>
#include <math.h>
#include <stdio.h>
#include <sys/time.h>
static inline double get_wall_time()
{
    struct timeval time;
    if (gettimeofday(&time,NULL))
    {
        return 0;
    }
    return (double)time.tv_sec + (double)time.tv_usec * .000001;
}
static double perf_start_time;
static double perf_end_time;
static inline double get_cpu_time()
{
    return (double)clock() / CLOCKS_PER_SEC;
}

static inline void perf_timer_start()
{
    perf_start_time = get_wall_time();
}

static inline void perf_timer_end()
{
    perf_end_time = get_wall_time();
}

static inline double perf_timer_get()
{
    return perf_end_time - perf_start_time;
}

static inline void perf_timer_print()
{
    printf("%.3g\n", perf_timer_get());
}

#define DIST_CACHE


struct Point
{
    short dx, dy;
#ifdef DIST_CACHE
    int f;
#endif
};

struct Grid
{
    int w, h;
    Point *grid;
    Grid(int width, int height) : w(width), h(height)
    {
        grid = new Point[(w + 2) * (h + 2)];
    }
    ~Grid()
    {
        delete[] grid;
    }
};
#ifdef DIST_CACHE
static const Point pointInside = { 0, 0, 0 };
static const Point pointEmpty = { 9999, 9999,9999*9999 };
#else
static const Point pointInside = { 0, 0 };
static const Point pointEmpty = { 9999, 9999 };
#endif

static inline Point Get(Grid &g, int x, int y)
{
    return g.grid[y * (g.w + 2) + x];
}

static inline void Put(Grid &g, int x, int y, const Point &p)
{
    g.grid[y * (g.w + 2) + x] = p;
}
#if 1

#ifndef DIST_CACHE
static inline void Compare(Grid &g, int x, int y, int offsetx, int offsety)
{
    Point p = Get(g, x, y);
    Point other = Get(g, x + offsetx, y + offsety);
    other.dx += offsetx;
    other.dy += offsety;
    int newdist = other.dx*other.dx+other.dy*other.dy;
    if (newdist <= p.dx*p.dx+p.dy*p.dy)
    {
        p = other;
    }
    Put(g, x, y, p);
}
#else
static inline void Compare(Grid &g, int x, int y, int offsetx, int offsety)
{
    int add;
    Point p = Get(g, x, y);
    Point other = Get(g, x + offsetx, y + offsety);
    if(offsety == 0) {
        add = 2 * other.dx + 1;
    }
    else if(offsetx == 0) {
        add = 2 * other.dy + 1;
    }
    else {
        add = 2 * (other.dy + other.dx + 1);
    }
    other.f += add;
    if (other.f < p.f)
    {
        p.f = other.f;
        if(offsety == 0) {
            p.dx = other.dx + 1;
            p.dy = other.dy;
        }
        else if(offsetx == 0) {
            p.dy = other.dy + 1;
            p.dx = other.dx;
        }
        else {
            p.dy = other.dy + 1;
            p.dx = other.dx + 1;
        }
    }
    Put(g, x, y, p);
}
#endif
static void GenerateSDF(Grid &g)
{
    /* forward scan */
    // #pragma omp parallel for
    for(int y = 1; y <= g.h; y++)
    {
        for(int x = 0; x <= g.w; x++)
        {
            Compare(g, x, y,  0, -1);
        }
        for(int x = 1; x <= g.w; x++)
        {
            Compare(g, x, y, -1,  0);
            Compare(g, x, y, -1, -1);
        }
        for(int x = g.w; x >= 0; x--)
        {
            Compare(g, x, y,  1,  0);
            Compare(g, x, y,  1, -1);
        }
    }

    /* backward scan */
    // #pragma omp parallel for
    for(int y = g.h-1; y > 0; y--)
    {
        for(int x = 0; x <= g.w; x++)
        {
            Compare(g, x, y,  0,  1);
        }
        for(int x = 1; x <= g.w; x++)
        {
            Compare(g, x, y, -1,  0);
            Compare(g, x, y, -1,  1);
        }
        for(int x = g.w-1; x >= 0; x--)
        {
            Compare(g, x, y,  1,  0);
            Compare(g, x, y,  1,  1);
        }
    }
}
#else
/* macro is a way faster than inline */
#define Compare(offsetx, offsety)                                              \
do {                                                                           \
    int add;                                                                   \
    Point other = Get(g, x + offsetx, y + offsety);                            \
    if(offsety == 0) {                                                         \
        add = 2 * other.dx + 1;                                                \
    }                                                                          \
    else if(offsetx == 0) {                                                    \
        add = 2 * other.dy + 1;                                                \
    }                                                                          \
    else {                                                                     \
        add = 2 * (other.dy + other.dx + 1);                                   \
    }                                                                          \
    other.f += add;                                                            \
    if (other.f < p.f)                                                         \
    {                                                                          \
        p.f = other.f;                                                         \
        if(offsety == 0) {                                                     \
            p.dx = other.dx + 1;                                               \
            p.dy = other.dy;                                                   \
        }                                                                      \
        else if(offsetx == 0) {                                                \
            p.dy = other.dy + 1;                                               \
            p.dx = other.dx;                                                   \
        }                                                                      \
        else {                                                                 \
            p.dy = other.dy + 1;                                               \
            p.dx = other.dx + 1;                                               \
        }                                                                      \
    }                                                                          \
} while(0)

static void GenerateSDF(Grid &g)
{
    for(int y = 1; y <= g.h; y++)
    {
        for(int x = 1; x <= g.w; x++)
        {
            Point p = Get(g, x, y);
            Compare(-1,  0);
            Compare( 0, -1);
            Compare(-1, -1);
            Compare( 1, -1);
            Put(g, x, y, p);
        }
        for(int x=g.w-1;x>=0;x--)
        {
            Point p = Get( g, x, y );
            Compare(1, 0 );
            Put( g, x, y, p );
        }
    }

    for(int y = g.h; y > 0; y--)
    {
        for(int x = g.w; x > 0; x--)
        {
            Point p = Get(g, x, y);
            Compare( 1,  0);
            Compare( 0,  1);
            Compare(-1,  1);
            Compare( 1,  1);
            Put(g, x, y, p);
        }
        for(int x=0;x<g.w;x++)
        {
            Point p = Get( g, x, y );
            Compare(-1, 0 );
            Put( g, x, y, p );
        }
    }
}
#endif



QImage dfcalculate(QImage *img, int distanceFieldScale, bool transparent)
{
    int x, y;
    int w = img->width(), h = img->height();
    QImage result(w, h, QImage::Format_RGB32);
    Grid grid1(w, h);
    Grid grid2(w, h);
    /* create 1-pixel gap */
    for(x = 0; x < w + 2; x++)
    {
        Put(grid1, x, 0, pointInside);
        Put(grid2, x, 0, pointEmpty);
        Put(grid1, x, h + 1, pointInside);
        Put(grid2, x, h + 1, pointEmpty);
    }
    for(y = 1; y <= h; y++)
    {
        Put(grid1, 0, y, pointInside);
        Put(grid2, 0, y, pointEmpty);
        // #pragma omp parallel for
        for(x = 1; x <= w; x++)
        {
            if(qGreen(img->pixel(x - 1, y - 1)) > 128)
            {
                Put(grid1, x, y, pointEmpty);
                Put(grid2, x, y, pointInside);
            }
            else
            {
                Put(grid1, x, y, pointInside);
                Put(grid2, x, y, pointEmpty);
            }
        }
        Put(grid1, w + 1, y, pointInside);
        Put(grid2, w + 1, y, pointEmpty);
    }
    perf_timer_start();
    /* calculation */
    #pragma omp parallel sections
    {
        #pragma omp section
        {
            GenerateSDF(grid1);
        }
        #pragma omp section
        {
            GenerateSDF(grid2);
        }
    }
    perf_timer_end();
    perf_timer_print();
    for(y = 1; y <= h; y++)
        for(x = 1; x <= w; x++)
        {
            #ifdef DIST_CACHE
            double dist1 = sqrt((double)Get(grid1, x, y).f);
            double dist2 = sqrt((double)Get(grid2, x, y).f);
            #else
            double dist1 = sqrt((double)(Get(grid1, x, y).dx*Get(grid1, x, y).dx+Get(grid1, x, y).dy*Get(grid1, x, y).dy));
            double dist2 = sqrt((double)(Get(grid2, x, y).dx*Get(grid2, x, y).dx+Get(grid2, x, y).dy*Get(grid2, x, y).dy));
            #endif
            double dist = dist1 - dist2;
            // Clamp and scale
            int c = fmod(dist2/16, 1) * 255;
            if(transparent)
                result.setPixel(x - 1, y - 1, qRgba(255,255,255,c));
            else
                result.setPixel(x - 1, y - 1, qRgb(c,c,c));
        }
    // printf("%d %d\n", w, h);
    return result;
}


QImage dfcalculate_(QImage *img, int distanceFieldScale, bool transparent)
{
    const int window = 300;
    int x, y;
    int w = img->width(), h = img->height();
    perf_timer_start();
    QImage result(w, h, QImage::Format_ARGB32);
    bool *ar = new bool[w*h];
    for(y = 0; y < h; y++)
    {
        for(x = 0; x < w; x++)
        {
            ar[y*w+x] = qGreen(img->pixel(x, y)) > 128;
        }
    }
    for(y = 0; y < h; y++)
    {
        for(x = 0; x < w; x++)
        {
            int min = 999999;
            for(int i = -window; i < window; i++)
            {
                for(int j = -window; j < window; j++)
                {
                    if(x+i >=0 && x+i < w && y+j >=0 && y+j < h)
                    {
                        int d = i*i+j*j;
                        if(d < min && ar[(y+j)*w+i+x])
                        {
                            min = d;
                        }
                    }
                }
            }
            int c = fmod(sqrt(min)/16, 1) * 255;
            result.setPixel(x, y, qRgb(c,c,c));
        }
    }
    perf_timer_end();
    perf_timer_print();
    printf("%d %d\n", w, h);
    return result;
}