[Patch] Add GeometryPrecisionReducer to C API

[smani]

The attached patch adds methods for using the GeometryPrecisionReducer to the C API. I'm using these methods in QGIS to achieve better results when performing operations on geometries which occasionally suffer from precision issues.

Related to the subject: I noticed that geos::precision::GeometryPrecisionReducer::targetPM is a const-reference to a geos::geom::PrecisionModel. Unless there are specific reasons for this, API-wise it would probably be cleaner to just store a copy of the PrecisionModel (since it is a tiny class) - this is what one would most likely expect to happen when seeing the "const geom::PrecisionModel &pm" in the constructor signature.

As far as the patch is concerned, storing the targetPM as a copy would remove the need to keeping the GEOSPrecisionModel around.

[smani]

Patch

[strk]

If I'm not wrong the geometries resulting from this kind of reduction would retain their own GeometryFactory which would point to the default (floating) precision model. Can you confirm ?

There's no direct connection between a Geometry and a PrecisionModel, rather the Geometry points to a GeometryFactory which points to a PrecisionModel.

What we're looking for is an API to allow for storing the precision within the geometries, so that operations on fixed-precision geometries will be carried with that precision. There have been some discussion about this for a GSoC project: http://trac.osgeo.org/geos/wiki/GSoC/CAPI_PrecisionModel and see also needs from Shapely here: ​https://github.com/Toblerity/Shapely/issues/110

[smani]

Yes indeed the geometry returned by GeometryPrecisionReducer::reduce retains the GeometryFactory of the input geometry.

The API you describe basically means exposing the GeometryFactory in the C API, along with methods to create a GeometryFactory with a given PrecisionModel, and then adding methods for all geometry types to create them with the specified GeometryFactory, right?

But considering a situation where the geometry is already given, exposing the GeometryPrecisionReducer would still be needed, since AFAICS it is the only way to reduce the precision of an already given geometry?

[strk]

Yes the GeometryPrecisionReducer would still need to be exposed, but would take a GeometryFactory rather than a PrecisionModel. We've been discussing a big about possible interfaces to manage these factories, or ways to hide them behind the ContextHandle.

May be better to discuss this on the list, best if with a proposal.

[smani]

I'll try and come up with a proposal soonish.

[smani]

Posted to mailing list: ​http://lists.osgeo.org/pipermail/geos-devel/2014-December/007076.html

[strk]

Sandro: is there an updated patch, or pull request, with the outcome of the discussion on the list ? I understand you're now using a C wrapper in QGIS to use that functionality, are there tests showing how using the precision reduction code worked with overlay operations ?

[smani]

The patch attached to this ticket is pretty much what I'm using in QGIS. I can create a pull request with that patch if that helps. The discussion on the mailing list deviated away from that patch towards a more implicit exposure of the precision reduction functionality, but the discussion stalled away with various issues with that idea still open. The patch attached to this ticket (and what is used in QGIS) is basically a straight 1:1 port of the C++ api to C.

[strk]

From what I can tell the C++ API is wider than the one you exposed to C. For example it has two kind of construction for the GeometryPrecisionReducer class. One takes a PrecisionModel, and the other takes a GeometryFactory. The one taking a GeometryFactory will store the GeometryFactory reference in the output Geometry, effectively making it "tagged" as having the target "Precision", the other (the only you expose to C) will just use the given PrecisionModel as a grid to reduce coordinates to, but the output geometry will not be different from before.

If what I described above is confirmed, the same thing could be exposed as a simple "GEOSGeometry_snapToGrid(<some_grid_config>)" function, with no need to expose a PrecisionModel _at_all_. This is surely not what I had in mind when thinking about exposing the PrecisionModel.

And I'd also argue that adding a dependency on GEOS-C++ library to QGIS for just snapping to a grid is probably overkill.

[smani]

Right, 1:1 was more like 1:1 of the functionality I needed.

It's been a while since I did the original experiments, but I think the important thing was that the actual precision model of the geometries changed in order to have the analysis operations run with that precision model.

Ultimately what I needed was a way to reduce the chance of getting these tiny spikes caused by imprecisions when for instance intersecting geometries, or cases where unions of neighboring geometries failed due to geos thinking they were disjoint (again due to precision issues).

I agree that having QGIS link agains the entire C++ GEOS library is an overkill, but at present it was the only way to get it merged into master. As soon as GEOS offers such functionality in the C library I'm happy to port the QGIS code to whatever API ultimately is implemented and remove the dependency again.

[strk]

What I'm saying is that from what I can read in code, I don't think you are having analisys operation run with the custom precision model. A test showing they are would be useful, but I found no test under QGIS.

Such a test would use input geometries that produce the _same_ coordinate values (same HEXWKB) after coming out from PrecisionReducer but produce spikes (or other issues) when passed directly to those analisys operation and do not produce them when passed after the PrecisionReduced threatment.

[smani]

If I recall correctly the reduce operation returns a geometry with a corresponding precision model. And I thought this to be sufficient to have the analysis operation run in reduced precision mode?

I'll see whether I can dig up my test geometries again and produce a test. But I definitely recall seeing a difference with and without the reduce operation.

[strk]

Perturbating the input (such as rounding its coordinate values) is often enough to change the outcome of operations. A good test would be computing the intersection between two segments long 10 times the grid cell size, crossing and having the endpoints at "scale" distance each other. If the computation is performed at the target precision model the intersection point will on the grid, otherwise it'll be in between.

[strk]

Here's a simple test (just checked with XMLTester):

A: LINESTRING(0 -10, 2 10) 
B: LINESTRING(2 -10, 0 10)
Operation: intersection
Result with FLOATING precision: POINT(1 0)
Result with FIXED precision (scale="0.5"): POINT(2 0)

Please see if you can reproduce that via the C wrapper of your patch, which would need to be accompanied by a unit test to be accepted, in any case.

[smani]

Ok, bad news:

​https://github.com/manisandro/QGIS/commit/1738337ddea6d3242abd22bd8fc3b9d6b2ce3b2b

    ********* Start testing of TestQgsGeos *********
    Config: Using QTest library 4.8.7, Qt 4.8.7
    PASS   : TestQgsGeos::initTestCase()
    FAIL!  : TestQgsGeos::lineIntersection(fixed) Compared values are not the same
       Actual (QString(resultwkt)): POINT (1.0000000000000000 0.0000000000000000)
       Expected (wkt_result): POINT (2.0000000000000000 0.0000000000000000)
       Loc: [/home/sandro/Documents/Devel/QGIS/qgis-master-2/tests/src/core/testqgsgeos.cpp(94)]
    PASS   : TestQgsGeos::cleanupTestCase()
    Totals: 2 passed, 1 failed, 0 skipped
    ********* Finished testing of TestQgsGeos *********

So intersection is running in full precision mode.

From ​http://lists.osgeo.org/pipermail/geos-devel/2015-January/007086.html I had interpreted that GEOSGeometryPrecisionReducer_reduce would change the geometry factory of the reduced geometry. Did I get that wrong, or is more needed to have the operations run in reduced precision mode? Because that is what I'm ultimately interested in, the simple snap-to-grid is a start but not the full solution (and that could indeed be implemented stand-alone, without pulling in the GEOS C++ API as you pointed out on the mailing list).

[strk]

GeometryPrecisionReducer can change or not change the geometry factory of the reduced geometry, depending on how it is constructed. If you pass a new GeometryFactory to the constructor, it will use it for the output geometry, otherwise it will not. Your implementation is passing a PrecisionModel rather than a pointer to a GeometryFactory, that's the problem.

This brings us to the need of exposing a GeometryFactory, more than a PrecisionModel.

[strk]

From what I can tell the discussion about GeomtryFactory management stalled here: ​https://lists.osgeo.org/pipermail/geos-devel/2015-January/007097.html

Should we re-take that on the mailing list ?

[smani]

First thing to decide IMO is how to proceed for QGIS 2.10: I can fix my geos_extra to construct the precision reducer correctly, get the fix commited, and then as far as QGIS 2.10 is concerned things are ok and then we can take the time to design a proper C API for GEOS. If this is okay with you, I'll do that and then yes, we can continue where we left of on the mailing list.

[smani]

• QGIS 2.12 clearly

[strk]

QGIS matters are better discussed on the qgis mailing list, let's use it: ​https://lists.osgeo.org/pipermail/qgis-developer/2015-October/039418.html

[strk]

A GEOSGeom_setPrecision function is in as of r4092. Tests are welcome, and if you find any issue or have enhancement proposals, please file a new ticket.

[strk]

Ticket retargeted after milestone deleted