[postgis-devel] new constrained delauney triangulation

[Nicklas Avén]

Hello


I have been working on a delauney-triangulation code.

Link to the code in the bottom

I need it booth for packaging maps to TilelessMap, but I also need it to 
be able to implement editing in the app.

For packaging I will use it in PostGIS.


The results so far looks so promising so I might suggest that it maybe 
could go into PostGIS core. Now it is in an extension.


All core logic comes from this paper

https://www.newcastle.edu.au/__data/assets/pdf_file/0019/22519/23_A-fast-algortithm-for-generating-constrained-Delaunay-triangulations.pdf


What I need is constrained triangulation, so I get only triangles in the 
polygons, not in between.

 From what I know PostGIS then only have ST_Tesselate, which is horribly 
slow and very, very picky with what it gets. Since it is more picky than 
ST_Isvalid it is very difficult to handle.


ST_DelaunayTriangles is quite fast, but cannot constrain (if I haven't 
missed something).


Anyway, what I need is also the possibility to get the result as indexes 
into the original point arrays. That is what the GPU wants. I haven't 
implemented that yet, just because I couldn't wrap my head around how to 
create arrays in PostgreSQL.


I also have a plan to suggest an extension to the twkb-format so the 
triangel indexes can be a part of a special geometry-type.

Now in TilelessMap I encode the triangle indexes as 3D pointsand keep 
them in a separate column. But with everything in a twkb geometry we 
have a format ready for pushing into a gpu. Just decode the points and 
tri-index into float arrays and serve to the GPU and it will be happy.

That is what I do in Tilelessmap.



Back to the triangulation code.

It seems pretty fast.

This is from the largest polygon I have got with 249667 vertex points in 
one single polygon (st_dump gives 1 row):

With ST_DelaunayTriangles unconstrained

test=# EXPLAIN ANALYZE select ST_DelaunayTriangles(geom) geom from 
my_riks WHERE ID = 14681;
                                                          QUERY PLAN
-----------------------------------------------------------------------------------------------------------------------------
  Index Scan using my_riks_pkey on my_riks  (cost=0.29..8.30 rows=1 
width=32) (actual time=2381.912..2386.171 rows=1 loops=1)
    Index Cond: (id = 14681)
  Planning time: 0.056 ms
  Execution time: 2386.185 ms
(4 rows)

Time: 2386.537 ms (00:02.387)


With the new code constrained

test=# EXPLAIN ANALYZE select triangulate(geom) geom from my_riks WHERE 
ID = 14681;
                                                         QUERY PLAN
---------------------------------------------------------------------------------------------------------------------------
  Index Scan using my_riks_pkey on my_riks  (cost=0.29..8.30 rows=1 
width=32) (actual time=522.360..527.080 rows=1 loops=1)
    Index Cond: (id = 14681)
  Planning time: 0.036 ms
  Execution time: 527.095 ms
(4 rows)

Time: 527.471 ms

Getting unconstrained from the new code is about the same speed, which 
is a little strange since the job to constrain is all done after the 
first round of Delauney-triangulation

test=# EXPLAIN ANALYZE select triangulate(geom,0) geom from my_riks 
WHERE ID = 14681;
                                                         QUERY PLAN
---------------------------------------------------------------------------------------------------------------------------
  Index Scan using my_riks_pkey on my_riks  (cost=0.29..8.30 rows=1 
width=32) (actual time=516.255..522.467 rows=1 loops=1)
    Index Cond: (id = 14681)
  Planning time: 0.039 ms
  Execution time: 522.482 ms
(4 rows)

Time: 522.793 ms

For smaller geometries the new code is about twice as fast as 
ST_DelaunayTriangles

St_Tessalate cannot be compared on larger polygons, it takes minutes.


The new code now handles some invalidities. There are also possibilities 
to handle more.

Right now the code doesn't always stop gracefully if you find a polygon 
it cannot handle.

I have even had the code unstoppable, so I had to shut down the machine.

So, be warned if yo want to test the code, it is a lot of error handling 
to do.

But if you can spend a reboot at worst, I would be glad if you can test 
it. I have tested on a few million polygons, which seems to work great.

The code also needs a lot of commenting to be easy to follow (and clean 
up and refactoring). I will do that especially if there is interest in 
the code.

To build as extension run make in folder postgis_extension

The code is here:

https://gitlab.com/nicklasaven/triangulation


Thanks


Nicklas Avén