IndexedPointInAreaLocator: Support all CoordinateSequence dimensions

The previous implementation relied on coordinates in all indexed segments
being offset by 24 bytes in order to optimize index construction and
query time. This commit retains the same performance performance while
allowing offsets of 16, 24, or 32 bytes to accommodate the different
coordinate types.

include/geos/algorithm/locate/IndexedPointInAreaLocator.h

@@ -54,22 +54,32 @@ namespace locate { // geos::algorithm::locate
 class GEOS_DLL IndexedPointInAreaLocator : public PointOnGeometryLocator {
 private:
     struct SegmentView {
-        SegmentView(const geom::Coordinate* p_p0, const geom::Coordinate* p_p1) : m_p0(p_p0) {
-            // All GEOS CoordinateSequences store their coordinates sequentially.
-            // Should that ever change, this assert will fail.
-            (void) p_p1;
-            assert(p_p0 + 1 == p_p1);
+        SegmentView(const geom::CoordinateXY* p0, const geom::CoordinateXY* p1) {
+            // There is a significant performance benefit in fitting our
+            // line segment into 8 bytes (about 15-20%). Because we know that
+            // p1 follows p0 in a CoordinateSequence, we know that the address
+            // of p1 is 16, 24, or 32 bytes greater than the address of p0.
+            // By packing this offset into the least significant bits of p0,
+            // we can retrieve both p0 and p1 while only using 8 byytes.
+            std::size_t os = static_cast<std::size_t>(reinterpret_cast<const double*>(p1) - reinterpret_cast<const double*>(p0)) - 2u;
+            m_p0 = reinterpret_cast<std::size_t>(p0) | os;
+
+            assert(&this->p0() == p0);
+            assert(&this->p1() == p1);
         }
 
-        const geom::Coordinate& p0() const {
-            return *m_p0;
+        const geom::CoordinateXY& p0() const {
+            auto ret = reinterpret_cast<const geom::CoordinateXY*>(m_p0 >> 2 << 2);
+            return *ret;
         }
 
-        const geom::Coordinate& p1() const {
-            return *(m_p0 + 1);
+        const geom::CoordinateXY& p1() const {
+            auto offset = (m_p0 & 0x03) + 2;
+            auto ret = reinterpret_cast<const geom::CoordinateXY*>(reinterpret_cast<double*>(m_p0 >> 2 << 2) + offset);
+            return *ret;
         }
 
-        const geom::Coordinate* m_p0;
+        std::size_t m_p0;
     };
 
     class IntervalIndexedGeometry {

src/algorithm/locate/IndexedPointInAreaLocator.cpp

@@ -32,6 +32,8 @@
 #include <algorithm>
 #include <typeinfo>
 
+using geos::geom::CoordinateXY;
+
 namespace geos {
 namespace algorithm {
 namespace locate {
@@ -73,7 +75,7 @@ void
 IndexedPointInAreaLocator::IntervalIndexedGeometry::addLine(const geom::CoordinateSequence* pts)
 {
     for(std::size_t i = 1, ni = pts->size(); i < ni; i++) {
-        SegmentView seg(&pts->getAt(i-1), &pts->getAt(i));
+        SegmentView seg(&pts->getAt<CoordinateXY>(i-1), &pts->getAt<CoordinateXY>(i));
         auto r = std::minmax(seg.p0().y, seg.p1().y);
 
         index.insert(index::strtree::Interval(r.first, r.second), seg);

tests/unit/algorithm/locate/IndexedPointInAreaLocatorTest.cpp

@@ -0,0 +1,89 @@
+#include <tut/tut.hpp>
+// geos
+#include <geos/geom/CoordinateSequence.h>
+#include <geos/geom/GeometryFactory.h>
+#include <geos/algorithm/locate/IndexedPointInAreaLocator.h>
+
+// std
+#include <sstream>
+#include <string>
+#include <memory>
+
+using geos::geom::CoordinateXY;
+using geos::geom::CoordinateSequence;
+using geos::geom::GeometryFactory;
+using geos::geom::Location;
+using geos::algorithm::locate::IndexedPointInAreaLocator;
+
+namespace tut {
+//
+// Test Group
+//
+
+// dummy data, not used
+struct test_indexedpointinarealocator_data {
+    const GeometryFactory& factory_ = *GeometryFactory::getDefaultInstance();
+};
+
+typedef test_group<test_indexedpointinarealocator_data> group;
+typedef group::object object;
+
+group test_indexedpointinarealocator_group("geos::algorithm::locate::IndexedPointInAreaLocator");
+
+// Test all CoordinateSequence dimensions
+template<>
+template<>
+void object::test<1>
+()
+{
+    // XY
+    CoordinateSequence seq_xy = CoordinateSequence::XY(0);
+    seq_xy.add(0.0, 0.0);
+    seq_xy.add(1.0, 0.0);
+    seq_xy.add(1.0, 1.0);
+    seq_xy.add(0.0, 1.0);
+    seq_xy.add(0.0, 0.0);
+    auto ls_xy = factory_.createLineString(seq_xy.clone());
+    IndexedPointInAreaLocator ipa_xy(*ls_xy);
+
+    CoordinateXY pt_boundary(0.5, 0);
+    CoordinateXY pt_interior(0.5, 0.5);
+    CoordinateXY pt_exterior(1.5, 0.5);
+
+    ensure_equals(ipa_xy.locate(&pt_boundary), Location::BOUNDARY);
+    ensure_equals(ipa_xy.locate(&pt_interior), Location::INTERIOR);
+    ensure_equals(ipa_xy.locate(&pt_exterior), Location::EXTERIOR);
+
+    // XYZ
+    CoordinateSequence seq_xyz = CoordinateSequence::XYZ(0);
+    seq_xyz.add(seq_xy);
+    auto ls_xyz = factory_.createLineString(seq_xyz.clone());
+
+    IndexedPointInAreaLocator ipa_xyz(*ls_xy);
+    ensure_equals(ipa_xyz.locate(&pt_boundary), Location::BOUNDARY);
+    ensure_equals(ipa_xyz.locate(&pt_interior), Location::INTERIOR);
+    ensure_equals(ipa_xyz.locate(&pt_exterior), Location::EXTERIOR);
+
+    // XYM
+    CoordinateSequence seq_xym = CoordinateSequence::XYM(0);
+    seq_xym.add(seq_xy);
+    auto ls_xym = factory_.createLineString(seq_xym.clone());
+
+    IndexedPointInAreaLocator ipa_xym(*ls_xy);
+    ensure_equals(ipa_xym.locate(&pt_boundary), Location::BOUNDARY);
+    ensure_equals(ipa_xym.locate(&pt_interior), Location::INTERIOR);
+    ensure_equals(ipa_xym.locate(&pt_exterior), Location::EXTERIOR);
+
+    // XYZM
+    CoordinateSequence seq_xyzm = CoordinateSequence::XYZM(0);
+    seq_xyzm.add(seq_xy);
+    auto ls_xyzm = factory_.createLineString(seq_xyzm.clone());
+
+    IndexedPointInAreaLocator ipa_xyzm(*ls_xy);
+    ensure_equals(ipa_xyzm.locate(&pt_boundary), Location::BOUNDARY);
+    ensure_equals(ipa_xyzm.locate(&pt_interior), Location::INTERIOR);
+    ensure_equals(ipa_xyzm.locate(&pt_exterior), Location::EXTERIOR);
+}
+
+}
+