多个 3D 阵列（曲线）的采样/插值

Question

我的目标是在 2D 和 3D 空间中以预定义的连续距离对曲线进行插值，以对多条曲线执行 PCA。

假设一个数据框包含多个 3D 数组（每个数组的大小不同）：

>>> df.curves
0    [[0.0, 0.0, 0.91452991453, 0.91452991453, 1.0]...
1    [[0.0, 0.0, 0.734693877551, 0.734693877551, 1....
2    [[0.0, 0.0, 1.0, 1.0, 1.0], [0.0, 0.6435643564...
3    [[0.0, 0.0, 0.551020408163, 0.551020408163, 1....
4    [[0.0, 0.0, 1.0, 1.0, 1.0], [0.0, 0.4389027431...
5    [[0.0, 0.0, 0.734693877551, 0.734693877551, 1....
Name: curves, dtype: object

>>> df.curves[0]
array([[ 0.        ,  0.        ,  0.73469388,  0.73469388,  1.        ],
       [ 0.        ,  0.1097561 ,  0.47560976,  0.5       ,  1.        ],
       [ 1.        ,  0.65036675,  0.08801956,  0.06845966,  0.        ]])

我们将维度命名为 x,y,z，其中所有维度的长度相同，x 和 y 维度单调递增：

3D 绘图

我尝试对数据进行等距采样，并允许具有统一采样率的曲线之间具有可比性。

一个简单的二维曲线采样函数（没有y dim）将是每个数据帧行：

def sample2DCurve(row, res=10, method='linear'):    
    # coords of interpolation
    xnew = np.linspace(0, 1, res)

    # call scipy interpolator interp1d
    # create interpolation function for 2D data
    sample2D = interpolate.interp1d(row[0], row[1], kind=method)

    # sample data points based on xnew
    znew = sample2D(xnew)

    return np.array([xnew, znew])

对于 3D 数据，我沿路径使用插值：

def sample3DCurves(row, res=10, method='linear'):
    #npts = row[0].size
    #p = np.zeros(npts, dtype=float)
    #for i in range(1, npts):
    #    dx = row[0][i] - row[0][i-1]
    #    dy = row[1][i] - row[1][i-1]
    #    dz = row[2][i] - row[2][i-1]
    #    v = np.array([dx, dy, dz])
    #    p[i] = p[i-1] + np.linalg.norm(v)
#==============================================================================
    # edit: cleaner algebra
    x, *y, z = row

    # vecs between subsequently measured points
    vecs = np.diff(row)

    # path: cum distance along points (norm from first to ith point)
    path = np.cumsum(np.linalg.norm(vecs, axis=0))
    path = np.insert(path, 0, 0)
#==============================================================================

    ## coords of interpolation
    coords = np.linspace(p[0], p[-1], res) #p[0]=0 p[-1]=max(p)

    # interpolation func for each axis with the path
    sampleX = interpolate.interp1d(p, row[0], kind=method)
    sampleY = interpolate.interp1d(p, row[1], kind=method)
    sampleZ = interpolate.interp1d(p, row[2], kind=method)

    # sample each dim
    xnew = sampleX(coords)
    ynew = sampleY(coords)
    znew = sampleZ(coords)

    return np.array([xnew, ynew, znew])

作为 3D 中的另一种方法，我想沿等值线执行插值，在 x,y-plane 中形成具有均匀半径的圆：

x,y-plane 中 [0,0,0] 周围的圆形等值线与 3D 相交

然后基于等值线与投影在x,y 平面中的（线性）内插曲线的交点对 z 值进行内插。

但我很难定义圆形线并与x,y-plane 中的曲线/路径矢量的投影相交。

非常感谢任何建议！ （也有其他语言 - R/Matlab 等）

Answer 1

对于我的后人来说，我的粗略（一个更简单和更 Pythonic 的代码非常感谢）解决方案。

如果这实际上是主成分分析研究曲线形状的有用预处理步骤，则保留问题/分析。

def seq_sampling(row, res=10, method='linear'):
    #3D sequential along x and y (isocircles):
    x, y, z = row

    #  distance to origin for each point (support vectors lengths)
    point_distance = np.linalg.norm(row[(0,2),], axis=0)

    # isocircle radii
    max_radius = math.sqrt(x[-1]+y[-1])
    radii = np.linspace(0, max_radius, res)

    # last (distance to origin) inner data points per circle (start point of segments)
    start_per_radius = [np.max(np.where(point_distance <= radius)) for radius in radii]

    # initialize coords
    new_x = np.zeros_like(radii)
    new_y = np.zeros_like(radii)
    new_z = np.zeros_like(radii)

    # assign first an last known coordinates
    new_x[0], new_y[0] = x[0], y[0] # 0, 0
    new_x[-1], new_y[-1] = x[-1], y[-1] # 1, 1

    for radius, startpoint in enumerate(start_per_radius[1:-1]):
        # intersect circles of radius with corresponding intersecting vectors
        # based on https://math.stackexchange.com/questions/311921/get-location-of-vector-circle-intersection
        # fix index count (starts with radius > 0)
        radius += 1

        # span line segment with point O outside and point I inside of iso-circle 
        endpoint = startpoint+1

        O_x = x[endpoint]
        O_y = y[endpoint]

        I_x  = x[startpoint]
        I_y  = y[startpoint]

        # coefficients
        a = (O_x-I_x)**2 + (O_y-I_y)**2
        b = 2*((O_x-I_x)*(I_y) + (O_y-I_y)*(I_y))
        c = (I_x)**2 + (I_y)**2 - radii[radius]**2

        # !radicant cannot be zero given:
        # each segment is defined by max point lying inside or on iso-circle and the next point
        # as both axis are monotonically (strict monotonically y) increasing the next point lies outside of the ico-circle
        # thus (in 2D) a segment is intersecting a circle by definition.
        t = 2*c / (- b - math.sqrt(b**2 - 4*a*c))

        #check if intersection lies on line segment / within boundaries of t [0,1]
        if (t >= 0) and (t <= 1):
        new_x[radius], new_y[radius] = (O_x - I_x)*t + I_x, (O_y-I_y)*t + I_y

    # interpolate new_y based on projected new_y 
    new_z = interpolate.interp1d(y, z, kind='linear')(new_y[1:-1])

    # assign first an last known coordinates
    new_z = np.insert(new_z,0,z[0])
    new_z = np.append(new_z,z[-1])

    return  np.array([new_x, new_y, new_z])

---