I think the code used to determine which triangle contains a certain point
should be factored out into its own TriFinder class,
+1 – this is a generally useful feature. In fact, it would be nice if
a lot of this were in a pacakge that deals with triangular meshes,
apart from MPL altogether (a scikit maybe?)
I hadn’t considered any wider interest in it. From a selfish matplotlib point of view, even if it was in say a scikit, we wouldn’t want to add an external dependency to the scikit so we would still need a copy of the source code within matplotlib anyway, just to do our interpolation for plotting purposes. The situation would be just like the delaunay code which exists in a scikit but we include it in the mpl source code to avoid the external dependency.
Once it is in mpl it is available for other projects to use if they wish, subject to the BSD-style license. There would be a question of who is responsible for maintaining it as I don’t particularly want to spread myself thinly on other open source projects when there are a lot of additions to mpl that I would like to do. There would be the danger for the other project of the reverse situation of our delaunay example, where we have code that needs improvement but it is essentially unmaintained, causing problems for users and embarrassment for developers.
I have a C++ TriFinder class
that I could modify to work within matplotlib, and it is O(log N) so should
be faster than your version for typical use cases.
What algorithm does this use? Is the code open source and/or availabel
for other projects?
I’m working on a package for working with unstructured grids in
general, and also have a use for “what triangle is this point in” code
for other purposes – and I havne’t found a fast, robust code for this
yet.
It is code I have written recently based on the trapezoidal map algorithm from the book ‘Computational Geometry: Algorithms and Applications’ by M. de Berg et al. It currently exists only on my main linux box, but it will be open source within mpl for others to use as mentioned above (subject to acceptance by the mpl developers of course). It is an excellent book that I wholeheartedly recommend to anyone with a passing interest in computational geometry. The algorithm itself looks painful initially (as do many of the algorithms in the book) but it reduces to a surprisingly small amount of code.
For well-formed triangulations (i.e. no duplicate points, no zero area triangles and no overlapping triangles) the algorithm is ‘sufficiently’ robust. It is not completely robust in the sense of Shewchuk’s robust predicates, but is designed quite cleverly (or luckily) to avoid many of the pitfalls that occur in naive point-in-triangle tests. For example, a naive implementation of a point-in-triangle test for a point on or very near the common edge of two adjacent triangles might return true for both triangles (which is usually fine), or false for both (which is catastrophic). The trapezoidal map avoids these problems by reducing the test to a single point-line test which is therefore guaranteed to return just one of the two triangles. It is possible to think of a situation that causes the algorithm to fail at a triangle which has such a small area that the slopes of two of the sides are within the machine precision of each other, but it is also possible to use the triangle connectivity to check for and resolve this problem. I haven’t done so yet and need to consider the tradeoff of effort required vs potential gain.
Someone who required guaranteed robustness could use Shewchuk’s point-line test with the algorithm. I would not do this with mpl however, as the correctness of the point-line test depends a lot on computer architecture and compiler flags. This would get out of date quickly and would need keen monitoring by someone with knowledge of and interest in the area, plus access to a lot of different computer architectures and OSes. I fail on all of those counts!
particularly as only a few days ago I committed to writing a triangular grid
interpolator for quad grids
what is a triangular interpolator for quad grids? sounds useful, too.
That was poor English from me. I meant interpolating from a triangular grid to a quad grid, typically to make use of the wide range of quad grid plotting functions like contour, pcolor, etc.
Ian
···
On 15 November 2012 21:25, Chris Barker <chris.barker@…236…> wrote:
On Wed, Nov 14, 2012 at 1:50 AM, Ian Thomas <ianthomas23@…149…> wrote: