Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Using qhull for Delaunay triangulation #2504

Merged
merged 5 commits into from Dec 20, 2013

Conversation

ianthomas23
Copy link
Member

This PR adds (some of) the Qhull library to matplotlib to perform Delaunay triangulation. Currently Delaunay triangulations are performed by the matplotlib.delaunay module, but this is not particularly robust. Qhull is much more robust, is open-source, and is used in many other projects including scipy.spatial, octave and matlab.

The qhull source code and license file are in extern/qhull/, and the Python/C wrapper in src/qhull_wrap.c. The setup scripts attempt to build the qhull extension in three ways (as other extensions): (1) using system-installed version via pkg-config, (2) using system-installed version without pkg-config, (3) using the local version in extern/qhull. I've tested all three methods using python 2 and 3 on linux, but not macos or windows as I don't use them. Any help here would be appreciated - you can confirm that everything is OK by running test_triangulation and test_mlab.

Qhull is used as a single atomic operation, i.e. when a triangulation is needed qhull is started up, the triangulation is calculated and returned, then qhull is closed down. There is no storage of qhull objects in memory in case they are needed later.

I've deprecated the matplotlib.delaunay module for future removal. Note to self: when it is time to remove it, the following need to be deleted: lib/matplotlib/delaunay, lib/matplotlib/tests/test_delaunay and lib/matplotlib/tests/baseline_images/test_delaunay, and the following modified: setup.py and setupext.py.

The delaunay code is used by matplotlib.tri.Triangulation, and hence by all of the pyplot.tri* functions. It is also used by mlab.griddata. I don't think that matplotlib.delaunay was specifically intended to be used directly by users, but it has been. This will now be harder to do which is by intent. I will be encouraging everyone to access the delaunay functionality via matplotlib.tri.Triangulation as this provides a standard interface, does consistency and error checking, and can be passed to any of the pyplot.tri* functions for plotting.

The PR is in two commits for easier review; the first is the addition of the qhull files, and the second contains my changes. I have added a number of new tests, in particular cases that used to fail using matplotlib.delaunay but now pass OK.

There are two decisions I have made that others may wish to question/alter:

(1) qhull expects a valid FILE* for writing errors to. scipy and octave use stderr. I haven't done this as I think writing to stderr is bad practice for a C extension in python, plus the error messages are too technical for our audience. So instead I create a temporary file using the ANSI function tmpfile() and write errors to it that are ignored. Instead I raise a python RuntimeError with a simpler explanation of the error. I have used tmpfile() rather than say writing to /dev/null as it is cross-platform. However, in case a user really does want the full details of the qhull error, e.g. to help with a bug report, if verbose mode is enabled then I do use stderr instead of tmpfile. This isn't elegant but I can't think of a better solution.

(2) mlab.griddata is a matlab-similar utility function to interpolate from an unstructured (i.e. triangular) grid to some other grid. It provides the option for linear or natural neighbor ('nn') interpolation, the latter being specific to delaunay triangulations. You could use either linear or nn interpolation with matplotlib.delaunay, but because this wasn't very robust if you had installed the natgrid mpl_toolkit you could use nn interpolation with that. Now the combination of qhull and matplotlib.tri.Triangulation does not support nn interpolation and I don't wish to do so as that would involve keeping the qhull data structures in memory in case they were subsequently needed for an nn interpolation, and that is a significant change that I can't justify for what is a matlab utility function. So now griddata has two options: for linear interpolation qhull is used, which is always available; for nn interpolation natgrid must be installed and if it isn't an error is raised. Unfortunately the historic default is nn interpolation so that a standard installation without natgrid will by default raise an error. Again this is not elegant but as I don't know why the default is nn I don't feel I can change it.

There is an argument that now the delaunay triangulations are robust, there is no need at all for natgrid and it could just be removed and the only interpolation option being linear, but I think this would be too big a change in one go. Finally I should add that I will be discouraging use of griddata and instead encouraging use of the matplotlib.tri.TriInterpolator classes which are more flexible and more powerful, so griddata will become less important than it is now.

@pelson
Copy link
Member

pelson commented Oct 9, 2013

This looks great. I'm not the biggest fan of shipping the dependency inside our own repo, but I guess it is consistent with our use of Agg and CXX, so I'm not going to make a song and dance about that.

Which platforms have you tested this on? I'd like to get some verification that this is working on windows and OSX at the least.

I've never got around to trying out the NADGRIDS extension myself, I think we should at least port this over to github and get some travis testing of it either here (it isn't being installed atm.) or in its own repo.

This would be a good subject for a discussion on our next developer hangout (2 weeks tomorrow).

@ianthomas23
Copy link
Member Author

@pelson: I've tested it comprehensively on linux (3 different distributions, 32 and 64-bit, with and without natgrid), but not at all on macos or windows as I don't use them.

I haven't looked at natgrid much except to confirm that it was generated using Cython which I have no experience of. For this PR I've opted for minimal change with respect to how it used to be, but I suspect that the more I look into it the more I would like to remove natgrid completely.

It may be a good subject for the next developer hangout but I can't make it (Thursdays, UK office hours are always bad for me). I can answer questions before and after though.

@mdboom
Copy link
Member

mdboom commented Nov 14, 2013

It looks like we have some PEP8 failures here. Once those are fixed, I think we'll just look over this for docstring and style issues, but I suspect the actual work here is probably fine.

@ianthomas23
Copy link
Member Author

The pep8 failures are due to the bootstrap code in _qhull.py that is autogenerated by distutils. It looks like this is a known problem as all the other bootstrap files for C/C++ extensions are listed as EXCLUDE_FILES in test_coding_standards.py. I will therefore do the same with the _qhull.py.

However, there are problems with the C code. @GBillotey has kindly tested the code on windows and identified a problem that I need to look at before this PR can proceed.

@ianthomas23
Copy link
Member Author

I have rebased this PR, dealt with the pep8 failure and fixed the problem on windows. For the latter, rather than creating a temporary file to write to which causes ownership problems on windows, I am instead telling qhull to write to /dev/null on POSIX-compliant OSes and nul on windows. And rather than doing some nasty #ifdef...#else...#endif in the C code to deal with the different OSes, I am instead using the handy python os.devnull as I think this is more future-proof.

I have fully tested the PR on linux and windows. There should be no problems with it working on macos.

@@ -102,6 +90,13 @@ def calculate_plane_coefficients(self, z):

@property
def edges(self):
"""
Return integer array of shape (?,2) containing all edges of
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should this be (N, 2)?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

No, this is intended. In the class documentation there are other arrays with lengths of npoints and ntri, which are fundamental quantities for a triangulation. You could argue that it should be (nedges, 2), but then you would need to explain that nedges is the number of edges. I prefer the more succinct use of the question mark.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Then say n_edges. The question mark makes it look like it was a mistake or
an oversight.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

OK, have changed it to nedges.

@pelson
Copy link
Member

pelson commented Dec 12, 2013

I think this is a most excellent (if admittedly huge) PR. I've gone over the code for smells and areas of obscurity, but I'm satisfied that what has been implemented has been implemented well and thoroughly.

I'm going to leave this up for 4 days (until Monday the 16th of December), but if there are no further objections, I will happily merge this.

It is almost inevitable that this will have some impact on the build for some users, but it is unrealistic IMHO to test every variant of OS/numpy etc. and it is for that reason we have release candidates 😉. Thank you @ianthomas23 for testing this in Windows and Linux.

👍 - great job!

pelson added a commit that referenced this pull request Dec 20, 2013
Using qhull for Delaunay triangulation
@pelson pelson merged commit 6afba8d into matplotlib:master Dec 20, 2013
@pelson
Copy link
Member

pelson commented Dec 20, 2013

Good job @ianthomas23! I've posted on the mpl-devel mailing list to try to get some real life testers before the release candidate phase.

@iNand
Copy link

iNand commented Feb 20, 2014

Thank you very much for this Patch! On my current (irregular) datasets, I encountered at least 3 different errors, using matplotlib 1.3.x calling matplotlib.mlab.griddata, while the developer branch and Qhull, solves (almost) all problems! I tried different datasets, and I found it really stable - I hope you can release matplotlib 1.4.x soon, so I could switch back to an official version.

@tacaswell
Copy link
Member

@ngoldbaum Reading this, it looks like we still have the old triangulation code, it just got moved.

@ngoldbaum
Copy link
Contributor

tacaswell added a commit to tacaswell/matplotlib that referenced this pull request Jan 10, 2015
tacaswell added a commit to tacaswell/matplotlib that referenced this pull request Feb 19, 2015
Removed in PR matplotlib#2504

merged to master as 6afba8d

Conflicts:
	doc/api/api_changes/code_removal.rst
tacaswell added a commit to tacaswell/matplotlib that referenced this pull request Apr 19, 2015
Removed in PR matplotlib#2504

merged to master as 6afba8d

Conflicts:
	doc/api/api_changes/code_removal.rst
          - documentation conflicts
tacaswell added a commit to tacaswell/matplotlib that referenced this pull request Apr 19, 2015
Removed in PR matplotlib#2504

merged to master as 6afba8d

Conflicts:
	doc/api/api_changes/code_removal.rst
          - documentation conflicts
tacaswell added a commit to tacaswell/matplotlib that referenced this pull request Sep 16, 2016
Removed in PR matplotlib#2504

merged to master as 6afba8d

Conflicts:
	doc/api/api_changes/code_removal.rst
          - documentation conflicts
tacaswell added a commit to tacaswell/matplotlib that referenced this pull request Sep 17, 2016
tacaswell added a commit to tacaswell/matplotlib that referenced this pull request Oct 19, 2016
@WombatBill
Copy link

I've been having trouble exporting the results of Delaunay triangulation in either Matplotlib or now through QHull in a form that I can use.
I eventually want to do some statistical analysis on the shape of the triangles (something done with an old archaeology program, BASP) to see whether the distribution of a whole lot of points are clustered, randomly distributed or patterned.

Matplotlib seems to give me sets of 4 points (the "simplices") that I've had trouble saving in anything like a useful format, but which I first assumed but have since found don't correspond to the old STL (Standardized Tessalation Language) structure: all 3 points in a triangular facet, plus the first point repeated, but rather to some two-facet polygon.
Just what QHull is giving me... I don't really know...

Anyway: could you give me some assistance with either exporting the simplices into a numpy array or a csv file or something, or what I can or cannot do with what comes out of QHull???

@dopplershift
Copy link
Contributor

dopplershift commented Mar 12, 2019

I think you're better off looking at SciPy's Delaunay code which also uses QHull, but is designed to give you numpy arrays.

@ianthomas23
Copy link
Member Author

@WombatBill A merged PR is not the place to put a usage question. It should be posted to the matplotlib users mailing list.

@dopplershift Please don't advise users to use scipy when it is not required.

If you google "matplotlib triangulation" the first result is https://matplotlib.org/api/tri_api.html. A cursory glance at the Triangulation class will show you what to do. Here is a minimal example:

>>> import matplotlib.tri as mtri
>>> x = [0, 0, 1, 1]
>>> y = [0, 1, 0, 1]
>>> triang = mtri.Triangulation(x, y)
>>> print(triang.triangles)
>>> print(type(triang.triangles))
[[3 1 0]
 [0 2 3]]
<class 'numpy.ndarray'> (2, 3)

which shows that 'triang.triangles' is a numpy array of shape (number of triangles, 3).

@WombatBill
Copy link

WombatBill commented Mar 12, 2019 via email

@WombatBill
Copy link

I've been using Scipy, but couldn't export it into any useful format.
This is why I went to QHull, trying to bypass the export issues, and figure out what the two programs are actually calculating.
And that is not at all clear.

Ideally, somehow, from a list of coordinates (x and y, or x, y and z) I should be able to get a list of the coordinates for individual triangles (or facets), and triangles have 3 coordinates, not the 4 I get from Scipy.

@ianthomas23
Copy link
Member Author

@WombatBill Read the link I gave you...

@WombatBill
Copy link

I am reading it; again...
I have read it before...

@ianthomas23
Copy link
Member Author

triang.triangles is an array of the indices of the points that comprise the triangles. If you want the x and y coordinates of those triangles, here is a simple example:

import matplotlib.tri as mtri
import numpy as np
x = np.asarray([0, 0, 1, 1])
y = np.asarray([0, 1, 0, 1])
triang = mtri.Triangulation(x, y)
triangles = triang.triangles
print('triangle indices', triangles)
print('triangle x-values', x[triangles])
print('triangle y-values', y[triangles])

The output is

triangle indices [[3 1 0]
 [0 2 3]]
triangle x-values [[1 0 0]
 [0 1 1]]
triangle y-values [[1 1 0]
 [0 0 1]]

@QuLogic QuLogic added this to the v1.4.0 milestone Mar 24, 2020
@ianthomas23 ianthomas23 deleted the qhull_delaunay branch July 8, 2021 18:21
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

10 participants