euclidian_distance_transform.cpp

/*
    Original implementation written in JS: https://github.com/mapbox/tiny-sdf/blob/master/index.js
    Paper on distance transform algorithm: http://cs.brown.edu/people/pfelzens/papers/dt-final.pdf
*/

#define INF 1e+20

void _edt1d(unsigned char * grid, unsigned int offset, unsigned int stride, unsigned int length, unsigned char * f, unsigned char * v, float * z)
{
	unsigned int q, k, r;
	float s;

	v[0] = 0;
	z[0] = -INF;
	z[1] = INF;

	for (q = 0; q < length; q++)
		f[q] = grid[offset + q * stride]; // copy the 1d data we want from the 2d data

	for (q = 1, k = 0, s = 0; q < length; q++)
	{
		do
		{
			r = v[k];
			s = (f[q] - f[r] + q * q - r * r) / (q - r) / 2;
		} while (s <= z[k] && --k > -1);

		k++;
		v[k] = q;
		z[k] = s;
		z[k + 1] = INF;
	}

	for (q = 0, k = 0; q < length; q++)
	{
		while (z[k + 1] < q) k++;
			r = v[k];
		grid[offset + q * stride] = f[r] + (q - r) * (q - r);
	}
}

void euclidian_distance_transform(unsigned char * data, unsigned int width, unsigned int height)
{
	int size = width > height ? width : height;
	unsigned char * f = new unsigned char[size];
	unsigned char * v = new unsigned char[size];
	float * z = new float[size + 1];

	for (int x = 0; x < width; x++)
		_edt1d(data, x, width, height, f, v, z);
	for (int y = 0; y < height; y++)
		_edt1d(data, y * width, 1, width, f, v, z);
}