linear interpolation of 4D waves

Share icon

Is there an equivalent to interp2D and interp3D for 4D waves? 

More specifically:

A function that returns an interpolated value for location P=(x, y, z, t) in a 4D scalar distribution srcWave.

I currently use code that extracts two adjacent chunks of a 4D wave, used interp3d on both chunks followed by linear interpolation between the two interp3d results. However, this is slow (the majority of time is spent in the chunk extraction).

NB: I'm only interested in linear interpolation of a 4D wave containing a scalar distribution sampled on a regular lattice

September 28, 2021 at 09:45 am - Permalink

Hello Fabrizio,

This is the first time in my recollection that this feature request has come up.  I will look into implementing it internally.

I'm not sure what you mean by "chunk extraction".  If you actually copy two chunks into some local arrays and apply 3D interpolation to each one you be running inefficiently.  I would consider direct element addressing in the formula for 4D interpolation.

A.G.

September 29, 2021 at 03:43 pm - Permalink

Here is an Igor procedure to perform the interpolation.  Note comments about improving efficiency if you need to call this function many times in a loop.

Function interp_4d(xx,yy,zz,tt,inWave)
    Variable xx,yy,zz,tt
    Wave inWave
   
    // May skip the following for efficiency if called from a loop:
    if(WaveExists(inWave)==0 || WaveType(inWave)==0 || WaveType(inWave)&1 || DimSize(inWave,3)<=0)
        Print "Bad input wave."
        return NaN
    endif
   
    Variable ix1,ix2,iy1,iy2,iz1,iz2,it1,it2
    Variable q1,q2,q3,q4
    Variable p1,p2,p3,p4
    Variable fx111,fx211,fx121,fx221,fxy11,fxy21,fxyz1
    Variable fx112,fx212,fx122,fx222,fxy12,fxy22,fxyz2

    // The following can be made more efficient by looping over the dimensions
    // and using waves for the various variables.  If this function is called from a loop
    // it would be useful to use a (read-only) wave containing the variables that gets passed
    // to each call.
   
    // dim0
    Variable x0=DimOffset(inWave,0)
    Variable dx=DimDelta(inWave,0)
    Variable nx=DimSize(inWave,0)
   
    ix1=floor((xx-x0)/dx)

    if(ix1<0 || ix1>=nx-1)
        return NaN
    endif
    ix2=ix1+1
    q1=(xx-(x0+ix1*dx))/dx
    p1=1-q1

    // dim1
    Variable y0=DimOffset(inWave,1)
    Variable dy=DimDelta(inWave,1)
    Variable ny=DimSize(inWave,1)
   
    iy1=floor((yy-y0)/dy)
    if(iy1<0 || iy1>=ny-1)
        return NaN
    endif
    iy2=iy1+1
    q2=(yy-(y0+iy1*dy))/dy
    p2=1-q2

    // dim2
    Variable z0=DimOffset(inWave,2)
    Variable dz=DimDelta(inWave,2)
    Variable nz=DimSize(inWave,2)
   
    iz1=floor((zz-z0)/dz)
    if(iz1<0 || iz1>=nz-1)
        return NaN
    endif
    iz2=iz1+1
    q3=(zz-(z0+iz1*dz))/dz
    p3=1-q3
   
    // dim3
    Variable t0=DimOffset(inWave,3)
    Variable dt=DimDelta(inWave,3)
    Variable nt=DimSize(inWave,3)
   
    it1=floor((tt-t0)/dt)
    if(it1<0 || it1>=nt-1)
        return NaN
    endif
    it2=it1+1
    q4=(tt-(t0+it1*dt))/dt
    p4=1-q4

    fx111=p1*inWave[ix1][iy1][iz1][it1]+q1*inWave[ix2][iy1][iz1][it1]
    fx211=p1*inWave[ix1][iy2][iz1][it1]+q1*inWave[ix2][iy2][iz1][it1]
    fx121=p1*inWave[ix1][iy1][iz2][it1]+q1*inWave[ix2][iy1][iz2][it1]
    fx221=p1*inWave[ix1][iy2][iz2][it1]+q1*inWave[ix2][iy2][iz2][it1]
    fx112=p1*inWave[ix1][iy1][iz1][it2]+q1*inWave[ix2][iy1][iz1][it2]
    fx212=p1*inWave[ix1][iy2][iz1][it2]+q1*inWave[ix2][iy2][iz1][it2]
    fx122=p1*inWave[ix1][iy1][iz2][it2]+q1*inWave[ix2][iy1][iz2][it2]
    fx222=p1*inWave[ix1][iy2][iz2][it2]+q1*inWave[ix2][iy2][iz2][it2]
    fxy11=p2*fx111+q2*fx211
    fxy21=p2*fx121+q2*fx221
    fxy12=p2*fx112+q2*fx212
    fxy22=p2*fx122+q2*fx222
    fxyz1=p3*fxy11+q3*fxy21
    fxyz2=p3*fxy12+q3*fxy22
    return p4*fxyz1+q4*fxyz2
End

 

September 30, 2021 at 12:52 pm - Permalink