3D Interpolation

This post is the fourth part of a multiple part series covering triangulation and interpolation using Igor Pro 7.

Gridded data

When the data are sampled on a rectangular grid you can use variations on linear interpolation depending on your application.  In the following examples we use the same data set:

Make/B/U/N=(11,13,15) ddd=gauss(x,5,3,y,6,3,z,7,3)*1e5

If we sample the XY plane at z=10 and display it as an image:

MatrixOp p10=ddd[][][10]
NewImage p10

we get an image that is pixelated because it is a low resolution sample of 8-bit integer data.

Figure 1: Low resolution pixelated image of p10.

Figure 1: Low resolution pixelated image of p10.

Interpolation Examples for gridded data

If you want to get a single interpolated value for an arbitrary XYZ point inside the volume defined by ddd, e.g., (x=1.1, y=2.25, z=3.7) you can use:

Print Interp3D(ddd,1.1, 2.25, 3.7)
25.48

If you need to interpolate a whole set of points that are not sampled on a rectangular grid, load their XYZ positions on a triplet wave and use Interp3DPath.  For example, we create a path containing 4 XYZ triplets:

Make/N=(4,3) myPath
myPath[0][0]= {1.7,3.2,5.7,2.75}
myPath[0][1]= {3.4,4.1,3.9,5.45}
myPath[0][2]= {2.8,6.4,4.9,7.89}

Perform the interpolation:

Interp3DPath ddd,myPath
Print W_Interpolated
  W_Interpolated[0]= {33.564,154.208,138.373,163.559}

If you need to interpolate the data to a new 3D wave of dimensions (100,150,200):

Create the wave:

Make/O/N=(100,150,200) new3DWave

Set the wave scaling to match the range of the original 3D wave:

CopyScales/i ddd,new3DWave

and finally compute the interpolated values:

new3DWave=Interp3D(ddd,x,y,z)

To find the interpolated plane corresponding to ddd[][][10] we first compute the new index using:

Print scaletoindex(new3DWave,10,2)
  142

MatrixOP/O p10Interp=new3DWave[][][142]
NewImage p10Interp
Figure 2: Image of layer 142 of new3DWave which is the interpolated form of layer 10 of the original data.

Figure 2: Image of layer 142 of new3DWave which is the interpolated form of layer 10 of the original data.

Note that the last row of p10Interp contains NaN values (shown in white above).  Normally interpolation routines return a NaN when the interpolation point is outside the domain of the data. In this case, the last row is 99 which corresponds to x value:

Print IndexToScale(new3DWave,99,0)
 10

But x=10 should be part of the domain of the data, e.g.,

Print Interp3D(ddd,10,0,10)
 5

The reason for the NaN becomes apparent when executing instead

Printf "%.16f" IndexToScale(new3DWave,99,0)
 10.0000000000000018

In other words, since Igor stores the wave scaling as two parameters, start and offset, copying the scales adjusts the start and offset but they are subject to floating point roundoff with the result that row 99 of new3DWave corresponds to x that is slightly larger to be out of bounds for the interpolation.

A pragmatic solution to this problem is to reduce the size of the offset parameter, e.g.,

SetScale/P x 0,0.1010101010101,"", new3DWave

Re-executing the interpolation:

new3DWave=Interp3D(ddd,x,y,z)
MatrixOP/O p10Interp=new3DWave[][][142]
Figure 3: Corrected interpolation to avoid some NaNs on the boundary.

Figure 3: Corrected interpolation to avoid some NaNs on the boundary.

3D Scatter data interpolation

Scatter data in this context takes the form of a 4-column wave where the first three columns contain the XYZ positions where the data were sampled and the last column contains the scalar values associated with each position.

The interpolation of 3D scatter data provides a continuous scalar value for any point inside the convex domain defined by the set of XYZ locations.  Igor returns NaN for points outside the convex domain.

In the following example I create a random sample of 500 points from the same data set we used above:

Make/o/n=(500,4) eee
eee[][0]=10*(enoise(1)+1)/2
eee[][1]=12*(enoise(1)+1)/2
eee[][2]=14*(enoise(1)+1)/2
eee[][3]=gauss(eee[p][0],5,3,eee[p][1],6,3,eee[p][2],7,3)*1e5

Next I compute the 3D triangulation of the data.  In this case the volume is subdivided into tetrahedra (not triangles).

Triangulate3d eee

The default result of Triangulate3D is stored in the wave M_3dVertexList.

If you want to interpolate only the plane corresponding to z=10 you can execute:

Interpolate3d/rngx={0,0.1,100}/rngy={0,0.1,120}/rngz={10,0.1,1} /Dest=sInterp triangulationwave=M_3dVertexList ,srcwave=eee
NewImage sInterp
ModifyImage sInterp ctab= {*,*,Rainbow256,0}
Figure 4: 3D interpolation of scatter data.

Figure 4: 3D interpolation of scatter data.

Note that you have NaNs in the corners suggesting that the random sampling of 500 points did not include edge points in this plane.  The interpolation, as in the 2D case, works by identifying the tetrahedron in which the interpolated point falls and then computing the interpolated value using the barycentric coordinates.

If you are interested in examining the details of the calculation you can execute instead:

Triangulate3d/out=2 eee

and then display the tetrahedra in Gizmo using a regular path object.

Figure 5: The tetrahedral decomposition of scatter data.  Here the density of points (500) makes it difficult to observe interior details but the absence of data from the corners is evident.

Figure 5: The tetrahedral decomposition of scatter data.  Here the density of points (500) makes it difficult to observe interior details but the absence of data from the corners is evident.