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

- Part 1: Planar Triangulations
- Part 2: Degeneracy and Sampling on a Rectangular Grid
- Part 3: Triangulation on the Surface of a Sphere and Triangulations in 3D
- Part 4: 3D Interpolation

## 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.

## 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

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]

## 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}

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.