Overview

S2ClosestEdgeQuery is a tool for finding the closest edge(s) between two geometries. By using the appropriate options, this class can answer questions such as:

  • Find the minimum distance between two geometries A and B.

  • Find all edges of geometry A that are within a distance d of geometry B.

  • Find the k edges of geometry A that are closest to a given point P.

The input geometries may consist of any number of points, polylines, and polygons (collectively known as shapes). Shapes do not need to be disjoint; they may overlap or intersect arbitrarily. The implementation is designed to be fast for both simple and complex geometries (e.g., one edge vs. 100 million edges).

The Index

S2ClosestEdgeQuery operates on two geometries, known as the index and the target. Each query returns the subset of edges in the index that have a specified relationship to the target (e.g., the closest edge to the target, or all edges within a distance d of the target).

The indexed geometry is provided as an S2ShapeIndex object, which is simply a collection of points, polylines, and/or polygons (possibly intersecting, as mentioned above). Here is a simple example showing how to build an index containing several polylines:

vector<S2Polyline*> polylines = ...;
S2ShapeIndex index;
for (S2Polyline* polyline : polylines) {
  index.Add(new S2Polyline::Shape(polyline));
}

An S2ClosestEdgeQuery object can then be constructed like so:

S2ClosestEdgeQuery query(&index);

Note that index is passed as a const pointer rather than a reference to reflect the fact that the index must persist for the lifetime of the query object.

The Target

The target represents the geometry that distances are measured to. It is specified as an object of type S2ClosestEdgeQuery::Target. There are Target subtypes for measuring the distance to a point, an edge, an S2Cell, or an entire S2ShapeIndex (which allows finding closest edges or measuring distances between two arbitrary collections of geometry).

For example, the following code tests whether any of the polylines in the example above is closer than dist_limit to a given point test_point:

S2Point test_point = ...;
S1ChordAngle dist_limit = ...;
S2ClosestEdgeQuery::PointTarget target(test_point);
if (query.IsDistanceLess(&target, dist_limit)) { ... }

Note that target is passed as a pointer; this is because Target objects may contain temporary data and can therefore only be used by one query at a time. (The target and query object are intended for use only within a single thread, as opposed to the S2ShapeIndex object which can be used in many different threads simultaneously.)

Query Methods

The fastest and simplest query method is IsDistanceLess, which returns true if the distance to the target is less than a given limit:

bool IsDistanceLess(Target* target, S1ChordAngle limit);

If you want to test whether the distance is less than or equal to “limit” instead, you can use IsDistanceLessOrEqual instead:

if (query.IsDistanceLessOrEqual(&target, limit)) { ... }

To compute the actual distance, you can use the GetDistance method:

S1ChordAngle GetDistance(Target* target);

However note that if you only need to compare the distance against some threshold value, it is much faster to call IsDistanceLess instead.

The return type of this method (S1ChordAngle) is an efficient representation of the distance between two points on the unit sphere. It can easily be converted to an S1Angle if desired (which measures angles in radians or degrees), and from there to an actual distance on the Earth spheroid. (See also Modeling Accuracy.)

The most general query method is FindClosestEdges:

std::vector<Result> FindClosestEdges(Target* target);

This method finds the closest edge(s) to the target that satisfy the query options (described in the following section). Each edge is returned in the form of a Result object that identifies the edge and specifies the distance to that edge:

class Result {
  Distance distance() const;  // The distance from the target to this edge.
  int32 shape_id() const;     // Identifies an indexed shape.
  int32 edge_id() const;      // Identifies an edge within the shape.
};

The Distance type is simply a thin wrapper around S1ChordAngle that provides some additional methods needed by the S2ClosestEdgeQuery implementation. The other two fields identify an edge in the S2ShapeIndex that was passed to the S2ClosestEdgeQuery constructor. The query object has two convenience methods that may be useful for processing result edges:

// Returns the endpoints of the given result edge.
S2Shape::Edge GetEdge(const Result& result) const;

// Returns the point on given result edge that is closest to "point".
S2Point Project(const S2Point& point, const Result& result) const;

For example, to retrieve the endpoints of each edge and also find the closest point on each result edge to a given target point, you could do this:

for (const auto& result : query.FindClosestEdges(&target)) {
  S2Shape::Edge edge = query.GetEdge(result);
  S2Point closest_point = query.Project(point, result);
}

There is also a convenience method for the common case where only one edge (the closest edge is needed). This method returns a single Result object rather than a vector:

Result FindClosestEdge(Target* target);

Options

S2ClosestEdgeQuery supports a number of options that control which edges are returned by each query. These options can be specified as a second argument to the constructor:

S2ClosestEdgeQuery::Options options;
S2ClosestEdgeQuery query(&index, options);

Unlike most classes in the library, S2ClosestEdgeQuery also allows options to be modified between calls to query methods via the mutable_options() accessor method. For example:

S1ChordAngle distance1 = query.GetDistance(&target1);
query.mutable_options()->set_max_distance(dist_limit);
S1ChordAngle distance2 = query.GetDistance(&target2);

This can be significantly more efficient compared to changing the options by constructing a new S2ClosestEdgeQuery object, because the query object contains a significant amount of internal state that is reused from one query to the next.

The most important options are the following:

// Specifies that at most "max_results" edges should be returned.
//
// REQUIRES: max_results >= 1
// DEFAULT: numeric_limits<int>::max()
int max_results() const;
void set_max_results(int max_results);

// Specifies that only edges whose distance to the target is less than
// "max_distance" should be returned.
//
// Note that edges whose distance is exactly equal to "max_distance" are
// not returned.  Normally this doesn't matter, because distances are not
// computed exactly in the first place, but if such edges are needed then
// see set_inclusive_max_distance() below.
//
// DEFAULT: Distance::Infinity()
void set_max_distance(S1ChordAngle max_distance);

// Like set_max_distance(), except that edges whose distance is exactly
// equal to "max_distance" are also returned.  Equivalent to calling
// set_max_distance(max_distance.Successor()).
void set_inclusive_max_distance(S1ChordAngle max_distance);

Note that you will always want to set one of these options before calling FindClosestEdges, because otherwise all the edges in the entire S2ShapeIndex will be returned. (This is not necessary when calling FindClosestEdge, GetDistance, or IsDistanceLess, because these methods all implicitly limit max_results() to 1.)

Another important option is include_interiors(), which determines whether distances are measured to the boundary and interior of polygons or only to their boundaries:

// Specifies that polygon interiors should be included when measuring
// distances.
//
// DEFAULT: true
bool include_interiors() const;
void set_include_interiors(bool include_interiors);

By default this option is true, which means for example that a point inside a polygon has a distance of zero. Note that in this situation there is no “closest edge” to return (since the minimum distance is attained at a point interior to the polygon), and therefore the FindClosestEdge(s) methods indicate this in the Result object(s) by setting edge_id to -1. (Such results can be identified using the is_interior() method.) For example:

for (const auto& result : query.FindClosestEdges(&target)) {
  if (result.is_interior()) {
    ProcessInterior(result.shape_id());
  } else {
    ProcessEdge(result.shape_id, result.edge_id);
  }
}

The following option can be useful when performance is critical:

// Specifies that edges up to max_error() further away than the true
// closest edges may be substituted in the result set, as long as such
// edges satisfy all the remaining search criteria (such as max_distance).
//
// DEFAULT: Distance::Zero()
Distance max_error() const;
void set_max_error(Distance max_error);

This option give the algorithm permission to stop the search early when the best possible improvement in distance drops below max_error(). For example, this option is used internally by the IsDistanceLess predicate to stop the search as soon as it finds any edge whose distance is less than the given limit, rather than continuing to search for an edge that is even closer (which would be wasted effort).

Note that this options only has an effect if max_results() is also specified; otherwise all edges closer than max_distance() will always be returned.

Target Types

The following target types exist in addition to the PointTarget type described earlier.

EdgeTarget computes the closest distance to a spherical geodesic edge. For example, to find the closest indexed edge to a target edge ab, you could do this:

S2ClosestEdgeQuery::EdgeTarget target(a, b);
auto result = query.FindClosestEdge(&target);

CellTarget computes the distance to an S2Cell, including the interior of the cell. This target type is used extensively internally to implement other algorithms. For example, to test whether the distance to cell is less than limit:

S2ClosestEdgeQuery::CellTarget target(cell);
if (query.IsDistanceLess(&target, limit)) { ... }

Finally, ShapeIndexTarget computes the distance to an S2ShapeIndex (containing an arbitrary collection of points, polylines, and polygons, possibly overlapping). For example, to calculate the distance between two S2ShapeIndexes, you can do this:

S2ClosestEdgeQuery query(&index1);
S2ClosestEdgeQuery::ShapeIndexTarget target(&index2);
S1ChordAngle distance = query.GetDistance(&target);

(As always, if you only want to compare the distance against a threshold value then it is much more efficient to call IsDistanceLess instead.)

The ShapeIndexTarget type also supports its own Options. In particular, if you want to measure the distance to the boundaries of polygons in the target index (ignoring the polygon interiors), you can call

target.set_include_interiors(false);

before executing the distance query. Note that the include_interiors() option on the query and the target are independent, so for example you can measure distance from the boundary of one polygon to the boundary and interior of another.

Modeling Accuracy

This class models the Earth as a sphere, which can lead to distance errors of up to 0.56% compared with modeling the Earth as an ellipsoid. While this is perfectly adequate for many purposes (e.g., determining nearby restaurants), some applications need higher accuracy.

This can be achieved by recomputing distances using using an external geodesy library such as geographiclib. For exact results, all distances would need to be computed this way (and in fact there is a templatized base class S2ClosestEdgeQueryBase that allows this), but in most cases it is sufficient to use spherical distances to determine which pair of points is closest, and then recalculate the distance accurately for that pair of points only. (This is because the spherical distance error varies relatively slowly over the Earth’s surface.) For example, if the target is a point then you could do the following:

S2ClosestEdgeQuery::PointTarget target(target_point);
auto result = query.FindClosestEdge(&target);
double meters;
if (result.is_empty()) {
  // The index is empty.
  meters = std::numeric_limits<double>::infinity();
} else if (result.is_interior()) {
  // The target point is inside an indexed shape.
  meters = 0.0;
} else {
  // Find the point on the result edge that is closest to the target.
  S2Point closest_point = query.Project(target_point, result);
  meters = GeoidDistance(target_point, closest_point);
}

The same technique can be used for more complex targets (ShapeIndexTarget) by executing a second query to find the closest pair of edges between the two geometry collections, finding the points on that edge pair that achieve the minimum spherical distance (see S2::GetEdgePairClosestPoints), and then calculating the ellipsoidal distance between those points:

S2ClosestEdgeQuery query1(&index1);
S2ClosestEdgeQuery::ShapeIndexTarget target2(&index2);
auto result1 = query1.FindClosestEdge(&target2);
double meters;
... like the code above ...
} else {
  // Get the edge from index1 (edge1) that is closest to index2.
  S2Shape::Edge edge1 = query1.GetEdge(result1);

  // Now find the edge from index2 (edge2) that is closest to edge1.
  S2ClosestEdgeQuery query2(&index2);
  S2ClosestEdgeQuery::EdgeTarget target1(edge1.v0, edge1.v1);
  auto result2 = query2.FindClosestEdge(&target1);
  S2Shape::Edge edge2 = query2.GetEdge(result2);

  // Find the closest point pair on edge1 and edge2.
  auto closest = S2::GetEdgePairClosestPoints(edge1.v0, edge1.v1,
                                              edge2.v0, edge2.v1);
  meters = GeoidDistance(closest.first, closest.second);
}

Numerical Accuracy

The spherical geodesic distance calculation used by this class is fast and accurate, but not exact. Expressing the errors in terms of distances on the Earth’s surface, for distances up to 10,000 km the error is less than 6 nanometers. However for points that are nearly antipodal, the error can be as large as 50 cm (dropping off to less than 1 micrometer for points that are 50 km away from being antipodal).

If you would like to measure distances conservatively by including these error bounds, you can use the following option:

options.set_conservative_max_distance(limit);

This will automatically increase limit by the maximum error in the distance calculation, ensuring that all edges whose true, exact distance is less than or equal to limit are returned (along with some edges whose true distance is slightly greater).

Algorithms that need to do exact distance comparisons can use this option to compute a set of candidate edges that can then be filtered further using exact predicates (such as s2pred::CompareEdgeDistance).

There is also a conservative version of the IsDistanceLessOrEqual predicate,

bool IsConservativeDistanceLessOrEqual(Target* target, S1ChordAngle limit);

which automatically increases the given limit by the maximum error in the distance calculation.