Motivation

Oslandia and IGN are involved in developping 3D spatial operations for PostGIS, with a use case focus on city modeling (CIM), and with a 3-year partial fundings from the European Union (FEDER).

CGAL

After PostGIS Paris Code Sprint, we really considered using the CGAL library, which already implements several 3D algorithms.
CGAL is also able to guarantee exact computations through the use of arbitrary precision numbers (if we choose the right CGAL kernel)

The question became: can we use CGAL as an underlying geometric computation framework for 3D operations with the additional benefit of exact computations ?

SFCGAL

We started by building an OGC Simple Feature compliant library on top of CGAL, called SFCGAL.

It is not yet complete, but main design elements have been already fixed. The class hierarchy is pretty standard and lots of commonalities may be found with other libs, like GEOS.

However, SFCGAL stores geometry coordinates using arbitrary precision numbers, as CGAL does.

On top of this class hierarchy, some 2D and 3D operations have already been implemented :

• 2D and 3D Intersection
• 2D and 3D Intersection test (e.g ST_Intersects)
• 2D and 3D Convex Hull
• 2D and 3D Triangulation
• 2D and 3D distances
• 3D Extrusion (from a 2D geometry)
• 2D Minkowski Sum

PostGIS-sfcgal

A PostGIS branch is also available including SFCGAL handling. It exposes SFCGAL functions to PostGIS.

PostGIS - SFCGAL architecture

All the new functions provided are isolated in a ’sfcgal’ schema.

Here is a list of already available functions :

• ST_Intersection / ST_3DIntersection (supporting all types of geometries, including solids)
• ST_Intersects / ST_3DIntersects (supporting all types of geometries, including solids)
• ST_ConvexHull / ST_3DConvexHull
• ST_Distance / ST_3DDistance
• ST_DelaunayTriangles / ST_3DDelaunayTriangles
• ST_Extrude
• (in progress) ST_Buffer

Comparison to GEOS

A simple benchmark helps us to compare GEOS-based spatial operations with SFCGAL-based ones. Comparison can only be made for 2D operations, since GEOS is mainly 2D. (Python benchmark script)

SFCGAL’s speed is comparable to GEOS’ speed (sometimes even a bit better).
However, for few non-native CGAL operations (area for instance), we rely on our own (quick) implementation, which still lacks optimizations, they are then not (yet) as fast as GEOS-based functions.

So using CGAL with maximum precision for native spatial operations could be done with keeping the same kind of performances than GEOS.
Put differently, with CGAL we gain arbitrary precision without performance loss.

Support for exact geometries

But we do not take yet the full benefit of an exact computation framework. For instance, consider this post from Paul on OpenGeo blog :

SELECT
  ST_Intersects(
    ST_Intersection(
      'LINESTRING(0 0,2 1)'::geometry,
      'LINESTRING(1 0,0 1)'::geometry),
    'LINESTRING(0 0,2 1)'::geometry);

It should answers ‘true’. However it is not the case with current PostGIS implementation, because coordinates are represented using floating points between two consecutive calls.

The SFCGAL library does not enforce floating coordinates and keep their exact representations.
But obviously we loose this precision when switching back and forth to double on each (de)serialization calls (to and from GSERIALIZED).

In order to take benefit of the exact representation within nested functions, we added a new (de)serialization format that supports arbitrary precision called ‘exact_geometry’.

Reusing the previous example, we could now write :

db=# SELECT
  sfcgal.ST_Intersects(
    sfcgal.ST_Intersection(
       'LINESTRING(0 0,2 1)'::exact_geometry,
       'LINESTRING(1 0,0 1)'::exact_geometry),
    'LINESTRING(0 0,2 1)'::exact_geometry);
st_intersection
-----------------
POINT(2/3 1/3)
(1 row)

Notice the new ‘a/b’ rational format that figures internal rational representation. SFCGAL is also able to parse this extended WKT.

Referenced geometries

When calling nested functions, we spend time serializing and deserializing on each function call. Even if the geometry output from the first function is really the same than the input of the second one.

Why do we need to serialize and then deserialize the same geometry ?

If you consider that (de)serialization from liblwgeom to CGAL exact representation really takes time, nested calls are performance evil (yeah no less).

So it implies a new geometry type called ‘ref_geometry’.
Ref_geometry only copy their memory address (converted to some integers) from functions to functions.

This implies:

• destruction of referenced geometries is postponed,
• they are not designed to live outside the query they first appeared in (and have to be converted back to a regular geometry type in order to be stored, using for instance sfcgal.ST_Geometry()).

Implementation considerations

The first point means, in PostgreSQL terms, that such objects must be attached to the correct memory context in order to be deleted at the ‘right’ time.

If the context is deleted too early, terminal functions might not be able to access our previously allocated geometries.
On the other hand, choosing a too static memory context will let objects dangling until late and will explode memory consumption.
The current choice (attaching to the current ExprContext when possible, to MessageContext otherwise) gives good results, but this would have to pass advanced tests where different kinds of evaluation contexts are produced by the postgresql execution engine.

The other difficulty at this point is that we are dealing with C++ objects.
Allocating something with palloc() allows it to be automatically deleted on memory context destruction. However no C++ destructors are called in that case.

But, memory contexts have support for user-supplied reset and deletion functions. Referenced geometries are then created on a custom child memory context. When this custom context is deleted (by its parent), C++ objects can be safely deleted.

Still with the same example, we get:

db=# SELECT sfcgal.ST_Intersection(
    'LINESTRING(0 0,2 1)'::ref_geometry,
    'LINESTRING(1 0,0 1)'::ref_geometry);
st_intersection
-----------------
POINT(2/3 1/3)
(1 row)

And, in order to illustrate the temporary nature of referenced geometries:

db=# CREATE TABLE temp_geo AS SELECT 'POINT(0 0)'::ref_geometry;
SELECT 1
db=# SELECT * FROM temp_geo;
NOTICE:  Referenced geometries must not be stored
ref_geometry
--------------
-deleted-
(1 row)

Referenced geometries are here made to be used with the SFCGAL::Geometry class. However, the mechanism is generic enough to be applied to other (C++) types.

Serialization comparison

A quick bench (4 nested calls to a function that does no processing, i.e sfcgal.ST_Copy), is useful to show performance differences between geometry representations.

Referenced geometries are way faster than (SFCGAL) serialized geometries and even still faster than native GSERIALIZED geometries.

Notice also that serialization of exact geometries is still slower than native serialization through GSERIALIZED. The overhead might be due to a lack of optimization on this part, but we feel this is due to the storage of arbitrary precision numbers.

Serialization benchmark: native vs. ref_geometry

Serialization benchmark: native vs. exact_geometry

Regarding memory consumption, the three approaches consume a comparable amount of memory.

How to test ?

To test SFCGAL and our postgis branch:

• Get and install CGAL 4.1 from cgal.org
• Get and install SFCGAL (setting ‘SFCGAL_USE_STATIC_LIBS’ to OFF)
• Get and install our postgis-sfcgal branch, (run ./configure with the ‘–with-sfcgal’ flag, refer to the README.md file for the other flags)
• Alternatively, you can apply a patch against PostGIS trunk (revision 11054 at the time of writing)
• Run postgis.sql on a database
• Try a 3D intersection between two cubes with (you can display the resulting WKT with the ‘DemoPlugin’ from SFCGAL’s viewer)

  WITH unit_cube AS (
      SELECT
        sfcgal.ST_MakeSolid(sfcgal.ST_Extrude(
          sfcgal.ST_Extrude(
            sfcgal.ST_Extrude('POINT(0 0)'::geometry, 1, 0, 0),
            0, 1, 0),
          0, 0, 1)
        )
      AS solid
    )
  SELECT ST_AsText(sfcgal.ST_3DIntersection( solid, ST_Translate( solid, 0.5, 0.5, 0 )))
  FROM unit_cube;
  

• Try using ‘exact_geometry’ and ‘ref_geometry’ types through ‘::’ cast or through conversion functions (sfcgal.ST_RefGeometry() and sfcgal.ST_ExactGeometry())

3D viewer

To render 3D operations, SFCGAL is shipped with a really basic 3D viewer based on Open Scene Graph and QT. It is able to connect to a PostgreSQL base. Usage is only for dev/debug purpose, consider it also as alpha.

SFCGAL viewer demo

We are also investigating the use of QGIS (with the Globe plugin) to serve as a 3D viewer of spatial queries.

PostGIS 3D demo

Discussion

Integration

We would like, as you could guess, to see PostGIS SFCGAL functions included into PostGIS.

Points to discuss/validate:

• As a first step, PostGIS-sfcgal support could be added in a quite near future, into PostGIS trunk, since its compilation is optional, it does not interfere with existing PostGIS functions through the use of a separate schema.
• The ‘exact_geometry’ type is important to keep the ability to nest several function calls without loosing geometry precision (and so 3D topology).
  On the long term however, what about a single geometry type that would optionally support arbitrary precision numbers ? This would probably need some efforts around GSERIALIZED functions. To be discussed.
• The ‘ref_geometry’ may first be considered optional. However it allows to keep performance with nested functions using SFCGAL. Current implementation is still to be considered a bit experimental and would have to pass more tests to be fully qualified, before a PostGIS integration. Any feedback related to C++ objects handling within the PostgreSQL framework is welcome.
• A quite minor point, is related to SFCGAL wrapper code written in C++. Some may find that it breaks the consistency of pure C PostGIS sources, but it was easier for early design and tests. Is keeping direct C++ calls from PostGIS not a big deal or, on the contrary, a C API is a need, for the sake of consistency ?

Roadmap

We still are involved in this project till (at least) end of 2014, so the current roadmap, for the next months, is:

• Several spatial processing are still missing. Especially the complete set of boolean predicates and constructions (which are present in CGAL, so more an interface task)
• There is no mechanism of geometry caching like GEOS has. This would be of interest regarding performances
• Fine tuning of some CGAL algorithms used are still needed.
• Hopefully, progress with PostGIS trunk integration

Feel free to share your thoughts, remarks and critics.