algorithm - 我的薄板样条插值实现的结果取决于自变量

Question

我实现了 thin plate spline算法(另请参见 this description)，以便使用 Python 对分散的数据进行插值。

当初始分散数据的边界框的纵横比接近 1 时，我的算法似乎可以正常工作。但是，缩放其中一个数据点坐标会改变插值结果。我创建了一个最小的工作示例，它代表了我要完成的工作。下面两幅图显示了 50 个随机点的插值结果。

首先，z = x^2在域x = [0, 3], y = [0, 120]上的插值:

如您所见，插值失败。现在，执行相同的过程，但在将 x 值缩放 40 倍后，我得到:

这一次，结果看起来更好。选择略有不同的比例因子会导致略有不同的插值。这表明我的算法有问题，但我找不到确切的问题。这是算法:

import numpy as np
import numba as nb

# pts1 = Mx2 matrix (original coordinates)
# z1   = Mx1 column vector (original values)
# pts2 = Nx2 matrix (interpolation coordinates)

def gen_K(n, pts1):
    K = np.zeros((n,n))

    for i in range(0,n):
        for j in range(0,n):
            if i != j:
                r = ( (pts1[i,0] - pts1[j,0])**2.0 + (pts1[i,1] - pts1[j,1])**2.0 )**0.5
                K[i,j] = r**2.0*np.log(r)

    return K

def compute_z2(m, n, pts1, pts2, coeffs):   
    z2 = np.zeros((m,1))

    x_min = np.min(pts1[:,0])
    x_max = np.max(pts1[:,0])
    y_min = np.min(pts1[:,1])
    y_max = np.max(pts1[:,1])

    for k in range(0,m):
        pt = pts2[k,:]

        # If point is located inside bounding box of pts1
        if (pt[0] >= x_min and pt[0] <= x_max and pt[1] >= y_min and pt[1] <= y_max):
            z2[k,0] = coeffs[-3,0] + coeffs[-2,0]*pts2[k,0] + coeffs[-1,0]*pts2[k,1]

            for i in range(0,n):
                r2 = ( (pts1[i,0] - pts2[k,0])**2.0 + (pts1[i,1] - pts2[k,1])**2.0 )**0.5

                if r2 != 0:
                    z2[k,0] += coeffs[i,0]*( r2**2.0*np.log(r2) )

        else:
            z2[k,0] = np.nan

    return z2

gen_K_nb = nb.jit(nb.float64[:,:](nb.int64, nb.float64[:,:]), nopython = True)(gen_K)
compute_z2_nb = nb.jit(nb.float64[:,:](nb.int64, nb.int64, nb.float64[:,:], nb.float64[:,:], nb.float64[:,:]), nopython = True)(compute_z2)

def TPS(pts1, z1, pts2, factor):
    n, m = pts1.shape[0], pts2.shape[0]

    P = np.hstack((np.ones((n,1)),pts1))
    Y = np.vstack((z1, np.zeros((3,1))))

    K = gen_K_nb(n, pts1)
    K += factor*np.identity(n)

    L = np.zeros((n+3,n+3))
    L[0:n, 0:n] = K
    L[0:n, n:n+3] = P
    L[n:n+3, 0:n] = P.T

    L_inv = np.linalg.inv(L)
    coeffs = L_inv.dot(Y)

    return compute_z2_nb(m, n, pts1, pts2, coeffs)

最后，这是我用来创建两个图的代码片段:

import matplotlib.pyplot as plt
import numpy as np

N = 50 # Number of random points

pts = np.random.rand(N,2)
pts[:,0] *= 3.0 # initial x values
pts[:,1] *= 120.0 # initial y values

z1 = (pts[:,0])**2.0

for scale in [1.0, 40.0]:

    pts1 = pts.copy()
    pts1[:,0] *= scale

    x2 = np.linspace(np.min(pts1[:,0]), np.max(pts1[:,0]), 40)
    y2 = np.linspace(np.min(pts1[:,1]), np.max(pts1[:,1]), 40)
    x2, y2 = np.meshgrid(x2, y2)

    pts2 = np.vstack((x2.flatten(), y2.flatten())).T
    z2 = TPS(pts1, z1.reshape(z1.shape[0], 1), pts2, 0.0)

    # Display
    fig = plt.figure(figsize=(4,3))
    ax = fig.add_subplot(111)
    C = ax.contourf(x2, y2, z2.reshape(x2.shape), np.linspace(0,9,10), extend='both')
    ax.plot(pts1[:,0], pts1[:,1], 'ok')
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    plt.colorbar(C, extendfrac=0)
    plt.tight_layout()

plt.show()

Answer 1

Thin Plate Spline 是标量不变的，这意味着如果按相同的因子缩放 x 和 y，结果应该相同。但是，如果您以不同方式缩放 x 和 y，则结果将不同。这是径向基函数的共同特征。一些径向基函数甚至不是标量不变的。