德莱尼三角剖分的欧几里德距离 - Scipy

jon*_*boy 5 python delaunay scipy pandas

spatial导入的包Scipy可以测量指定点之间的欧几里德距离。是否可以使用Delaunay包装返回相同的测量值?使用df下面的方法,按 分组测量所有点之间的平均距离Time。但是,我希望使用 Delaunay 三角测量来测量平均距离。

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial import Delaunay

df = pd.DataFrame({
    'Time' : [1,1,1,1,2,2,2,2],                  
    'A_X' : [5, 5, 6, 6, 4, 3, 3, 4], 
    'A_Y' : [5, 6, 6, 5, 5, 6, 5, 6],                         
        })

def make_points(x):
    return np.array(list(zip(x['A_X'], x['A_Y'])))

points = df.groupby("Time").apply(make_points)

for p in points:
    tri = Delaunay(p)
    ax.triplot(*p.T, tri.simplices)
Run Code Online (Sandbox Code Playgroud)

可以使用下面的方法测量所有点之间的平均距离,但我希望包含 Delaunay。

 avg_dist = (df.groupby(['Time'])
             .apply(lambda x: spatial.distance.pdist
             (np.array(list(zip(x['A_X'], x['A_Y']))))
             .mean() if len(x) > 1 else 0)
             .reset_index()
             )
Run Code Online (Sandbox Code Playgroud)

预期输出:

   Time         0
0     1  1.082842
1     2  1.082842
Run Code Online (Sandbox Code Playgroud)

meT*_*sky 5

你可以试试这个功能

from itertools import combinations
import numpy as np
    
def edges_with_no_replacement(points):
    
    # get the unique coordinates
    points = np.unique(points.loc[:,['A_X','A_Y']].values,return_index=False,axis=0)
    if len(points) <= 1: return 0
    # for two points, no triangle
    # I think return the distance between the two points make more sense? You can change the return value to zero.
    if len(points) == 2: return np.linalg.norm(points[0]-points[1])
    
    tri = Delaunay(points)
    triangles = tri.simplices
    # get all the unique edges 
    all_edges = set([tuple(sorted(edge)) for item in triangles for edge in combinations(item,2)])
    # compute the average dist 
    return np.mean([np.linalg.norm(points[edge[0]]-points[edge[1]]) for edge in all_edges])
Run Code Online (Sandbox Code Playgroud)

该函数将首先找到给定三角形的所有唯一边,然后返回三角形边的平均长度。应用此功能

avg_dist = (df.groupby(['Time']).apply(edges_with_no_replacement).reset_index())
Run Code Online (Sandbox Code Playgroud)

输出是

    Time    0
0   1   1.082843
1   2   1.082843
Run Code Online (Sandbox Code Playgroud)

请注意,如果点在同一条线上,该函数edges_with_no_replacement仍然会抛出QhullError,例如

Delaunay(np.array([[1,2],[1,3],[1,4]]))
Run Code Online (Sandbox Code Playgroud)

所以,你必须确保这些点不在同一条线上。