2017-07-29 105 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