pythonscipycomputer-visionshapelyconcave-hull

Estimating an area of an image generated by a set of points (Alpha shapes??)


I have a set of points in an example ASCII file showing a 2D image. enter image description here 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. enter image description here

Update
This is how my real data looks like: enter image description here

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.


Solution

  • 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:

    A fuzzy Delaunay triangulation

    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:

    Buffered polygons

    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:

    Fuzzy Delaunay triangulation with coloured polygons

    Finally, cull off all of the polygons which are too large:

    Fuzzy Delaunay triangulation with large polygons removed

    #!/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])