Here is a C# code for 1D squared euclidean distance transform according to the Felzenszwald & Huttenlocher's paper:
private static void DistanceTransform(double[] dataInput, ref double[] dataOutput)
{
int n = dataInput.Length;
int k = 0;
int[] v = new int[n];
double[] z = new double[n + 1];
v[0] = 0;
z[0] = Double.NegativeInfinity;
z[1] = Double.PositiveInfinity;
double s;
for (int q = 1; q < n; q++)
{
while (true)
{
s = (((dataInput[q] + q * q) - (dataInput[v[k]] + v[k] * v[k])) / (2.0 * q - 2.0 * v[k]));
if (s <= z[k])
{
k--;
}
else
{
break;
}
}
k++;
v[k] = q;
z[k] = s;
z[k + 1] = Double.PositiveInfinity;
}
k = 0;
for (int q = 0; q < n; q++)
{
while (z[k + 1] < q)
{
k++;
}
dataOutput[q] = ((q - v[k]) * (q - v[k]) + dataInput[v[k]]);
}
}
This can be readily used for binary and grayscale images by applying it first on image columns and then rows (or vice versa, of course).
The transform is indeed very fast.
Here are the source and output images:


The black pixels have value 0 and the white have some large value (have to be larger than largest possible squared distance in the images but not infinity) so that the transform returns distance from the black pixels and the white ones are ommited.
To get true euclidean distance transform, simply take a square root of each pixel from the output image.