Linear tri interpolator #1582

Merged
merged 5 commits into from Jan 10, 2013

Projects

None yet

5 participants

@ianthomas23
Member

This is the first PR for issue #1521. It adds a TrapezoidMapTriFinder that is used to determine the triangles in which (x,y) points lie and a LinearTriInterpolator to interpolate scalar fields defined on triangular grids. Have included tests for both classes, a pylab_example for the interpolator and an whizzy interactive event_handling example for the trifinder, and have updated sphinx docs, whats_new and CHANGELOG.

I'll repeat the observation from #1521 that the delaunay module already has a linear interpolator, but it cannot be used with user-specified triangulations unlike the one in this PR which can be used with user-specified and delaunay-calculated triangulations. The new one is much more flexible in accepting arrays of points to iterate over as they can be any shape, e.g. 1D for a transect through a triangulation, 2D to create a quad mesh to pass to contourf, etc.

It is possible that the TrapezoidMapTriFinder is overkill for what we need. It is O(log N) but may be outperformed by a simpler 'walking across the triangulation' algorithm if most users are passing in hi-res regular 2D grids to interpolate to. My preference is to release this code out into the wild as it is and get feedback from users before adding/modifying the trifinder.

Member

In the absence of comments from anyone else, I will check this in myself. I'll leave it a few days more in case there are any objections.

Member

Just gave this a visual inspection, and the code looks very well commented and fairly easy to follow. Good job, @ianthomas23

@mdboom mdboom commented on an outdated diff Jan 7, 2013
lib/matplotlib/tri/_tri.cpp
+ {
+ // The 3 triangle points have identical x and y! The best
+ // we can do here is take normal = (0,0,1) and for the
+ // constant p take the mean of the 3 points' z-values.
+ normal = XYZ(0.0, 0.0, 1.0);
+ point0.z = (point0.z + point1.z + point2.z) / 3.0;
+ }
+ }
+
+ *planes++ = -normal.x / normal.z; // x
+ *planes++ = -normal.y / normal.z; // y
+ *planes++ = normal.dot(point0) / normal.z; // constant
+ }
+ }
+
+ return Py::asObject((PyObject*)planes_array);
mdboom
mdboom Jan 7, 2013 Owner

This function should probably be wrapped in a try ... catch to ensure that Py_DECREF(z) is called. It seems unlikely an exception would get thrown here, but a future change may introduce something that could.

Owner
mdboom commented Jan 7, 2013

Looks good on quick experimentation. 👍 from me.

@NelleV NelleV and 2 others commented on an outdated diff Jan 8, 2013
examples/pylab_examples/triinterp_demo.py
@@ -0,0 +1,36 @@
+"""
+Interpolation from triangular grid to quad grid.
+"""
+import matplotlib.pyplot as plt
+import matplotlib.tri as mtri
+import numpy as np
+import math
+
+
+# Create triangulation.
+x = np.asarray([0, 1, 2, 3, 0.5, 1.5, 2.5, 1, 2, 1.5])
+y = np.asarray([0, 0, 0, 0, 1, 1, 1, 2, 2, 3])
+triangles = [[0,1,4], [1,2,5], [2,3,6], [1,5,4], [2,6,5], [4,5,7], [5,6,8],
NelleV
NelleV Jan 8, 2013 Contributor

Tiny nitpicks: this file is not completely pep8 compliant: there should be spaces after ,, and spaces around * (l. 18)

ianthomas23
ianthomas23 Jan 8, 2013 Member

I think this line is more readable without the spaces after commas, but as I am not particularly bothered either way I am happy to change it.

No extras spaces are needed around any of the * in the line
z = np.cos(1.5*x)*np.cos(1.5*y)
From PEP8: "If operators with different priorities are used, consider adding whitespace around the operators with the lowest priority(ies). Use your own judgement". There are no operator priorities causing confusion here; the line is fine as it is. I assume you are using the pep8 tool which is sometimes a little too enthusiastic - it's recommendation of
z = np.cos(1.5 * x) * np.cos(1.5 * y)
is very poor.

NelleV
NelleV Jan 8, 2013 Contributor

I don't have to use the pep8 tool to see this isn't pep8. I prefer z = np.cos(1.5 * x) * np.cos(1.5 * y) above z = np.cos(1.5*x)*np.cos(1.5*y)

This is a matter of choosing whether matplotlib follows strict pep8 convention or not. I always follow the strict pep8 conventation, and many scientific python projects do (with the pep8 tool or flake8 to check the compliance). I believe matplotlib should do the same (it avoids the long discussions on whether X is better than Y and the numerous commits to reformat the code such as the ones I've been doing recently on matplotlib's codebase), but I don't think there has been a clear decision on that.

ianthomas23
ianthomas23 Jan 8, 2013 Member

I think you are missing the point. I quoted the line from PEP8 that specifically states "use your own judgement". It is impossible to strictly follow the PEP8 document without exercising judgement, i.e. there is no definite right or wrong answer here.

The coding guidelines state that I should follow PEP8. I have. If you think the coding guidelines should be changed to follow pep8 (the tool, not the document) then you should start such a discussion on the matplotlib-devel mailing list.

efiring
efiring Jan 9, 2013 Owner

On 2013/01/08 6:47 AM, Varoquaux wrote:

In examples/pylab_examples/triinterp_demo.py:

@@ -0,0 +1,36 @@
+"""
+Interpolation from triangular grid to quad grid.
+"""
+import matplotlib.pyplot as plt
+import matplotlib.tri as mtri
+import numpy as np
+import math
+
+
+# Create triangulation.
+x = np.asarray([0, 1, 2, 3, 0.5, 1.5, 2.5, 1, 2, 1.5])
+y = np.asarray([0, 0, 0, 0, 1, 1, 1, 2, 2, 3])
+triangles = [[0,1,4], [1,2,5], [2,3,6], [1,5,4], [2,6,5], [4,5,7], [5,6,8],

I don't have to use the pep8 tool to see this isn't pep8. I prefer |z =
np.cos(1.5 * x) * np.cos(1.5 * y)| above |z = np.cos(1.5_x)_np.cos(1.5*y)|

Nelle, please look at the examples below the lines in PEP 8 that Ian
quoted. They make it perfectly clear that PEP 8 does not mandate spaces
around binary operators. Personally, in this case I think the most
readable option would be to put spaces only around the middle operator.
The goal is clarity and readability, which sometimes requires
judgement. There is no such thing as "strict PEP 8 compliance"; it is
an oxymoron.

Eric

NelleV
NelleV Jan 9, 2013 Contributor

This is a very recent change in PEP8 I was not aware of.

Else, yes, there are stricter ways to consider the document than others (I don't consider Twisted as following PEP8, and others might). I also believe the less the developper need to think and the more tools are here to do the work and check whether the code is valid, the better it is. When it comes to styling, it's not a matter of judgement, but of culture and personal preferences, which one should generally avoid in a project.

A fix has been integrated to pep8 (the tool) to reflect the change in the PEP.

@NelleV NelleV and 1 other commented on an outdated diff Jan 8, 2013
lib/matplotlib/tests/test_triangulation.py
@@ -131,3 +131,122 @@ def test_no_modify():
tri = mtri.Triangulation(points[:,0], points[:,1], triangles)
edges = tri.edges
assert_array_equal(old_triangles, triangles)
+
+def test_trifinder():
+ # Test points within triangles of masked triangulation.
+ x,y = np.meshgrid(np.arange(4), np.arange(4))
+ x = x.ravel()
+ y = y.ravel()
+ triangles = [[0,1,4], [1,5,4], [1,2,5], [2,6,5], [2,3,6], [3,7,6], [4,5,8],
NelleV
NelleV Jan 8, 2013 Contributor

Same as above.

ianthomas23
ianthomas23 Jan 8, 2013 Member

Will change.

@NelleV NelleV commented on the diff Jan 8, 2013
lib/matplotlib/tri/triinterpolate.py
+ if trifinder is not None and not isinstance(trifinder, TriFinder):
+ raise ValueError('Expected a TriFinder object')
+ self._trifinder = trifinder or self._triangulation.get_trifinder()
+
+
+class LinearTriInterpolator(TriInterpolator):
+ """
+ A LinearTriInterpolator performs linear interpolation on a triangular grid.
+
+ Each triangle is represented by a plane so that an interpolated value at
+ point (x,y) lies on the plane of the triangle containing (x,y).
+ Interpolated values are therefore continuous across the triangulation, but
+ their first derivatives are discontinuous at edges between triangles.
+ """
+ def __init__(self, triangulation, z, trifinder=None):
+ """
NelleV
NelleV Jan 8, 2013 Contributor

I've noticed matplotlib has it's own guidelines from docstrings (by that I mean that it doesn't uses numpy's guideline). Are these documented somewhere ?

mdboom
mdboom Jan 8, 2013 Owner

It uses "raw" Sphinx/rst without the Numpy extensions. To the extent that they are documented, they are documented in the Sphinx documentation. But there is no sense of structure as defined in the Numpy docstring standard. MEP10 hopes to address this in the future.

NelleV
NelleV Jan 8, 2013 Contributor

I can help on this MEP, both with the finalization of the document and the implementation. Is there anyone in charge I should ask?

mdboom
mdboom Jan 8, 2013 Owner

I wrote the original MEP -- I just haven't had much time to see it through. If you're game, I'd love to see you take this on. We can continue the discussion on matplotlib-devel if you'd like...

@ianthomas23 ianthomas23 merged commit e4d4984 into matplotlib:master Jan 10, 2013

1 check failed

default The Travis build failed
Details
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment