pgr_pointsaspolygon returns a self intersecting polygon which does not cover underlying points

[nilgundag]

I am using pgrouting 2.0 library to create an alpha shape around my set of points. It works nearly perfect but I have a case where pgr_pointsaspolygon returns a self intersecting polygon which does not really cover the underlying points in the table:

This is how i call the pgr_pointsaspolygon function :

SELECT 1 as a,* FROM  pgr_pointsaspolygon('SELECT id, ST_X(the_geom) AS x, ST_Y(the_geom) AS y FROM test')

[dkastl]

Thank you for the bug report.

Would you be able to provide sample data, that makes this case re-producable?

[nilgundag]

You can find dump of test2 database in the following file:

https://gist.github.com/nilgundag/868e4532a3780a88effd/raw/9865c398e569585f191f5e575f2a6ec9b8b70660/test2.sql

And the geometry would be: SELECT * FROM pgr_pointsaspolygon('SELECT * FROM test2')

The points have SRID 4326.

[sanak]

This issue's cause seems to be CGAL computation part in core alphashape function. After changing "src/driving_distance/src/alpha_drivedist.cpp" line 150 as follows, I got a valid result. (Please see attached image.) But I don't know whether removing "*6" is valid fix or not.

  //A.set_alpha(*A.find_optimal_alpha(1)*6);
  Alpha_iterator opt = A.find_optimal_alpha(1);
  A.set_alpha(*opt);

[dkastl]

@sanak Thanks! Some time ago I already tried to figure out what *6 was supposed to be for.

Just the link to code line: https://github.com/pgRouting/pgrouting/blob/develop/src/driving_distance/src/alpha_drivedist.cpp#L150

[nilgundag]

So, is it safe for me to just remove *6 there or should I write tests for different cases and compare the results to make sure all other cases still works fine?

[dkastl]

Well, if you could test some different cases and see if removing *6 causes any issues, that would be great.

[sanak]

Okay, I will ask it to CGAL ML in this weekend.

[nilgundag]

@sanak thank you very much. I can write random tests, but they probably wont fully cover all possible scenerios.

[sanak]

Sorry, the cause was not *6, but CGAL Alpha shape result edges ordering part (find_next_edge function) in case of existing holes in result. https://github.com/pgRouting/pgrouting/blob/master/src/driving_distance/src/alpha_drivedist.cpp#L88

I tried 2 patterns("Check used segment index" and "Join triangles") to consider holes in result with the following example code. ("Join triangles" pattern was very slow, so I also tried "Check used segment index" pattern.)

Example code: https://gist.github.com/sanak/9156747 Result images:

And it will be necessary to think about the followings.

• Return inner rings if those exist. (Need to change pgr_alphashape return type. - http://docs.pgrouting.org/dev/src/driving_distance/doc/dd_alphashape.html )
• Path an arbitrary alpha value to pgr_alphashape function.

[woodbri]

I'm not familiar with the CGAL algorithms but I thought we could eliminate CGAL by using the following algorithm.

1. compute the driving distance costs to each node in the network
2. use postgis to create a Delaney triangulation by inserting the nodes (note 1)
3. if we use the distance as the z-value, we basically create a bowl with the low spot at the start node and walls of the bowl climbing the further you go outward.
4. intersect the bowl with horizontal plans with the z value of the contour(s) you want (note 2)
5. pass the intersection results to postgis build area to create polygons

notes:

1. I can't remember if the Delaney triangulation works in 3D or 2D only. It 2D, you need to then take a second pass over the triangles and add the z value
2. st_intersection may be slow, in which case we can create a very fast plan-triangle intersector. this is trivial code and will be much faster then st_intersection

[sanak]

@woodbri Thanks for comment ! I still don't know Alpha shape algorithm internal, but if eliminating CGAL is possible, it will be good to maintain this project code. Implementing Alpha shape function to GEOS (or JTS) may be also interesting.

[sanak]

About returning holes or multiple polygon, ST_BuildArea may be available in current pgr_pointsAsPolygon function with returning NULL values' row as ring separator from pgr_alphashape. I will try to check it later. http://postgis.net/docs/ST_BuildArea.html

[woodbri]

This might need to get pushed out to the 2.2.0 Milestone.

[cvvergara]

@sanak I am revising all the issues that are to be handled in the 2.2 version, and I was wondering if this one was solved on the 2.1 version.

[dkastl]

Since it was merged it was probably resolved. We may want to add a test case though ... new issue for that and close this one?

[cvvergara]

ok. Marking as fixed on develop due to the fact that we milestone 2.1.0 has being closed. Closing.