ST_DelaunayTriangles

[strk]

Expose delaunay triangulation to postgis. Interface details to be decided. This ticket is mostly to gather use cases at this stage.

[woodbri]

Adding some use cases for discussion that I presented on the list.

I wrote code years ago in C that would take a cloud of 3D points insert them into a Delaunay triangular surface, then slice the surface triangles with a plane(s) into intersecting edges and then composite rings from the edges. It seems like we are very close to being able to do this in PostGIS. I'm not sure if we have the Delaunay triangularization yet.

I'm thinking of something like:

select b.z, a.* from
  st_buildarea(
    st_collect(
      intersection(
        delaunay_triangles(
          select the_geom from points),
          st_translate(
            st_expand(select extents(the_geom) from points, 0.1),
            0.0, 0.0, b.z
          )
        )
      )
    )
  ) as a,
  (select 1.0 * generate_series(2, 10, 2) as z) as b
 group by b.z;

So I'm thinking of something like this is the use case, where this would create contour rings for Z at 2, 4, 6, 8, and 10 based on the triangulated surface created from the points. My SQL is probably broken, but it is only to present the idea.

This use case is very broad in its applicability to real world cases. Think of the points as XY as in location of a sensor or some other point of measurement and the Z as the measurement value. Creating a 3D triangulated surface and cutting it with z-planes creates contours of equal value results. Create a vertical plane(s) and you are generating elevation profiles. and if you have 2d path like a bicycle route, a GPS track, etc and extrude that into a vertical panel the intersects the surface, you intersection becomes the measurement profile of the path. If the surface is elevation, then you have an elevation profile, if the measurement was chemical concentrations in a water body and the track is fish, then the profile because an exposure profile for the fish.

So some use cases of this are:

1. driving distance - compute the costs to the node using Dijkstra and the XY is the node location and Z is the cost.
2. weather air-presure - XY is sensor location Z is the bars, rings are isobars
3. XY location, Z=altitude, height contours one of the more obvious cases

And I think that being able to build a pipeline of functions like I offered above would be extremely useful. in all of these use cases.

While I do not need it specifically at the moment, anyone trying to render 3D landscapes might want to be able to extract the triangles or might want to be able to apply hillshading to the the triangulated surface and create a postGIS raster as a result. For the hillshading task I believe they would like to get the average normal at the each triangle node, this is needed for both Gouraud and Phong shading algorithms, being able to extract that might be of interest in the future.

[strk]

JTS 1.12 has Delaunay Triangulation package which should be ported to GEOS. See  http://trac.osgeo.org/geos/ticket/487

This ticket is in search of a sponsor

[strk]

It's actually from JTS-1.11:  http://lin-ear-th-inking.blogspot.com/2009/04/delaunay-triangulation-in-jts-111.html

[strk]

News from the GEOS side of things: the C++ port is ready in trunk (to become 3.4.0) and it's now time to expose it trough the C-API. Reading the C++ interface there are 2 possible outputs: a collection of lines or a collection of polygons. The former is smaller, and I guess you could form the latter by passing the former to ST_Polygonize.

So, what do you think would be more appropriate to return ?

[strk]

For the record, the ticket for the C-API side of GEOS is  http://trac.osgeo.org/geos/ticket/565

[strk]

So I ended up exposing an API which allows both collection-of-polygons and multlinestring outputs, so it only needs to be decided how to access the functionality from postgis. How should the function look like to request one or the other, and what should the default be (if any).

Note that the function also accept a tolerance.

[woodbri]

It seems to me that having something like:

select the_geom 
  from delaunay_triangles(sql_for_points text, tolerance double precision, multilinestring_out boolean);
select the_geom 
  from delaunay_triangles(sql_for_points text, tolerance double precision);

where the function returns a set of geometries based on multilinestring_out flag. If true you get multilinestrings if false or absent you get polygons. You mentioned a tolerance so I tossed that in, but what is it for?

I like the set returning function because it is easy to throw ST_Collection() around it if you want a collection but much harder to explode it is you do not want that and it support your idea of pipelining.

BTW. Nice work on this!

[strk]

Given that GEOS gives us the already-collected representation there's really no gain in a set-returning function.

The flag could be fine but maybe not a boolean but an integer, to keep the door open for more flags in the future (ala GML output functions)

r9993 adds the liblwgeom part (with an int flag)

[strk]

Have fun with r9994 -- adds SQL entry point, SQL level test and documentation.

[pracine]

This is definitely a nice feature but how can I build a Delaunay triangulation from a table of, say, 1 000 000 points?

[strk]

select ST_DelaunayTriangles(ST_Collect(geom)) FROM one_million_points;