Triangular grid interpolation and refinement #1521

ianthomas23 opened this Issue Nov 20, 2012 · 3 comments

3 participants

Matplotlib Developers member

As discussed on the matplotlib-devel mailing list, @GBillotey has recently written and donated code to perform interpolation and grid refinement for triangular grids. It is very useful functionality. As it needs some rewriting and I have an alternative algorithm to add to it, Geoffroy and I have decided to share the work. We intend to submit 2 or 3 pull requests to cover all the functionality.

In a spirit of collaboration, I am documenting more of the high-level design than usual so that we both know what we are doing and to allow other developers to comment - there are a couple of specific questions I have below.

Firstly I should point out that the existing delaunay module has a linear interpolator, but it is only useful for those using that module to perform the triangulation for them. Those specifying their own triangulation do not currently have such functionality. The additions here will provide linear and cubic interpolators for those specifying their own triangulations or using the delaunay module, plus a grid refiner to allow higher-resolution grids to be generated and plotted. Also, a class used to determine in which triangle a particular (x,y) point lies is required for the interpolation routines, and this may be of additional use in its own right.

For the interpolators, there will be an abstract base class TriInterpolator, and concrete derived classes LinearTriInterpolator and CubicTriInterpolator. They will define __call__(x,y) where x and y are arrays with the same shape and any dimension. 1D arrays would be used to calculate a transect through a triangulation, 2D arrays to obtain a quad grid for e.g. a contourf call, etc. Usage example:

triang = Triangulation(x, y, triangles)  # Existing code.
interp = TriLinearInterpolator(triang, z)  # z array of same shape as x & y.
xnew, ynew = np.meshgrid(np.linspace(0, 1, 10), np.linspace(0, 1, 10))
znew = interp(xnew, ynew)

QUESTION 1: what to do with interpolation points outside the triangulation? We can either return a default floating-point value such as np.nan, or used a fully-fledged masked array. My initial preference was the masked array, but now I think that is overkill and the np.nan option is better.

For the grid refiner there will be an abstract base class TriRefiner, and a derived concrete class TriRefiner (to be decided). This will allow other refiners to be added later.

To find in which triangles particular (x,y) points are located, there will be a separate TriFinder class. My preference, as above, is to have TriFinder as an abstract base class and the derived concrete class is TrapezoidMapTriFinder, named after the algorithm used. Again, this approach allows other TriFinder-derived classes to be added later if required.

Now TrapezoidMapTriFinder does some once-only initialisation of a search structure to make lookups fast. Once one has been created and initialised, it makes sense to reuse it if more than one interpolator is used. Hence my preference here is to allow a particular TriFinder to be specified in an interpolator constructor if desired, but by default each interpolator will call triangulation.get_tri_finder() to retrieve and reuse the triangulation's default TriFinder. Initially a Triangulation will not have a TriFinder defined but it will be created when it is first requested. This approach means that most users can completely ignore TriFinders and just deal with the required interpolators - a single TriFinder will be created behind the scenes and reused for speed. QUESTION 2: Is this Triangulation-TriFinder arrangement OK? They will have circular references to each other, but I think this is acceptable given that it makes the standard usage so easy.



I maybe wong but I make the following suggestion:

While it is clear that a TriFinder needs a reference to a Triangulation instance, why would Triangulation need to hold a reference to the full TriFinder ?
In a weaker way it could store the once-only computed data structure (I dont't know the algo in details, but it could store for instance a tree of triangle indexes). When computed once, this data structure would be reused by any Trifinder pointing to this specific Triangulation (but without explicit cross-references).

So the data structure would be more or less like this:
TriInterpolator -> TriFinder -> Triangulation -> once-only computed data

Also, not for the standard user but a close() or finalise() method of TriInterpolator could do the cleanup of Triangulation object or explicitely remove the cross-reference.
That way an 'advanced user' could decide that he will not need the Triangulation additionnal data structure any more.

Matplotlib Developers member

@GBillotey: sure, we could do it your suggested way, but I think the simple implementation with the circular references is OK. Python's garbage collector can deal with circular references in simple cases like this, and we have the option of using the weakref module if we want to be explicit.

I suppose that rather than asking if the implementation was acceptable I was really asking if the use of circular references was considered by others to be bad design, in case we had an unwritten policy on such things.

Matplotlib Developers member

I suppose that rather than asking if the implementation was acceptable I was really asking if the use of circular references was considered by others to be bad design, in case we had an unwritten policy on such things.

Personally: yes. But matplotlib uses circular references all over the place (did you know an Artist holds the axes it may [or may not] belong too...), so if it makes implementation easier, go for it!

@ianthomas23 ianthomas23 closed this Mar 7, 2013
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment