Given a convex polygon with integer coordinates vertices, how can you transform it so you can always find a point inside it with integer coordinates

convex-geometrygeometrypolygons

Given a convex polygon, with all its vertices having integer coordinates, you are allowed to perform one transformation T to all of its vertices (which are points), such that:

T transforms a point in another (possibly the same) point.
If P is a point with integer coordinates, T(P) will be a point with integer coordinates.
Let u, v be two consecutive vertices of the original polygon. If a point P was originally on the edge {u, v} of the original polygon, T(P) will be a point that is on the edge {T(u), T(v)} of the transformed polygon.
If a point P was originally in the interior of the original polygon, T(P) will be a point that is in the interior of the transformed polygon.
If a point P was originally outside the original polygon, T(P) will be a point that is outside the transformed polygon.

After applying the transformation T to the polygon, you must find a point Q in the interior the transformed polygon with integer coordinates.

Notice that Q must not be on an edge of the transformed polygon, but strictly in the interior of it.

You may assume that:

We are working in 2D.
No two vertices of the original polygon are the same.
No three consecutive vertices of the original polygon are collinear.

I'm interested in finding a transformation T that satisfies those conditions for all polygons with vertices with integer coordinates, and a way to generate a point with integer coordinates in the interior of the transformed polygon.

Some observations I've come to, but couldn't prove:

It appears that all convex polygons with integer coordinates vertices such that

The maximum Y of any of its points minus the minimum Y of any of its points is 1; or
The maximum X of any of its points minus the minimum X of any of its points is 1

will yield no points with integer coordinates in the interior of them.

One particular example is the triangle formed by {0, 0}, {0, 1}, {1, 0}.

Triangle of the example

For this case, notice that any of the following transformations will still yield a polygon with no integer coordinates points inside of it:
1. (x, y) -> (2x, y)
  Transformation 1
2. (x, y) -> (x, 2y)
  Transformation 2
3. (x, y) -> (2x, 2y)
  Transformation 3
But the following transformations will:
1. (x, y) -> (-3x, -3y)
  Transformation 1A
2. (x, y) -> (6x, 4y)
  Transformation 2A

Best Answer

The transformation $T(x,y) = (3x,3y)$ always works.

Here's why. Let $A = (x_0, y_0)$, $B = (x_1, y_1)$, and $C = (x_2, y_2)$ be three vertices of the polygon. Then their centroid $G = (\frac{x_0+x_1+x_2}{3}, \frac{y_0 + y_1 + y_2}{3})$ is in the interior of $\triangle ABC$, and therefore in the interior of the polygon.

The linear transformation has all the nice properties we want, and in particular it preserves interiors. So we know that $T(G)$ is in the interior of the transformed polygon.

But it's also true that $T(G)$ has integer coordinates: more precisely, the integer coordinates $(x_0 + x_1 + x_2, y_0 + y_1 + y_2)$.

In fact, as long as the polygon is not a triangle, the transformation $T(x,y) = (2x,2y)$ will work. For this, let $A = (x_0, y_0)$ and $B = (x_1, y_1)$ be two vertices of the polygon that are not adjacent. Then the midpoint of $AB$, the point $M = (\frac{x_0 + x_1}{2}, \frac{y_0 + y_1}{2})$, is in the interior of the polygon, and by a similar argument as for the previous transformation, $T(M) = (x_0 + x_1, y_0 + y_1)$ is an integer point inside the transformed polygon.

Related Solutions

How many vertices can the intersection of two convex polygons have

The edges of the new polygons are parts of edges of the old polygons. While each edge of polygon 0 my be parted into up to three line segments by polygon 1, at most one of these line segments can become an edge of the new polygon. The same holds for the edges of polygon 1. Hence polygon 2 cannot have more than $$n_0+n_1$$ edges.

Find closest point to polygon shape from any coordinate

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.

Best Answer

Related Solutions

How many vertices can the intersection of two convex polygons have

Find closest point to polygon shape from any coordinate

Related Question