問題描述

我在 我想估計這些點填充的總面積.該平面內的某些地方沒有被任何點填充，因為這些區域已被屏蔽.我猜想估計面積可能是應用凹殼 或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])


                        
這篇關于估計由一組點生成的圖像區域(Alpha 形狀??)的文章就介紹到這了，希望我們推薦的答案對大家有所幫助，也希望大家多多支持html5模板網！
                
【網站聲明】本站部分內容來源于互聯網,旨在幫助大家更快的解決問題，如果有圖片或者內容侵犯了您的權益，請聯系我們刪除處理，感謝您的支持！