This issue is due to a limitation associated with the Delaunay-based linear interpolation technique, which is not necessarily unique. In other words, certain datasets can be triangulated in different ways and still respect the Delaunay criterion.
To illustrate, a rectangle can be triangulated by two triangles. However, the triangles can be arranged in two ways: with the diagonal downwards \ or upwards /. Both are perfectly valid Delaunay triangulations. Linear interpolation is calculated across these triangular surfaces, and thus interpolations can differ depending on the triangulation.
Furthermore, the Qhull-based "griddata" adds noise to perturb the data set each time the triangulation is performed. This means that each time a degenerate dataset is evaluated using the Qhull-based "griddata" the resulting value can potentially be different. In other words the result can potentially be different when calling "griddata" repeatedly in the same session.
We have recognized this issue and we have provided a more robust and more stable solution in R2013a, namely the "scatteredInterpolant" class. Users are encouraged to use "scatteredInterpolant" instead of "griddata".
Also, if you want to interpolate over a data set that is already on a regular grid, you could use "interp2" rather than "griddata" because the sample data is not scattered. "interp2" is not susceptible to inconsistencies because the technique is not triangulation-based.
For more information, you may find the following blog entry useful:
Additional information:
Both "griddata" and "scatteredInterpolant" can only interpolate data representing a single-valued function. That is, a given sample point (x,y) must correspond to a unique value z. For example, "griddata" cannot interpolate points on the surface of a sphere, but it can interpolate points on a hemisphere that is properly oriented to satisfy the single-valued constraint. "griddata" can handle duplicate values in close proximity; for example due to sampling overlap or noise, so long as the data approximates a single-valued function. The problem arises when we have data lying on opposite sides of the X-Y plane.
If a least-squares curve-fit approach to the grid data (similar to a surface fit) is needed instead of a nearest-neighbor direct interpolation, so as to make it more robust to "noisy" data - we recommend Lowess smoothing:
Note: Lowess smoothing fits a surface to the data, the name is misleading.
