4

How can I implement 2-dimensional interpolation in FORTRAN where the data looks like as shown below. x and y are two co ordinates and z is a value dependent on them x is spaced evenly but y is not uniformly spaced and the maximum value of y corresponding to uniform values of x keep on increasing. Without losing much of accuracy-

  • What is the simplest method to obtain a value of z based on a given x and y?
  • What is the fastest method to obtain a value of z based on a given x and y?

Thanks SM

x    y    z
-----------
0   0     -
0   0.014 -
0   0.02  -

.....
....

0.1 0     -
0.1 0.02  -
0.1 0.03  - 

.......
.....

1.0  0     -
1.0  0.05  -
1.0  0.08  -

.......
.......
Floris
  • 45,857
  • 6
  • 70
  • 122
user1117812
  • 121
  • 2
  • 7

3 Answers3

8

I am going to assume you have already read your data into an array N x 3, following the format you gave. I am assuming you don't know ahead of time what the X spacing is - you definitely don't know the Y spacing as it varies. Thus I would recommend the following strategy:

  • Figure out the X spacing: start at the first row, and go through the X elements until you see a change in value. You now know XSTART and XSTEP - you will need those later.
  • Do a binary search in your array for value X, until you find a value XFOUND such that XFOUND < X < XFOUND + XSTART
  • Assuming you are pointing "somewhere in the list", you find the corresponding Y value - depending on whether it's bigger or smaller than the value you need, you step up or down the array until you find the first entry < Y. The corresponding values are X11, Y11, Z11. The next line in the array has X12 Y12 and Z12.
  • You need two more points before you can do the interpolation - repeat this process, looking for the "next larger value of X". This will give you XYZ21 and XYZ22
  • Now you can think about calculating the interpolated Z value. In general there are different techniques, with different accuracies:
  • "Nearest neighbor": find the closest point, and use its value of Z (simplest, least accurate)
  • "linear interpolation": find the three closest points, and do a linear interpolation of values based on the relative distances
  • "higher order estimates": to do this, you typically need to create a full connectivity mapping of the grid points, so you can do spline interpolation and get a smooth interpolation that will typically be more accurate at points between the grid points (assuming that the function being described by the samples is in fact a smooth function!)

My FORTRAN is a bit rusty - hope this is some help.

PS - potentially a simpler approach is to use the fact that the X values are already evenly spaced. This allows you to make a better interpolation. See this picture:

enter image description here

Floris
  • 45,857
  • 6
  • 70
  • 122
  • Thanks for the answer. I have been trying to implement this. But the question is how to convert to rectangular grid? The data that I have is almost parabolic where x is uniformly spaced but y varies from 0 to max then 0. – user1117812 Jul 06 '13 at 16:11
  • 1
    Are you saying the y data is not sorted? – Floris Jul 06 '13 at 23:02
  • Yes the Y data is not sorted and I don't think it can be as the value of y is dependent on x like y^2 = 4*a*x (not exactly though). – user1117812 Jul 09 '13 at 12:50
  • 3
    The method I describe only requires you to find the four "non rectangular" corner points, where `y11 < Y, y21 < Y, y12 > Y, y22 > Y`. Even for a non-sorted array you can find those points; after that, use the equations I gave to convert to rectangular grid and interpolate. Do look at @Simon's suggestion to look at GEOMPACK to obtain Delaunay triangulation - then you could use [barycentric interpolation](http://classes.soe.ucsc.edu/cmps160/Fall10/resources/barycentricInterpolation.pdf) . Doing this once to obtain values on a rectangular grid would be the way to go - then use fast interpolation. – Floris Jul 09 '13 at 13:51
  • 1
    +1: If speed is a critical issue, creating a rectangular grid and then interpolating in that will very likely be quicker than interpolating in a triangular grid. Constructing a triangulation once and interpolating in it to generate the rectangular grid would be a good approach, I think. I'd note that there is the potential to introduce additional error by interpolating in a grid that was, itself, originally constructed by interpolation. Therefore, for a case where accuracy is critical, I'd suggest retaining the triangular grid for the interpolation. – Simon Jul 10 '13 at 06:49
  • I've done this in spherical coordinates. (I have code I can share, if anyone's interested!) This will definitely be MUCH faster than @Simon's answer, but relies on assumptions about the problem geometry. Determining the indices of the four corners in this approach is likely the biggest bottleneck. If the y-values were sorted — perhaps in advance, if they stay the same? — the binary search algorithm would be an optimal way to do this. – jvriesem Jan 30 '19 at 21:04
6

Before finding the fastest way to solve this problem, I'd suggest finding a way to solve the problem. In brief, I'd suggest:

1) Find a Delaunay triangulation of the (x,y) points. Fortran code to do this is in GEOMPACK, for instance.

2) To interpolate, given the Delaunay triangulation, find the triangle that contains the point where the interpolation is to be done, then interpolate the z value, based on the location of the point relative to each of the triangle vertices.

(EDIT) I had forgotten the name of the method (if I ever knew it) but, thanks to @Floris, a good approach to interpolate in a triangle is called barycentric interpolation, which finds the interpolated value based on the ratios of the areas in the three smaller triangles into which a large triangle can be split by drawing lines from a point inside the large triangle to each of the corners of the triangle. The area of each small triangle can be found from the length of each side of the triangle using Heron's formula.

If speed improvements were needed, I think they could be obtained mainly from finding a quick way to find the triangle that contains the interpolation point but I'd want to profile some test runs before choosing which bits of the code to optimise.

Simon
  • 10,679
  • 1
  • 30
  • 44
  • For [x,y,z] points in 3-D space, barycentric interpolation of a point at [x,y,zi] within a 2-D [x1,y1] [x2,y2] [x3,y3] Delaunay facet is equivalent to the exact position on the 3-D Delaunay facet if the z values are elevations, yes? This is assuming there are no vertically overlapping facets. – Adam Erickson Jan 12 '16 at 22:31
  • 1
    An updated link for GEOMPACK: https://people.sc.fsu.edu/~jburkardt/f77_src/geompack/geompack.html. – jvriesem Jan 30 '19 at 21:07
0

I suggest using one of the TOMS interpolation algorithms: https://people.sc.fsu.edu/~jburkardt/f77_src/toms526/toms526.html https://people.sc.fsu.edu/~jburkardt/f77_src/toms660/toms660.html https://people.sc.fsu.edu/~jburkardt/f77_src/toms790/toms790.html

Eli Lansey
  • 1,130
  • 2
  • 15
  • 28