How do you find the distance distribution of a polygon?

Question

As part of something I'm coding for work, I need to find the distribution of distance from a point to a polygonal object (described by a set of coordinates $(x_1,y_1)...(x_n,y_n)$). In other words, a function/algorithm/process for "How much of the object is at a given distance?"

After a bit of scribbling, I haven't been able to find a satisfactory solution, and if this is a known technique elsewhere, it is named something that I have not thought to search for.

Answer 1

Let $f(x,y)$ be the mass distribution of the object (the density of the object at point $x,y$).

Change the coordinates to polar coordinates: $$g(\rho,\theta) = \rho f(\rho \cos\theta, \rho \sin\theta).$$

Then the function $g(\rho,\theta)$ tells you the amount of mass of the object at distance $\rho$ and angle $\theta$.

Finally, integrate over $\theta$ to find the density of rhe object ar distance $\rho$: $$m(\rho) =\int_0^{2\pi} g(\rho,\theta) \mathrm d \theta.$$

If the density if uniform, then you are effectively looking for the lengths of the arcs of circles of different radii that fall inside the polygon.

One algorithm to do this easily would be

1. lay out all the points in a $\rho, \theta$ plane, so that vertical lines have constant distance from the origin.
2. Select a fine enough grid of points along the $\rho$ dimension and draw vertical lines at these points.
3. Find the intersection of the vertical lines with the sides of the polygon.
4. Sum the distances of all the pairs of points that correspond to lines inside the polygon.

Here is the result of running the algorithm

and here is the code to generate the picture

import matplotlib.pyplot as plt
import numpy as np


def getpoly():
    """Create polynomial as set of vertices"""
    x = [1.2, 1.4, 1.35, 1.45, 1.6, 1.7, 1.4, 1.6, 1.3]
    y = [0.7, 0.8, 0.83, 0.85, 0.95, 1.1, 0.95, 1.1, 1.0]
    x.append(x[0])
    y.append(y[0])
    return x, y


def plotpoly(ax, xs, ys):
    """Plot the polynomial"""
    for i in range(1, len(xs)):
        ax.plot(xs[i-1:i+1], ys[i-1:i+1], color='C0')
        xm = 0.5*(xs[i-1] + xs[i])
        ym = 0.5*(ys[i-1] + ys[i])
        ax.text(xm, ym, f"{i}")
    ax.set_xlim(0, max(xs))
    ax.set_ylim(0, max(ys))
    ax.set_xlabel("x")
    ax.set_ylabel("y")


def isbetween(x, a, b):
    if a < b:
        return a <= x <= b
    else:
        return b <= x <= a


def intersectpoly(q, xs, ys):
    """Find intersections of vertical line at q with polygon"""
    p = []
    for i in range(len(xs)-1):
        if isbetween(q, xs[i], xs[i+1]):
            y = ys[i] + (ys[i+1] - ys[i])/(xs[i+1] - xs[i]) * (q - xs[i])
            p.append(y)
    return p


def getpolar(xs, ys):
    xs = np.asarray(xs)
    ys = np.asarray(ys)
    return np.sqrt(xs**2+ ys**2), np.arctan2(ys, xs)


def sumarcsegments(r, xs, ys):
    p = intersectpoly(r, xs, ys)
    s = 0.0
    for i in range(0, len(p), 2):
        s += r*(p[i+1] - p[i])
    return np.abs(s)


if __name__ == "__main__":
    x, y = getpoly()
    r, th = getpolar(x, y)

    fig, (ax1, ax2) = plt.subplots(1, 2)

    plotpoly(ax1, x, y)

    rs = np.linspace(min(r), max(r), 100)
    ds = np.zeros(100)
    for i, rho in enumerate(rs):
        ds[i] = sumarcsegments(rho, r, th)

    ax2.plot(rs, ds)
    ax2.set_xlim(0.0, max(rs))
    ax2.set_xlabel("distance from origin")
    ax2.set_ylabel("density")
    plt.show()