Interpolation in a triangular mesh (tri.Triangulation)

Hi!

I had recently to develop interpolators for a function defined at the nodes
of a user-specified triangular mesh.
(Beside interpolation, it can help producing higher-quality tricontour
plots, using interpolation on a refined mesh and matplotlib tricontour
function.)

Being a regular user of matplotlib, I would be happy if it can be useful to
others...
The code is hosted here:
https://github.com/GBillotey/trimesh-interpolator.git

Please let me know if it this dev. can be useful and if the code needs some
cleaning (I am not a computer scientist, only a mechanical engineer)

Cheers,
Geoffroy.

···

--
View this message in context: http://matplotlib.1069221.n5.nabble.com/Interpolation-in-a-triangular-mesh-tri-Triangulation-tp39586.html
Sent from the matplotlib - devel mailing list archive at Nabble.com.

Hi Geoffroy

This will certainly be very useful. I need to spend some time looking at it and seeing how it would best fit within the matplotlib framework, particularly as only a few days ago I committed to writing a triangular grid interpolator for quad grids and it would be sensible to group these interpolators together in some way.

I’ll get back to you when I’ve had time to look at it.

Thanks for your efforts!
Ian

···

On 28 October 2012 20:17, GBillotey <geoffroy.billotey@…149…> wrote:

Hi!

I had recently to develop interpolators for a function defined at the nodes

of a user-specified triangular mesh.

(Beside interpolation, it can help producing higher-quality tricontour

plots, using interpolation on a refined mesh and matplotlib tricontour

function.)

Being a regular user of matplotlib, I would be happy if it can be useful to

others…

The code is hosted here:

https://github.com/GBillotey/trimesh-interpolator.git

Please let me know if it this dev. can be useful and if the code needs some

cleaning (I am not a computer scientist, only a mechanical engineer)

Cheers,

Geoffroy.

Hi Geoffroy,

I have had some time to look at your TriLinearInterpolator in some detail (the other two files only briefly). I would indeed like to add something like this to matplotlib - the mesh refinement looks very nice and the interpolators would be useful to many people.

As you suspected, the code does need significant changes before we can include it. Some are merely cosmetic, as all code must adhere to PEP8 and the matplotlib coding guidelines, but there are also some functional and performance improvements. For example, your wavefront method for finding the triangle containing a certain point must be able to deal with masked triangulations and indeed triangulations that are discontinuous, for example two islands in a masked-out ocean, which is unusual but must be supported. In terms of performance, there is much explicit looping within numpy arrays that could be improved using other numpy array commands, and would also reduce the length of the source code. There is an argument for some of the performance-critical code to be in C/C++.

I think the code used to determine which triangle contains a certain point should be factored out into its own TriFinder class, so that (1) it does not need to be replicated in the two interpolator classes, and (2) different algorithms can be easily swapped if necessary. 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.

I expect that this is probably more work than you anticipated when you asked if the code needed any improvement! I propose the following: if you are happy to give matplotlib your source code as it stands and for us to include it under our BSD-style license, then I will take on the responsibility of getting it into a form that will be accepted by the other developers. I will acknowledge your contribution in both the source code and on the web site, something like “based on code contributed by Geoffroy Billotey”.

Alternatively, if you would like to use this as an excuse to learn how to contribute to matplotlib more actively but don’t want to take on everything, then we could divide up the work so that first I write my C++ log(N) TriFinder class and the linear interpolator that uses it, and then you could modify the cubic interpolator following the format of the linear interpolator and using my guidance as and when you need it.

Let me know your preference,
Ian

P.S. Never apologise for not being a computer scientist! Many of our developers, myself included, are proper scientists or engineers!!!

···

On 29 October 2012 09:37, Ian Thomas <ianthomas23@…149…> wrote:

Hi Geoffroy

This will certainly be very useful. I need to spend some time looking at it and seeing how it would best fit within the matplotlib framework, particularly as only a few days ago I committed to writing a triangular grid interpolator for quad grids and it would be sensible to group these interpolators together in some way.

I’ll get back to you when I’ve had time to look at it.

Thanks for your efforts!
Ian

On 28 October 2012 20:17, GBillotey <geoffroy.billotey@…149…> wrote:

Hi!

I had recently to develop interpolators for a function defined at the nodes

of a user-specified triangular mesh.

(Beside interpolation, it can help producing higher-quality tricontour

plots, using interpolation on a refined mesh and matplotlib tricontour

function.)

Being a regular user of matplotlib, I would be happy if it can be useful to

others…

The code is hosted here:

https://github.com/GBillotey/trimesh-interpolator.git

Please let me know if it this dev. can be useful and if the code needs some

cleaning (I am not a computer scientist, only a mechanical engineer)

Cheers,

Geoffroy.

Hello Ian,

Thank you for your second proposition ; I find it very interesting in fact.

(But beware that, as I do not feel able to do this on my own, your code exemple / guidance would be needed for sure… )

I will have a look at http://matplotlib.org/devel/index.html

···

I think your idea of having a separate TriFinder class is quite good. My search algo. is not optimised, only avoiding a O(N) for interpolations along a path or line.

Some other parts of the code (especially the loop over the (x,y) points in call() ) may also be performance-critical, at least for the ‘cubic’ interpolator. But I dont have a clear picture on how to embed C++ code in python so I would need your example to figure out what is possible.

Regards
Geoffroy.

2012/11/14 Ian Thomas <ianthomas23@…149…>

Hi Geoffroy,

I have had some time to look at your TriLinearInterpolator in some detail (the other two files only briefly). I would indeed like to add something like this to matplotlib - the mesh refinement looks very nice and the interpolators would be useful to many people.

As you suspected, the code does need significant changes before we can include it. Some are merely cosmetic, as all code must adhere to PEP8 and the matplotlib coding guidelines, but there are also some functional and performance improvements. For example, your wavefront method for finding the triangle containing a certain point must be able to deal with masked triangulations and indeed triangulations that are discontinuous, for example two islands in a masked-out ocean, which is unusual but must be supported. In terms of performance, there is much explicit looping within numpy arrays that could be improved using other numpy array commands, and would also reduce the length of the source code. There is an argument for some of the performance-critical code to be in C/C++.

I think the code used to determine which triangle contains a certain point should be factored out into its own TriFinder class, so that (1) it does not need to be replicated in the two interpolator classes, and (2) different algorithms can be easily swapped if necessary. 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.

I expect that this is probably more work than you anticipated when you asked if the code needed any improvement! I propose the following: if you are happy to give matplotlib your source code as it stands and for us to include it under our BSD-style license, then I will take on the responsibility of getting it into a form that will be accepted by the other developers. I will acknowledge your contribution in both the source code and on the web site, something like “based on code contributed by Geoffroy Billotey”.

Alternatively, if you would like to use this as an excuse to learn how to contribute to matplotlib more actively but don’t want to take on everything, then we could divide up the work so that first I write my C++ log(N) TriFinder class and the linear interpolator that uses it, and then you could modify the cubic interpolator following the format of the linear interpolator and using my guidance as and when you need it.

Let me know your preference,
Ian

P.S. Never apologise for not being a computer scientist! Many of our developers, myself included, are proper scientists or engineers!!!

On 29 October 2012 09:37, Ian Thomas <ianthomas23@…149…> wrote:

Hi Geoffroy

This will certainly be very useful. I need to spend some time looking at it and seeing how it would best fit within the matplotlib framework, particularly as only a few days ago I committed to writing a triangular grid interpolator for quad grids and it would be sensible to group these interpolators together in some way.

I’ll get back to you when I’ve had time to look at it.

Thanks for your efforts!
Ian

On 28 October 2012 20:17, GBillotey <geoffroy.billotey@…149…> wrote:

Hi!

I had recently to develop interpolators for a function defined at the nodes

of a user-specified triangular mesh.

(Beside interpolation, it can help producing higher-quality tricontour

plots, using interpolation on a refined mesh and matplotlib tricontour

function.)

Being a regular user of matplotlib, I would be happy if it can be useful to

others…

The code is hosted here:

https://github.com/GBillotey/trimesh-interpolator.git

Please let me know if it this dev. can be useful and if the code needs some

cleaning (I am not a computer scientist, only a mechanical engineer)

Cheers,

Geoffroy.

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

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.

-Chris

···

On Wed, Nov 14, 2012 at 1:50 AM, Ian Thomas <ianthomas23@...149...> wrote:

--

Christopher Barker, Ph.D.
Oceanographer

Emergency Response Division
NOAA/NOS/OR&R (206) 526-6959 voice
7600 Sand Point Way NE (206) 526-6329 fax
Seattle, WA 98115 (206) 526-6317 main reception

Chris.Barker@...236...

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 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 pretty sure there is an O(log n) algorithm in the Numerical
Recipes book. It requires you to construct the triangulation in a
specific way (this allows one to set up a tree data structure of
triangles nicely). There may be others that I am not aware of though.

···

On Thu, Nov 15, 2012 at 3:25 PM, Chris Barker <chris.barker@...236...> wrote:

On Wed, Nov 14, 2012 at 1:50 AM, Ian Thomas <ianthomas23@...149...> wrote:

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.

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.

-Chris

--

Christopher Barker, Ph.D.
Oceanographer

Emergency Response Division
NOAA/NOS/OR&R (206) 526-6959 voice
7600 Sand Point Way NE (206) 526-6329 fax
Seattle, WA 98115 (206) 526-6317 main reception

Chris.Barker@...236...

------------------------------------------------------------------------------
Monitor your physical, virtual and cloud infrastructure from a single
web console. Get in-depth insight into apps, servers, databases, vmware,
SAP, cloud infrastructure, etc. Download 30-day Free Trial.
Pricing starts from $795 for 25 servers or applications!
http://p.sf.net/sfu/zoho_dev2dev_nov
_______________________________________________
Matplotlib-devel mailing list
Matplotlib-devel@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/matplotlib-devel

--
Damon McDougall
http://www.damon-is-a-geek.com
Institute for Computational Engineering Sciences
201 E. 24th St.
Stop C0200
The University of Texas at Austin
Austin, TX 78712-1229

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:

I think this is the standard method used for point-in-triangulation tests for delaunay triangulations, as the search structure is created as the triangulation is built. We need to support non-delaunay triangulations that are specified by the user, which requires a different approach.

Ian

···

On 16 November 2012 05:14, Damon McDougall <damon.mcdougall@…149…> wrote:

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 pretty sure there is an O(log n) algorithm in the Numerical

Recipes book. It requires you to construct the triangulation in a

specific way (this allows one to set up a tree data structure of

triangles nicely). There may be others that I am not aware of though.