#!/usr/bin/python3
"""Test a linear-time discrete approximate SDF (signed distance function).

This algorithm approximates the Euclidean distance transform (for
either the outside or inside of a set of pixels) with a worst-case
error of ±0.65% in 24 add-and-pairwise-min operations per pixel,
producing level sets that are the Minkowski sum of the seed set with
an irregular icositetragon.  This is a confirmation of previously
calculated results.

Here we are only testing it on a single-pixel-set seed, relying on the
structure of the algorithm to claim that it works for all possible
subsets of pixels.

Faster, cruder versions (like ±4.1% with 8 operations per pixel) and
slower, more precise versions are also possible.

This implementation does many passes over the data and has inner loops
in interpreted Python, so it is going to be very slow.

"""
import numpy


def assess(test, cx, cy):
    "Measure worst-case error of an approximation of Euclidean distance from (cx, cy)."
    xsize, ysize = test.shape
    y, x = numpy.meshgrid(range(ysize), range(xsize)) # nasty confusing Matlab order
    assert y.shape == x.shape == test.shape
    euclidean = ((x - cx)**2 + (y - cy)**2)**(1/2)
    assert euclidean.shape == test.shape
    mask = (euclidean != 0)
    error = test[mask] / euclidean[mask]
    return error.min(), error.max()


def seed(cx, cy, xsize, ysize):
    "Initialize a conservative approximation of distance from (cx, cy)."
    vals = numpy.full((xsize, ysize), numpy.inf)
    vals[cx, cy] = 0
    return vals


def downright(arr, dx, dy):
    "Apply a recursive IIR dilation down and right."
    assert dy >= 0
    assert dy > 0 or dx >= 0
    dist = (dx**2 + dy**2)**(1/2)

    xsize, ysize = arr.shape
    xrange = range(dx, xsize) if dx >= 0 else range(0, xsize + dx)
    for y in range(dy, ysize):
        for x in xrange:
            arr[x, y] = min(arr[x, y], arr[x - dx, y - dy] + dist)


def left(arr, dx):
    "Apply a recursive IIR dilation left."
    assert dx > 0
    xsize, _ = arr.shape
    for x in range(xsize - dx - 1, -1, -1):
        arr[x, ...] = numpy.minimum(arr[x, ...], arr[x + dx, ...] + dx)


def up(arr, dx, dy):
    "Apply a recursive IIR dilation up."
    assert dy > 0
    dist = (dx**2 + dy**2)**(1/2)

    xsize, ysize = arr.shape
    xrange = range(dx, xsize) if dx >= 0 else range(0, xsize + dx)
    for y in range(ysize - dy - 1, -1, -1):
        for x in xrange:
            arr[x, y] = min(arr[x, y], arr[x - dx, y + dy] + dist)


def main():
    print("Initial distance function drawn from {0, ∞}:")
    cx, cy = 3, 4
    a = seed(cx=cx, cy=cy, xsize=24, ysize=10)
    print(a)

    print()
    print("Four cartesian-direction (4-neighbor) recursive relaxation passes")
    print("compute exact Manhattan distances.")
    downright(a, 1, 0)
    downright(a, 0, 1)
    left(a, 1)
    up(a, 0, 1)
    print(a)

    print()
    print("These are off by at most √2 from Euclidean distance.")
    print(assess(a, cx=cx, cy=cy))

    print()
    print("Diagonal (8-neighbor) passes compute octagonal approximation.")
    downright(a, 1, 1)
    downright(a, -1, 1)
    up(a, -1, 1)
    up(a, 1, 1)
    print(a.round(2))

    print("Error is reduced to [-0, +8.2%]:")
    print(assess(a, cx=cx, cy=cy))

    print()
    print("Expanding neighborhood to include (1, 2) and (1, 3) in each octant")
    print("reduces error to 1.3% (irregular icositetragon approximation):")
    for n in [2, 3]:
        for direction in [downright, up]:
            direction(a, n, 1)
            direction(a, -n, 1)
            direction(a, 1, n)
            direction(a, -1, n)

    print(a.round(2))

    print(assess(a, cx=cx, cy=cy))

    print()
    print("Scaling down estimates by 0.65% reduces error to ±0.65%:")
    a /= 1.0130722**(1/2)
    print(a.round(3))
    print(assess(a, cx=cx, cy=cy))
    

if __name__ == '__main__':
    main()