测量连续点的曲率

Question

我有一个点列表（数万数量级），我需要使用python识别两件事：

1- 这些点之间的连续点组（abs(x2-x1)

2-各组的曲率/半径

这是一组样本点：

[[331, 400], [331, 1200], [332, 400], [332, 486], [332, 522], [332, 655], [332, 1200], [332, 3800], [332, 3877], [332, 3944], [332, 3963], [332, 3992], [332, 4050], [333, 400], [333, 486], [333, 522], [333, 560], [333, 588], [333, 655], [333, 700], [333, 1200], [333, 3800], [333, 3877], [333, 3944], [333, 3963], [333, 3992], [333, 4050], [334, 400], [334, 486], [334, 522], [334, 558], [334, 586], [334, 654], [334, 697], [334, 1200], [334, 3800], [334, 3877], [334, 3944], [334, 3963], [334, 3992], [334, 4050], [335, 400], [335, 486], [335, 521], [335, 556], [335, 585], [335, 653], [335, 695], [335, 1200], [335, 3800], [335, 3877], [335, 3944], [335, 3963], [335, 3992], [335, 4050], [336, 400], [336, 486], [336, 520], [336, 555], [336, 584], [336, 651], [336, 693], [336, 1200], [336, 3800], [336, 3877], [336, 3944], [336, 3963], [336, 3992], [336, 4050], [337, 400], [337, 486], [337, 554], [337, 583], [337, 649], [337, 692], [337, 1200], [337, 3800], [337, 3877], [337, 3944], [337, 3963], [337, 3992], [337, 4050], [338, 377], [338, 400], [338, 486], [338, 553], [338, 582], [338, 647], [338, 691], [338, 1200], [338, 3800], [338, 3877], [338, 3944], [338, 3963], [338, 3992], [338, 4050], [339, 377], [339, 400], [339, 486], [339, 553], [339, 581], [339, 585], [339, 644], [339, 654], [339, 690], [339, 706], [339, 1200], [339, 3800], [339, 3877], [339, 3944], [339, 3963], [339, 3992], [339, 4050], [340, 376], [340, 400], [340, 486], [340, 552], [340, 580], [340, 585], [340, 641], [340, 655], [340, 689], [340, 713], [340, 1200], [340, 3800], [340, 3877], [340, 3944], [340, 3963], [340, 3992], [340, 4050], [341, 376], [341, 400], [341, 486], [341, 552], [341, 579], [341, 585], [341, 639], [341, 655], [341, 688], [341, 715], [341, 1200], [341, 3800], [341, 3877], [341, 3944], [341, 3963], [341, 3992], [341, 4050], [342, 375]、[342、400]、[342、486]、[342、552]、[342、578]、[342、585]、 [342, 637], [342, 655], [342, 688], [342, 717], [342, 1200], [342, 3800], [342, 3858], [342, 3925], [342, 3954], [342, 4011], [342, 4050], [342, 4107], [343, 374], [343, 400], [343, 486], [343, 521], [343, 552], [343, 577], [343, 585], [343, 635], [343, 642], [343, 687], [343, 718], [343, 1200], [343, 3800], [343, 3858], [343, 3925], [343, 3954], [343, 4011], [343, 4050], [343, 4107], [344, 373], [344, 400], [344, 486], [344, 521], [344, 552], [344, 576], [344, 585], [344, 633], [344, 642], [344, 687], [344, 719], [344, 1200], [344, 3800], [344, 3858], [344, 3925], [344, 3954], [344, 4011], [344, 4050], [344, 4107], [345, 372], [345, 400], [345, 486], [345, 521], [345, 552], [345, 575], [345, 585], [345, 630], [345, 642], [345, 687], [345, 720], [345, 1200], [345, 3800], [345, 3858], [345, 3925], [345, 3954], [345, 4011], [345, 4050], [345, 4107], [346, 370], [346, 400], [346, 486], [346, 521], [346, 552], [346, 574], [346, 585], [346, 628], [346, 642], [346, 686], [346, 721], [346, 1200], [346, 3800], [346, 3858], [346, 3925], [346, 3954], [346, 4011], [346, 4050], [346, 4107], [347, 368], [347, 400], [347, 486], [347, 521], [347, 552], [347, 572], [347, 585], [347, 626], [347, 642], [347, 686], [347, 721], [347, 1200], [347, 3800], [347, 3858], [347, 3925], [347, 3954], [347, 4011], [347, 4050], [347, 4107], [348, 366], [348, 400], [348, 487], [348, 521], [348, 552], [348, 570], [348, 585], [348, 624], [348, 642], [348, 686], [348, 721], [348, 1200], [348, 3800], [348, 3858], [348, 3925], [348, 3954], [348, 4011], [348, 4050]、[348、4107]、[349、364]、[349、400]、[349、487]、[349、521]、 [349, 553], [349, 568], [349, 585], [349, 622], [349, 642], [349, 686]、[349、722]、[349、1200]、[349、3800]、[349、3858]、[349、3925]、 [349, 3954], [349, 4011], [349, 4050], [349, 4107], [350, 362], [350, 400], [350, 487], [350, 521], [350, 553], [350, 585], [350, 619], [350, 642], [350, 686], [350, 722], [350, 1200], [350, 3800], [350, 3858], [350, 3925], [350, 3954], [350, 4011], [350, 4050], [350, 4107], [351, 357], [351, 400], [351, 487], [351, 521], [351, 554], [351, 585], [351, 619], [351, 642], [351, 686], [351, 722], [351, 1200], [351, 3800], [351, 3819], [351, 3858], [351, 3877], [351, 3915], [351, 3934], [351, 3963], [351, 3992], [351, 4050], [351, 4069]、[351、4107]、[352、355]、[352、373]、[352、400]、[352、487]、 [352, 520], [352, 555], [352, 585], [352, 621], [352, 642], [352, 686]、[352、722]、[352、1200]、[352、3800]、[352、3819]、[352、3858]、 [352, 3877], [352, 3915], [352, 3934], [352, 3963], [352, 3992], [352, 4050], [352, 4069], [352, 4107], [353, 353], [353, 375], [353, 400], [353, 487], [353, 520], [353, 556], [353, 585], [353, 623], [353, 642], [353, 686], [353, 722], [353, 1200], [353, 3800], [353, 3819], [353, 3858], [353, 3877], [353, 3915], [353, 3934], [353, 3963], [353, 3992], [353, 4050], [353, 4069], [353, 4107], [354, 351], [354, 376], [354, 400], [354, 487], [354, 520], [354, 558], [354, 584], [354, 625], [354, 642], [354, 686], [354, 721], [354, 1200], [354, 3800], [354, 3819], [354, 3858], [354, 3877]]

Answer 1

这将为您提供集群和list of angles：

from sklearn.cluster import DBSCAN
from scipy.spatial import distance
from scipy.optimize import curve_fit
import numpy as np, math
data = [[331, 400], [331, 1200], [332, 400], [332, 486], [332, 522]] #....

def angle(pt1, pt2):
    x1, y1 = pt1
    x2, y2 = pt2
    inner_product = x1*x2 + y1*y2
    len1 = math.hypot(x1, y1)
    len2 = math.hypot(x2, y2)
    return math.acos(inner_product/(len1*len2))

db=DBSCAN(eps=1,min_samples=2,metric='precomputed').fit(
  distance.squareform(distance.pdist(data)))
core_samples = db.core_sample_indices_
labels = db.labels_
n_clusters_ = len(set(labels)) - (1 if -1 in labels else 0)
unique_labels = set(labels)

for k in unique_labels:
    class_members = [index[0] for index in np.argwhere(labels == k)]
    cluster_core_samples = [index for index in core_samples
                            if labels[index] == k]
    curve = np.array([data[index] for index in class_members])
    print k, curve, [angle(p1,p2) for p1,p2 in zip(curve,curve[1:])]

Answer 2

得到簇后，你必须得到角度。

现在获取角度有点棘手，因为正如我在评论中所说，它不太可能保持不变。

所以我打算做的是将每个集群转换为多项式并得到它的导数 这样你就可以得到给定点的角度，所以首先polyfit。

polyfit 将获取点列表并将其转换为给定顺序的多项式。本质上，订单越高，您获得的准确性就越高，但它有它的问题。您可能需要进行试验。

所以

z = np.polyfit(ClusterXCoordinates, ClusterYCoordinates, 5)// last part are the degree, this is where you experiment 

没有一些ploy1d魔法

p = np.poly1d

我们取导数 现在得到曲线中两点之间的角度

point = p(x_coordinate)//will give you the y of this cord after polyfit so it's continues and derivable  

所以请选择derivative

// remember p computes the value of the function
    def df(x, h=0.1e-5):
//All above code to create p
       return ( p(x+h/2) - p(x-h/2) )/h
//Get slope
slope = df(pointToCheckFromCluster)
//calculate angle in rad by defintion

math.atan(slope);

http://docs.scipy.org/doc/numpy/reference/generated/numpy.poly1d.html#numpy.poly1d