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 <geo...@gm...> 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 <ian...@gm...> 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 <geo...@gm...> 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 <ian...@gm...> > 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 <ian...@gm...> 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 <geo...@gm...> 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. >>> >> >
On Wed, Nov 14, 2012 at 1:50 AM, Ian Thomas <ian...@gm...> wrote: > 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 -- 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 Chr...@no...
On 15 November 2012 21:25, Chris Barker <chr...@no...> wrote: > On Wed, Nov 14, 2012 at 1:50 AM, Ian Thomas <ian...@gm...> wrote: > > > 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 Thu, Nov 15, 2012 at 3:25 PM, Chris Barker <chr...@no...> wrote: > On Wed, Nov 14, 2012 at 1:50 AM, Ian Thomas <ian...@gm...> wrote: > >> 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. > > 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 > > Chr...@no... > > ------------------------------------------------------------------------------ > 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 > Mat...@li... > 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
On 16 November 2012 05:14, Damon McDougall <dam...@gm...>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. > 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