Fwd: Increasing the maximum number of points in the spatial.Delaunay triangulation
I originally wrote this to the scipy-user group, but I think it might be more appropriate here. I am working on an image processing algorithm that relies on the Delaunay triangulation of a 3d image. Unfortunately, I need to work on images with more than 2**24 points that go into the Delaunay triangulation. This is problematic because the Qhull routine labels the input vertices with a 24-bit label. I have tried editing the source code for qhull here in libqhull.h, lines 352, 380 and 663 - 664: unsigned id:32; /* unique identifier, =>room for 8 flags, bit field matches qh.ridge_id */ ... unsigned id:32; /* unique identifier, bit field matches qh.vertex_id */ ... unsigned ridge_id:32; /* ID of next, new ridge from newridge() */ unsigned vertex_id:32; /* ID of next, new vertex from newvertex() */ poly.c, lines 569, 1019 - 1023: long ridgeid; ... if (qh ridge_id == 0xFFFFFFFF) { qh_fprintf(qh ferr, 7074, "\ qhull warning: more than %d ridges. ID field overflows and two ridges\n\ may have the same identifier. Otherwise output ok.\n", 0xFFFFFFFF); and poly2.c, lines 2277 - 2279: if (qh vertex_id == 0xFFFFFFFF) { qh_fprintf(qh ferr, 6159, "qhull error: more than %d vertices. ID field overflows and two vertices\n\ may have the same identifier. Vertices will not be sorted correctly.\n", 0xFFFFFFFF); This compiles fine, but so far, this produces a memory leak (as far as I can tell) when I put more than 2**24 points into the algorithm. I have emailed Brad Barber, who indicated that this changing the ridgeid field was possible. I am also working on a divide and conquer approach, but am having trouble deleting the old Delaunay object before instantiating the next one. I can send that code along if there is no good way to increase the number of points in the algorithm. I think it would be useful to have an algorithm that can accept more points, given the increased speed and processing power of modern computers. Thanks for your help.
22.06.2015, 19:48, Cantwell Carson kirjoitti: [clip]
This compiles fine, but so far, this produces a memory leak (as far as I can tell) when I put more than 2**24 points into the algorithm. I have emailed Brad Barber, who indicated that this changing the ridgeid field was possible. I am also working on a divide and conquer approach, but am having trouble deleting the old Delaunay object before instantiating the next one. I can send that code along if there is no good way to increase the number of points in the algorithm.
I think it would be useful to have an algorithm that can accept more points, given the increased speed and processing power of modern computers. Thanks for your help.
I can't immediately help you with the Qhull internals, but other comment: Your best bet with that would probably be to try to get your patch included in Qhull itself. The Qhull code is shipped with Scipy only for convenience, and is unmodified from the 2012.1 version --- shipping a patched version might not be a very good thing if the Qhull upstream is still responsive. Of course, Qhull is not the only game in town, at least for low-dimensional triangulations, and it could be interesting to have more featureful alternatives for those cases.
Are there other ideas for what could ship with scipy? I have tried CGAL and mayavi and they function, but no better than Qhull (IMHO). Matplotlib also has a routine, but I doubt the redundancy in algorithm would be beneficial. For example, here is a CGAL wrapper that I tried, along with associated comments about the speed: def Delaunay_wrapper(points): ## Designed for integer point coordinates. from CGAL.CGAL_Kernel import Point_3 from CGAL.CGAL_Triangulation_3 import Delaunay_triangulation_3 ## Carry out the triangulation (very fast) L=[] for coordinates in points: L.append(Point_3(*coordinates)) T=Delaunay_triangulation_3(L) ## Read the information at the pointers T into pts_dict (very slow) pts_dict = {} for idx in range(4): pts_dict[idx] = [] for cell in T.finite_cells(): pts_dict[idx].append((cell.vertex(idx).point())) pts_dict[idx] = [map(float, str(x).split(' ')) for x in pts_dict[idx]] ## reslice the list to group the cells together cell_list = zip(pts_dict[0],pts_dict[1],pts_dict[2],pts_dict[3]) ## return a list of n cells comprised of 4 points with 3 coords return cell_list I was using the SWIG wrapper to expose the CGAL methods to Python. The triangulation is very fast, but writing them into a Python object is very slow. If this could be sped up, there might be a benefit to making it an option in Scipy. As for the number of points, I recommend that we insert a warning or error when the number of points in Qhull algorithm exceeds 2**24 so that the method can exit without crashing, instead of killing the whole process as it does now. On Mon, Jun 22, 2015 at 2:54 PM, Pauli Virtanen <pav@iki.fi> wrote:
22.06.2015, 19:48, Cantwell Carson kirjoitti: [clip]
This compiles fine, but so far, this produces a memory leak (as far as I can tell) when I put more than 2**24 points into the algorithm. I have emailed Brad Barber, who indicated that this changing the ridgeid field was possible. I am also working on a divide and conquer approach, but am having trouble deleting the old Delaunay object before instantiating the next one. I can send that code along if there is no good way to increase the number of points in the algorithm.
I think it would be useful to have an algorithm that can accept more points, given the increased speed and processing power of modern computers. Thanks for your help.
I can't immediately help you with the Qhull internals, but other comment:
Your best bet with that would probably be to try to get your patch included in Qhull itself. The Qhull code is shipped with Scipy only for convenience, and is unmodified from the 2012.1 version --- shipping a patched version might not be a very good thing if the Qhull upstream is still responsive.
Of course, Qhull is not the only game in town, at least for low-dimensional triangulations, and it could be interesting to have more featureful alternatives for those cases.
_______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
22.06.2015, 22:07, Cantwell Carson kirjoitti:
Are there other ideas for what could ship with scipy? I have tried CGAL and mayavi and they function, but no better than Qhull (IMHO). Matplotlib also has a routine, but I doubt the redundancy in algorithm would be beneficial.
CGAL might be more robust as it AFAIK uses exact predicates. It is however GPL-licensed, so it cannot be used in Scipy without making the whole bundle GPL. I know there are a number of other libraries doing Delaunay in 2d/3d, but can't name any from top of my head. Last time I looked, I don't remember seeing much license-compatible code however, and performance/scaling issues such as here probably aren't immediately apparent on first sight.
I recommend that we insert a warning or error when the number of points in Qhull algorithm exceeds 2**24 so that the method can exit without crashing, instead of killing the whole process as it does now.
Seems like a good idea, patches accepted. Ticket: https://github.com/scipy/scipy/issues/4985
participants (2)
-
Cantwell Carson -
Pauli Virtanen