Finding the largest area axis-parallel square with a known center in a polygon

Finding the maximum volume shape that will fit inside of another shape is a common class of geometric optimization problem that occurs in many applications (computer graphics, collision detection, GIS, etc.). General problem formulations, such as finding the largest area axis-parallel rectangle in a polygon, involve complex analysis and do not have simple algorithmic solutions that are easy to implement.

Finding the largest area axis-parallel square with a known center in a polygon is a sufficiently restricted problem to admit a simple $O(n)$ time and $O(1)$ space solution. The analysis, algorithm, and implementation are very straightforward. I have found this problem formulation useful in computing cropping regions for eye images stenciled by hidden area meshes in virtual reality systems, but it should be general enough to be useful in other contexts.

Previous Work

Finding the largest axis-parallel square with a known center is much easier than the more general problem of finding the largest area axis-parallel rectangle in a polygon. Daniels et al. provide an $O(n\log^2n)$ algorithm for finding the largest area axis-parallel rectangle in a general polygon and a $\Omega(n\log n)$ lower bound for both self-intersecting polygons and general polygons [2]. DePano et al. give algorithms for computing largest squares and equilateral triangles in convex polygons [3]. Alt et al. give $\Theta(\log n)$ algorithm for computing the largest inscribed isothetic rectangle in a convex polygon [1].

Preliminaries

The general form of a line equation is written as $ax + by + c=0$ . Given two points $(x_1, y_1)$ and $(x_2, y_2)$ , we can compute coefficients $a$ , $b$ , and $c$ as follows. First, compute the slope as rise over run $m = (y_2 - y_1)/(x_2 - x_1)$ . Second, plug this value $m$ and one of the points into the point-slope form, given by $y - y_1 = m(x - x_1)$ , and gather all terms on the left-hand side. This reveals that $a = y_1 - y_2$ , $b = x_2 - x_1$ , and $c = y_2x_1 - y_1x_2$ .

The general form of a line equation can be used as an orientation test by plugging in the coordinates of a point and observing the sign. Negative values indicate the point is to the left of the line through the two points, 0 indicates that the point is on the line, and positive values indicate the point is to the right of the line. We use this fact to structure our line intersection code similar to Prasad [4], using the coefficients to subsequently compute intersection points.

For full derivations and interactive demos, see Computing General Form Line Coefficients From Two Points and Computing the Intersection of Two Lines.

Algorithm

Given a polygon $P$ and point $c$ , we compute a single value: the size $s$ , or halfwidth of a square $S$ contained in $P$ and centered at $c$ . We observe that as $S$ expands, its expansion will be stopped either by one of its sides colliding with a vertex of $P$ , or by one of its corners colliding with an edge of $P$ . Thus, for each edge $p_ip_{i+1}$ of $P$ , we can see if either $p_i$ impedes the relevant side of $S$ , or if $p_ip_{i+1}$ impedes the expansion a corner of $S$ . By doing this for all edges, we can compute the maximum value of $s$ .

Our algorithm maintains the loop invariant that $s$ is as large as possible, given the subset of $P$ that we have examined so far. We initialize $s = \infty$ , satisfying the invariant on line 5. We maintain the invariant each time through the loop by examining all possible values that could reduce the value of $s$ : the distance between $p_i$ and $c$ along the relevant axis, and the intersection of $p_ip_j$ with two diagonal lines passing through $c$ . At termination, i == P.size(), implying that every edge of $P$ has been examined, and $s$ is as large as possible with respect to all of $P$ .

Expanding a square. Two cases stop the expansion, side collision and corner collision. We account for side collisions by examining the distance of $P$ ‘s vertices to the center $c$ . We account for corner collisions by computing intersections between two diagonal lines and each edge of $P$ .

void MaxSquare(const std::vector<Vector2>& P, const Vector2& c, float& s) {
    if (P.size() < 2) return;

    // Descend over all vertices/edges to find maximum halfwidth.
    s = std::numeric_limits<float>::infinity();
    size_t i = 0;
    size_t j = 1;
    Vector2 pi{ P[i].x - c.x, P[i].y - c.y };
    while (i < P.size()) {
        // Center input on c so computations have more precision.
        const Vector2 pj{ P[j].x - c.x, P[j].y - c.y };

        // For N and S quadrants, minimize on Y, for E and W, minimize on X.
        const float absiX = fabsf(pi.x);
        const float absiY = fabsf(pi.y);
        s = std::min(absiX < absiY ? absiY : absiX, s);

        // Compute coefficients for line through I and J.
        const float lIJA = pj.y - pi.y;
        const float lIJB = pi.x - pj.x;
        const float lIJC = pj.x * pi.y - pi.x * pj.y;

        // Compute denominator values for each line pair (diagonals x edge IJ).
        const float denIJ045 = lIJB + lIJA;
        const float denIJ135 = lIJB - lIJA;

        // If edge IJ straddles a diagonal, consider x coord of intersection.
        // (No need to look at both coords since x ~= y on a diagonal.)
        if ((pi.x < pi.y != pj.x < pj.y) && denIJ045 != 0)
            s = std::min(fabsf(-lIJC / denIJ045), s);

        if ((pi.x < -pi.y != pj.x < -pj.y) && denIJ135 != 0)
            s = std::min(fabsf(lIJC / denIJ135), s);

        pi = pj;
        ++i;
        j = (j + 1) % P.size();
    }
}

Since we examine each edge only once and perform a constant amount of work per edge, our running time is $O(n)$ . We compute a temporary, constant amount of information each time through the loop, using $O(1)$ space.

Edge Cases and Robustness

We require that all input polygons have at least three vertices on line 2, although the code can be easily adapted to handle degenerate polygons. While we are concerned with computing the maximum square inside a polygon, we note that our algorithm will compute a useful square even if the specified center lies outside of $P$ . Our algorithm also handles polygons with edge crossings.

We center $P$ on $c$ on lines 8 and 11 in order to increase the precision of subsequent calculations and to simplify many expressions, since this zeroes out the $C$ coefficient for both diagonal line equations. For polygon vertex coordinates within a factor of 2 of $c$ this subtraction is exact [5]. However, all branching decisions in the main loop may be affected by floating-point roundoff error in their operands. This can result in two problems: misclassifying $p_i$ into the wrong quadrant or a diagonal intersecting the wrong edge.

First, we may encounter errors when choosing which axis we will minimize along on line 16. However, incorrect branching behavior will only occur when vertex $p_i$ is nearly on one of the diagonals. In these cases, $p_i.x$ is approximately equal to $p_i.y$ , so choosing the wrong branch will have minimal impact.

Second, we may encounter errors when choosing whether edge $p_ip_j$ intersects either diagonal on lines 29 and 32. However, we compute orientation values for each vertex only once, so if a diagonal misses an intersection with an edge, it will instead intersect the neighboring edge. Although it is possible that intersecting the wrong edge will result in a significantly different intersection point, our algorithm avoids entirely missing intersections as is possible through naively shooting rays at each edge.

References

Helmut Alt, David Hsu, and Jack Snoeyink. Computing the largest inscribed isothetic rectangle. In CCCG, pages 67-72. Citeseer, 1995.
Karen Daniels, Victor Milenkovic, and Dan Roth. Finding the largest area axis-parallel rectangle in a polygon. Computational Geometry, 7(1-2):125-148, 1997.
A DePano, Yan Ke, and J O’Rourke. Finding largest inscribed equilateral triangles and squares. In Proc. 25th Allerton Conf. Commun. Control Comput, pages 869-878, 1987.
Mukesh Prasad. Intersection of line segments. In Graphics Gems II, pages 7-9. Elsevier, 1991.
Pat H Sterbenz. Floating-point computation. page 138, 1974.

Clinton Freeman

Computer Science, Design