Triangulation and Interpolation on the Surface of a Sphere

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

In the previous two posts of these series I looked at triangulation and interpolation in the plane.  In this post I consider triangulation and interpolation when the data are sampled on the surface of a sphere.  Without loss of generality, we may consider a set of points {x, y, z} such that x^2+y^2+z^2=1.  

Figure 1 shows 25 random points on the surface of a sphere together with the Delaunay triangulation and the corresponding Voronoi tessellation.

Figure 1: 25 random points on the surface of a sphere.  Delaunay triangulation shows in green and Voronoi polygons are in blue.

Figure 1: 25 random points on the surface of a sphere.  Delaunay triangulation shows in green and Voronoi polygons are in blue.

The random set was generated using the function:

//  R=1, theta is [-pi/2,pi/2] and phi [0,2pi].
Function randomSphereData(num)
    Variable num
    
    Make/O/N=(num,3) regionalScatter
    Variable i,theta,phi
    
    for(i=0;i<num;i+=1)
        theta=(pi/2)*enoise(1)
        phi=pi*enoise(1)+pi
        regionalScatter[i][0]=1*sin(theta)*cos(phi)
        regionalScatter[i][1]=1*sin(theta)*sin(phi)
        regionalScatter[i][2]=1*cos(theta)
    endfor
End

You can perform the triangulation on the surface of the sphere using the command:

SphericalTriangulate regionalScatter

The operation creates the wave M_TriangulationData that contains all the relevant information for drawing Figure 1.  In most situations you would not be interested in the triangulation details. Instead, you would use the wave M_TriangulationData to interpolate scalar data on the surface of the sphere given scalar values sampled at the regionalScatter locations say {vi}.  To compute the interpolated data you need three waves:  the first wave is M_TriangulationData that you got from SphericalTriangulate.  The second wave (dataPointsWave) is a 4 column wave that combines the scatter XYZ locations (regionalScatter wave) in the first 3 columns and {vi} as the 4th column.  You can create this wave using the command:

Concatenate {regionalScatter,viWave}, dataPointsWave

The third wave is a triplet wave that specifies the XYZ target locations on the sphere where you want to compute the interpolated values.  There are really no restrictions on the target locations wave as long as the XYZ triplets lie on the surface of the sphere.  If your target is a flattened image, e.g., some projection of the sphere, you can use the required XY resolution of the image and the projection mapping to create a target triplet wave of points on the surface of the sphere.

If you are using Gizmo to display the result of the interpolation, practical consideration suggest that the target wave should match the XYZ points of a parametric surface representing the sphere.  Note that in general you can convert a parametric wave into a triplet wave using:

// Function to convert a Gizmo's parametric wave into a triplet wave.
Function parametricToTriplet(inWave)
    Wave inWave
    
    Variable np=DimSize(inWave,0)*DimSize(inWave,1)
    Make/O/N=(np,3) paramTripletWave
    MatrixOP/o/free aa=Redimension(inWave[][][0],np,1)
    paramTripletWave[][0]=aa[p]
    MatrixOP/o/free aa=Redimension(inWave[][][1],np,1)
    paramTripletWave[][1]=aa[p]
    MatrixOP/o/free aa=Redimension(inWave[][][2],np,1)
    paramTripletWave[][2]=aa[p]
End

You can perform the spherical interpolation using the command: 

SphericalInterpolate M_TriangulationData, dataPointsWave, paramTripletWave

The result is stored in the wave W_SphericalInterpolation and contains sequential scalar values corresponding to rows in paramTripletWave.  One way to display the interpolated values is to create a color wave for the parametric surface representing the sphere.  This is accomplished using the function:

// pWave is the 3D parametric wave used to draw the sphere.
// Pass W_SphericalInterpolation for interpWave.
Function createParametricColorWave(pWave,interpWave)
    Wave pWave,interpWave
    
    Variable rows=dimsize(pwave,0)
    Variable cols=dimSize(pWave,1)
    Make/O/N=(rows,cols,4) pWaveColor=1 
    Variable mmax=WaveMax(interpWave)
    Variable mmin=WaveMin(interpWave)
    ColorTab2Wave Rainbow256
    Wave M_Colors
    MatrixOP/O M_Colors=FP32(M_Colors)/65535         
    Variable nor=(DimSize(M_Colors,0)-1)/(mmax-mmin)
    Variable i, np=rows*cols,index,rr,cc
    for(i=0;i<np;i+=1)
        index=trunc((interpWave[i]-mmin)*nor)
        rr=mod(i,rows)
        cc=trunc(i/rows)
        pWaveColor[rr][cc][0]=M_Colors[index][0]
        pWaveColor[rr][cc][1]=M_Colors[index][1]
        pWaveColor[rr][cc][2]=M_Colors[index][2]
    endfor
End

The following example of interpolation on the surface of a sphere uses the procedure described above to analyze EEG data.  In this application EEG signals are measured by a standard net of sensors placed on a human head which, for the purpose of the analysis, is approximated by a spherical shape.  The known locations of the sensors replace the regionalScatter data in the example above.  The EEG signals are measured simultaneously at all sensors and their values at some time t0 correspond to the scalar set {vi} above.  The result of the interpolation is shown in Figure 2 which consists of three Gizmo subwindows displaying the same data at different orientations.

Figure 2: Three views of interpolated EEG signals

Figure 2: Three views of interpolated EEG signals