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 polygonpartitioner, and it’s primary task is breaking up orthogonal polygons into nonoverlapping 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^{1}. 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) linesegment that lies within the polygon and has at least one endpoint that is a concave vertex. The other endpoint of this linesegment 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 counterclockwise relative to the positive xaxis.
An outline of the algorithm
The main challenge with partitioning orthogonal polygons is the presence of chords. The overall algorithm for partitioning a chordfree 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 chordfree case. Here is the algorithm for partitioning chordfree polygons:
Input: An orthogonal polygon $P$ without any chords.
 Setup; extracting linesegments. Let $E$ be the set of linesegments along the boundary of $P$.
 Create interior line. For each concave vertex $v_i$:
 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$.
 Add $l_i$ to $E$.
 Extract rectangles. The linesegments in $E$ (implicitly) partition the polygon. Use them to find the nonoverlapping rectangles.
I will go through each step in detail for the remainder of this post.
Step 1: Extracting linesegments
Let $E$ be the set of linesegments 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 xcoordinate, or have the same ycoordinate.
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$^{2}. 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_{i1}, v_i$) makes relative to the xaxis, 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 xaxis. 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 xaxis 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$:
 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$.
 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.
 Sorting. Sort the boundary points by their xcoordinate. Call this sorted collection $W$, and let $w_i$ denote the $i$th vertex of $W$ ($w_i$ has xcoordinate $x_i$ and ycoordinate $y_i$).
 Sweeping. For $i = 1, \ldots, n$:
 If $y_i$ is not in $T$ then add it $T$, otherwise remove $y_i$ from $T$.
 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)$.
 Add the edge $[w_i, \bar{w}_i]$ to $C$.
 Return $C$.
The sweeping step
You should think of linesweeping as moving an imaginary vertical^{3} line from left to right across a group of line segments. This image (taken from topcoder.com) shows the linesweeping algorithm mid sweep.
Notes on adding and removing points in your tree. The ycoordinates 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 ycoordinate 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 ycoordinate. 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 linesegments in $E$ (implicitly) partition the polygon. Use them to find the nonoverlapping rectangles.
Well that’s vague. Here is how I did it.
Input: The set of edges $C$ created in step2, and the set of convex vertices denoted by $U$.

Setup. Create three empty lists, $C_{ul}$, $C_{ll}$, and $C_{lr}$, which will represent the upperleft, lowerleft, and lowerright corners of all the rectangles respectively. Let $R$ denote an (empty) collection of rectangles.
 Rectangle corner identification. For $z \in C \cup U$:
 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}$
 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}$
 If $z$ is a convex vertex;
 Sort $C_{ul}$ by the points’ xcoordinate, using the ycoordinate to break any ties. Sort $C_{lr}$ by the points’ ycoordinate, using the xcoordinate to break ties.
 Rectangle construction. For $p \in C_{ll}$:
 Find the point in $C_{ul}$ that is closest to $p$ but greater than it, prioritizing $p$’s xcoordinate over its ycoordinate. Call this point $u_1$.
 Find the point in $C_{lr}$ that is closest to $p$ but greater than it, prioritizing $p$’s ycoordinate over its xcoordinate. Call this point $u_2$.
 Add the rectangle with upperleft point $u_1$ and lowerright point $u_2$ to $R$.
 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$)^{4}. 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 lowerleft point of some rectangle. The upperleft point of this rectangle, $u_p$, will have the same xcoordinate as $p$ but will have a ycoordinate that is greater than $p$’s ycoordinate. In fact, out of all the points in $C_{ul}$ with the same xcoordinate as $p, u_p$ will be the point whose ycoordinate is nearest to $p$’s but is still greater than it.
Similarly, the lowerright point of this rectangle, $l_p$, will have the same ycoordinate as $p$ but will have an xcoordinate that is greater than $p$’s xcoordinate. And out of all the points in $C_{lr}$ with the same ycoordinate as $p, l_p$ will be the point whose xcoordinate 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 step1, we loop through all the boundary points and interior lines doing a constant amount of work per element. Thus, step1 has $\Theta(n)$ time complexity. In step2, 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^{5}. 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 chordfree version above. Don’t be fooled, it’s way more complicated.
Anyway, the methods implemented in the polygonpartitioner 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.

This problem (partitioning polygons) originally came up as part of a work related project. The algorithms I developed for my employer were completely different from the ones available publicly on my github page (it was also written in a different language). Also, the code I wrote for my employer could do more (as in handle more edge cases) than what I have in the public project. ↩

In order for three points along the boundary to have an angle of zero, they all must lie along the same linesegment. We have eliminated this possibility during the preprocessing step so this can’t happen. ↩

We would be moving an imaginary horizontal line for bottom to top if we were linesweeping horizontally. ↩

The upper right point on a rectangle points to the right, so it would have an angle if 0$^\circ$. ↩

The minimum number of rectangles for an orthogonal polygon is $n$ / 2 + $h$ − $c$ − 1, where $n$ is the number of vertices, $h$ is the number of holes, and $c$ is the number of chords. See section 3 of Eppstein et. al (here is the arXiv preprint of that paper). ↩