The partitioning problem

I recently completed a computational geometry related project. Its focus was on orthogonal polygons, which are polygons with edges that intersect at right angles. The project is called polygon-partitioner, and it’s primary task is breaking up orthogonal polygons into non-overlapping rectangles. It’s secondary task is approximating arbitrary polygons by orthogonal polygons. It was really challenging to find algorithms that do these things and I had a lot of fun doing it.

Here is a picture of what the output looks like. The input to the partitioning function is the exterior (in red), and the output is a list of rectangles.

I can’t take complete credit for the algorithm that partitions a polygon. The overall algorithm was (probably) first introduced in [Ferrari et. al, 1984]. The author of this blog post writes out the high level algorithm.

Motivation

Apparently, partitioning polygons is useful for VLSI design problems. VLSI, or Very Large Scale Integration, is used for designing computer chips. That’s probably how this problem became interesting to researchers in the days of yore, but I don’t care about all that¹. I was just having fun.

When I went about implementing the algorithm from [Ferrari et. al, 1984], there were a number of issues that came up that didn’t seem to have obvious solutions. I also couldn’t find much related to extracting the rectangles once you have successfully drawn the lines that make them (that is, once you’ve drawn the gray lines from the above image, return a list of rectangles). I’m writing this post to lay out my solutions. If these things are discussed elsewhere feel free to reach out and let me know.

A high level partitioning algorithm

Before we write down the high level algorithm, we need to define a few terms related to vertices (or corner points) on the boundary of our polygon.

Some definitions

• A convex vertex is a point on the boundary that makes a 90 degree angle with the interior of the polygon. The corner points of rectangles are always convex vertices.
• A concave vertex is a point that makes a 270 degree angle with the interior of the polygon.
• A vertical (horizontal) interior line is any vertical (horizontal) line-segment that lies within the polygon and has at least one endpoint that is a concave vertex. The other endpoint of this line-segment doesn’t need to be a vertex on the boundary. Each gray line in the above figure is an interior line.
• And lastly, a chord is an interior line whose endpoints are both concave vertices.

Note that all angles are measured counter-clockwise relative to the positive x-axis.

An outline of the algorithm

The main challenge with partitioning orthogonal polygons is the presence of chords. The overall algorithm for partitioning a chord-free polygon is nicely laid out in this blog post by Aaron Sterling. That post also talks about what to do to handle chords, and so does this stackoverflow answer. The algorithm I implemented handles orthogonal polygons with chords but its overall structure is similar to the chord-free case. Here is the algorithm for partitioning chord-free polygons:

Input: An orthogonal polygon $P$ without any chords.

1. Setup; extracting line-segments. Let $E$ be the set of line-segments along the boundary of $P$.
2. Create interior line. For each concave vertex $v_i$:
   1. Create a vertical (or horizontal) interior line that starts at $v_i$ and ends at the closest horizontal (vertical) line segment in $E$. Call this line segment $l_i$.
   2. Add $l_i$ to $E$.
3. Extract rectangles. The line-segments in $E$ (implicitly) partition the polygon. Use them to find the non-overlapping rectangles.

I will go through each step in detail for the remainder of this post.

Step 1: Extracting line-segments

Let $E$ be the set of line-segments along the boundary of $P$.

Although this seems self explanatory, you can do a little extra work here to make your life easier during the later steps. In particular, you can identify which vertices along the boundary are concave ones, and choose how you want to create your interior lines when the time comes.

After this step, each vertex should have two pieces of data associated with it: a tag that specifies whether the vertex is concave, and an angle (that is either 0, 90, 180, or 270) which specifies the “direction of the point”.

Preprocessing the input

First make sure the input polygon doesn’t have any extraneous vertices along the boundary. These are points that can be removed without changing the shape of the polygon. For example, if you are given a square then this square should only have 4 points. If it has more than 4 points then you know that at least one point can be removed without changing the shape. I accomplished this by looking at triplets of adjacent vertices and making sure that all three vertices don’t have the same x-coordinate, or have the same y-coordinate.

Identifying concave vertices

We will cycle through the vertices along the boundary in the clockwise direction. For each vertex $v_i$, we also consider the vertices immediately before and after $v_i$. Look at the angle these three points make relative to the interior of the polygon, it will be either 270$^\circ$ or 90$^\circ$². If it is 90 then $v_i$ is convex, otherwise $v_i$ is concave.

Setup for interior lines

When making interior lines, we need to choose a direction to “extend” concave vertices. To do so, look at the angle that the vector ($v_{i-1}, v_i$) makes relative to the x-axis, and use this angle as the direction to extend concave vertices. Because this is an orthogonal polygon, this angle will be either 0, 90, 180, or 270. Use this angle to help you decide, when the time comes, which direction to extend $v_i$ when creating the interior lines.

Time complexity

Calculating these angles takes a constant amount of computational time per vertex, so this step is $\Theta(n)$ time complexity, where $n$ is the number of vertices in the polygon.

An example

Let’s take the polygon below as an example. We will be traversing the boundary clockwise.

When trying to figure out whether $v_1$ is concave or not, we need to look at the angle the line path $(v_6, v_1, v_2)$ makes relative to the interior of the polygon. Turns out it is 90$^\circ$ so $v_1$ is convex. When calculating its “direction”, we look at the angle the vector $(v_6, v_1)$ makes relative to the x-axis. This is the same thing as saying “what direction are we heading in if we start at $v_6$ and end at $v_1$”. We can tell that this is going west, so the angle will be 180 degrees.

Similarly, when it comes time to consider the concavity of $v_3$, we will also consider $v_2$ and $v_4$. The angle that the path $(v_2, v_3, v_4)$ makes relative to the interior is 270$^\circ$, so $v_3$ is concave. Also, the angle that $(v_2, v_3)$ makes relative to the x-axis is 0 degrees, so when it comes time to make our interior lines, we will “extend” $v_3$ to the right (which we have shown in the above picture).

Step 2: Creating interior lines

For each concave vertex $v_i$:

1. Create a vertical (or horizontal) interior line that starts at $v_i$ and ends at the closest horizontal (vertical) line segment in $E$. Call this line segment $l_i$.
2. Add $l_i$ to $E$.

I opted to use a line sweeping based approach for this part of the algorithm. I’ve read that you can create these interior lines in $O(n)$, such as in [Liou et. al, 1990], but their algorithm seemed too complicated so I used line sweeping instead. My approach has $\Theta(n \log n)$ time complexity.

The line sweeping algorithm

Suppose we want to create vertical interior lines. The line sweeping algorithm I used works as follows:

Input: Let $P$ be a simple orthogonal polygon without any chords, and $B$ be its boundary points (its vertices). Each vertex along the boundary is assumed to have additional metadata about its concavity and direction. Let $T$ be an empty binary search tree, and $C$ be an empty container of lines.

1. Sorting. Sort the boundary points by their x-coordinate. Call this sorted collection $W$, and let $w_i$ denote the $i$-th vertex of $W$ ($w_i$ has x-coordinate $x_i$ and y-coordinate $y_i$).
2. Sweeping. For $i = 1, \ldots, n$:
   1. If $y_i$ is not in $T$ then add it $T$, otherwise remove $y_i$ from $T$.
   2. If $w_i$ is a concave vertex with associated angle that is 90 (or 270), then find the closest $\bar{y}_i$ in $T$ that is greater (or less) than $y_i$. Create a point $\bar{w}_i = (x_i, \bar{y}_i)$.
   3. Add the edge $[w_i, \bar{w}_i]$ to $C$.
3. Return $C$.

The sweeping step

You should think of line-sweeping as moving an imaginary vertical³ line from left to right across a group of line segments. This image (taken from topcoder.com) shows the line-sweeping algorithm mid sweep.

Notes on adding and removing points in your tree. The y-coordinates in the tree $T$ represent line segments that you are currently sweeping across. When you add $y_i$ to $T$, you are taking note that, at this moment of the line sweeping process, there is a horizontal line segment that is open and its y-coordinate is $y_i$. When you see the same $y_i$ again, you are closing that line segment.

Notes on creating $\bar{w}_i$. If the vertex is concave with associated angle 90, then using the binary search tree to find the closest $\bar{y}$ is really efficient. Each search of the tree takes $\Theta(\log m)$ operations to complete, where $m$ is the number of elements in the tree.

Horizontal lines

After we create the vertical lines, we need to create the horizontal lines. For this part, add the new points from $C$ to the set of boundary points $B$ and call this set $\bar{B}$. Then use $\bar{B}$ as input into the above line sweeping algorithm. Of course, you will now need to sort $\bar{B}$ by each point’s y-coordinate. For part 1 of the sweeping step, replace $y_i$ with $x_i$. For part 2 of the sweeping step, $\bar{x}_i$ will be the closest value in $T$ that is greater (or less) than ${x}_i$ where $w_i$ has associated angle 0 (or 180).

Time complexity

In this step, we sort the points on our boundary, a $\Theta(n \log n)$ operation. We also do $n$ operations with our binary search tree. Each operation on the tree — checking if the tree already contains a value, inserting or removing a value, finding a closest point greater than (less than) some value — has a time complexity of $O(\log n)$. In total, this part of the algorithm has a total time complexity of $\Theta(n \log n)$.

Step 3: Extracting the rectangles

The line-segments in $E$ (implicitly) partition the polygon. Use them to find the non-overlapping rectangles.

Well that’s vague. Here is how I did it.

Input: The set of edges $C$ created in step-2, and the set of convex vertices denoted by $U$.

1. Setup. Create three empty lists, $C_{ul}$, $C_{ll}$, and $C_{lr}$, which will represent the upper-left, lower-left, and lower-right corners of all the rectangles respectively. Let $R$ denote an (empty) collection of rectangles.

2. Rectangle corner identification. For $z \in C \cup U$:
   1. If $z$ is a convex vertex;
      • If $z$’s angle is 90 then add it to $C_{ul}$
      • If $z$’s angle is 180 then add it to $C_{ll}$
      • If $z$’s angle is 270 then add it to $C_{lr}$
   2. If $z = (w, \bar{w})$ is an interior line;
      • If $w$’s angle is 0 then add $w$ to $C_{ll}$ and add $\bar{w}$ to $C_{lr}$
      • If $w$’s angle is 90 then add $w$ to $C_{lr}$ and add $\bar{w}$ to $C_{ul}$
      • If $w$’s angle is 180 then add $\bar{w}$ to $C_{ll}$ and add $\bar{w}$ to $C_{ul}$
      • If $w$’s angle is 270 then add $w$ to $C_{ul}$ and add $\bar{w}$ to $C_{lr}$ and $C_{ll}$
3. Sort $C_{ul}$ by the points’ x-coordinate, using the y-coordinate to break any ties. Sort $C_{lr}$ by the points’ y-coordinate, using the x-coordinate to break ties.
4. Rectangle construction. For $p \in C_{ll}$:
   1. Find the point in $C_{ul}$ that is closest to $p$ but greater than it, prioritizing $p$’s x-coordinate over its y-coordinate. Call this point $u_1$.
   2. Find the point in $C_{lr}$ that is closest to $p$ but greater than it, prioritizing $p$’s y-coordinate over its x-coordinate. Call this point $u_2$.
   3. Add the rectangle with upper-left point $u_1$ and lower-right point $u_2$ to $R$.
5. Return $R$.

Rectangle corner identification

The main thing to realize is that you can identify points on a rectangle using the point’s “direction” (or angle), and knowing whether the point is convex.

Let’s use an arbitrary rectangle to see how the logic works. The lower left point in a rectangle is pointing left (so has an angle of 180$^\circ$), the upper left point on a rectangle points up (angle of 90$^\circ$), and the lower right point on a rectangle points down (an angle of 270$^\circ$)⁴. All of these points are convex points.

The only trick is realizing that the same implication above goes in both directions. If you are convex and: point left, then you are the lower left point of some rectangle; point up, then you are the upper left point of some rectangle; and point down, then you are the lower right point of some rectangle.

For interior lines, you work through similar logic as in the convex point case.

Rectangle construction

Let $p \in C_{ll}$. This point will be the lower-left point of some rectangle. The upper-left point of this rectangle, $u_p$, will have the same x-coordinate as $p$ but will have a y-coordinate that is greater than $p$’s y-coordinate. In fact, out of all the points in $C_{ul}$ with the same x-coordinate as $p, u_p$ will be the point whose y-coordinate is nearest to $p$’s but is still greater than it.

Similarly, the lower-right point of this rectangle, $l_p$, will have the same y-coordinate as $p$ but will have an x-coordinate that is greater than $p$’s x-coordinate. And out of all the points in $C_{lr}$ with the same y-coordinate as $p, l_p$ will be the point whose x-coordinate is nearest to $p$’s but is still greater than it.

Finding these points is relatively straightforward since we can define a total ordering on the set $C_{ul}$ (sort by x and then by y if there are ties) and $C_{lr}$ (sort by y and then by x if there are ties).

In python, the code would be

import bisect
import collections

Point = collections.namedtuple('Point', ['x', 'y'])
Rectangle = collections.namedtuple('Rectangle', ['ul', 'lr'])

def extract(ul, ll, lr):
    cul = sorted(ul, key=lambda p: (p.x, p.y))
    clr = sorted(lr, key=lambda p: (p.y, p.x))
    slr = list(map(lambda p: Point(p.y, p.x), clr))

    for point in ll:
        i1 = bisect.bisect_right(cul, point)
        i2 = bisect.bisect_right(slr, Point(point.y, point.x))
        yield Rectangle(cul[i1], clr[i2])

Time complexity

In this step, we do a lot of sorting and searching within sorted lists. In step-1, we loop through all the boundary points and interior lines doing a constant amount of work per element. Thus, step-1 has $\Theta(n)$ time complexity. In step-2, we sort two subsets of our boundary points. Since there are fewer rectangles than there are points on the boundary, this is an $\Theta(n \log n)$ operation⁵. We also do $O(n)$ search operations, each having $\Theta(\log n)$ time complexity, yielding a total time complexity of $\Theta(n \log n)$ for this step of the partitioning algorithm.

Closing

I just want to stress that the above logic only works when you know there aren’t any chords in the input polygon. Also, implementing a method that handles polygons with chords is really annoying. The overview algorithm from that blog post I keep linking to looks slightly more complicated than the chord-free version above. Don’t be fooled, it’s way more complicated.

Anyway, the methods implemented in the polygon-partitioner package handle polygons with chords. Just note that when there are chords in the input polygon, I don’t guarantee to return a partitioning with the minimum number of rectangles.