问题描述
我在
我想估算这些点填充的总面积。该平面内有一些地方未被任何点填充，因为这些区域已被掩盖。我认为估算该区域可能是实用的，方法是使用凹面船体或 alpha形状。
我尝试
更新
实际数据如下：
我的问题是估算面积的最佳方法是什么前面提到的形状？我无法弄清楚该代码无法正常工作是怎么了？！任何帮助将不胜感激。
好的，这是个主意。 Delaunay三角剖分将生成不加选择的大三角形。这也将是有问题的，因为将仅生成三角形。
因此，我们将生成所谓的模糊Delaunay三角剖分。我们将所有点放入kd树中，对于每个点 p ，请查看其最近的 k 邻居。 kd-tree可以使速度更快。
对于每个 k 个邻居，找到到焦点的距离 p 。使用该距离生成权重。我们希望附近的点比更远的点更受青睐，因此此处使用指数函数 exp（-alpha * dist）是合适的。使用加权距离来构建描述每个点的概率的概率密度函数。
现在，从该分布中进行大量绘制。经常选择附近的点，而较远的点则较少。对于绘制点，记下为焦点绘制了多少次。结果是一个加权图，其中图中的每个边都连接附近的点，并按选择对的频率进行加权。
现在，从图中删除其所有边重量太小。这些是可能没有连接的点。结果如下所示：
现在，让我们将所有剩余的边缘扔到
用大的多边形覆盖物来区分多边形整个区域将产生用于三角剖分的多边形。可能还要等一下。结果看起来像这样：
最后，剔除所有太大的多边形：
＃！/ usr / bin / env python
将numpy导入为np
将matplotlib.pyplot导入为plt
import random
import scipy
import scipy.spatial
import networkx as nx
importly
import shape.geometry
import matplotlib
dat = np.loadtxt（'test.asc'）
xycoors = dat [：，0：2]
xcoors = xycoors [：，0]＃便捷别名
ycoors = xycoors [：，1]＃便利别名
npts = len（dat [：，0]）＃点数
dist = s cipy.spatial.distance.euclidean
def GetGraph（xycoors，alpha = 0.0035）：
kdt = scipy.spatial.KDTree（xycoors）＃构建kd树以快速查找邻居
G = nx.Graph（）
npts = np.max（xycoors.shape）
对于x范围（npts）：
G.add_node（x）
dist ，idx = kdt.query（xycoors [x ,:]，k = 10）＃获取到邻居的距离，不包括中心点
dist = dist [1：]＃丢弃中心点
idx = idx [1：]＃删除中心点
pq = np.exp（-alpha * dist）＃附近点的指数加权
pq = pq / np.sum（pq）＃转换为PDF
choices = np.random.choice（idx，p = pq，size = 50）＃根据PDF
为选择c的邻居选择＃如果G.has_edge（x，则将邻居插入图
，c）：＃已经见过的邻居
G [x] [c] ['weight'] + = 1＃强度连接
其他：
G.add_edge（x，c，weight = 1）＃新邻居；建立连接
返回G
def PruneGraph（G，cutoff）：
newg = G.copy（）
bad_edges = set（）
for x在newg中：
for k，v在newg [x] .items（）中：
如果v ['weight'] bad_edges.add（（x，k））
在bad_edges中的b：
尝试：
newg.remove_edge（* b）
除了nx.exception.NetworkXError：
传递
返回newg
def PlotGraph（xycoors，G，cutoff = 6）：
xcoors = xycoors [：，0]
ycoors = xycoors [：，1]
G = PruneGraph（G，cutoff）
plt.plot（xcoors，ycoors， o）
对于范围x中的x（npts）：
对于G [x]中的k，v。 items（）：
plt.plot（（xcoors [x]，xcoors [k]），（ycoors [x]，ycoors [k]），'k-'，lw = 1）
plt .show（）
def GetPolys（xycoors，G）：b $ b $＃获取连接图中所有点的线
xcoors = xycoors [：，0]
ycoors = xycoors [：，1]
行= []
对于x在范围内（npts）：
对于k，v在G [x] .items（）： b $ b lines.append（ （（（xcoors [x]，ycoors [x]），（xcoors [k]，ycoors [k]）））
＃获取区域
的边界xmin = np.min（xycoors [：，0 ]）
xmax = np.max（xycoors [：，0]）
ymin = np.min（xycoors [：，1]）$ b $ b ymax = np.max（xycoors [：,, 1]）$ b $ b mls = shapely.geometry.MultiLineString（lines）＃将线束
mlsb = mls.buffer（2）＃将线变成狭窄的多边形
bbox = shapely.geometry.box （xmin，ymin，xmax，ymax）＃生成背景多边形
polys = bbox.difference（mlsb）＃减去以生成多边形
返回多边形
def PlotPolys（polys，area_cutoff ）：
图，ax = plt.subplots（figsize =（8，8））
用于polys中的多边形：
如果polygon.area mpl_poly = matplotlib patch.Polygon（np.array（polygon.exterior），alpha = 0.4，facecolor = np.random.rand（3,1））
ax.add_patch（mpl_poly）
ax.autoscale（）
fig.show（）
＃功能开始于此
G = GetGr aph（xycoors，alpha = 0.0035）
＃选择一个值，该值会剥去该直方图左侧的适当量
weights = sorted（[x的[v ['weight']]在G中表示k，v在G [x] .items（）]）
plt.hist（weights，bins = 20）; plt.show（）
PlotGraph（xycoors，G ，cutoff = 6）＃绘制图表以确保我们的截止点正常。可能需要一段时间
prunedg = PruneGraph（G，cutoff = 6）＃修剪图形
polys = GetPolys（xycoors，prunedg）＃从图形
中获取多边形=已排序（多边形中p的p。面积）
plt.plot（面积）
plt.hist（areas，bins = 20）; plt.show（）
area_cutoff = 150000
PlotPolys（polys，area_cutoff = area_cutoff）
good_polys =（[如果p.area ＜area_cutoff，则为p中的p中的p的p））
total_area = sum（[good_polys中p为p的区域中的p] ）
I have a set of points in an example ASCII file showing a 2D image.I would like to estimate the total area that these points are filling. There are some places inside this plane that are not filled by any point because these regions have been masked out. What I guess might be practical for estimating the area would be applying a concave hull or alpha shapes.I tried this approach to find an appropriate alpha value, and consequently estimate the area.
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')
I get this result but I would like that this method could detect the hole in the middle.
Update
This is how my real data looks like:
My question is what is the best way to estimate an area of the aforementioned shape? I can not figure out what has gone wrong that this code doesn't work properly?!! Any help will be appreciated.
Okay, here's the idea. A Delaunay triangulation is going to generate triangles which are indiscriminately large. It's also going to be problematic because only triangles will be generated.
Therefore, we'll generate what you might call a "fuzzy Delaunay triangulation". We'll put all the points into a kd-tree and, for each point p, look at its k nearest neighbors. The kd-tree makes this fast.
For each of those k neighbors, find the distance to the focal point p. Use this distance to generate a weighting. We want nearby points to be favored over more distant points, so an exponential function exp(-alpha*dist) is appropriate here. Use the weighted distances to build a probability density function describing the probability of drawing each point.
Now, draw from that distribution a large number of times. Nearby points will be chosen often while farther away points will be chosen less often. For point drawn, make a note of how many times it was drawn for the focal point. The result is a weighted graph where each edge in the graph connects nearby points and is weighted by how often the pairs were chosen.
Now, cull all edges from the graph whose weights are too small. These are the points which are probably not connected. The result looks like this:
Now, let's throw all of the remaining edges into shapely. We can then convert the edges into very small polygons by buffering them. Like so:
Differencing the polygons with a large polygon covering the entire region will yield polygons for the triangulation. THIS MAY TAKE A WHILE. The result looks like this:
Finally, cull off all of the polygons which are too large:
#!/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])
这篇关于估计由一组点（“ Alpha形状” ??）生成的图像的面积的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！