PythonAlphaShape：基于点集估算图像区域的Alpha形状算法解析

好的,这是个主意. Delaunay三角剖分将产生不加区别的大三角形.它也会有问题,因为只会生成三角形.

因此,我们将生成您可能称之为“模糊Delaunay三角剖分”的东西.我们将所有点都放入kd树中,对于每个点p,查看它的k个最近邻点. kd树让这个快.

对于这些k个邻居中的每一个,找到到焦点p的距离.使用此距离生成加权.我们希望附近的点比更远的点更受青睐,所以这里指数函数exp(-alpha * dist)是合适的.使用加权距离建立概率密度函数,描述绘制每个点的概率.

现在,从该分布中抽出很多次.将经常选择附近的积分,而不太经常选择更远的积分.对于绘制的点,记下为焦点绘制的次数.结果是加权图,其中图中的每条边连接附近的点,并根据选择对的频率进行加权.

现在,从权重太小的图表中剔除所有边缘.这些是可能没有连接的点.结果如下&＃xff1a;

现在,让我们将所有剩余的边缘扔到shapely.然后我们可以通过缓冲它们将边缘转换为非常小的多边形.像这样&＃xff1a;

使用覆盖整个区域的大多边形来区分多边形将产生用于三角测量的多边形.可能还要等一下.结果如下&＃xff1a;

最后,剔除所有太大的多边形&＃xff1a;

#!/usr/bin/env python

import numpy as np

import matplotlib.pyplot as plt

import random

import scipy

import scipy.spatial

import networkx as nx

import shapely

import shapely.geometry

import matplotlib

dat &＃61; np.loadtxt(&＃39;test.asc&＃39;)

xycoors &＃61; dat[:,0:2]

xcoors &＃61; xycoors[:,0] #Convenience alias

ycoors &＃61; xycoors[:,1] #Convenience alias

npts &＃61; len(dat[:,0]) #Number of points

dist &＃61; scipy.spatial.distance.euclidean

def GetGraph(xycoors, alpha&＃61;0.0035):

kdt &＃61; scipy.spatial.KDTree(xycoors) #Build kd-tree for quick neighbor lookups

G &＃61; nx.Graph()

npts &＃61; np.max(xycoors.shape)

for x in range(npts):

G.add_node(x)

dist, idx &＃61; kdt.query(xycoors[x,:], k&＃61;10) #Get distances to neighbours, excluding the cenral point

dist &＃61; dist[1:] #Drop central point

idx &＃61; idx[1:] #Drop central point

pq &＃61; np.exp(-alpha*dist) #Exponential weighting of nearby points

pq &＃61; pq/np.sum(pq) #Convert to a PDF

choices &＃61; np.random.choice(idx, p&＃61;pq, size&＃61;50) #Choose neighbors based on PDF

for c in choices: #Insert neighbors into graph

if G.has_edge(x, c): #Already seen neighbor

G[x][c][&＃39;weight&＃39;] &＃43;&＃61; 1 #Strengthen connection

else:

G.add_edge(x, c, weight&＃61;1) #New neighbor; build connection

return G

def PruneGraph(G,cutoff):

newg &＃61; G.copy()

bad_edges &＃61; set()

for x in newg:

for k,v in newg[x].items():

if v[&＃39;weight&＃39;]

bad_edges.add((x,k))

for b in bad_edges:

try:

newg.remove_edge(*b)

except nx.exception.NetworkXError:

pass

return newg

def PlotGraph(xycoors,G,cutoff&＃61;6):

xcoors &＃61; xycoors[:,0]

ycoors &＃61; xycoors[:,1]

G &＃61; PruneGraph(G,cutoff)

plt.plot(xcoors, ycoors, "o")

for x in range(npts):

for k,v in G[x].items():

plt.plot((xcoors[x],xcoors[k]),(ycoors[x],ycoors[k]), &＃39;k-&＃39;, lw&＃61;1)

plt.show()

def GetPolys(xycoors,G):

#Get lines connecting all points in the graph

xcoors &＃61; xycoors[:,0]

ycoors &＃61; xycoors[:,1]

lines &＃61; []

for x in range(npts):

for k,v in G[x].items():

lines.append(((xcoors[x],ycoors[x]),(xcoors[k],ycoors[k])))

#Get bounds of region

xmin &＃61; np.min(xycoors[:,0])

xmax &＃61; np.max(xycoors[:,0])

ymin &＃61; np.min(xycoors[:,1])

ymax &＃61; np.max(xycoors[:,1])

mls &＃61; shapely.geometry.MultiLineString(lines) #Bundle the lines

mlsb &＃61; mls.buffer(2) #Turn lines into narrow polygons

bbox &＃61; shapely.geometry.box(xmin,ymin,xmax,ymax) #Generate background polygon

polys &＃61; bbox.difference(mlsb) #Subtract to generate polygons

return polys

def PlotPolys(polys,area_cutoff):

fig, ax &＃61; plt.subplots(figsize&＃61;(8, 8))

for polygon in polys:

if polygon.area

mpl_poly &＃61; matplotlib.patches.Polygon(np.array(polygon.exterior), alpha&＃61;0.4, facecolor&＃61;np.random.rand(3,1))

ax.add_patch(mpl_poly)

ax.autoscale()

fig.show()

#Functional stuff starts here

G &＃61; GetGraph(xycoors, alpha&＃61;0.0035)

#Choose a value that rips off an appropriate amount of the left side of this histogram

weights &＃61; sorted([v[&＃39;weight&＃39;] for x in G for k,v in G[x].items()])

plt.hist(weights, bins&＃61;20);plt.show()

PlotGraph(xycoors,G,cutoff&＃61;6) #Plot the graph to ensure our cut-offs were okay. May take a while

prunedg &＃61; PruneGraph(G,cutoff&＃61;6) #Prune the graph

polys &＃61; GetPolys(xycoors,prunedg) #Get polygons from graph

areas &＃61; sorted(p.area for p in polys)

plt.plot(areas)

plt.hist(areas,bins&＃61;20);plt.show()

area_cutoff &＃61; 150000

PlotPolys(polys,area_cutoff&＃61;area_cutoff)

good_polys &＃61; ([p for p in polys if p.area

total_area &＃61; sum([p.area for p in good_polys])