估计由一组点（Alpha 形状？？）生成的图像区域
Estimating an area of an image generated by a set of points (Alpha shapes??)
我在显示二维图像的 example ASCII file 中有一组点。
我想估计这些点填充的总面积。这个平面内有一些地方没有被任何点填充，因为这些区域已被屏蔽掉。我想估计面积可能实用的方法是应用 凹包 或 alpha 形状 。
我尝试 this approach 找到合适的 alpha 值，然后估计面积。
from shapely.ops import cascaded_union, polygonize
import shapely.geometry as geometry
from scipy.spatial import Delaunay
import numpy as np
import pylab as pl
from descartes import PolygonPatch
from matplotlib.collections import LineCollection
def plot_polygon(polygon):
fig = pl.figure(figsize=(10,10))
ax = fig.add_subplot(111)
margin = .3
x_min, y_min, x_max, y_max = polygon.bounds
ax.set_xlim([x_min-margin, x_max+margin])
ax.set_ylim([y_min-margin, y_max+margin])
patch = PolygonPatch(polygon, fc='#999999',
ec='#000000', fill=True,
zorder=-1)
ax.add_patch(patch)
return fig
def alpha_shape(points, alpha):
if len(points) < 4:
# When you have a triangle, there is no sense
# in computing an alpha shape.
return geometry.MultiPoint(list(points)).convex_hull
def add_edge(edges, edge_points, coords, i, j):
"""
Add a line between the i-th and j-th points,
if not in the list already
"""
if (i, j) in edges or (j, i) in edges:
# already added
return
edges.add( (i, j) )
edge_points.append(coords[ [i, j] ])
coords = np.array([point.coords[0]
for point in points])
tri = Delaunay(coords)
edges = set()
edge_points = []
# loop over triangles:
# ia, ib, ic = indices of corner points of the
# triangle
for ia, ib, ic in tri.vertices:
pa = coords[ia]
pb = coords[ib]
pc = coords[ic]
# Lengths of sides of triangle
a = np.sqrt((pa[0]-pb[0])**2 + (pa[1]-pb[1])**2)
b = np.sqrt((pb[0]-pc[0])**2 + (pb[1]-pc[1])**2)
c = np.sqrt((pc[0]-pa[0])**2 + (pc[1]-pa[1])**2)
# Semiperimeter of triangle
s = (a + b + c)/2.0
# Area of triangle by Heron's formula
area = np.sqrt(s*(s-a)*(s-b)*(s-c))
circum_r = a*b*c/(4.0*area)
# Here's the radius filter.
#print circum_r
if circum_r < 1.0/alpha:
add_edge(edges, edge_points, coords, ia, ib)
add_edge(edges, edge_points, coords, ib, ic)
add_edge(edges, edge_points, coords, ic, ia)
m = geometry.MultiLineString(edge_points)
triangles = list(polygonize(m))
return cascaded_union(triangles), edge_points
points=[]
with open("test.asc") as f:
for line in f:
coords=map(float,line.split(" "))
points.append(geometry.shape(geometry.Point(coords[0],coords[1])))
print geometry.Point(coords[0],coords[1])
x = [p.x for p in points]
y = [p.y for p in points]
pl.figure(figsize=(10,10))
point_collection = geometry.MultiPoint(list(points))
point_collection.envelope
convex_hull_polygon = point_collection.convex_hull
_ = plot_polygon(convex_hull_polygon)
_ = pl.plot(x,y,'o', color='#f16824')
concave_hull, edge_points = alpha_shape(points, alpha=0.001)
lines = LineCollection(edge_points)
_ = plot_polygon(concave_hull)
_ = pl.plot(x,y,'o', color='#f16824')
我得到了这个结果，但我希望这种方法可以检测到中间的洞。
更新
这是我的真实数据的样子：
我的问题是，估计上述形状面积的最佳方法是什么？我无法弄清楚这段代码无法正常工作出了什么问题？！！任何帮助将不胜感激。
编辑：
我注意到您有自己的代码来计算 alpha 形状，
Delaunay三角形的面积就在那里，所以计算形状的面积就更容易了...
如果要将三角形添加到 alpha 形多边形中，只需添加三角形的面积。
如果你想检测孔...添加一个辅助阈值以避免添加面积大于阈值的三角形。对于此示例，值 max_area = 99999 将删除孔。
唯一的问题是您创建图形输出的方式，因为您看不到孔。
def alpha_shape(points, alpha, max_area):
if len(points) < 4:
# When you have a triangle, there is no sense
# in computing an alpha shape.
return geometry.MultiPoint(list(points)).convex_hull , 0
def add_edge(edges, edge_points, coords, i, j):
"""
Add a line between the i-th and j-th points,
if not in the list already
"""
if (i, j) in edges or (j, i) in edges:
# already added
return
edges.add( (i, j) )
edge_points.append(coords[ [i, j] ])
coords = np.array([point.coords[0]
for point in points])
tri = Delaunay(coords)
total_area = 0
edges = set()
edge_points = []
# loop over triangles:
# ia, ib, ic = indices of corner points of the
# triangle
for ia, ib, ic in tri.vertices:
pa = coords[ia]
pb = coords[ib]
pc = coords[ic]
# Lengths of sides of triangle
a = np.sqrt((pa[0]-pb[0])**2 + (pa[1]-pb[1])**2)
b = np.sqrt((pb[0]-pc[0])**2 + (pb[1]-pc[1])**2)
c = np.sqrt((pc[0]-pa[0])**2 + (pc[1]-pa[1])**2)
# Semiperimeter of triangle
s = (a + b + c)/2.0
# Area of triangle by Heron's formula
area = np.sqrt(s*(s-a)*(s-b)*(s-c))
circum_r = a*b*c/(4.0*area)
# Here's the radius filter.
# print("radius", circum_r)
if circum_r < 1.0/alpha and area < max_area:
add_edge(edges, edge_points, coords, ia, ib)
add_edge(edges, edge_points, coords, ib, ic)
add_edge(edges, edge_points, coords, ic, ia)
total_area += area
m = geometry.MultiLineString(edge_points)
triangles = list(polygonize(m))
return cascaded_union(triangles), edge_points, total_area
旧答案：
要计算不规则简单多边形的面积，您可以使用 Shoelace formula 和边界的 CCW 坐标作为输入。
如果要检测云内部的孔，则必须删除外接半径大于次级阈值的 Delaunay 三角形。
理想情况是：计算 Delaunay 三角剖分并使用当前的 alpha 形状进行过滤。然后，计算每个三角形的外接半径，去掉外接半径远大于平均外接半径的三角形。
要计算带孔的不规则多边形的面积，请对每个孔边界使用鞋带公式。以CCW（正向）顺序输入外边界，得到面积。然后按CW（负）顺序输入每个孔的边界，得到面积的（负）值。
这里有一个想法：使用 k-means clustering.
您可以在 Python 中按如下方式完成此操作：
from sklearn.cluster import KMeans
import numpy as np
import matplotlib.pyplot as plt
dat = np.loadtxt('test.asc')
xycoors = dat[:,0:2]
fit = KMeans(n_clusters=2).fit(xycoors)
plt.scatter(dat[:,0],dat[:,1], c=fit.labels_)
plt.axes().set_aspect('equal', 'datalim')
plt.gray()
plt.show()
使用您的数据，得出以下结果：
现在，您可以获取顶部簇和底部簇的凸包，并分别计算每个簇的面积。然后将这些区域相加成为它们合并区域的估计量，但巧妙地避免了中间的洞。
要微调您的结果，您可以使用聚类的数量和算法的不同启动次数（该算法是随机的，通常 运行 不止一次）。
你问，比如，两个簇是否总是在中间留下一个洞。我使用以下代码对此进行了试验。我生成一个均匀分布的点，然后切出一个随机大小和方向的椭圆来模拟一个洞。
#!/usr/bin/env python3
import sklearn
import sklearn.cluster
import numpy as np
import matplotlib.pyplot as plt
PWIDTH = 6
PHEIGHT = 6
def GetPoints(num):
return np.random.rand(num,2)*300-150 #Centered about zero
def MakeHole(pts): #Chop out a randomly orientated and sized ellipse
a = np.random.uniform(10,150) #Semi-major axis
b = np.random.uniform(10,150) #Semi-minor axis
h = np.random.uniform(-150,150) #X-center
k = np.random.uniform(-150,150) #Y-center
A = np.random.uniform(0,2*np.pi) #Angle of rotation
surviving_points = []
for pt in range(pts.shape[0]):
x = pts[pt,0]
y = pts[pt,1]
if ((x-h)*np.cos(A)+(y-k)*np.sin(A))**2/a/a+((x-h)*np.sin(A)-(y-k)*np.cos(A))**2/b/b>1:
surviving_points.append(pt)
return pts[surviving_points,:]
def ShowManyClusters(pts,fitter,clusters,title):
colors = np.array([x for x in 'bgrcmykbgrcmykbgrcmykbgrcmyk'])
fig,axs = plt.subplots(PWIDTH,PHEIGHT)
axs = axs.ravel()
for i in range(PWIDTH*PHEIGHT):
lbls = fitter(pts[i],clusters)
axs[i].scatter(pts[i][:,0],pts[i][:,1], c=colors[lbls])
axs[i].get_xaxis().set_ticks([])
axs[i].get_yaxis().set_ticks([])
plt.suptitle(title)
#plt.show()
plt.savefig('/z/'+title+'.png')
fitters = {
'SpectralClustering': lambda x,clusters: sklearn.cluster.SpectralClustering(n_clusters=clusters,affinity='nearest_neighbors').fit(x).labels_,
'KMeans': lambda x,clusters: sklearn.cluster.KMeans(n_clusters=clusters).fit(x).labels_,
'AffinityPropagation': lambda x,clusters: sklearn.cluster.AffinityPropagation().fit(x).labels_,
}
np.random.seed(1)
pts = []
for i in range(PWIDTH*PHEIGHT):
temp = GetPoints(300)
temp = MakeHole(temp)
pts.append(temp)
for name,fitter in fitters.items():
for clusters in [2,3]:
np.random.seed(1)
ShowManyClusters(pts,fitter,clusters,"{0}: {1} clusters".format(name,clusters))
考虑 K-Means 的结果：
至少在我看来，当 "hole" 将数据分成两个单独的 blob 时，使用两个集群似乎效果最差。 （在这种情况下，当椭圆的方向与包含样本点的矩形区域的两条边重叠时会发生这种情况。）使用三个聚类可以解决大部分困难。
您还会注意到 K-means 在第 1 列第 3 行以及第 3 列第 4 行上产生了一些违反直觉的结果。查看 sklearn 的聚类方法 here 显示以下比较图像：
由此看来，似乎 SpectralClustering 产生的结果与我们想要的一致。在上面的相同数据上尝试此操作可以解决提到的问题（请参阅第 1 列，第 3 行和第 3 列，第 4 行）。
上述内容表明，具有三个聚类的谱聚类应该足以应对此类大多数情况。
虽然您似乎打算做一个凹形，但这里有一条替代路线，速度非常快，我认为会给您非常稳定的读数：
创建一个以参数 (int radiusOfInfluence) 作为参数的函数。在函数内 运行 一个以那个为半径的体素过滤器。然后只需将该圆的面积 (pi*AOI^2) 乘以云中剩余点的数量。这应该会给你一个相对稳健的面积估计，并且对孔洞和奇怪的边缘非常有弹性。
需要考虑的一些事项：
-这会给你一个积极的区域超调，因为边缘超出了一个半径。对此进行调整的修改可能是 运行 统计离群值去除过滤器（在反向模式下）以获取统计边缘点。然后可以假设每个边缘点的大约一半位于形状之外，在乘以面积之前从总计数中减去找到的点数的一半。
-影响半径在很大程度上决定了此函数的孔洞检测，因为较大的影响半径将允许单个点覆盖更大的区域，而且通过调整统计离群值过滤器上的标准截止值，您可以更积极地检测内部孔洞和以这种方式调整您的区域。
它确实回避了您所追求的问题，因为这更像是假设一组合理分布的样本的射击精度/射击分组类型评估。您的方法有点假设您的外边缘点是可能的绝对限制（根据情况，这可能是一个合理的假设）
编辑------------------------
我没有时间写出示例代码，但我可以进一步解释以帮助理解。
其核心是voxel filter。非常简单，它在 x、y 坐标中设置一个种子点，然后在整个 space 上创建一个网格，该网格在用户指定的过滤器半径的两个轴上都有单位（网格间距）。在每个网格框中，它会将所有点平均为一个点。这对于这个概念非常重要，因为它几乎完全消除了重叠问题。
第二部分（逆stat outlier removal）只是一点小聪明来加强你的边缘配合。基本上，stat outlier 是通过查看每个点到其 (k) 个最近邻点的距离来消除噪声的。在为每个点生成到 k 个最近邻居的平均距离后，它会建立一个直方图和一个用户定义的参数作为保留或删除点的二进制阈值。当反转并设置为合理的截止值（~0.75 std 应该有效）时，它会删除对象主体中的所有点（即只留下边缘点）。这一点很重要的原因是，从技术上讲，这些点超出了对象边界 1 个半径。虽然有些边缘角度是锐角，有些边缘角度是钝角（即大于或小于半个溢出圆），但每个点去掉 1/2 的圆面积应该在整个对象上给你一个很好的边缘拟合改进.
请记住，归根结底，这只会给您一个数字。就压力测试而言，我建议创建已知区域的人为点云和/或创建一个图形输出来显示您放置圆圈和半圆圈的位置（如果您喜欢的话，可以朝向对象的内部）。
您想要改进此方法的旋钮是：
体素过滤器半径，每个点的影响区域（实际上可以与体素过滤器半径分开控制，尽管它们应该保持非常接近彼此），标准截止。
希望这有助于澄清，干杯！
好的，这就是我的想法。 Delaunay 三角剖分将生成不分青红皂白的大三角形。这也会有问题，因为只会生成三角形。
因此，我们将生成您可能称之为 "fuzzy Delaunay triangulation" 的内容。我们将把所有点放入 kd-tree 中，并且对于每个点 p，查看它的 k 最近邻点。 kd-tree 让这个变得很快。
对于每个 k 个邻居，找出到焦点 p 的距离。使用此距离生成权重。我们希望附近的点比更远的点更受青睐，因此指数函数 exp(-alpha*dist) 在这里是合适的。使用加权距离构建描述绘制每个点的概率的概率密度函数。
现在，从那个分布中抽取很多次。附近的点将经常被选择，而较远的点将被较少地选择。对于绘制的点，记下它为焦点绘制了多少次。结果是一个加权图，其中图中的每条边都连接附近的点，并根据选择对的频率进行加权。
现在，剔除图中权重太小的所有边。这些是可能没有连接的点。结果如下所示：
现在，让我们把所有剩余的边都扔进shapely。然后我们可以通过缓冲将边转换为非常小的多边形。像这样：
将多边形与覆盖整个区域的大多边形进行差分会产生用于三角剖分的多边形。可能还要等一下。结果如下所示：
最后，剔除所有过大的多边形：
#!/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 = np.loadtxt('test.asc')
xycoors = dat[:,0:2]
xcoors = xycoors[:,0] #Convenience alias
ycoors = xycoors[:,1] #Convenience alias
npts = len(dat[:,0]) #Number of points
dist = scipy.spatial.distance.euclidean
def GetGraph(xycoors, alpha=0.0035):
kdt = scipy.spatial.KDTree(xycoors) #Build kd-tree for quick neighbor lookups
G = nx.Graph()
npts = np.max(xycoors.shape)
for x in range(npts):
G.add_node(x)
dist, idx = kdt.query(xycoors[x,:], k=10) #Get distances to neighbours, excluding the cenral point
dist = dist[1:] #Drop central point
idx = idx[1:] #Drop central point
pq = np.exp(-alpha*dist) #Exponential weighting of nearby points
pq = pq/np.sum(pq) #Convert to a PDF
choices = np.random.choice(idx, p=pq, size=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]['weight'] += 1 #Strengthen connection
else:
G.add_edge(x, c, weight=1) #New neighbor; build connection
return G
def PruneGraph(G,cutoff):
newg = G.copy()
bad_edges = set()
for x in newg:
for k,v in newg[x].items():
if v['weight']<cutoff:
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=6):
xcoors = xycoors[:,0]
ycoors = xycoors[:,1]
G = 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]), 'k-', lw=1)
plt.show()
def GetPolys(xycoors,G):
#Get lines connecting all points in the graph
xcoors = xycoors[:,0]
ycoors = xycoors[:,1]
lines = []
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 = np.min(xycoors[:,0])
xmax = np.max(xycoors[:,0])
ymin = np.min(xycoors[:,1])
ymax = np.max(xycoors[:,1])
mls = shapely.geometry.MultiLineString(lines) #Bundle the lines
mlsb = mls.buffer(2) #Turn lines into narrow polygons
bbox = shapely.geometry.box(xmin,ymin,xmax,ymax) #Generate background polygon
polys = bbox.difference(mlsb) #Subtract to generate polygons
return polys
def PlotPolys(polys,area_cutoff):
fig, ax = plt.subplots(figsize=(8, 8))
for polygon in polys:
if polygon.area<area_cutoff:
mpl_poly = matplotlib.patches.Polygon(np.array(polygon.exterior), alpha=0.4, facecolor=np.random.rand(3,1))
ax.add_patch(mpl_poly)
ax.autoscale()
fig.show()
#Functional stuff starts here
G = GetGraph(xycoors, alpha=0.0035)
#Choose a value that rips off an appropriate amount of the left side of this histogram
weights = sorted([v['weight'] for x in G for k,v in G[x].items()])
plt.hist(weights, bins=20);plt.show()
PlotGraph(xycoors,G,cutoff=6) #Plot the graph to ensure our cut-offs were okay. May take a while
prunedg = PruneGraph(G,cutoff=6) #Prune the graph
polys = GetPolys(xycoors,prunedg) #Get polygons from graph
areas = sorted(p.area for p in polys)
plt.plot(areas)
plt.hist(areas,bins=20);plt.show()
area_cutoff = 150000
PlotPolys(polys,area_cutoff=area_cutoff)
good_polys = ([p for p in polys if p.area<area_cutoff])
total_area = sum([p.area for p in good_polys])
我在显示二维图像的 example ASCII file 中有一组点。
alpha 值，然后估计面积。
from shapely.ops import cascaded_union, polygonize
import shapely.geometry as geometry
from scipy.spatial import Delaunay
import numpy as np
import pylab as pl
from descartes import PolygonPatch
from matplotlib.collections import LineCollection
def plot_polygon(polygon):
fig = pl.figure(figsize=(10,10))
ax = fig.add_subplot(111)
margin = .3
x_min, y_min, x_max, y_max = polygon.bounds
ax.set_xlim([x_min-margin, x_max+margin])
ax.set_ylim([y_min-margin, y_max+margin])
patch = PolygonPatch(polygon, fc='#999999',
ec='#000000', fill=True,
zorder=-1)
ax.add_patch(patch)
return fig
def alpha_shape(points, alpha):
if len(points) < 4:
# When you have a triangle, there is no sense
# in computing an alpha shape.
return geometry.MultiPoint(list(points)).convex_hull
def add_edge(edges, edge_points, coords, i, j):
"""
Add a line between the i-th and j-th points,
if not in the list already
"""
if (i, j) in edges or (j, i) in edges:
# already added
return
edges.add( (i, j) )
edge_points.append(coords[ [i, j] ])
coords = np.array([point.coords[0]
for point in points])
tri = Delaunay(coords)
edges = set()
edge_points = []
# loop over triangles:
# ia, ib, ic = indices of corner points of the
# triangle
for ia, ib, ic in tri.vertices:
pa = coords[ia]
pb = coords[ib]
pc = coords[ic]
# Lengths of sides of triangle
a = np.sqrt((pa[0]-pb[0])**2 + (pa[1]-pb[1])**2)
b = np.sqrt((pb[0]-pc[0])**2 + (pb[1]-pc[1])**2)
c = np.sqrt((pc[0]-pa[0])**2 + (pc[1]-pa[1])**2)
# Semiperimeter of triangle
s = (a + b + c)/2.0
# Area of triangle by Heron's formula
area = np.sqrt(s*(s-a)*(s-b)*(s-c))
circum_r = a*b*c/(4.0*area)
# Here's the radius filter.
#print circum_r
if circum_r < 1.0/alpha:
add_edge(edges, edge_points, coords, ia, ib)
add_edge(edges, edge_points, coords, ib, ic)
add_edge(edges, edge_points, coords, ic, ia)
m = geometry.MultiLineString(edge_points)
triangles = list(polygonize(m))
return cascaded_union(triangles), edge_points
points=[]
with open("test.asc") as f:
for line in f:
coords=map(float,line.split(" "))
points.append(geometry.shape(geometry.Point(coords[0],coords[1])))
print geometry.Point(coords[0],coords[1])
x = [p.x for p in points]
y = [p.y for p in points]
pl.figure(figsize=(10,10))
point_collection = geometry.MultiPoint(list(points))
point_collection.envelope
convex_hull_polygon = point_collection.convex_hull
_ = plot_polygon(convex_hull_polygon)
_ = pl.plot(x,y,'o', color='#f16824')
concave_hull, edge_points = alpha_shape(points, alpha=0.001)
lines = LineCollection(edge_points)
_ = plot_polygon(concave_hull)
_ = pl.plot(x,y,'o', color='#f16824')
我得到了这个结果，但我希望这种方法可以检测到中间的洞。
更新
这是我的真实数据的样子：
我的问题是，估计上述形状面积的最佳方法是什么？我无法弄清楚这段代码无法正常工作出了什么问题？！！任何帮助将不胜感激。
编辑：
我注意到您有自己的代码来计算 alpha 形状， Delaunay三角形的面积就在那里，所以计算形状的面积就更容易了...
如果要将三角形添加到 alpha 形多边形中，只需添加三角形的面积。
如果你想检测孔...添加一个辅助阈值以避免添加面积大于阈值的三角形。对于此示例，值 max_area = 99999 将删除孔。
唯一的问题是您创建图形输出的方式，因为您看不到孔。
def alpha_shape(points, alpha, max_area):
if len(points) < 4:
# When you have a triangle, there is no sense
# in computing an alpha shape.
return geometry.MultiPoint(list(points)).convex_hull , 0
def add_edge(edges, edge_points, coords, i, j):
"""
Add a line between the i-th and j-th points,
if not in the list already
"""
if (i, j) in edges or (j, i) in edges:
# already added
return
edges.add( (i, j) )
edge_points.append(coords[ [i, j] ])
coords = np.array([point.coords[0]
for point in points])
tri = Delaunay(coords)
total_area = 0
edges = set()
edge_points = []
# loop over triangles:
# ia, ib, ic = indices of corner points of the
# triangle
for ia, ib, ic in tri.vertices:
pa = coords[ia]
pb = coords[ib]
pc = coords[ic]
# Lengths of sides of triangle
a = np.sqrt((pa[0]-pb[0])**2 + (pa[1]-pb[1])**2)
b = np.sqrt((pb[0]-pc[0])**2 + (pb[1]-pc[1])**2)
c = np.sqrt((pc[0]-pa[0])**2 + (pc[1]-pa[1])**2)
# Semiperimeter of triangle
s = (a + b + c)/2.0
# Area of triangle by Heron's formula
area = np.sqrt(s*(s-a)*(s-b)*(s-c))
circum_r = a*b*c/(4.0*area)
# Here's the radius filter.
# print("radius", circum_r)
if circum_r < 1.0/alpha and area < max_area:
add_edge(edges, edge_points, coords, ia, ib)
add_edge(edges, edge_points, coords, ib, ic)
add_edge(edges, edge_points, coords, ic, ia)
total_area += area
m = geometry.MultiLineString(edge_points)
triangles = list(polygonize(m))
return cascaded_union(triangles), edge_points, total_area
旧答案：
要计算不规则简单多边形的面积，您可以使用 Shoelace formula 和边界的 CCW 坐标作为输入。
如果要检测云内部的孔，则必须删除外接半径大于次级阈值的 Delaunay 三角形。 理想情况是：计算 Delaunay 三角剖分并使用当前的 alpha 形状进行过滤。然后，计算每个三角形的外接半径，去掉外接半径远大于平均外接半径的三角形。
要计算带孔的不规则多边形的面积，请对每个孔边界使用鞋带公式。以CCW（正向）顺序输入外边界，得到面积。然后按CW（负）顺序输入每个孔的边界，得到面积的（负）值。
这里有一个想法：使用 k-means clustering.
您可以在 Python 中按如下方式完成此操作：
from sklearn.cluster import KMeans
import numpy as np
import matplotlib.pyplot as plt
dat = np.loadtxt('test.asc')
xycoors = dat[:,0:2]
fit = KMeans(n_clusters=2).fit(xycoors)
plt.scatter(dat[:,0],dat[:,1], c=fit.labels_)
plt.axes().set_aspect('equal', 'datalim')
plt.gray()
plt.show()
使用您的数据，得出以下结果：
现在，您可以获取顶部簇和底部簇的凸包，并分别计算每个簇的面积。然后将这些区域相加成为它们合并区域的估计量，但巧妙地避免了中间的洞。
要微调您的结果，您可以使用聚类的数量和算法的不同启动次数（该算法是随机的，通常 运行 不止一次）。
你问，比如，两个簇是否总是在中间留下一个洞。我使用以下代码对此进行了试验。我生成一个均匀分布的点，然后切出一个随机大小和方向的椭圆来模拟一个洞。
#!/usr/bin/env python3
import sklearn
import sklearn.cluster
import numpy as np
import matplotlib.pyplot as plt
PWIDTH = 6
PHEIGHT = 6
def GetPoints(num):
return np.random.rand(num,2)*300-150 #Centered about zero
def MakeHole(pts): #Chop out a randomly orientated and sized ellipse
a = np.random.uniform(10,150) #Semi-major axis
b = np.random.uniform(10,150) #Semi-minor axis
h = np.random.uniform(-150,150) #X-center
k = np.random.uniform(-150,150) #Y-center
A = np.random.uniform(0,2*np.pi) #Angle of rotation
surviving_points = []
for pt in range(pts.shape[0]):
x = pts[pt,0]
y = pts[pt,1]
if ((x-h)*np.cos(A)+(y-k)*np.sin(A))**2/a/a+((x-h)*np.sin(A)-(y-k)*np.cos(A))**2/b/b>1:
surviving_points.append(pt)
return pts[surviving_points,:]
def ShowManyClusters(pts,fitter,clusters,title):
colors = np.array([x for x in 'bgrcmykbgrcmykbgrcmykbgrcmyk'])
fig,axs = plt.subplots(PWIDTH,PHEIGHT)
axs = axs.ravel()
for i in range(PWIDTH*PHEIGHT):
lbls = fitter(pts[i],clusters)
axs[i].scatter(pts[i][:,0],pts[i][:,1], c=colors[lbls])
axs[i].get_xaxis().set_ticks([])
axs[i].get_yaxis().set_ticks([])
plt.suptitle(title)
#plt.show()
plt.savefig('/z/'+title+'.png')
fitters = {
'SpectralClustering': lambda x,clusters: sklearn.cluster.SpectralClustering(n_clusters=clusters,affinity='nearest_neighbors').fit(x).labels_,
'KMeans': lambda x,clusters: sklearn.cluster.KMeans(n_clusters=clusters).fit(x).labels_,
'AffinityPropagation': lambda x,clusters: sklearn.cluster.AffinityPropagation().fit(x).labels_,
}
np.random.seed(1)
pts = []
for i in range(PWIDTH*PHEIGHT):
temp = GetPoints(300)
temp = MakeHole(temp)
pts.append(temp)
for name,fitter in fitters.items():
for clusters in [2,3]:
np.random.seed(1)
ShowManyClusters(pts,fitter,clusters,"{0}: {1} clusters".format(name,clusters))
考虑 K-Means 的结果：
至少在我看来，当 "hole" 将数据分成两个单独的 blob 时，使用两个集群似乎效果最差。 （在这种情况下，当椭圆的方向与包含样本点的矩形区域的两条边重叠时会发生这种情况。）使用三个聚类可以解决大部分困难。
您还会注意到 K-means 在第 1 列第 3 行以及第 3 列第 4 行上产生了一些违反直觉的结果。查看 sklearn 的聚类方法 here 显示以下比较图像：
由此看来，似乎 SpectralClustering 产生的结果与我们想要的一致。在上面的相同数据上尝试此操作可以解决提到的问题（请参阅第 1 列，第 3 行和第 3 列，第 4 行）。
上述内容表明，具有三个聚类的谱聚类应该足以应对此类大多数情况。
虽然您似乎打算做一个凹形，但这里有一条替代路线，速度非常快，我认为会给您非常稳定的读数：
创建一个以参数 (int radiusOfInfluence) 作为参数的函数。在函数内 运行 一个以那个为半径的体素过滤器。然后只需将该圆的面积 (pi*AOI^2) 乘以云中剩余点的数量。这应该会给你一个相对稳健的面积估计，并且对孔洞和奇怪的边缘非常有弹性。
需要考虑的一些事项：
-这会给你一个积极的区域超调，因为边缘超出了一个半径。对此进行调整的修改可能是 运行 统计离群值去除过滤器（在反向模式下）以获取统计边缘点。然后可以假设每个边缘点的大约一半位于形状之外，在乘以面积之前从总计数中减去找到的点数的一半。
-影响半径在很大程度上决定了此函数的孔洞检测，因为较大的影响半径将允许单个点覆盖更大的区域，而且通过调整统计离群值过滤器上的标准截止值，您可以更积极地检测内部孔洞和以这种方式调整您的区域。
它确实回避了您所追求的问题，因为这更像是假设一组合理分布的样本的射击精度/射击分组类型评估。您的方法有点假设您的外边缘点是可能的绝对限制（根据情况，这可能是一个合理的假设）
编辑------------------------
我没有时间写出示例代码，但我可以进一步解释以帮助理解。
其核心是voxel filter。非常简单，它在 x、y 坐标中设置一个种子点，然后在整个 space 上创建一个网格，该网格在用户指定的过滤器半径的两个轴上都有单位（网格间距）。在每个网格框中，它会将所有点平均为一个点。这对于这个概念非常重要，因为它几乎完全消除了重叠问题。
第二部分（逆stat outlier removal）只是一点小聪明来加强你的边缘配合。基本上，stat outlier 是通过查看每个点到其 (k) 个最近邻点的距离来消除噪声的。在为每个点生成到 k 个最近邻居的平均距离后，它会建立一个直方图和一个用户定义的参数作为保留或删除点的二进制阈值。当反转并设置为合理的截止值（~0.75 std 应该有效）时，它会删除对象主体中的所有点（即只留下边缘点）。这一点很重要的原因是，从技术上讲，这些点超出了对象边界 1 个半径。虽然有些边缘角度是锐角，有些边缘角度是钝角（即大于或小于半个溢出圆），但每个点去掉 1/2 的圆面积应该在整个对象上给你一个很好的边缘拟合改进.
请记住，归根结底，这只会给您一个数字。就压力测试而言，我建议创建已知区域的人为点云和/或创建一个图形输出来显示您放置圆圈和半圆圈的位置（如果您喜欢的话，可以朝向对象的内部）。
您想要改进此方法的旋钮是： 体素过滤器半径，每个点的影响区域（实际上可以与体素过滤器半径分开控制，尽管它们应该保持非常接近彼此），标准截止。
希望这有助于澄清，干杯！
好的，这就是我的想法。 Delaunay 三角剖分将生成不分青红皂白的大三角形。这也会有问题，因为只会生成三角形。
因此，我们将生成您可能称之为 "fuzzy Delaunay triangulation" 的内容。我们将把所有点放入 kd-tree 中，并且对于每个点 p，查看它的 k 最近邻点。 kd-tree 让这个变得很快。
对于每个 k 个邻居，找出到焦点 p 的距离。使用此距离生成权重。我们希望附近的点比更远的点更受青睐，因此指数函数 exp(-alpha*dist) 在这里是合适的。使用加权距离构建描述绘制每个点的概率的概率密度函数。
现在，从那个分布中抽取很多次。附近的点将经常被选择，而较远的点将被较少地选择。对于绘制的点，记下它为焦点绘制了多少次。结果是一个加权图，其中图中的每条边都连接附近的点，并根据选择对的频率进行加权。
现在，剔除图中权重太小的所有边。这些是可能没有连接的点。结果如下所示：
现在，让我们把所有剩余的边都扔进shapely。然后我们可以通过缓冲将边转换为非常小的多边形。像这样：
将多边形与覆盖整个区域的大多边形进行差分会产生用于三角剖分的多边形。可能还要等一下。结果如下所示：
最后，剔除所有过大的多边形：
#!/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 = np.loadtxt('test.asc')
xycoors = dat[:,0:2]
xcoors = xycoors[:,0] #Convenience alias
ycoors = xycoors[:,1] #Convenience alias
npts = len(dat[:,0]) #Number of points
dist = scipy.spatial.distance.euclidean
def GetGraph(xycoors, alpha=0.0035):
kdt = scipy.spatial.KDTree(xycoors) #Build kd-tree for quick neighbor lookups
G = nx.Graph()
npts = np.max(xycoors.shape)
for x in range(npts):
G.add_node(x)
dist, idx = kdt.query(xycoors[x,:], k=10) #Get distances to neighbours, excluding the cenral point
dist = dist[1:] #Drop central point
idx = idx[1:] #Drop central point
pq = np.exp(-alpha*dist) #Exponential weighting of nearby points
pq = pq/np.sum(pq) #Convert to a PDF
choices = np.random.choice(idx, p=pq, size=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]['weight'] += 1 #Strengthen connection
else:
G.add_edge(x, c, weight=1) #New neighbor; build connection
return G
def PruneGraph(G,cutoff):
newg = G.copy()
bad_edges = set()
for x in newg:
for k,v in newg[x].items():
if v['weight']<cutoff:
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=6):
xcoors = xycoors[:,0]
ycoors = xycoors[:,1]
G = 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]), 'k-', lw=1)
plt.show()
def GetPolys(xycoors,G):
#Get lines connecting all points in the graph
xcoors = xycoors[:,0]
ycoors = xycoors[:,1]
lines = []
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 = np.min(xycoors[:,0])
xmax = np.max(xycoors[:,0])
ymin = np.min(xycoors[:,1])
ymax = np.max(xycoors[:,1])
mls = shapely.geometry.MultiLineString(lines) #Bundle the lines
mlsb = mls.buffer(2) #Turn lines into narrow polygons
bbox = shapely.geometry.box(xmin,ymin,xmax,ymax) #Generate background polygon
polys = bbox.difference(mlsb) #Subtract to generate polygons
return polys
def PlotPolys(polys,area_cutoff):
fig, ax = plt.subplots(figsize=(8, 8))
for polygon in polys:
if polygon.area<area_cutoff:
mpl_poly = matplotlib.patches.Polygon(np.array(polygon.exterior), alpha=0.4, facecolor=np.random.rand(3,1))
ax.add_patch(mpl_poly)
ax.autoscale()
fig.show()
#Functional stuff starts here
G = GetGraph(xycoors, alpha=0.0035)
#Choose a value that rips off an appropriate amount of the left side of this histogram
weights = sorted([v['weight'] for x in G for k,v in G[x].items()])
plt.hist(weights, bins=20);plt.show()
PlotGraph(xycoors,G,cutoff=6) #Plot the graph to ensure our cut-offs were okay. May take a while
prunedg = PruneGraph(G,cutoff=6) #Prune the graph
polys = GetPolys(xycoors,prunedg) #Get polygons from graph
areas = sorted(p.area for p in polys)
plt.plot(areas)
plt.hist(areas,bins=20);plt.show()
area_cutoff = 150000
PlotPolys(polys,area_cutoff=area_cutoff)
good_polys = ([p for p in polys if p.area<area_cutoff])
total_area = sum([p.area for p in good_polys])