within

[mapsam]

@artemp I'm wondering if I can use intersection.hpp to perform a "within"-esque operation?I'd like to basically get a boolean response if two geometries intersect - in my case a point-in-polygon operation "is this point within this polygon?"

Looks like boost has a within function - perhaps I'll use this for now?

refs: mapbox/vtquery#36

[artemp]

@mapsam - Absolutely ^ , using boost::geometry::within is a good plan for now. Just be aware that within is used as an optimisation step for polygons in closest_point, so perhaps we should avoid calling it multiple times.
Also, within with default strategy will return false if query point is exactly on a boundary. See boost docs for more info.

[mapsam]

@artemp thanks for noting that closet_point uses boost::geometry::within internally - seeing that it sets the distance to 0.0 exactly - so I'm now just checking if (cpinfo.distance == 0.0) and assuming it is a truthy PIP hit.

[springmeyer]

so I'm now just checking if (cpinfo.distance == 0.0) and assuming it is a truthy PIP hit.

@mapsam that feels like a good solution to me. One thing I noticed about the within usage inside of closest_point_impl.hpp is that the multipolygon case (

spatial-algorithms/include/mapbox/geometry/algorithms/closest_point_impl.hpp

Lines 64 to 86 in 5dfca06

	template <typename MultiGeometry>
	result_type operator() (MultiGeometry const& multi_geom) const
	{
	result_type result;
	bool first = true;
	for (auto const& geom : multi_geom)
	{
	if (first)
	{
	first = false;
	result = std::move(operator()(geom));
	}
	else
	{
	auto sub_result = operator()(geom);
	if (sub_result.distance < result.distance)
	{
	result = std::move(sub_result);
	}
	}
	}
	return result;
	}

) calls out to the single polygon case (

spatial-algorithms/include/mapbox/geometry/algorithms/closest_point_impl.hpp

Lines 33 to 56 in 5dfca06

	result_type operator() (mapbox::geometry::polygon<coordinate_type> const& poly) const
	{
	info_type info;
	if (boost::geometry::within(pt_, poly))
	{
	return result_type(static_cast<double>(pt_.x), static_cast<double>(pt_.y), 0.0);
	}
	bool first = true;
	for (auto const& ring : poly)
	{
	info_type ring_info;
	boost::geometry::closest_point(pt_ ,ring, ring_info);
	if (first)
	{
	first = false;
	info = std::move(ring_info);
	}
	else if (ring_info.distance < info.distance)
	{
	info = std::move(ring_info);
	}
	}
	return result_type(info.closest_point.x, info.closest_point.y, info.distance);
	}

). If I'm reading this correctly I think it means you'll get back cpinfo.distance == 0.0 if the query point is within any one part of a multipolygon and it will check all polygon parts before returning. Does this seem like the correct behavior for your use case? One performance issue I see is that if the first part of a multipolygon is directly hit (cpinfo.distance == 0.0), we still process the others, which is unlikely needed.

[mapsam]

@springmeyer looks like you're reading that correctly to me. I wonder if closest_point_impl can break out of the loop if it gets a hit? Although, this means it wouldn't take into account any overlapping geometries in a multipolygon (is that allowed?).

[springmeyer]

I wonder if closest_point_impl can break out of the loop if it gets a hit?

Yes, that could work - @artemp thoughts? Something like this perhaps?

diff --git a/include/mapbox/geometry/algorithms/closest_point_impl.hpp b/include/mapbox/geometry/algorithms/closest_point_impl.hpp
index 2103c77..94411cb 100644
--- a/include/mapbox/geometry/algorithms/closest_point_impl.hpp
+++ b/include/mapbox/geometry/algorithms/closest_point_impl.hpp
@@ -72,6 +72,10 @@ struct closest_point
             {
                 first = false;
                 result = std::move(operator()(geom));
+                // early return if exact hit
+                if (result.distance == 0.0) {
+                    return result;
+                }
             }
             else
             {

Note: probably need to use an epsilon since comparing floating point numbers is not ideal.

Although, this means it wouldn't take into account any overlapping geometries in a multipolygon (is that allowed?).

From a generic spatial-algorithms perspective I assume our docs would say this is invalid, but would need @artemp's confirmation.

From a vector tile perspective overlapping geometries would be invalid, so I don't think this would need supported (and breaking early would be fine).

[artemp]

2.2.12 MultiPolygon

A MultiPolygon is a MultiSurface whose elements are Polygons..
The assertions for MultiPolygons are :
1.
The interiors of 2 Polygons that are elements of a MultiPolygon may not intersect.
∀ M ∈ MultiPolygon, ∀ Pi, Pj ∈ M.Geometries(), i ≠ j, Interior(Pi) ∩ Interior(Pj) = ∅
2.
The Boundaries of any 2 Polygons that are elements of a MultiPolygon may not ‘cross’ and may touch
at only a finite number of points. (Note that crossing is prevented by assertion 1 above).
∀ M ∈ MultiPolygon, ∀ Pi, Pj ∈ M.Geometries(), ∀ ci ∈ Pi.Boundaries(), cj ∈ Pj.Boundaries()
ci ∩ cj = {p1, ....., pk | pi ∈ Point, 1 <= i <= k}
3. A MultiPolygon is defined as topologically closed.
4. A MultiPolygon may not have cut lines, spikes or punctures, a MultiPolygon is a Regular, Closed
point set:
∀ M ∈ MultiPolygon, M = Closure(Interior(M))
5.
The interior of a MultiPolygon with more than 1 Polygon is not connected, the number of connected
components of the interior of a MultiPolygon is equal to the number of Polygons in the MultiPolygon.
The boundary of a MultiPolygon is a set of closed curves (LineStrings) corresponding to the boundaries of
its element Polygons. Each curve in the boundary of the MultiPolygon is in the boundary of exactly 1
element Polygon, and every curve in the boundary of an element Polygon is in the boundary of the
MultiPolygon.
The reader is referred to work by Worboys, et. al (7, 8) and Clementini, et. al (5, 6) for work on the
definition and specification of MultiPolygons.

Yep, great catch! Because we're assuming that geometries are valid and simple in OGC sense - ^ early return is totally justified.

/cc @mapsam @springmeyer

[artemp]

the same applies to MultiLineStrings ^ and in-fact this optimisation makes sense even for invalid geometries.

[artemp]

@springmeyer -
comparing to 0.0 is fine in this case because we're setting distance to 0.0 in the first place.
my bad, always use boost::geometry::math::equal for equality test of floating point numbers.
fixed in 37a43c4