alphahull.py

import fiona
import shapely.geometry as geometry
from shapely.ops import cascaded_union, polygonize
from scipy.spatial import Delaunay
import numpy as np
import math
        

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):
        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])
    print 'starting delauney'
    tri = Delaunay(coords)
    print 'finished delauney'
    edges = set()
    edge_points = []
    # loop over triangles:
    # ia, ib, ic = indices of corner points of the
    # triangle
    print 'iterating'
    for ia, ib, ic in tri.vertices:
        pa = coords[ia]
        pb = coords[ib]
        pc = coords[ic]
        # Lengths of sides of triangle
        a = math.sqrt((pa[0]-pb[0])**2 + (pa[1]-pb[1])**2)
        b = math.sqrt((pb[0]-pc[0])**2 + (pb[1]-pc[1])**2)
        c = math.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 = math.sqrt(s*(s-a)*(s-b)*(s-c))
        circum_r = a*b*c/(4.0*area)
        #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)
    print 'polygonizing'
    m = geometry.MultiLineString(edge_points)
    triangles = list(polygonize(m))
    return cascaded_union(triangles), edge_points

## INPUT
areaa = str(raw_input('area: '))
input_file = 'points/{0}.csv'.format(areaa)

points = []
with open(input_file) as kf:
    for line in kf.readlines()[1:]:
        x,y = float(line.split(';')[0]),float(line.split(';')[1])
        points.append(geometry.Point(x,y))
print 'points read'
    
concave_hull, edge_points = alpha_shape(points,
                                        alpha=0.008)
print 'func finished'
                                        
filename = 'hulls/alphahullWKT_{0}.wkt'.format(areaa)
with open(filename,'w') as wkt:
    wkt.write('wkt\n')
    for item in concave_hull:
        wkt.write(str(item) + '\n')
print 'ready'