The Explorer

Human fascination with exploring the unknown has persisted across time and space. From ancient texts like The Epic of Gilgamesh to childhood classics like Dora The Explorer, the archetype is firmly rooted in human experience.

Our history is rife with individuals embodying the explorer ethos to extraordinary degree. In fact, my DBQ for AP US History concerned itself with this exact topic. I'll include a lightly adapted snippet from the introduction for reference.

...After securing independence from Britain, the United States began looking westward. Under the Articles of Confederation, the US established a system of creating new states that could join the union under the Land Ordinances of 1785 and 1787. Improved relationships with foreign powers led to deals such as Jay's Treaty, the Louisiana Purchase, and the Treaty of 1818, which fueled America's westward expansion. Despite the rise of sectionalist tensions that would drive the United States down a path of disunion, it is clear that Manifest Destiny played a key role in making the United States stretch from coast to coast.

Interested in the archetype, I eventually found myself at Point Nemo.

Point Nemo

I want you to imagine that you are an explorer. Cast amongst the ranks of those legendary icons from history, your task is to sail atop the farthest point from land. A brief search on the National Ocean Service lands you at 48° 52.6′S, 123° 23.6′W better known as Point Nemo. A staggering 2,688 miles from Dunice Island, it is the most inaccessible point in the ocean.

Your browser does not support SVG
Point Nemo and the Surrounding Area

The natural question, of course, is how one might go about finding Point Nemo for themselves.

Theory

At this point, we are going to take a brief sidestep into theory. Our first task to define what we are looking for.

Consider a set of points in a metric space. The set divides the space it lives in into three other sets: the Interior, Exterior, and Boundary.

Your browser does not support SVG
The Interior, Exterior, and Boundary

A polygon \( P \) is a set of points whose boundary forms a sequence of edges \( E \). We define the pole of inaccessiblility as the point within the polygon that is farthest from its edges.

\[ \text{pole} = \arg\max_{p \in P} \left( \min_{e \in E} \mathrm{d}(p, e) \right) \]

Solving this optimization problem is not easy in general. One approach starts by partitioning the plane into regions closest to each vertex. This tesselation is called a Voronoi Diagram and may be constructed in \( O(n \log{n}) \) time via Fortune's Algorithm.

Your browser does not support SVG
Generalized Voronoi Diagram

With the Voronoi Diagram on hand, we search through every node (intersection of three cells) and record the best so far. This is guranteed to yield the optimal solution.

PolyLabel

While the procedure above produces the pole, we might ask ourselves if we can trade precision for speed? As it so happens, we can!

Description

PolyLabel is an iterative grid-based algorithm, which starts by covering the polygon with square cells and iteratively subdividing them in the order of the most promising ones. The procedure applies to any complex polygon (e.g., with holes) and is guranteed to produce the optimal solution to a given precision.

I'll leave it as an excersize for the motivated reader to check out the MapBox README for detail.

Proof

Consider a Polygon \( P \) and Cell \( C \) with center \( c \) and half-side \( h \). Let \( \mathrm{d} \) be a function that yields the minimum distance from a Point to a Polygon. We make the following claim: the maximum distance \( m \) from any point \( p \in C \) to \( P \) is bounded by \( \mathrm{d}(P, p) \le \mathrm{d}(P, c) + h \sqrt{2} \). We justify this claim by way of contradiction.

Assume, for the sake of contradiction, that there exists a point \( p \in C \) such that \( \mathrm{d}(P, p) > m \). The distance from \( p \) to the edge closest to \( c \) must be greater than than \( m \) because \( \mathrm{d}(P, p) \) is the minimum distance. This implies that \( p \) lies outside of \( C \). This contradicts the assumption that \( p \in C \). Therefore, we reject the premise that \( p \) exists.

Your browser does not support SVG
Complex Polygon Covered by a Grid of Cells

The algorithm may now be understood as follows:

  1. Cover the polygon with Cells. The pole must be in a Cell.
  2. If the maximum distance of a Cell exceeds the current best distance, it is promising and we should subdivide.

By continually subdividing, we eventually hone in on the optimal solution. The only wrinkle is that this process will continue indefinitely unless we specify a tolerance.

Implementation

To implement this algorithm, I created a Polygon and Cell Class. The Polygon Class precomputes data (e.g., centroid, boundary, etc) and provides a member function to compute the signed distance from the Polygon. If the sign is positive, the Point is inside. If the sign is negative, the Point is outside. The Cell Class maintains the center, half-side, and potential. Direct comparison between Cells is supported.

Upon profiling, it became clear that the key optimization for this algorithm is efficiently computing the distance from a given Point. A vectorized solution is shared below:

                
                    def compute_distance(self, point):
                        """Compute the distance from a point to the polygon."""
                        scalars = np.einsum("ij,ij->i", self.norm, point - self.start)
                        clips = np.clip(scalars, 0, 1).reshape(-1, 1)
                        d = np.min(
                            np.norm(
                                self.start + self.diff * clips - point, 
                                axis=1
                            )
                        )
                        return d if self.inside(point) else -d
                
            

The boundary of the Polygon is defined by start, stop, and diff vectors. The distance computation has two parts. First, we compute the projection of the point onto each edge and clip the result to lie on the edge. This scalar tells us how far to move in the direction of diff. Then, we compute the distance from each edge to the point and return the minimum. Checking whether a point is inside a polygon is done via Ray-Casting. You can find the complete implementation on GitHub.

Resources and References

If you are interested in further reading, I recommend checking out Paul Bourke's notes on Polygons and Meshes and Stefan Huber disertation on the Vroni Project.

Comment Section