2017-07-29 72 views
3

我的目标是在预定义的连续距离内插入2D和3D空间中的曲线,以在多条曲线上执行PCA。多个3D阵列(曲线)的采样/插值

假设多个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.  ]]) 

让我们命名尺寸xyz其中所有尺寸具有相同的长度和xy尺寸是单调增加:

3D图

我尝试将数据采样为等距,并允许具有均匀采样率的曲线之间的可比性。

用于2D曲线(无y暗淡)一个简单的采样函数将每个数据帧的行:

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 - 具有统一半径的平面上沿着等值线形成圆进行插值:

Circular等值线周围[0,0,0]在x,y - 平面与3D路口

然后z值是基于与在x,y投影的(线性)插值曲线的等值线的交点内插 - 平面。

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

任何建议非常感谢! (也用其他语言 - R/Matlab等)

+0

既然你知道异圆的半径可能试图找到该曲线穿过与z轴的距离飞机?只是一个想法,有趣的问题。 – DrBwts

回答

0

对于后代我的粗糙(更简单和更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]) 

3D Plot with circular iso-lines

+0

认为你需要编辑你的缩进 – DrBwts

+1

@DrBwts谢谢。我修好了它。 – hard

+0

另外我修正了一个bug,因为当'x'没有改变时,点似乎错误地插入了段。插值现在在'new_y'上,它严格**单调(增加)。 – hard