问题描述

我在 我想估计这些点填充的总面积.该平面内的某些地方没有被任何点填充，因为这些区域已被屏蔽.我猜想估计面积可能是应用凹壳 或alpha 形状.我尝试了

更新
这就是我的真实数据的样子:

我的问题是估计上述形状面积的最佳方法是什么?我无法弄清楚这段代码不能正常工作出了什么问题?！！任何帮助将不胜感激.

好的，这就是想法.Delaunay 三角剖分将生成不分青红皂白的大三角形.这也会有问题，因为只会生成三角形.

因此，我们将生成您可能称之为模糊 Delaunay 三角剖分"的内容.我们将把所有点放入一个 kd-tree 中，并且对于每个点 p，查看它的 k 最近邻居.kd-tree 使这个速度更快.

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

现在，从该分布中提取大量数据.将经常选择附近的点，而不太经常选择较远的点.对于绘制的点，请记下为焦点绘制了多少次.结果是一个加权图，图中的每条边都连接附近的点，并根据选择对的频率进行加权.

现在，从图中剔除权重太小的所有边.这些是可能没有连接的点.结果如下所示:

现在，让我们将所有剩余的边放入

用覆盖整个区域的大多边形来区分多边形将产生用于三角剖分的多边形.可能还要等一下.结果如下所示:

最后，剔除所有过大的多边形:

#!/usr/bin/env python将 numpy 导入为 np将 matplotlib.pyplot 导入为 plt随机导入导入 scipy导入 scipy.spatial将 networkx 导入为 nx匀称地进口导入 shapely.geometry导入 matplotlibdat = np.loadtxt('test.asc')xycoors = dat[:,0:2]xcoors = xycoors[:,0] #方便别名ycoors = xycoors[:,1] #方便别名npts = len(dat[:,0]) #点数dist = scipy.spatial.distance.euclideandef GetGraph(xycoors, alpha=0.0035):kdt = scipy.spatial.KDTree(xycoors) #构建kd-tree用于快速邻居查找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) #转换为PDFChoices = np.random.choice(idx, p=pq, size=50) #根据PDF选择邻居for c in selection: #将邻居插入图中if G.has_edge(x, c): #已经见过的邻居G[x][c]['weight'] += 1 #加强连接别的:G.add_edge(x, c, weight=1) #新邻居；建立连接返回 Gdef PruneGraph(G,cutoff):newg = G.copy()bad_edges = 设置()对于 newg 中的 x:对于 newg[x].items() 中的 k,v:如果 v['weight']
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])