Skip to content

fd_interpolate

fd_interpolate(P, gs, h, corner=None)

Bi/Trilinear interpolation matrix

Given a regular finite-difference grid described by the number of nodes on each side, the grid spacing, and the location of the bottom-left-front-most corner node, and a list of points, construct a sparse matrix of bilinear interpolation weights so that P = W @ x

Parameters:

Name Type Description Default
P numpy double array

Matrix of interpolated point coordinates

 required
gs numpy int array

Grid size [nx,ny(,nz)]

 required
h numpy double array

Spacing between grid points [hx,hy(,hz)]

 required
corner

Location of the bottom-left-front-most corner node

 None

Returns:

Name Type Description
W scipy sparse.csr_matrix

Sparse matrix such that if x are the grid nodes, P = W @ x
See Also

regular_square_mesh.

Notes

The ordering in the output is consistent with the mesh built in regular_square_mesh

Examples:

# Grid parameters
gs = np.array([19,15])
h = 1.0/(gs-1)
# Build a grid
x, y = np.meshgrid(np.linspace(0,1,gs[0]),np.linspace(0,1,gs[1]))
V = np.concatenate((np.reshape(x,(-1, 1)),np.reshape(y,(-1, 1))),axis=1)
# Random set of points
P = np.random.rand(10,2)

# Random grid corner
corner = np.random.rand(2)

# Displace by corner
P = P + np.tile(corner,(P.shape[0],1))
V = V + np.tile(corner,(V.shape[0],1))
# Interpolation matrix
from gpytoolbox import fd_interpolate
W = fd_interpolate(P,gs=gs,h=h,corner=corner)
# Now, P = W @ V
Source code in src/gpytoolbox/fd_interpolate.py 
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
def fd_interpolate(P,gs,h,corner=None):
    """Bi/Trilinear interpolation matrix

    Given a regular finite-difference grid described by the number of nodes on each side, the grid spacing, and the location of the bottom-left-front-most corner node, and a list of points, construct a sparse matrix of bilinear interpolation weights so that P = W @ x

    Parameters
    ----------
    P : numpy double array 
        Matrix of interpolated point coordinates
    gs : numpy int array
        Grid size [nx,ny(,nz)]
    h : numpy double array
        Spacing between grid points [hx,hy(,hz)]
    corner: numpy double array (optional, default None)
        Location of the bottom-left-front-most corner node

    Returns
    -------
    W  : scipy sparse.csr_matrix
        Sparse matrix such that if x are the grid nodes, P = W @ x

    See Also
    --------
    regular_square_mesh.

    Notes
    -----
    The ordering in the output is consistent with the mesh built in regular_square_mesh

    Examples
    --------
    ```python
    # Grid parameters
    gs = np.array([19,15])
    h = 1.0/(gs-1)
    # Build a grid
    x, y = np.meshgrid(np.linspace(0,1,gs[0]),np.linspace(0,1,gs[1]))
    V = np.concatenate((np.reshape(x,(-1, 1)),np.reshape(y,(-1, 1))),axis=1)
    # Random set of points
    P = np.random.rand(10,2)

    # Random grid corner
    corner = np.random.rand(2)

    # Displace by corner
    P = P + np.tile(corner,(P.shape[0],1))
    V = V + np.tile(corner,(V.shape[0],1))
    # Interpolation matrix
    from gpytoolbox import fd_interpolate
    W = fd_interpolate(P,gs=gs,h=h,corner=corner)
    # Now, P = W @ V
    ```
    """

    dim = P.shape[1]

    if corner is None:
        corner = np.zeros(dim)

    indeces = np.floor( (P - np.tile(corner,(P.shape[0],1)))/np.tile(h,(P.shape[0],1)) ).astype(int)
    I = linspace(0,P.shape[0]-1,P.shape[0],dtype=int)
    if dim==2:
        indeces_vectorized = indeces[:,0] + gs[0]*indeces[:,1]
        I = np.concatenate((I,I,I,I))
        J = np.concatenate(( indeces_vectorized,indeces_vectorized+gs[0],indeces_vectorized+1,indeces_vectorized+1+gs[0] ))
    else:
        indeces_vectorized = (indeces[:,1] + gs[1]*indeces[:,2])*gs[0] + indeces[:,0]
        I = np.concatenate((I,I,I,I,I,I,I,I))
        J = np.concatenate(( indeces_vectorized,indeces_vectorized+gs[0],indeces_vectorized+1,indeces_vectorized+1+gs[0], indeces_vectorized+(gs[1]*gs[0]),indeces_vectorized+gs[0]+(gs[1]*gs[0]),indeces_vectorized+1+(gs[1]*gs[0]),indeces_vectorized+1+gs[0]+(gs[1]*gs[0]) ))

    # Position in the bottom left corner
    vij = np.tile(corner,(P.shape[0],1)) + indeces*np.tile(h,(P.shape[0],1))
    # Normalized position inside cell
    vij = (P - vij)/np.tile(h,(P.shape[0],1))
    # Coefficients wrt to each corner of cell
    if dim==2:
        coeff_00 = (1-vij[:,1])*(1-vij[:,0]) # bottom left
        coeff_01 = (1-vij[:,1])*vij[:,0] # bottom right
        coeff_10 = vij[:,1]*(1-vij[:,0]) # top left
        coeff_11 = vij[:,1]*vij[:,0] # top right
        # concatenate (watch that order is consistent with J!!)
        vals = np.concatenate((coeff_00,coeff_10,coeff_01,coeff_11))
        mat_dim = gs[0]*gs[1]
    else:
        coeff_000 = (1-vij[:,1])*(1-vij[:,0])*(1-vij[:,2]) # bottom left
        coeff_010 = (1-vij[:,1])*vij[:,0]*(1-vij[:,2]) # bottom right
        coeff_100 = vij[:,1]*(1-vij[:,0])*(1-vij[:,2]) # top left
        coeff_110 = vij[:,1]*vij[:,0]*(1-vij[:,2]) # top right
        coeff_001 = (1-vij[:,1])*(1-vij[:,0])*vij[:,2] # bottom left
        coeff_011 = (1-vij[:,1])*vij[:,0]*vij[:,2] # bottom right
        coeff_101 = vij[:,1]*(1-vij[:,0])*vij[:,2] # top left
        coeff_111 = vij[:,1]*vij[:,0]*vij[:,2] # top right
        # concatenate (watch that order is consistent with J!!)
        vals = np.concatenate((coeff_000,coeff_100,coeff_010,coeff_110,coeff_001,coeff_101,coeff_011,coeff_111))
        mat_dim = gs[0]*gs[1]*gs[2]
    # Build scipy matrix
    # to-do: maybe add warning if this happens?
    I = I[(J>=0)*(J<mat_dim)]
    vals = vals[(J>=0)*(J<mat_dim)]
    J = J[(J>=0)*(J<mat_dim)]
    W = csr_matrix((vals,(I,J)),shape=(P.shape[0],mat_dim))
    return W