Geos ccw analysis

[BuonOmo]

Summary

We needed to avoid performance loss when computing ring direction for rgeo/rgeo-geojson#47. This led me to this implementation of a ccw? method for GEOS C API, which is way faster than the ruby implementation (see rgeo/rgeo-geojson#47 (comment))

Commit details

Add clockwise analysis for GEOS_CAPI 8b640b0

Add a RGeo::Geos::Analysis.ccw? method which itself uses the
GEOSCoordSeq_isCCW_r to determine if a geometry is counter-clockwise.

Also add an error.h header to help throwing errors when there is an
issue in the c extension. And add a RGeo::GeosError for GEOS related
errors.

Change from require to require_relative 163a68f

This helps a lot a developer which has a global version installed and
tries to contribute to the library, since it will strictly look for
related files...

For review

I did some opiniated changes, I'd be glad to reverse some if you do not agree on those :)

[keithdoggett]

Other c_api tests add the if RGeo::Geos.capi_supported? statement after the closure to the test class. I like this approach, so at some point, we could change all of them to this format. We can save that for another PR though.

[BuonOmo]

Yeah I saw the if RGeo::Geos.capi_supported? and the WARNING, but this is hard to see IMHO. I'm all for changing that. I won't be able to work on it until next week. If you have time and you want to add a commit to this PR I'd be OK. Or I'll open an issue and we'll add it to our low priority queue :)

[keithdoggett]

This looks great! I just have 2 other things to note about it.

1. Does it make sense to add a LinearRing#ccw? method to the LinearRing feature class? This is not how JTS and libgeos handle it, but it could be convenient to include it.
2. Your build failed because this is only supported in libgeos version 3.7+. One thought I have on how to handle this is to edit the extconf.rb file to test if the isCCW function exists and define a macro like RGEO_IS_CCW_DEFINED and put an #ifdef statement for it at the top of the analysis.c/h files. Then we could check if the RGeo::Geos::Analysis module is defined inside the ccw? function. These steps may not work perfectly right but something along these lines I believe could fix the issue.

Edit: I forgot that have_func will define the macro automatically.

[BuonOmo]

Does it make sense to add a LinearRing#ccw? method to the LinearRing feature class? This is not how JTS and libgeos handle it, but it could be convenient to include it.

Definitely, we're in ruby, convenience should always be the main focus!

Your build failed because this is only supported in libgeos version 3.7+. One thought I have on how to handle this is to edit the extconf.rb file to test if the isCCW function exists and define a macro like RGEO_IS_CCW_DEFINED and put an #ifdef statement for it at the top of the analysis.c/h files. Then we could check if the RGeo::Geos::Analysis module is defined inside the ccw? function. These steps may not work perfectly right but something along these lines I believe could fix the issue.

Yes, adding a have_func to the makefile is the move!

[BuonOmo]

Nice ! I've added one last commit, which unfortunately is a test that fails. When a ring has only colinear points, the result is not the same whether we're using geos or a simple factory...

What's your opinion on that one ? IMHO, we could remove the test and add a note about colinear points, but I don't know if there is a spec on that somewhere.

[keithdoggett]

We should either remove it or add some sort of handling for invalid rings. The JTS description of the isCCW method says that it will only work with valid rings and may return incorrect results otherwise. The colinear ring that you made in the test is not valid, so there should be no expectation the result is correct.

I think we should just leave it as is and leave a note about it. It could be annoying to raise an error everytime there's an issue and I think the handling of invalid geoms needs to be re-addressed in v3 anyways.

[BuonOmo]

@keithdoggett After digging a bit more into that, I think this is an issue. Since uses_lenient_assertions has not been set to true, the geometry validation should not have been skipped. I've made a commit on that, and fixing also a potential memory leak. I'm not sure if it should be part of this pr, or maybe another one since it happens on different methods anyway.

require_relative "lib/rgeo"

factory = RGeo::Geos.factory
p1 = factory.point(1, 1)
p2 = factory.point(2, 1)
p3 = factory.point(5, 1)
l = factory.linear_ring([p1, p2, p3, p1])

# Both those assertions should not be true:
l.is_ring? == false
factory.property(:uses_lenient_assertions) == nil

IMO, we should add a notice saying that only if you use lenient_assertions, this is an issue. If this flag is not on, I don't think users would expect an unknown behaviour in case of invalid ring. What do you think ? The commit is a ~100 liner:

Force linear_string validation

diff --git a/ext/geos_c_impl/factory.h b/ext/geos_c_impl/factory.h
index 91ab26c..96b5b04 100644
--- a/ext/geos_c_impl/factory.h
+++ b/ext/geos_c_impl/factory.h
@@ -82,12 +82,12 @@ typedef struct {
   int buffer_resolution;
 } RGeo_FactoryData;
 
-#define RGEO_FACTORYFLAGS_LENIENT_MULTIPOLYGON 1
-#define RGEO_FACTORYFLAGS_SUPPORTS_Z 2
-#define RGEO_FACTORYFLAGS_SUPPORTS_M 4
-#define RGEO_FACTORYFLAGS_SUPPORTS_Z_OR_M 6
-#define RGEO_FACTORYFLAGS_PREPARE_HEURISTIC 8
-
+#define RGEO_FACTORYFLAGS_LENIENT_MULTIPOLYGON 0b00001
+#define RGEO_FACTORYFLAGS_SUPPORTS_Z           0b00010
+#define RGEO_FACTORYFLAGS_SUPPORTS_M           0b00100
+#define RGEO_FACTORYFLAGS_SUPPORTS_Z_OR_M      0b00110
+#define RGEO_FACTORYFLAGS_PREPARE_HEURISTIC    0b01000
+#define RGEO_FACTORYFLAGS_LENIENT_ASSERTIONS   0b10000
 
 /*
   Wrapped structure for Geometry objects.
diff --git a/ext/geos_c_impl/geometry.c b/ext/geos_c_impl/geometry.c
index 3a34f24..dda69b0 100644
--- a/ext/geos_c_impl/geometry.c
+++ b/ext/geos_c_impl/geometry.c
@@ -1047,6 +1047,7 @@ static VALUE method_geometry_invalid_reason(VALUE self)
     str = GEOSisValidReason_r(self_data->geos_context, self_geom);
     if (str) {
       result = rb_str_new2(str);
+      GEOSFree_r(self_data->geos_context, str);
     }
     else {
       result = rb_str_new2("Exception");
diff --git a/ext/geos_c_impl/line_string.c b/ext/geos_c_impl/line_string.c
index 0f5bc64..1d491d6 100644
--- a/ext/geos_c_impl/line_string.c
+++ b/ext/geos_c_impl/line_string.c
@@ -17,6 +17,7 @@
 #include "line_string.h"
 
 #include "coordinates.h"
+#include "errors.h"
 
 RGEO_BEGIN_C
 
@@ -500,21 +501,29 @@ static VALUE cmethod_create_line_string(VALUE module, VALUE factory, VALUE array
 
 static VALUE cmethod_create_linear_ring(VALUE module, VALUE factory, VALUE array)
 {
-  VALUE result;
   GEOSCoordSequence* coord_seq;
   RGeo_FactoryData* factory_data;
   GEOSGeometry* geom;
+  char* invalid_reason;
 
-  result = Qnil;
   coord_seq = coord_seq_from_array(factory, array, 1);
-  if (coord_seq) {
-    factory_data = RGEO_FACTORY_DATA_PTR(factory);
-    geom = GEOSGeom_createLinearRing_r(factory_data->geos_context, coord_seq);
-    if (geom) {
-      result = rgeo_wrap_geos_geometry(factory, geom, factory_data->globals->geos_linear_ring);
-    }
+  if (!coord_seq) { return Qnil; }
+
+  factory_data = RGEO_FACTORY_DATA_PTR(factory);
+  geom = GEOSGeom_createLinearRing_r(factory_data->geos_context, coord_seq);
+  if (!geom) { return Qnil; }
+
+  if ((factory_data->flags & RGEO_FACTORYFLAGS_LENIENT_ASSERTIONS) || GEOSisValid_r(factory_data->geos_context, geom)) {
+    return rgeo_wrap_geos_geometry(factory, geom, factory_data->globals->geos_linear_ring);
+  }
+  invalid_reason = GEOSisValidReason_r(factory_data->geos_context, geom);
+  if (invalid_reason) {
+    // TODO: This is bad, looking for a canonical solution https://stackoverflow.com/questions/64834575
+    GEOSFree_r(factory_data->geos_context, invalid_reason);
+    rb_raise(geos_error, "%s", invalid_reason);
+  } else {
+    rb_raise(rgeo_error, "Unknown error.");
   }
-  return result;
 }
 
 
diff --git a/ext/geos_c_impl/preface.h b/ext/geos_c_impl/preface.h
index 002c924..e59115b 100644
--- a/ext/geos_c_impl/preface.h
+++ b/ext/geos_c_impl/preface.h
@@ -2,6 +2,8 @@
   Preface header for GEOS wrapper
 */
 
+// GEOS Header to ensure we can only use the thread-safe API.
+#define GEOS_USE_ONLY_R_API
 
 #ifdef HAVE_GEOS_C_H
 #ifdef HAVE_GEOSSETSRID_R
diff --git a/lib/rgeo/geos/capi_factory.rb b/lib/rgeo/geos/capi_factory.rb
index a9c36c8..431f8cd 100644
--- a/lib/rgeo/geos/capi_factory.rb
+++ b/lib/rgeo/geos/capi_factory.rb
@@ -33,6 +33,7 @@ module RGeo
             raise Error::UnsupportedOperation, "GEOS cannot support both Z and M coordinates at the same time."
           end
           flags |= 8 unless opts_[:auto_prepare] == :disabled
+          flags |= 16 if opts_[:uses_lenient_assertions]
 
           # Buffer resolution
           buffer_resolution_ = opts_[:buffer_resolution].to_i
@@ -256,8 +257,12 @@ module RGeo
         _buffer_resolution
       end
 
-      # Returns true if this factory is lenient with MultiPolygon assertions
+      def lenient_assertions?
+        _flags & 0x10 != 0
+      end
 
+      # Returns true if this factory is lenient with MultiPolygon assertions
+      # TODO: This should be defined in the c code, using the header RGEO_FACTORYFLAGS_LENIENT_MULTIPOLYGON
       def lenient_multi_polygon_assertions?
         _flags & 0x1 != 0
       end
diff --git a/test/cartesian_analysis_test.rb b/test/cartesian_analysis_test.rb
index 1cdbc11..086fd0d 100644
--- a/test/cartesian_analysis_test.rb
+++ b/test/cartesian_analysis_test.rb
@@ -15,10 +15,15 @@ class CartesianAnalysisTest < Minitest::Test # :nodoc:
     end
 
     def line
-      p1 = @factory.point(1, 1)
-      p2 = @factory.point(2, 1)
-      p3 = @factory.point(5, 1)
-      @factory.line_string([p1, p2, p3, p1])
+require_relative "lib/rgeo"
+
+factory = RGeo::Geos.factory(uses_lenient_assertions: true)
+p1 = factory.point(1, 1)
+p2 = factory.point(2, 1)
+p3 = factory.point(5, 1)
+l = factory.linear_ring([p1, p2, p3, p1])
+l.is_ring? == false
+factory.property(:uses_lenient_assertions) == nil
     end
 
     def clockwise_triangle

[keithdoggett]

@BuonOmo you raise a good point here. The CAPIFactory and FFIFactory do not raise errors on invalid geometry creation like the Cartesian/Geographic/Projection factories. I don't think that this should be included in this PR since this would be a breaking change for anyone using the GEOS factories.

We need to figure out how we want the validation method to work for v3, but I'm of the opinion that allowing instantiation of invalid geometries is fine, but a validity check should be done before certain operations (area, ccw?, etc.) to prevent weird results unless uses_lenient_assertions is set. I think this is closer to how GEOS/JTS behaves, but we can discuss that later.

[BuonOmo]

@keithdoggett, ok correct this should not be part of this PR. I'll open an issue soon.

I think we should discuss about planning for v3 in our next monthly :)

I've added a note to the PR, I'll let you review and I think we'll be good to go !

[keithdoggett]

@BuonOmo this looks good now. I think the last thing we should do is remove the test_ccw_p_line test from the CartesianAnalysisTest module since it fails when GEOS is version 3.7+ and we've determined that using this method with invalid geoms will have weird behavior.

[BuonOmo]

@keithdoggett merged ! I think we could try and merge #230 and then make a minor. Does it look or for you ?