[postgis-devel] new constrained delauney triangulation

Mon Mar 11 10:45:59 PDT 2019

Hi,

On Mon, Mar 11, 2019 at 1:21 AM Nicklas Avén <nicklas.aven at jordogskog.no>
wrote:

> 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.
>

PostGIS is GPL2, your code is GPL3. Please change before I look at it so I
can still do something for PostGIS in Delaunay :)

> 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.
>

This is cool. There's a hidden function in SFCGAL that does the CDT:
https://trac.osgeo.org/postgis/ticket/4198 - but in the end it's the same
algo as Tesselate?

CDT is needed for properly generating isochrones:
https://www.patreon.com/posts/isochrones-are-20933638 - so I would love to
see CDT in core replacing ST_DelaunayTriangles and helping to iron the
issues!

CDT will also help replacing ST_ConcaveHull with a lot faster thing.

 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 believe this is also useful for compressed versions of TINs, that should
be stored as half edges.

> 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.
>

TRIANGLE is already part of geometry. What if it was just supported in
serialization as is?

>
> 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.
>

What do you use for inCircle predicate?

>
> 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.
>

Would be cool if it can be made a part of PostGIS. Then we can utilize the
coverage infrastructure making sure we can generate tests.

Also, it is a single-geometry-input thing? Then we can make a fuzzer for it
and let Fuzzie look for the crashers.

>
> 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
>
> _______________________________________________
> postgis-devel mailing list
> postgis-devel at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/postgis-devel

-- 
Darafei Praliaskouski
Support me: http://patreon.com/komzpa
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-devel/attachments/20190311/e353fd65/attachment-0001.html>