- c - 在位数组中找到第一个零
- linux - Unix 显示有关匹配两种模式之一的文件的信息
- 正则表达式替换多个文件
- linux - 隐藏来自 xtrace 的命令
我在 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')
我的问题是,估计上述形状面积的最佳方法是什么?我无法弄清楚这段代码无法正常工作出了什么问题?!!任何帮助将不胜感激。
最佳答案
好的,这就是想法。 Delaunay 三角剖分将生成不分青红皂白的大三角形。这也会有问题,因为只会生成三角形。
因此,我们将生成您所谓的“模糊 Delaunay 三角剖分”。我们会将所有点放入 kd 树中,并针对每个点 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])
关于python - 估计由一组点(Alpha 形状??)生成的图像区域,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/41268547/
我正在尝试使用以下 keytool 命令为我的应用程序生成 keystore : keytool -genkey -alias tomcat -keystore tomcat.keystore -ke
编辑:在西里尔正确解决问题后,我注意到只需将生成轴的函数放在用于生成标签的函数下面就可以解决问题。 我几乎读完了 O'Reilly 书中关于 D3.js 的教程,并在倒数第二页上制作了散点图,但是当添
虽然使用 GraphiQL 效果很好,但我的老板要求我实现一个用户界面,用户可以在其中通过 UI 元素(例如复选框、映射关系)检查呈现给他们的元素并获取数据,这样做将为该人生成 graphql 输入,
我尝试在 Netbean 6.8 中使用 ws-import 生成 Java 类。我想重新生成 jax-ws,因为在 ebay.api.paypalapi 包中发现了一个错误(我认为该错误是由于 Pa
我有一个 perl 脚本,它获取系统日期并将该日期写入文件名。 系统日期被分配给 TRH1 变量,然后它被设置为一个文件名。 $TRH1 =`date + %Y%m%d%H%M`; print "TR
我是 Haskell 的新手,需要帮助。我正在尝试构建一种必须具有某种唯一性的新数据类型,因此我决定使用 UUID 作为唯一标识符: data MyType = MyType { uuid ::
我制作了一个脚本,它可以根据 Mysql 数据库中的一些表生成 XML。 该脚本在 PHP 中运行。 public function getRawMaterials($apiKey, $format
所以这是我的项目中的一个问题。 In this task, we will use OpenSSL to generate digital signatures. Please prepare a f
我在 SAS LIFEREG 中有一个加速故障时间模型,我想绘制它。因为 SAS 在绘图方面非常糟糕,我想实际重新生成 R 中曲线的数据并将它们绘制在那里。 SAS 提出了一个尺度(在指数分布固定为
我正在为 Django 后端制作一个样板,并且我需要能够使它到达下一个下载它的人显然无法访问我的 secret key 的地方,或者拥有不同的 key 。我一直在研究一些选项,并在这个过程中进行了实验
我正在创建一个生成采购订单的应用程序。我可以根据用户输入的详细信息创建文本文件。我想生成一个看起来比普通文本文件好得多的 Excel。有没有可以在我的应用程序中使用的开源库? 最佳答案 目前还没有任何
我正在尝试使用 ScalaCheck 为 BST 创建一个 Gen,但是当我调用 .sample 方法时,它给了我 java.lang.NullPointerException。我哪里错了? seal
已关闭。此问题需要 debugging details 。目前不接受答案。 编辑问题以包含 desired behavior, a specific problem or error, and the
我尝试编写一些代码,例如(在verilog中): parameter N = 128; if (encoder_in[0] == 1) begin 23 binary_out = 1;
我正忙于在 Grails 项目中进行从 MySQL 到 Postgres 的相当复杂的数据迁移。 我正在使用 GORM 在 PostGres 中生成模式,然后执行 MySQL -> mysqldump
如何使用纯 XSLT 生成 UUID?基本上是寻找一种使用 XSLT 创建独特序列的方法。该序列可以是任意长度。 我正在使用 XSLT 2.0。 最佳答案 这是一个good example 。基本上,
我尝试安装.app文件,但是当我安装并单击“同步”(在iTunes中)时,我开始在设备上开始安装,然后停止,这是一个问题,我不知道在哪里,但我看到了我无法解决的奇怪的事情: 最佳答案 似乎您没有在Xc
自从我生成 JavaDocs 以来已经有一段时间了,我确信这些选项在过去 10 年左右的时间里已经得到了改进。 我能否得到一些有关生成器的建议,该生成器将输出类似于 .Net 文档结构的 JavaDo
我想学习如何生成 PDF,我不想使用任何第三方工具,我想自己用代码创建它。到目前为止,我所看到的唯一示例是我通过在第 3 方 dll 上打开反射器查看的代码,以查看发生了什么。不幸的是,到目前为止我看
我正在从 Epplus 库生成 excel 条形图。 这是我成功生成的。 我的 table 是这样的 Mumbai Delhi Financial D
我是一名优秀的程序员,十分优秀!