How to find closest point to polygon shape from any coordinate?

Question

I have a 2D polygon with known coordinates $ABCDEFG\,(x, y)$, and also some points inside and outside the shape, $K(x)$.

I want to figure out a formula, how to calculate the nearest point $C(x)$ and its coordinates $(x, y)$ from $K(x)$ on the shape. The polygon can be different, so the points so the formula needs to be generalized. Sadly I have no idea how to do this, so any help appreciated.

Answer 1

Let's assume that the $N$-sided polygon has vertices $(x_i, y_i)$ in clockwise or counterclockwise order, $i = 1 \dots N$, with $(x_0, y_0) = (x_N, y_N)$, and we wish to find the distance from point $(x_p, y_p)$ to the perimeter of the polygon.

As Théophile mentioned in a comment to the question, the simple (but perhaps not optimal) method is to find the minimum distance between the point and each edge line segment.

The way I prefer to find the distance between the line segment $(x_A, y_A) - (x_B, y_B)$ and point $(x_p, y_p)$, is parametrise the line segment with $0 \le t \le 1$, $$\left\lbrace \begin{aligned} x(t) & = (1 - t) x_A + t x_B \\ y(t) & = (1 - t) y_A + t y_B \\ \end{aligned} \right. \tag{1}\label{G1}$$ The line segment corresponds to $0 \le t \le 1$, so if $t \le 0$, the closest point to $(x_p, y_p)$ is $(x_A, y_A)$; if $t \ge 1$ the closest point to $(x_p, y_p)$ is $(x_B, y_B)$; otherwise, the closest point to $(x_p, y_p)$ is $(x(t), y(t))$.

There are many ways to find $t$ that corresponds to the minimum distance between $(x_p, y_p)$ and the line extending the edge line segment. Suffice it to say, they all yield a single solution, $$t = \frac{ ( x_p - x_A ) ( x_B - x_A ) + ( y_p - y_A ) ( y_B - y_A ) }{ (x_B - x_A)^2 + (y_B - y_A)^2 }$$

(If you precalculate $\vec{m}_i = \left(\frac{x_B - x_A}{(x_B-x_A)^2 + (y_B-y_A)^2} ,\, \frac{y_B - y_A}{(x_B-x_A)^2 + (y_B-y_A)^2}\right)$ for each edge $i$, you can calculate $t_i = (\vec{p} - \vec{v}_i) \cdot \vec{m}_i$ very efficiently, with just two scalar multiplications, two scalar subtractions, and one scalar addition per edge. To calculate the point when $0 \lt t_i \lt 1$, additional four scalar multiplications, two additions, and two subtractions are needed. The distance squared is another two scalar multiplications and an addition. Therefore, this way the worst case cost per edge is eight scalar multiplications and nine additions or multiplications – and that is quite acceptable and efficient, in computer geometry terms.)

Because the distance squared is a monotonic function of the distance, we don't need the square root for each edge line segment; we can find the minimum squared distance instead. This minimum distance $d_i^2$ for the edge line segment ending at vertex $i$ is $$d_i^2 = \begin{cases} (x_{i-1} - x_p)^2 + (y_{i-1} - y_p)^2, & t \le 0 \\ (x(t) - x_p)^2 + (y_(t) - y_p)^2, & 0 \lt t \lt 1 \\ (x_i - x_p)^2 + (y_i - y_p)^2, & 1 \le t \\ \end{cases} \tag{2}\label{G2}$$

Here is a CC0-licensed Python3 example program, that implements a Point class to describe 2D points and vecors, a Polygon class to describe closed 2D polygons, and when run, generates a random 10-sided polygon and ten random points, and saves example.svg, an SVG image you can open in your browser, showing the polygon (in gray on white background), the random points (red dots), and blue lines from each random point to the closest point on the polygon perimeter.

Save the following as example.py, and run it using e.g. python3 example.py:

# SPDX-License-Identifier: CC0-1.0
# -*- coding: utf-8 -*-
from math import sqrt as _sqrt, sin as _sin, cos as _cos, atan2 as _atan2, pi as _pi, inf as _inf

class Point(tuple):
    """A simple 2D point or vector type

       Point(x, y)
       Point(point)
       Point(x=xcoord, y=ycoord)
       Point(r=radius, theta=angle_in_radians)"""

    def __new__(cls, *args, x=None, y=None, r=None, theta=None):
        if x is not None and y is not None:
            return tuple.__new__(cls, (float(x), float(y)))
        if r is not None and theta is not None:
            return tuple.__new__(cls, (float(r)*_cos(float(theta)), float(r)*_sin(float(theta))))
        if len(args) == 1 and isinstance(args[0], (list, tuple, range)) and len(args[0]) == 2:
            return tuple.__new__(cls, (float(args[0][0]), float(args[0][1])))
        if len(args) == 2:
            return tuple.__new__(cls, (float(args[0]), float(args[1])))
        raise ValueError("Invalid parameters to Point().")

    def __init__(self, *args, **kwargs):
        """Points are immutable"""
        pass

    @property
    def x(self):
        """Point.x: Point x coordinate"""
        return self[0]

    @property
    def y(self):
        """Point.y: Point y coordinate"""
        return self[1]

    @property
    def r(self):
        """Point.r: Point distance from origin"""
        return _sqrt(self[0]*self[0] + self[1]*self[1])

    @property
    def theta(self):
        """Point.theta: Angle in radians from point to positive x axis, with respect to origin"""
        return _atan2(self[1], self[0])

    @property
    def ccw(self):
        """Point.ccw: Point rotated 90 degrees counterclockwise around origin"""
        return tuple.__new__(self.__class__, (-self[1], self[0]))

    @property
    def cw(self):
        """Point.cw: Point rotated 90 degrees clockwise around origin"""
        return tuple.__new__(self.__class__, (self[1], -self[0]))

    @property
    def norm(self):
        """Point.norm: Euclidean norm (distance to origin)"""
        return _sqrt(self[0]*self[0] + self[1]*self[1])

    @property
    def normsqr(self):
        """Point.normsqr: Euclidean norm squared"""
        return self[0]*self[0] + self[1]*self[1]

    def dot(self, other):
        """A.dot(B): Dot product between two 2D vectors"""
        if not isinstance(other, Point):
            other = Point(other)
        return self[0]*other[0] + self[1]*other[1]

    def cross(self, other):
        """A.cross(B): The 2D analog of vector cross product A x B"""
        if not isinstance(other, Point):
            other = Point(other)
        return self[0]*other[1] - self[1]*other[0]

    def __bool__(self):
        """Point == True if it is not (0,0)"""
        return (self[0] != 0) or (self[1] != 0)

    def __abs__(self):
        """abs(Point): Euclidean distance from origin"""
        return _sqrt(self[0]*self[0] + self[1]*self[1])

    def __pos__(self):
        """+Point"""
        return tuple.__new__(self.__class__, (self[0], self[1]))

    def __neg__(self):
        """-Point"""
        return tuple.__new__(self.__class__, (-self[0], -self[1]))

    def __add__(self, other):
        """Point + Point"""
        if not isinstance(other, Point):
            other = Point(other)
        return tuple.__new__(self.__class__, (self[0]+other[0], self[1]+other[1]))
    def __radd__(self, other):
        """Point + Point"""
        if not isinstance(other, Point):
            other = Point(other)
        return tuple.__new__(self.__class__, (other[0]+self[0], other[1]+self[1]))
    def __iadd__(self, other):
        """Point += Point"""
        if not isinstance(other, Point):
            other = Point(other)
        return tuple.__new__(self.__class__, (self[0]+other[0], self[1]+other[1]))

    def __sub__(self, other):
        """Point + Point"""
        if not isinstance(other, Point):
            other = Point(other)
        return tuple.__new__(self.__class__, (self[0]-other[0], self[1]-other[1]))
    def __rsub__(self, other):
        """Point - Point"""
        if not isinstance(other, Point):
            other = Point(other)
        return tuple.__new__(self.__class__, (other[0]-self[0], other[1]-self[1]))
    def __isub__(self, other):
        """Point -= Point"""
        if not isinstance(other, Point):
            other = Point(other)
        return tuple.__new__(self.__class__, (self[0]-other[0], self[1]-other[1]))

    def __mul__(self, other):
        """Point * number"""
        if isinstance(other, (int, float)):
            return tuple.__new__(self.__class__, (self[0]*other, self[1]*other))
        return NotImplemented
    def __rmul__(self, other):
        """number * Point"""
        if isinstance(other, (int, float)):
            return tuple.__new__(self.__class__, (other*self[0], other*self[1]))
        return NotImplemented
    def __imul__(self, other):
        """Point *= number"""
        if isinstance(other, (int, float)):
            return tuple.__new__(self.__class__, (self[0]*other, self[1]*other))
        return NotImplemented

    def __truediv__(self, other):
        """Point / number"""
        if isinstance(other, (int, float)):
            return tuple.__new__(self.__class__, (self[0]/other, self[1]/other))
        return NotImplemented
    def __rtruediv__(self, other):
        """number / Point is invalid"""
        return NotImplemented
    def __itruediv__(self, other):
        """Point /= number"""
        if isinstance(other, (int, float)):
            return tuple.__new__(self.__class__, (self[0]/other, self[1]/other))
        return NotImplemented


class Polygon:
    """A 2D polygon class"""

    __slots__ = ['_points', '_iedges']

    def __init__(self, points=[]):
        self._points = []
        self._iedges = None
        self.append(*points)

    def __len__(self):
        """len(Polygon): Number of points in polygon"""
        return len(self._points)

    def __getitem__(self, i):
        """Polygon[i]: Vertex i of polygon as a Point"""
        if isinstance(i, int):
            return self._points[i]
        else:
            raise TypeError("Vertex numbers must be integers, not %s" % str(type(i)))

    def __iter__(self):
        """Iterator over polygon points"""
        for p in self._points:
            yield p

    def aabb(self):
        """Returns the diagonal of the axis-aligned bounding box for the polygon"""
        if len(self._points) < 1:
            return None, None

        minX, minY = _inf, _inf
        maxX, maxY =-_inf,-_inf
        for p in self._points:
            minX, maxX = min(minX, p.x), max(p.x, maxX)
            minY, maxY = min(minY, p.y), max(p.y, maxY)

        return Point(minX,minY), Point(maxX,maxY)

    def append(self, *points):
        """Add one or more vertices to the polygon"""
        for p in points:
            if isinstance(p, Point):
                self._points.append(p)
            else:
                self._points.append(Point(p))
        if len(points) > 0:
            self._iedges = None

    def precalculate(self):
        """Precalculate iedge vectors, for faster stable calculations"""
        n = len(self._points)
        if n > 1:
            self._iedges = []
            pB = self._points[-1]
            for i in range(0, n):
                pA = pB
                pB = self._points[i]
                pAB = pB - pA
                DD = pAB.normsqr
                if DD > 0:
                    self._iedges.append(pAB/DD)
                else:
                    self._iedges.append(Point(0,0))
        elif n == 1:
            self._iedges = [ Point(0,0) ]
        else:
            self._iedges = None

    def perimeter_nearest_to(self, p):
        """Returns the point on the perimeter of the polygon closest to the specified point."""
        if not isinstance(p, Point):
            p = Point(p)

        n = len(self._points)
        if n > 1:
            nearestNormsqr = _inf
            nearest = None

            # Make sure the iedge vectors have been precalculated.
            if self._iedges is None:
                self.precalculate()

            pB = self._points[-1]
            for i in range(0, n):
                pA = pB
                pB = self._points[i]

                t = (p - pA).dot(self._iedges[i])
                if t <= 0:
                    q = pA
                elif t < 1:
                    q = (1-t)*pA + t*pB
                else:
                    q = pB

                qq = (q - p).normsqr
                if qq < nearestNormsqr:
                    nearest = q
                    nearestNormsqr = qq

        elif n == 1:
            nearest = self._points[0]
            nearestNormsqr = (nearest - p).normsqr

        else:
            nearest = None
            nearestNormsqr = _inf

        return nearest


if __name__ == '__main__':
    from random import Random
    uniform = Random().uniform

    # For simplicity, the output is an SVG image with viewBox (0,0) - (1000, 1000)
    poly = Polygon([ Point(uniform(10,990),uniform(10,990)) for i in range(0, 10) ])

    with open("example.svg", "w", encoding="utf-8") as out:
        out.write("<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n")
        out.write("<svg xmlns=\"http://www.w3.org/2000/svg\" viewBox=\"0 0 1000 1000\">\n")

        # Fill the background with white.
        out.write("<rect x=\"0\" y=\"0\" width=\"1000\" height=\"1000\" stroke=\"none\" fill=\"#ffffff\" />\n")

        # Fill the polygon background with light gray.
        out.write("<path stroke=\"none\" fill=\"#cccccc\" fill-rule=\"nonzero\" d=\"M")
        for p in poly:
            out.write(" %.3f,%.3f" % p)
        out.write(" z\"/>\n")

        # Experiment with ten random points.
        points = [ Point(uniform(5,995),uniform(5,995)) for i in range(0, 10) ]

        # Draw red dots at each point.
        for p in points:
            out.write("<circle cx=\"%.3f\" cy=\"%.3f\" r=\"5\" stroke=\"#000000\" fill=\"#ff0000\" />\n" % p)

        # Draw a blue line from each point to the closest point on the perimeter.
        out.write("<path stroke=\"#0000ff\" fill=\"none\" d=\"")
        for p in points:
            n = poly.perimeter_nearest_to(p)
            out.write("M%.3f,%.3f L%.3f,%.3f " % (*p, *n))
        out.write("\" />\n")

        # Draw the polygon perimeter with thin black lines.
        out.write("<path stroke=\"#000000\" fill=\"none\" d=\"M")
        for p in poly:
            out.write(" %.3f,%.3f" % p)
        out.write(" z\"/>\n")

        out.write("</svg>\n")

The classes and methods have descriptions, so you can run pydoc3 example (after saving the above example.py) to see those.

Answer 2

I've taken the time to translate @Glärbo's fantastic answer into C# using Unity Vector2 for the points. Note that the original solution has a bunch of helper functions that are not used in the specific answer; I haven't translated these across.

using System.Collections.Generic;
using System.Linq;
using UnityEngine;

public class ClosePointPolygon {
  private Vector2[] points;
  private Vector2[] iEdges;
  private float nearestNormSqr = float.MaxValue;

  public ClosePointPolygon(Vector2[] points) {
    this.points = points;
  }

  public void Precalculate() {
    if (points.Length == 0) iEdges = new Vector2[0];
    else if (points.Length == 1) iEdges = new Vector2[1] { new Vector2(0, 0) };
    else {
      List<Vector2> list = new List<Vector2>();
      Vector2 pB = points.Last();
      for (int i = 0; i < points.Length; i++) {
        Vector2 pA = pB;
        pB = points[i];
        Vector2 pAB = pB - pA;
        float DD = pAB.NormSqr();
        if (DD > 0f)
          list.Add(pAB / DD);
        else
          list.Add(Vector2.zero);
      }
      iEdges = list.ToArray();
    }
  }

  public Vector2 NearestPointFrom(Vector2 p) {
    int len = points.Length;
    Vector2 nearest = Vector2.positiveInfinity;
    if (len > 1) {
      nearestNormSqr = float.MaxValue;
      nearest = Vector2.positiveInfinity;

      if (iEdges == null) Precalculate();

      Vector2 pB = points.Last();
      for (int i = 0; i < len; i++) {
        Vector2 pA = pB;
        pB = points[i];

        Vector2 q = Vector2.positiveInfinity;
        float t = Vector2.Dot(p - pA, iEdges[i]);
        if (t <= 0f)
          q = pA;
        else if (t < 1f)
          q = (1f - t) * pA + t * pB;
        else
          q = pB;

        float qq = (q - p).NormSqr();
        if (qq < nearestNormSqr) {
          nearest = q;
          nearestNormSqr = qq;
        }
      }
    } else if (len == 1) {
      nearest = points[0];
      nearestNormSqr = (nearest - p).NormSqr();
    } else {
      Debug.LogError("NearestPointFrom points is empty");
      nearestNormSqr = float.MaxValue;
    }
    return nearest;
  }
}

public static class Extensions {
  public static float NormSqr(this Vector2 v) {
    return (v.x * v.x) + (v.y * v.y);
  }
}

Example showing points (red) connected to closest points (black):