Ray to irregular concave poly intersection

Question

Is it possible to take a 2d vector representing position and find the closest point of intersection to an irregular polygon with concave sides? For example this shape.

I think I could do it if I had a point to line function, but the only one I have seems to turn the line into a ray, and all I get is the perpendicular intersection of the ray.

Answer 1

The solution is obviously the minimum distance between the point and the line segments corresponding to the sides of the polygon.

If $d( x_0 , y_0 , x_1 , y_1 , x_p , y_p )$ is a function that yields the distance between point $( x_p , y_p )$ and the line segment from $( x_0 , y_0 )$ to $( x_1 , y_1 )$, and the polygon has $N$ vertices, numbered from $1$ to $N$, then the distance between the point $( x_p , y_p )$ and the polygon is $$\begin{aligned} \min\bigr( & d( x_1 , y_1 , x_2 , y_2 , x_p , y_p ), \\ \, & d( x_2 , y_2 , x_3 , y_3 , x_p , y_p ), \\ \, & \dots \\ \, & d( x_{N-1} , y_{N-1} , x_N , y_N , x_p , y_p ) \\ \, & d( x_N , y_N , x_1 , y_1 , x_p , y_p ) \bigr) \end{aligned}$$

Note the line segment that corresponds to the last polygon edge, between the last vertex $(x_N , y_N )$ and first vertex $(x_1 , y_1 )$.

---

There are several questions and answers regarding the distance between a point and a line segment here, but for the purposes of education, I shall explain how one might go about discovering it on their own.

First, let's parametrise our line segment. Let the two endpoints (vertices in the polygon) be $(x_0 , y_0)$ and $(x_1 , y_1)$, and $0 \le t \le 1$ be the parameter: $$\left\lbrace\begin{aligned} x(t) &= x_0 + t ( x_1 - x_0 ) = (1 - t) x_0 + t x_1 \\ y(t) &= y_0 + t ( y_1 - y_0 ) = (1 - t) y_0 + t y_1 \end{aligned} \right. \tag{1}\label{NA1}$$ Although the two forms are mathematically equivalent, the rightmost form is most numerically stable when using floating-point math. (The difference is observable when the difference in magnitudes of the coordinates is large; i.e., when one coordinate is very close to zero, and the other coordinate is huge.)

Calculating distance involves a square root. Because square root is a monotonically increasing function, instead of distance, we can use the distance squared to find the minimum distance.

The distance squared between a point on the line segment corresponding to $t$, and point $(x_p, y_p)$, is $$s(t) = \left( x(t) - x_p \right)^2 + \left( y(t) - y_p \right)^2 \tag{2}\label{NA2}$$

When the distance and distance squared reach the minimum (or maximum), their derivative is zero. Therefore, the point $t$ on the line that is closest to $(x_p , y_p)$ is $$\frac{d\,s(t)}{d\,t} = 0 \tag{3}\label{NA3}$$ If we solve $\eqref{NA3}$ for $t$, we find $$t = \frac{ ( x_p - x_0 ) ( x_1 - x_0 ) + ( y_p - y_0 )( y_1 - y_0 ) }{ ( x_1 - x_0 )^2 + ( y_1 - y_0 )^2 } \tag{4}\label{NA4}$$

Remember that we parametrised the line segment with $0 \le t \le 1$. If $\eqref{NA4}$ says $t \lt 0$, then the closest point between the line and the point is on the line extending from $( x_0 , y_0 )$ to infinity, and $(x_0 , y_0)$ is the closest point on the line segment to $(x_p , y_p)$. Correspondingly, if $t \gt 1$, then the closest point is on the line extending from $( x_1 , y_1 )$ to infinity, and $(x_1 , y_1)$ is the closest point on the line segment to $(x_p, y_p)$.

If we define the distance between two points as $$L(x_a , y_a , x_b , y_b) = \sqrt{ (x_b - x_a)^2 + (y_b - y_a)^2} \tag{5}\label{NA5}$$ then we can define the minimum distance between point $(x_p , y_p)$ and line segment from $(x_0 , y_0)$ to $(x_1 , y_1)$ as follows:

First, we compute $t$ as in $\eqref{NA4}$: $$t = \frac{ ( x_p - x_0 ) ( x_1 - x_0 ) + ( y_p - y_0 )( y_1 - y_0 ) }{ ( x_1 - x_0 )^2 + ( y_1 - y_0 )^2 }$$ Then, we compute the minimum distance depending on the value of $t$, using $\eqref{NA5}$: $$d(x_0 , y_0 , x_1 , y_1) = \begin{cases} L(x_0, y_0, x_p, y_p), & t \lt 0 \\ L(x(t), y(t), x_p, y_p), & 0 \le t \le 1 \\ L(x_1, y_1, x_p, y_p), & t \gt 1 \end{cases}$$

The overall cost for computing the distance is one square root, one division, six multiplications, and seven additions or subtractions when $t \lt 0$ or $t \gt 1$ (and one uses temporary variables); and one square root, one division, ten multiplications, and ten additions or substractions otherwise (still using temporary variables). In other words, it is not as lightweight as calculating the distance between a point and a line, but only a few additions, subtractions, and multiplications more — which is quite acceptable.

---

In programming, you can optimize the calculation, seeing that in $\eqref{NA4}$, the denominator is never negative. For example, in pseudocode, one might write the function returning the minimum distance between point (xp,yp) and the line segment from (x0,y0) to (x1,y1) as

Function LineSegPointDistance(x0, y0, x1, y1, xp, yp):
    Let x01 = x1 - x0
    Let y01 = y1 - y0
    Let x0p = xp - x0
    Let y0p = yp - y0
    Let td = x01*x01 + y01*y01

    # Is the line segment degenerate, a point only?
    If td == 0, Then:
        Return sqrt(x0p*x0p + y0p*y0p)
    End If

    Let tn = x0p*x01 + y0p*y01

    If tn <= 0, Then:
        # t <= 0, so closest point to (xp, yp) is (x0, y0).
        Return Sqrt(x0p*x0p + y0p*y0p)
    End If

    If tn >= td, Then:
        # t >= 1, so closest point to (xp, yp) is (x1, y1).
        Let x1p = xp - x1
        Let y1p = yp - y1
        Return Sqrt(x1p*x1p + y1p*y1p)
    End If

    Let t1 = tn / td
    Let t0 = 1 - t1
    Let xd = t0*x0 + t1*x1 - xp
    Let yd = t0*y0 + t1*y1 - yp
    Return Sqrt(xd*xd + yd*yd)
End Function

Note that the above pseudocode function also handles degenerate line segments (where x0 == x1 and y0 == y1) correctly.

(This also happens to be one of the rare cases where floating-point comparison (of tn) to exact zero is okay. Personally, I like to write it as tn <= 0.0, but that's just habit; a reminder to myself and other programmers that there are no negative tn cases to consider and even if they did occur, they could be folded into this case. The reason why the comparison is okay, is because we only calculate tn/td when we have already verified that 0 <= tn/td <= 1; if td were to be very small, we already known tn is of similar or smaller magnitude, and there should be no problem in calculating the division, even with limited-precision floating-point numbers.)

The pseudocode to calculate the distance between point (xp, yp) and a polygon with vertices x[N], y[N], with N > 1, and indexing the vertices from 1 to N, inclusive, can be something like follows:

Function LineSegPointDistanceSquared(x0, y0, x1, y1, xp, yp):
    Let x01 = x1 - x0
    Let y01 = y1 - y0
    Let x0p = xp - x0
    Let y0p = yp - y0
    Let td = x01*x01 + y01*y01

    # Is the line segment degenerate, a point only?
    If td == 0, Then:
        Return (x0p*x0p + y0p*y0p)
    End If

    Let tn = x0p*x01 + y0p*y01

    If tn <= 0, Then:
        # t <= 0, so closest point to (xp, yp) is (x0, y0).
        Return (x0p*x0p + y0p*y0p)
    End If

    If tn >= td, Then:
        # t >= 1, so closest point to (xp, yp) is (x1, y1).
        Let x1p = xp - x1
        Let y1p = yp - y1
        Return (x1p*x1p + y1p*y1p)
    End If

    Let t1 = tn / td
    Let t0 = 1 - t1
    Let xd = t0*x0 + t1*x1 - xp
    Let yd = t0*y0 + t1*y1 - yp
    Return (xd*xd + yd*yd)
End Function

Function PointToPolygonEdgeDistance(x[], y[], N, xp, yp):
    Let minSquared = LineSegPointDistanceSquared(x[N], y[N], x[1], y[1], xp, yp)
    For i = 1 to N-1:
        Let squared = LineSegPointDistanceSquared(x[i], y[i], x[i+1], y[i+1], xp, yp)
        If squared < minSquared, Then:
            Let minSquared = squared
        End If
    End For
    return Sqrt(minSquared)
End Function

The above again utilises the fact that square root and square are monotonically increasing functions (for all nonnegative arguments). In the loop checking each edge of the polygon, we can compare distances squared to find the minimum, and only do one final square root to calculate the actual distance. (Square root is a pretty slow floating-point operation compared to e.g. multiplications and summation or subtraction, so avoiding doing it for each line segment gives a nice speedup.)

Do note that this approach yields the minimum distance, and not signed distance. That means you cannot use this function to determine if a point is inside or outside the polygon. This approach only returns the distance to the nearest edge.