A universal tiler and polygon decompositor

tiler.py

# To use this, you will want to change the location of the saved data
# at least and possibly the functions called below (at the end of the file).

# You first need to install the Python Shapely library for manipulating polygons:

# http://trac.gispython.org/lab/wiki/Shapely

from math import pi, sin, cos, atan2
import shapely.geometry
import sys, copy

def position_opposite(large_order,small_order):
    sym_lookup = {}
    s0 = large_order/small_order
    #opposites only make sense if s0 is even
    if ((s0 % 2) == 0):
        for i in range(0,s0):
            if i % 2 == 0:
                sym_lookup[i] = i
            else:
                sym_lookup[i] = (i+(s0/2)-1) % s0
    else:
        sym_lookup = position_sequential(large_order,small_order)
    return sym_lookup

def position_6(large_order,small_order):
    #return [3, 2, 1, 0, 5, 4]
    return [4,5,0,1,2,3]
        
def position_sequential(large_order,small_order):
    sym_lookup = {}
    s0 = large_order/small_order
    for i in range(0,s0):
        sym_lookup[i] = i
    return sym_lookup

def position_alternating(large_order,small_order):
    sym_lookup = {}
    s0 = large_order/small_order
    for i in range(0,s0/2):
        sym_lookup[i] = 2*i
    for i in range(0,s0/2):
        sym_lookup[i+s0/2] = 2*i+1
    return sym_lookup

class Tiling(object):
    """Stores and manipulates a tile patch"""
    def __init__(self):
        self.pc = 0
        self.polygons = []
        self.ec = 0
        self.edges = []
        self.vc = 0
        self.vertices = []
        self.spc = 0
        self.shapely_polygons = []
        self.stack = []

    def add_to_shapely_polygons(self,pe):
        poly = []
        for e in pe:
            poly.append(e[0])
        self.shapely_polygons.append(shapely.geometry.Polygon(poly))
        self.spc += 1
        
    def get_vertex_list(self,v0,e0,p0):
        error = False
        e = self.edges[e0]
        v = self.vertices[v0]
        poly = [[v.x,v.y]]
        v1 = e.v1
        if v1 == v0:
            v1 = e.v0
        i = 0
        # cycle around the vertices until we return to v0
        while v1 != v0:
            if i > len(self.vertices):
                error = True
                break
            v = self.vertices[v1]
            poly.append([v.x,v.y])
            e0 = self.get_other_edge2(v1,e0,p0)
            v1 = self.get_other_vertex(v1,e0)
            i += 1
        
        if error:
            print "cannot compute vertex list"
            sys.exit()
        else:
            return poly
    
    def get_edge_list(self,v0,e0,p0):
        error = False
        e = self.edges[e0]
        poly = [e0]
        v1 = e.v1
        if v0 == v1:
            v1 = e.v0
        i = 0
        # cycle around the vertices until we return to v0
        while v1 != v0:
            e0 = self.get_other_edge2(v1,e0,p0)
            if e0 == -1 or i > len(self.vertices):
                error = True
                break
            v1 = self.get_other_vertex(v1,e0)
            poly.append(e0)
            i += 1
        if error:
            print "cannot compute edge list"
            sys.exit()
        else:
            return poly
        
    def convert_to_shapely_polygon(self,p,r=True):
        e0 = self.polygons[p].edges[0]
        e = self.edges[e0]
        v0 = e.v0
        poly = self.get_vertex_list(v0,e0,p)                
        if poly == False or len(poly) < 3:
            #sanity checks
            return False
        else:
            return shapely.geometry.Polygon(poly)

    def add_to_shapely_polygons2(self,poly):
        self.shapely_polygons.append(shapely.geometry.Polygon(poly))
        self.spc += 1

    def polygon_to_svg(self,p):
        if self.shapely_polygons[p] == False:
            return ''
        points = []
        for x,y in list(self.shapely_polygons[p].exterior.coords):
            points.append(str(20*x+500)+","+str(20*y+500))
        r = '<polygon style="stroke:#000000;stroke-width:1;fill:#F0F0F0" points="'
        r += " ".join(points)
        r += '"/>'+"\n"
        return r

    def labels_to_svg(self):
        r = ''
        for v in self.vertices:
            r += '<text x="'+str(20*v.x+500)
            r += '" y="'+str(20*v.y+500)
            r += '">('+str(v.num)+")</text>\n"
        for e in self.edges:
            vx = (0.5)*(self.vertices[e.v0].x+self.vertices[e.v1].x)
            vy = (0.5)*(self.vertices[e.v0].y+self.vertices[e.v1].y)
            r += '<text x="'+str(20*vx+500)
            r += '" y="'+str(20*vy+500)
            r += '">'+str(e.num)+"</text>\n"
        return r

    # TODO: change the data structures to make this more efficient to compute?
    def get_other_edge2(self,v0,e0,p0):
        for e in self.polygons[p0].edges:
            if e != e0 and (self.edges[e].v0 == v0 or self.edges[e].v1 == v0):
                return e
        return -1

    def add_to_vertex_sweep(self,v,e,p):
        vertex = self.vertices[v]
        if v in []:
            print "@@@************ adding to sweeps", v, e, p
            print "current sweep_e", vertex.sweep_e
            print "current sweep_p", vertex.sweep_p
        e_prev = -1
        p_prev = -1
        e2 = -1
        p2 = -1

        if len(vertex.sweep_e) <= 1:
            # the sweep direction is not yet established
            if p not in vertex.sweep_p:
                vertex.sweep_p.append(p)
            if e not in vertex.sweep_e:
                vertex.sweep_e.append(e)
        else:
            # need to add only in the sweep direction
            if vertex.sweep_e[-1] == e:
                if p not in vertex.sweep_p:
                    vertex.sweep_p.append(p)
            elif vertex.sweep_e[0] == e:
                if p not in vertex.sweep_p:
                    vertex.sweep_p = [p] + vertex.sweep_p
            else:
                # try the polygon method
                if vertex.sweep_p[-1] == p:
                    if e not in vertex.sweep_e:
                        vertex.sweep_e.append(e)
                elif vertex.sweep_p[0] == p:
                    if e not in vertex.sweep_e:
                        vertex.sweep_e = [e] + vertex.sweep_e
                else:
                    #try the other edge method
                    #TODO: should the get_other_edge code be here?
                    e2 = self.get_other_edge2(v,e,p)
                    if vertex.sweep_e[-1] == e2:
                        if p not in vertex.sweep_p:
                            vertex.sweep_p.append(p)
                    elif vertex.sweep_e[0] == e2:
                        if p not in vertex.sweep_p:
                            vertex.sweep_p = [p] + vertex.sweep_p
                    else:
                        # try the other polygon method
                        p2 = self.get_other_polygon(e,p)
                        if vertex.sweep_p[-1] == p2:
                            if p not in vertex.sweep_p:
                                vertex.sweep_p.append(p)
                            if e not in vertex.sweep_e:
                                vertex.sweep_e.append(e)
                        elif vertex.sweep_p[0] == p2:
                            if p not in vertex.sweep_p:
                                vertex.sweep_p = [p] + vertex.sweep_p
                            if e not in vertex.sweep_e:
                                vertex.sweep_e = [e] + vertex.sweep_e
                        else:
                            # try the previous polygon method
                            p_prev = self.get_previous_polygon(v,e,p)
                            if vertex.sweep_p[-1] == p_prev:
                                if p not in vertex.sweep_p:
                                    vertex.sweep_p.append(p)
                                if e not in vertex.sweep_e:
                                    vertex.sweep_e.append(e)
                            elif vertex.sweep_p[0] == p_prev:
                                if p not in vertex.sweep_p:
                                    vertex.sweep_p = [p] + vertex.sweep_p
                                if e not in vertex.sweep_e:
                                    vertex.sweep_e = [e] + vertex.sweep_e
                            else:
                                # try the previous edge method
                                e_prev = self.get_other_edge2(v,e,p_prev)
                                p_prev = self.get_previous_polygon(v,e,p)
                                if vertex.sweep_p[-1] == p_prev:
                                    if p not in vertex.sweep_p:
                                        vertex.sweep_p.append(p)
                                    if e not in vertex.sweep_e:
                                        vertex.sweep_e.append(e)
                                elif vertex.sweep_p[0] == p_prev:
                                    if p not in vertex.sweep_p:
                                        vertex.sweep_p = [p] + vertex.sweep_p
                                    if e not in vertex.sweep_e:
                                        vertex.sweep_e = [e] + vertex.sweep_e
                                else:
                                    pass
                                    
    def add_to_sweeps(self,edge_ids,e0,p0):
        for e in edge_ids:
            if e != e0 and e != -1:
                self.add_to_vertex_sweep(self.edges[e].v0,e,p0)
                self.add_to_vertex_sweep(self.edges[e].v1,e,p0)

    def get_other_polygon(self,e,p):
        plist = self.polygons_for([e])
        if len(plist) < 2:
            #print "warning: cannot find both polygons for", e, p
            return plist[0]
        if plist[0] == p:
            return plist[1]
        else:
            return plist[0]

    def get_previous_polygon(self,v,e0,p):
        elist = self.polygons[p].edges
        for e in elist:
            if e != e0 and (self.edges[e].v0 == v or self.edges[e].v1 == v):
                return self.get_other_polygon(e,p)
        return -1

    def get_previous_vertex(self,v0,e0):
        for p in self.polygons:
            elist = p.edges
            if e0 in elist:
                for e in elist:
                    if e != e0:
                        if self.edges[e].v0 == v0:
                            return self.edges[e].v1
                        elif self.edges[e].v1 == v0:
                            return self.edges[e].v0
        return -1

    def create_polygon(self,e0,p,poly_edges):

        #TODO: try to simplify this
        #In particular, some of the new_vertex calls are redundant
        edge_ids = self.get_edge_ids(poly_edges)
        self.add_to_sweeps(edge_ids,e0,p)
        vl = []
        for i in range(0,len(edge_ids)):
            e = edge_ids[i]
            if e == -1:
                # need a new edge
                vfi = self.new_vertex(poly_edges[i][0])
                vsi = self.new_vertex(poly_edges[i][1])
                e = self.new_edge(vfi,vsi)
            else:
                vfi = self.edges[e].v0
                vsi = self.edges[e].v1
            if e != e0:
                self.add_to_polygon(p,e)
            self.add_to_vertex_sweep(vfi,e,p)
            self.add_to_vertex_sweep(vsi,e,p)

        for pe in poly_edges:
            vi = self.new_vertex(pe[0])
            if vi not in vl:
                vl.append(vi)
        return vl

    def get_edge_ids(self,poly_edges):
        edge_ids = []
        for pe in poly_edges:
            e = self.get_edge_from_poly_edge(pe)
            edge_ids.append(e)
        return edge_ids

    def add_to_polygon(self,p,e):
        self.polygons[p].edges.append(e)
        
    def remove_from_polygon(self,p,e):
        v0 = self.edges[e].v0
        v1 = self.edges[e].v1
        
        self.polygons[p].edges.remove(e)

    #search vertex list looking for duplicates
    # TODO: find something more efficient?
    def new_vertex(self,v):
        for vi in range(0,self.vc):
            if not different_vertex(v,self.vertices[vi]):
                return vi
        v = Vertex(v[0],v[1])
        v.num = self.vc
        self.vertices.append(v)
        self.vc += 1
        return v.num

    def new_edge(self,v0,v1):
        for ei in range(0,self.ec):
            if not different_edge(v0,v1,self.edges[ei]):
                return ei
        e = Edge(v0,v1)
        e.num = self.ec
        self.edges.append(e)
        self.ec += 1
        return e.num

    def new_polygon(self,e0,ptype,make_internal=True):
        p = Polygon(e0,ptype,make_internal)
        p.num = self.pc
        self.polygons.append(p)
        self.pc += 1
        return p.num

    def get_edge_from_poly_edge(self,pe):
        for ei in range(0,self.ec):
            e = self.edges[ei]
            v0 = [self.vertices[e.v0].x, self.vertices[e.v0].y]
            v1 = [self.vertices[e.v1].x, self.vertices[e.v1].y]
            t0 = equal_vertex(v0,pe[0]) and equal_vertex(v1,pe[1])
            t1 = equal_vertex(v0,pe[1]) and equal_vertex(v1,pe[0])
            if t0 or t1:
                return ei
        return -1

    def start_polygon(self,n):
        v0 = self.new_vertex([-0.5,0.0])
        v1 = self.new_vertex([0.5,0.0])
        e = self.new_edge(v1,v0)
        ab,external,le = self.get_leading_edge2(v0,0)
        self.complete_at_vertex(v0,[n],external,1,[],[])

    def star_around(self,a):
        v0 = self.new_vertex([-0.5,0.0])
        v1 = self.new_vertex([0.5,0.0])
        e = self.new_edge(v0,v1)
        n = int(2*pi/a + 0.5)

        for i in range(0,n):
            p = compute_rhombus(v1,e,a)
            e = get_other_edge(v1,e,p)

    def edges_for(self,v):
        ev = []
        for e in self.edges:
            if (e.v0 == v) or (e.v1 == v):
                ev.append(e.num)
        return ev

    # elist is a list of edges
    def polygons_for(self,elist):
        pe = []
        for p in self.polygons:
            for e in elist:
                if e in p.edges:
                    pe.append(p.num)
                    break
        return pe

    def compute_circle_gap_at_vertex(self,v):
        at = 0
        vertex = self.vertices[v]
        for p in vertex.sweep_p:
            g = vertex.sweep_lu[p]
            if g > 0:
                ng = (pi / g) * (g-2)
                #ng = d[g]['p']
            else:
                ng = -g
            at += ng
        return 2*pi - at

    def get_poly_edge_from_edge(self,e):
        v0 = self.vertices[e[0]]
        v1 = self.vertices[e[1]]
        return [[v0.x,v0.y],[v1.x,v1.y]]

    def validate_sweep(self,v,species,rspecies):
        if len(species) == 0:
            # an empty species turns off validation
            return True
        # search for a solution to the sweep
        sweep = self.get_sweep(v)
        # if the sweep is complete, backstep to make sure that this is a legal vertex
        if abs(self.compute_circle_gap_at_vertex(v)) < 0.0001:
            sweep = sweep[1:]
        #print "in validate_sweep, checking sweep", convert_sweep(sweep)
        for cspecies in [species, rspecies]:
            for i in range(0,len(cspecies)):
                #print "against vertex type", convert_sweep(cspecies[i])
                pl = get_plist_from_sweep(sweep,cspecies[i])
                if len(pl) != 0:
                    return True
        return False
    
    def find_overlapping_polygons(self,pl):
        overlaps = []
        p0 = []
        for e in pl:
            p0.append(e[0])
        poly = shapely.geometry.Polygon(p0)
        i = 0
        for p in self.shapely_polygons:
            try:
                if poly.intersects(p):
                    pintersect = poly.intersection(p)
                    if (type(pintersect) is shapely.geometry.polygon.Polygon) and pintersect.area > 0.01:
                        #print "overlap detected", i, p,poly,pintersect,pintersect.area
                        overlaps.append(i)
            except:
                print "overlap intersect function error"
                return []
            i += 1
        return overlaps
    
    def punch(self,v,e,n):
        external = [self.edges[e].v0,self.edges[e].v1]
        vraw = [self.vertices[v].x,self.vertices[v].y]
        e0 = self.get_poly_edge_from_edge(external)
        parity = -1
        if n > 0:
            r = compute_polygon(vraw,e0,n,parity)
        else:
            r = compute_rhombus(vraw,e0,-n,parity)
        
        # only do something if the new polygon overlaps exactly one
        # other polygon
        overlaps = self.find_overlapping_polygons(r[0])
        if len(overlaps) == 0:
            print "Error: no overlap in punch"
            sys.exit()
        elif len(overlaps) > 1:
            print "Error: multiple overlaps in punch"
            sys.exit()
        else :
            ps = t.polygons_for([e])
            if len(ps) > 0:
                p0 = ps[0]
            else:
                print "Error in punch: cannot find existing polygon for edge ", e
                sys.exit()
            poly_edges = r[0]
            root_edge = r[1]
            a = r[2]
            e0 = self.get_edge_from_poly_edge(root_edge)
            p = self.new_polygon(e0,n)
            vl = self.create_polygon(e0,p,poly_edges)        
            # add the sweep angles to the vertices
            # TODO: this code is wrong for punches
            if n > 0:
                for vi in vl:
                    vertex = self.vertices[vi]
                    vertex.sweep_lu[p] = n
                    vertex.sweep_plu[p] = n
            else:
                # rhombuses have two angles, so the sweep code is a bit more complex
                if v == vl[0] or v == vl[2]:
                    ang = [a,pi-a,a,pi-a]
                else:
                    ang = [pi-a,a,pi-a,a]
    
                for i in range(0,len(vl)):
                    vertex = self.vertices[vl[i]]
                    vertex.sweep_lu[p] = -ang[i]
                    if ang[i] > (pi/2):
                        vertex.sweep_plu[p] = int(2*pi/(pi-ang[i])+0.5)
                    else:
                        vertex.sweep_plu[p] = -int((2*pi/ang[i])+0.5)
            self.add_to_shapely_polygons(poly_edges)
            esp = self.polygons[p].edges
            esp0 = copy.deepcopy(self.polygons[p0].edges)
            for e0 in esp:
                if e0 in esp0:
                    self.remove_from_polygon(p0,e0)
                else:
                    self.polygons[p0].edges.append(e0)
            
            self.shapely_polygons[p0] = self.convert_to_shapely_polygon(p0)
        
        return True
    
    def split(self,v0,e0,p0):
        #take the edges from the end of e0 to the beginning of e1
        # copy it and translate it to the beginning of e0
        e1 = self.find_equal_slope_edge(v0,e0,p0)
        edges = self.get_edge_list(v0,e0,p0)
        if edges == False:
            return False
        es = []
        for e in edges:
            es.append(e)
            if e == e1:
                break
        v = self.vertices[v0]
        v0x = v.x
        v0y = v.y
        v1 = self.get_other_vertex(v0,e0)
        v = self.vertices[v1]
        v1x = v.x
        v1y = v.y
        vx = v0x - v1x
        vy = v0y - v1y
        es2 = copy.deepcopy(es)
        for e1 in es[1:-1]:
            e = self.edges[e1]
            pe1 = self.get_poly_edge_from_edge([e.v0,e.v1])
            nv0 = self.new_vertex([pe1[0][0]+vx,pe1[0][1]+vy])
            nv1 = self.new_vertex([pe1[1][0]+vx,pe1[1][1]+vy])
            e = self.new_edge(nv0,nv1)
            if e not in es2:
                es2.append(e)
                self.add_to_polygon(p0,e)
        # make a new polygon out of the edges
        # sets ptype to 1 to signify a worm
        p = self.new_polygon(es2[0],1) 
        for e in es2[1:]:
            self.add_to_polygon(p,e)
        s = self.convert_to_shapely_polygon(p)
        self.shapely_polygons.append(s)
        self.spc += 1
        #remove edges from original polygon
        for e in es:
            self.remove_from_polygon(p0,e)
        self.shapely_polygons[p0] = self.convert_to_shapely_polygon(p0)

    def clip_all(self,p0):
        """clip worm completely into rhombs"""
        n = len(self.polygons[p0].edges)
        while n > 4:
            self.clip(p0)
            n = len(self.polygons[p0].edges)
    
    def clip(self,p0):
        num_edges = len(self.edges)
        es = copy.deepcopy(self.polygons[p0].edges)
        v0 = self.edges[es[0]].v0
        vs = self.get_vertex_list(v0,es[0],p0)
        n = len(es)/2
        v0 = self.edges[es[1]].v0
        v1 = self.edges[es[n+1]].v0
        e3 = self.new_edge(v0,v1)
        if e3 < num_edges:
            #sanity check
            #oops, try the other vertices
            v0 = self.edges[es[1]].v1
            v1 = self.edges[es[n+1]].v1
            e3 = self.new_edge(v0,v1)
            
        # make a new worm out of the edges
        p = self.new_polygon(es[0],1)
        self.remove_from_polygon(p0,es[0])
        self.add_to_polygon(p,es[1])
        self.remove_from_polygon(p0,es[1])
        self.add_to_polygon(p,e3)
        self.add_to_polygon(p,es[n+1])
        self.remove_from_polygon(p0,es[n+1])
        self.polygons[p0].edges = [e3] + self.polygons[p0].edges
        s = self.convert_to_shapely_polygon(p)
        self.shapely_polygons.append(s)
        self.shapely_polygons[p0] = self.convert_to_shapely_polygon(p0)
    
    def find_equal_slope_edge(self,v0,e0,p0):
        es = self.get_edge_list(v0,e0,p0) 
        vs = self.get_vertex_list(v0,e0,p0)
        if vs == False or len(vs) < 3:
            # sanity check
            return None
        v0 = vs[0]
        v1 = vs[1]
        vax = v1[0]
        vay = v1[1]
        dx0 = vax - v0[0]
        dy0 = vay - v0[1]

        i = 1
        for v in vs[2:]:
            dx = vax - v[0]
            dy = vay - v[1]
            if abs(dx) < 0.0001 and abs(dx0) < 0.0001:
                return es[i]
            elif abs(dx-dx0) < 0.0001 and abs(dy-dy0) < 0.0001:
                return es[i]
            vax = v[0]
            vay = v[1]
            i += 1
        
        dx = v0[0] - vax
        dy = v0[1] - vay
        if abs(dx) < 0.0001 and abs(dx0) < 0.0001:
            return es[i]
        elif abs(dx-dx0) < 0.0001 and abs(dy-dy0) < 0.0001:
            return es[i]
        return None
    
    def rhombic_triangle(self,m):
        sa = 1.0/(3*m)
        self.start_polygon(-pi*sa)        
        stn = -1
        for i in range(2,m):
            stn = stn + i - 1
            stn2 = 2*stn
            a = -pi*i*sa
            self.extend_edge(stn+i-1,stn2+1,[a])
            for n in range(0,i-2):
                self.extend_edge(stn+n,stn2+2*(n+1),[a])
            self.extend_edge(stn+2*i-1,stn2+2*(i-1),[a])
        i = m
        stn = stn + i - 1
        stn2 = 2*stn
        self.extend_edge(stn+i-1,stn2+1,[3])
        for n in range(0,i-2):
            self.extend_edge(stn+n,stn2+2*(n+1),[3])
        self.extend_edge(stn+2*i-1,stn2+2*(i-1),[3])
        
    def fan(self,m):
        sa = 1.0/(m)
        self.start_polygon(-pi*sa)        
        stn = -1
        for i in range(2,m):
            stn = stn + i - 1
            stn2 = 2*stn
            a = -pi*i*sa
            self.extend_edge(stn+i-1,stn2+1,[a])
            for n in range(0,i-2):
                self.extend_edge(stn+n,stn2+2*(n+1),[a])
            self.extend_edge(stn+2*i-1,stn2+2*(i-1),[a])
        
    def universal_tiler_tu(self,n,decomp_polygon):
        
        self.rhombic_triangle(n)
        tri = copy.deepcopy(self)
        
        v0 = (n*(n-1) + 4*n)/2
        e0 = (2*n*(n-1) + 6*n)/2 -1
        inc = len(self.edges) - n
        
        for i in range(0,5):
            self.merge(v0,e0+i*inc,tri,v0,e0-1)
        incv = v0 - n
        for i in range(0,5):
            if i == 0:
                extra = 0
            else:
                extra = 1
            self.merge(v0-1+i*incv+extra, inc-1+i*inc+extra,tri,v0-1,inc-2)
        self.merge(0, inc-1+5*inc,tri,v0-1,inc-2)
        
        self.merge(0,2*(inc-1+5*inc),decomp_polygon,0,0)
        self.merge(v0-1,14*n+7*n*(n-1)-2,decomp_polygon,0,0)
        
    def get_decompose_data(self,v0,e0,p0):
        es0 = self.get_edge_list(v0,e0,p0)
        if len(es0) < 2:
            #sanity check
            print "unable to get decompose data"
            sys.exit()
        es0a = es0[1]
        v0a = self.get_other_vertex(v0,es0a)
        es0b = es0[-1]
        v0b =  self.get_other_vertex(v0,es0b)
        return [[v0a,es0a,-1],[v0b,es0b,1]]
    
    # given adjacent edges, find a common vertex
    def get_common_vertex(self,e0,e1):
        e0 = self.edges[e0]
        e1 = self.edges[e1]
        if e0.v0 == e1.v0:
            return e0.v0
        elif e0.v0 == e1.v1:
            return e0.v0
        elif e0.v1 == e1.v0:
            return e0.v1
        elif e0.v1 == e1.v1:
            return e0.v1
        else:
            return False
        
    def decompose_zonagon(self,p0,as_worms=True):
        n = len(self.polygons[p0].edges)
        if n % 2 == 1:
            #sanity check, this is not a zonagon!
            return
        for i in range(0,(n/2)-2):
            e0 = self.polygons[p0].edges[0]
            v0 = self.edges[e0].v0
            es = self.get_edge_list(v0,e0,p0)
            va = self.get_common_vertex(es[0],es[1])
            d0 = self.get_decompose_data(va, es[0],p0)

            self.split(d0[0][0],d0[0][1],p0)
            p1 = self.pc - 1
            p = self.polygons[p1]
            if as_worms:
                p.ptype = 1
            else:
                self.clip_all(p1)
            
    def get_sym_edge_index(self,i,symclasses):
        for j in range(0,len(symclasses)):
            if i == symclasses[j]:
                return j
        return False
    
    def get_split_data(self,v0,e0,p0,p1,do_extra_inc=False):
        es0 = self.get_edge_list(v0,e0,p0)
        n = len(es0)
        extra_inc = 0
        if n % 2 == 0:
            m = (n/2)-1
            if do_extra_inc:
                extra_inc = 1
        else:
            m = n/2
        es1 = self.polygons[p1].edges
        found = False
        # seek along the p0 edges starting with e0
        # until we find an edge that is in p1        
        for i in range(1,n):
            if es0[i] in es1:
                found = True
                break
        if not found:
            #sanity check
            print "Error: could not find edge in get_split_data"
            sys.exit()
        
        i0 = i
        # now try the same thing but in the other direction
        for i in range(1,n):
            if es0[n-i] in es1:
                found = True
                break
        if not found:
            #sanity check
            print "Error: could not find edge in get_split_data"
            sys.exit()        
        i1 = i
        
        # we have found our split edges in both directions, 
        # so return the results
        if i0 > 0:
            i0b = i0 - 1
        else:
            i0b = n - 1
        if i0 < n - 1:
            i0a = i0 + 1
        else:
            i0a = 0
            
        if i1 > 0:
            i1a = i1 - 1
        else:
            i1a = n - 1
        if i1 < n - 1:
            i1b = i1 + 1
        else:
            i1b = 0
        va = self.get_common_vertex(es0[i0a],es0[i0])
        r0 = [va,es0[i0],m-i0+1,es0[i0a]-es0[i0]]
        
        vb = self.get_common_vertex(es0[n-i1a-1],es0[n-i1-1])
        r1 = [vb,es0[n-i1],m-i1+1+extra_inc,es0[i1]-es0[i1b]]
        
        return [r0,r1]        
    
    def decompose_polygon(self,p0,n,small_order,position_function=position_opposite,eliminate_sym_classes=True):
        # By default, this function tries to arrange the
        # small order regular polygons in rings consisting of rhombs 
        # and two regular polygons of opposite symmetry classes.
        # If the large order polygon has an even number of sides, 
        # this ensures that the polygon inside each ring is 
        # a zonagon and so can optionally be decomposed into rhombs.
        
        if n == 0:
            return
        large_order = len(self.polygons[p0].edges)
        s0 = large_order/small_order
        sym_lookup = position_function(large_order,small_order)
        if small_order % 2 == 0:
            if eliminate_sym_classes and (large_order == n*small_order):
                m1 = n - 1
            else:         
                m1 = n
        else:
            if large_order == n*small_order:
                m1 = n - 1
            else:
                m1 = n
        print sym_lookup
        for i in range(0,m1):
            vs,es,symclasses = self.report(p0,large_order,small_order)
            qi0 = self.get_sym_edge_index(sym_lookup[i],symclasses)
            e0 = self.polygons[p0].edges[qi0]
            v0 = self.edges[e0].v0
            if qi0 == 0:
                qi1 = len(es)-1
            else:
                qi1 = qi0 - 1
            va = self.get_common_vertex(es[qi0],es[qi1])
            self.punch(va, es[qi0],small_order)
            p1 = self.pc - 1
            r = self.get_split_data(va,es[qi0],p1,p0,eliminate_sym_classes)
            for k in range(0,2):
                s,f,m1,inc = r[k]
                for p in range(0,m1):
                    self.split(s+p*inc,f+p*inc,p0)
        
        self.report(p0,large_order,small_order)

    def clip_all_worms(self):
        # clip the worms into rhombs
        num_polygons = self.pc
        for p in range(1,num_polygons):
            poly = self.polygons[p]
            if len(poly.edges) > 4 and poly.ptype == 1:
                self.clip_all(p)                    

# a specialized function that creates and decomposes a regular polygon
# with 6n sides where n is odd.
    def decompose6_odd(self,n):      
        #self.start_polygon(6*n)
        m = n/2
        self.punch(0,0,n)
        p0 = self.pc - 1
        self.punch(3*n+1,3*n,n)
        p1 = p0 + 1
        d0 = self.get_decompose_data(0,0,p0)
        d1 = self.get_decompose_data(3*n+1,3*n,p1)
        d = [d0[0],d1[0],d1[1],d0[1]]
        for i in range(0,4):
            s,f,inc = d[i]
            for p in range(0,m):
                self.split(s+p*inc,f+p*inc,0)
        
        self.report(0,6*n,n)
        e0 = self.polygons[0].edges[0]
        v0 = self.edges[e0].v0
        es = self.get_edge_list(v0,e0,0)
        va = self.get_common_vertex(es[n-1],es[n])
        vb = self.get_common_vertex(es[3*n-1],es[3*n])
        
        self.punch(va,es[n-1],n)
        p0 = self.pc - 1
        self.punch(vb,es[3*n-1],n)
        p1 = p0 + 1
        d0 = self.get_decompose_data(va,es[n-1],p0)
        d1 = self.get_decompose_data(vb,es[3*n-1],p1)
        
        d = [d0[1],d1[1],d0[0],d1[0]]
        
        for i in range(0,4):
            s,f,inc = d[i]
            for p in range(0,m):
                self.split(s+p*inc,f+p*inc,0)
                
        self.report(0,6*n,n)
        
        e0 = self.polygons[0].edges[0]
        v0 = self.edges[e0].v0
        es = self.get_edge_list(v0,e0,0)
        va = self.get_common_vertex(es[n+m],es[n+m+1])
        
        self.punch(va, es[n+m],n)
        p0 = self.pc - 1
        d0 = self.get_decompose_data(va, es[n+m],p0)
        
        d = [d0[0],d0[1]]
        
        for i in range(0,2):
            s,f,inc = d[i]
            for p in range(0,m):
                self.split(s+p*inc,f+p*inc,0)

    def complete_at_vertex(self,v,alist,external,parity,species,rspecies):
        success = True
        vraw = [self.vertices[v].x,self.vertices[v].y]
        e0 = self.get_poly_edge_from_edge(external)
        results = []
        new_polygons = []
        for n in alist:
            if self.insufficient_gap(v,n):
                success = False
                break
            if n > 0:
                r = compute_polygon(vraw,e0,n,parity)
            else:
                r = compute_rhombus(vraw,e0,-n,parity)
            if overlaps_polygons(r[0], self.shapely_polygons+new_polygons):
                if v not in []:
                    success = False
                    break
            e0 = get_other_edge(vraw,e0,r[0])
            if len(e0) <= 1:
                success = True
                break
            #so far, so good
            results.append([n,r])
            # add to temporary shapely list
            spoly = []
            for e in r[0]:
                spoly.append(e[0])
            new_polygons.append(shapely.geometry.Polygon(spoly))
        if success:
            # it is possible to add all the polygons in the list without overlap,
            # so now we will attempt to add them to the tiling
            #save state
            state = [self.vc,self.ec,self.pc]
            for poly in results:
                n = poly[0]
                r = poly[1]
                poly_edges = r[0]
                root_edge = r[1]
                a = r[2]
                e0 = self.get_edge_from_poly_edge(root_edge)
                p = self.new_polygon(e0,n)
                vl = self.create_polygon(e0,p,poly_edges)

                # now attempt to add the sweep angles to the vertices
                if n > 0:
                    for vi in vl:
                        vertex = self.vertices[vi]
                        vertex.sweep_lu[p] = n
                        vertex.sweep_plu[p] = n
                        if not self.validate_sweep(vi,species,rspecies):
                            success = False
                            vb = vi
                            break
                else:
                    # rhombuses have two angles, so the sweep code is a bit more complex
                    if v == vl[0] or v == vl[2]:
                        ang = [a,pi-a,a,pi-a]
                    else:
                        ang = [pi-a,a,pi-a,a]

                    for i in range(0,len(vl)):
                        vertex = self.vertices[vl[i]]
                        vertex.sweep_lu[p] = -ang[i]
                        if ang[i] > (pi/2):
                            vertex.sweep_plu[p] = int(2*pi/(pi-ang[i])+0.5)
                        else:
                            vertex.sweep_plu[p] = -int((2*pi/ang[i])+0.5)
                        if not self.validate_sweep(vl[i],species,rspecies):
                            success = False
                            vb = vl[i]
                            break
    
                if success:
                    self.add_to_shapely_polygons(poly_edges)
                else:
                    # backtrack
                    sb = self.get_sweep(vb)
                    self.vc,self.ec,self.pc = state
                    self.vertices = self.vertices[:self.vc]
                    self.edges = self.edges[:self.ec]
                    self.polygons = self.polygons[:self.pc]
                    self.shapely_polygons = self.shapely_polygons[:self.pc]
                    self.clean_sweeps()
                    break

        return success

    def get_external_edges(self,v):
        ev = []
        for e in self.edges:
            if (e.v0 == v) or (e.v1 == v):
                if not e.internal:
                    ev.append(e.num)
        if len(ev) < 3:
            return ev
        # We have too many edges so need to determine which
        # ones are really external. If an edge belongs to
        # two polygons, it is internal.
        ep = {}
        for p in self.polygons:
            for ei in ev:
                if ei in p.edges:
                    if ei not in ep:
                        ep[ei] = []
                    ep[ei].append(p.num)
        ee = []
        for ei in ev:
            if ei in ep:
                if len(ep[ei]) == 2:
                    # TODO: make the next line work with backtracking
                    #edges[ei].internal = True
                    pass
                else:
                    ee.append(ei)
        return ee

    def get_leading_edge2(self,v,gap):
        ee = self.get_external_edges(v)
        if len(ee) == 0:
            return [0,[],-1]
        elif len(ee) == 1:
            ne = self.edges[ee[0]]
            nv0 = self.vertices[ne.v0]
            nv1 = self.vertices[ne.v1]
            return [0,[ne.v1,ne.v0],ne.num]
        e0 = self.edges[ee[0]]
        e1 = self.edges[ee[1]]
        if e0.v0 == v:
            v00 = self.vertices[e0.v1]
            v01 = self.vertices[v]
        else:
            v00 = self.vertices[e0.v0]
            v01 = self.vertices[v]

        if e1.v0 == v:
            v11 = self.vertices[e1.v1]
            v10 = self.vertices[v]
        else:
            v11 = self.vertices[e1.v0]
            v10 = self.vertices[v]

        if abs(gap-pi) <= 0.001:
            # need to resolve ambiguity
            v2 = self.get_previous_vertex(v,ee[0])
            vx = self.vertices[v2].x
            vy = self.vertices[v2].y
            ab = angle_between([v00.x-v01.x,v00.y-v01.y],[vx-v10.x,vy-v10.y])
            if ab > 0:
                return [ab,[v00.num,v],e0.num]
            else:
                return [ab,[v,v00.num],e0.num]
        else:
            ab = angle_between([v00.x-v01.x,v00.y-v01.y],[v11.x-v10.x,v11.y-v10.y])
            if abs(abs(ab) - gap) < 0.001:
                if ab > 0:
                    return [ab,[v,v00.num],e0.num]
                else:
                    return [ab,[v00.num,v],e0.num]
            else:
                if ab < 0:
                    return [ab,[v,v00.num],e0.num]
                else:
                    return [ab,[v00.num,v],e0.num]

    def insufficient_gap(self,v,n):
        if n < 0:
            ang = -n
        else:
            ang = (pi / n) * (n-2)
            #ang = d[n]['p']
        gap = self.compute_circle_gap_at_vertex(v)
        return ang - gap > 0.001

    # make sure that the proposed tiling is legal
    # at both ends of an edge

    def verify_external(self,v,e):
        if self.edges[e].internal:
            return False
        v2 = self.get_other_vertex(v,e)
        gap2 = self.compute_circle_gap_at_vertex(v2)
        return gap2 > 0

    def dump_polygon(self,p):
        print p, str(len(self.polygons[p].edges)) + "-gon"
        for e in self.polygons[p].edges:
            v0 = self.vertices[edges[e].v0]
            v1 = self.vertices[edges[e].v1]
            print "Edge[",e,"]:",self.edges[e].internal,self.edges[e].v0,[v0.x,v0.y],self.edges[e].v1,[v1.x,v1.y]

    def get_other_vertex(self,v,e):
        edge = self.edges[e]
        if edge.v0 != v:
            return edge.v0
        else:
            return edge.v1

    def get_sweep(self,v):
        sweep = []
        for p in self.vertices[v].sweep_p:
            sweep.append(self.vertices[v].sweep_lu[p])
        return sweep

    def get_psweep(self,v):
        sweep = []
        for p in self.vertices[v].sweep_p:
            sweep.append(self.vertices[v].sweep_plu[p])
        return sweep

    def extend_edge(self,v,e,plist):
        external = [self.edges[e].v0,self.edges[e].v1]
        if not self.complete_at_vertex(v,plist,external,1,[],[]):
            print "unable to complete vertex",v,"in extend_edge"

    def add_to_vertices(self,sweep,plist):
        for vi in range(0,self.vc):
            self.add_to_vertex(vi,sweep,plist)

    def clean_sweeps(self):
        for v in self.vertices:
            new_e = []
            vertex = self.vertices[v.num]
            for e in vertex.sweep_e:
                if e < self.ec:
                    new_e.append(e)
            vertex.sweep_e = new_e
            new_p = []
            for p in vertex.sweep_p:
                if p < self.pc:
                    new_p.append(p)
                else:
                    if p in vertex.sweep_lu:
                        del vertex.sweep_lu[p]
            vertex.sweep_p = new_p

    def clean_sweeps2(self,current_vc,todo):
        for v in self.vertices:
            new_e = []
            vertex = self.vertices[v.num]
            for e in vertex.sweep_e:
                if e < self.ec:
                    new_e.append(e)
            vertex.sweep_e = new_e
            new_p = []
            for p in vertex.sweep_p:
                if p < self.pc:
                    new_p.append(p)
                else:
                    if p in vertex.sweep_lu:
                        del vertex.sweep_lu[p]
            if len(vertex.sweep_p) != len(new_p):
                todo.append(v.num)
            vertex.sweep_p = new_p
        todo2 = []
        for i in todo:
            if i < vc and i not in todo2:
                todo2.append(i)
        return todo2
    
    def insert_slope_class(self,n,dx,dy,ec):
        ec2 = []
        found = False
        for [dx1,dy1,nums] in ec:
            if abs(dx1) < 0.0001 and abs(dx) < 0.0001:
                ec2.append([dx1,dy1,nums + [n]])
                found = True
            elif abs(dx) >= 0.0001 and abs(dx1) >= 0.0001:
                if abs((dy/dx)-(dy1/dx1)) < 0.0001:
                    ec2.append([dx1,dy1,nums + [n]])
                    found = True
                else:
                    ec2.append([dx1,dy1,nums])
            else:
                ec2.append([dx1,dy1,nums])
        if not found:
            ec2.append([dx,dy,[n]])
        return ec2
    
    def get_angle_data(self,large_order,small_order,v0,vs):
        if large_order % 2 == 0:
            r = large_order/(2*pi)
        else:
            r = large_order/pi
        n = large_order/small_order
        v0 = self.vertices[v0]
        v0 = [v0.x,v0.y]
        va = v0
        ec = []
        angles = []
        symclasses = []
        i = 0
        for v in vs[1:]:
            dx = v[0] - va[0]
            dy = v[1] - va[1]
            a = r*atan2(dy,dx)
            angles.append(a)
            if a > 0:
                s = int(a+0.5) % n
            else:
                s = int(a-0.5) % n
            symclasses.append(s)
            ec = self.insert_slope_class(i,dx,dy,ec)
            va = v
            i += 1
        
        dx = v0[0] - va[0]
        dy = v0[1] - va[1]
        a = r*atan2(dy,dx)
        angles.append(a)
        if a > 0:
            s = int(a+0.5) % n
        else:
            s = int(a-0.5) % n
        symclasses.append(s)
        ec = self.insert_slope_class(i,dx,dy,ec)
        return [angles, symclasses,ec]
    
    def get_data(self,p0,large_order,small_order):
        e0 = self.polygons[p0].edges[0]
        v0 = self.edges[e0].v0
        es = self.get_edge_list(v0,e0,p0)
        vs = self.get_vertex_list(v0,e0,p0)
        [angles, symclasses,ec] = self.get_angle_data(large_order,small_order,v0,vs)
        return [vs,es,symclasses] 
    
    def report(self,p0,large_order,small_order):
        """A version of get_data which prints out the information as well"""
        print
        print "****** Polygon report ******"
        
        print "Polygon",p0,"has",len(self.polygons[p0].edges),"edges"
        e0 = self.polygons[p0].edges[0]
        v0 = self.edges[e0].v0
        es = self.get_edge_list(v0,e0,p0)
        print "edge list",es
        vs = self.get_vertex_list(v0,e0,p0)
        print "vertex list",vs
        [angles, symclasses,ec] = self.get_angle_data(large_order,small_order,v0,vs)
        print "with angles",angles
        print "with symclasses",symclasses
        print "and",len(ec),"slope classes",ec
        return [vs,es,symclasses]

    def merge(self,v0,e0,t1,v1,e1):
        theta = self.get_merge_rotate(v0,e0,t1,v1,e1)
        vxi,vyi = self.get_merge_translate(v0,e0,t1,v1,e1)
        for p in t1.polygons:
            self.merge_polygon(v0,t1,p,vxi,vyi,theta)

    def get_merge_translate(self,v0,e0,t1,v1,e1):
        v0x = self.vertices[v0].x
        v0y = self.vertices[v0].y
        v1x = t1.vertices[v1].x
        v1y = t1.vertices[v1].y
        return (v0x-v1x), (v0y-v1y)

    def get_merge_rotate(self,v0,e0,t1,v1,e1):
        vxi,vyi = self.get_merge_translate(v0,e0,t1,v1,e1)

        # get the other vertices
        v0other = self.get_other_vertex(v0,e0)
        v1other = t1.get_other_vertex(v1,e1)

        #translate the other vertices (mapping v0 onto (0,0))
        v0x = self.vertices[v0other].x - self.vertices[v0].x
        v0y = self.vertices[v0other].y - self.vertices[v0].y
        v1x = t1.vertices[v1other].x + vxi - self.vertices[v0].x
        v1y = t1.vertices[v1other].y + vyi - self.vertices[v0].y

        #return the angle between the other vertices
        return angle_between([v1x,v1y],[v0x,v0y])

    def merge_polygon(self,v0,t1,p,vxi,vyi,theta):
        # get the vertices and map them into the new coordinates
        # we first use vxi, vyi to translate t2 coordinates into self
        # coordinates, then rotate the new polygons around
        # v0x, voy by theta radians
        v0x = self.vertices[v0].x
        v0y = self.vertices[v0].y
        edges = []
        for ei in p.edges:
            v0 = t1.vertices[t1.edges[ei].v0]
            v1 = t1.vertices[t1.edges[ei].v1]
            x0,y0 = rotate([v0.x+vxi-v0x,v0.y+vyi-v0y],theta)
            x1,y1 = rotate([v1.x+vxi-v0x,v1.y+vyi-v0y],theta)
            v0n = self.new_vertex([x0+v0x,y0+v0y])
            v1n = self.new_vertex([x1+v0x,y1+v0y])
            edges.append(self.new_edge(v0n,v1n))

        pedges = []
        for x,y in list(t1.shapely_polygons[p.num].exterior.coords):
            x0,y0 = rotate([x+vxi-v0x,y+vyi-v0y],theta)
            pedges.append([x0+v0x,y0+v0y])

        # we have the edges, so create the new polygon
        p0 = self.new_polygon(edges[0],p.ptype,False)
        for ei in edges[1:]:
            self.add_to_polygon(p0,ei)
        self.add_to_shapely_polygons2(pedges)
        
    def toSVG(self,include_labels=False):
        xml = """<?xml version="1.0" encoding="UTF-8" standalone="no"?>
<svg width="1500" height="1500"
xmlns:dc="http://purl.org/dc/elements/1.1/"
   xmlns:cc="http://web.resource.org/cc/"
   xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#"
   xmlns:svg="http://www.w3.org/2000/svg"
   xmlns="http://www.w3.org/2000/svg">
>
<defs>
<style type="text/css">
<![CDATA[
text {
font-family: Arial, sans-serif;
font-size:2;
background-color:#ffff00;
fill:#ff0000;
}
]]>
</style>
</defs>
"""
        for p in range(0,self.pc):
            xml += self.polygon_to_svg(p)
        if include_labels:    
            xml += self.labels_to_svg()
            
        xml += "</svg>"
        
        return xml

class Vertex(object):
    """Vertex class with public x, y, num and sweep attributes """

    def __init__(self,x,y):
        self.x = x
        self.y = y
        self.num = 0
        self.sweep = []
        self.sweep_e = []
        self.sweep_p = []
        self.sweep_lu = {}
        self.sweep_plu = {}
        self.parity = 0

class Edge(object):
    """Edge class with public v0, v1, num and internal attributes """

    def __init__(self,v0,v1,internal=False):
        self.v0 = v0
        self.v1 = v1
        self.num = 0
        self.internal = False

class Polygon(object):
    """Polygon class with num, edges list and polygon/rhombus type """

    def __init__(self,e0,ptype,make_internal=True):
        self.edges = [e0]
        self.ptype = ptype
        self.num = 0

class LineSegment(object):
    """line segment is an object with four values x1, y1, x2,y2"""

    def __init__(self,x1,y1,x2,y2):
        self.x1 = x1
        self.y1 = y1
        self.x2 = x2
        self.y2 = y2

def compute_polygon(vs,e0,n,parity=1):
    if n < 0:
        return compute_rhombus(vs,e0,n,parity)

    if parity == 1:
        v0 = e0[1]
        v1 = e0[0]
    else:
        v0 = e0[0]
        v1 = e0[1]

    a = pi/n

    #compute the edges
    r = 1/(2*sin(a))
    r_internal = r*cos(a)
    # compute edge angle
    ea = atan2(v1[1]-v0[1],v1[0]-v0[0])
    phi = a - ea
    x = 0.5
    y = -r_internal
    xc, yc = rotate([x,y],ea)
    centre = [xc+v0[0], yc+v0[1]]
    poly_edges = get_polygon_edges(r,centre,phi,n)
    return [poly_edges,e0,a]

#compute a new rhombus given an edge and an angle
# and add it to the end of the polygon list
# rewrite and add correct sweep bit plus fix angles
def compute_rhombus(vs,e0,a,parity=1):
    if parity == 1:
        v0 = e0[0]
        v1 = e0[1]
    else:
        v0 = e0[1]
        v1 = e0[0]

    # the following trick generates the correct rhombus shape based
    # upon the vertex ordering (vertex vs should always have angle a)
    if equal_vertex(vs,v0):
        flip = True
        a2 = a
    else:
        flip = True
        a2 = pi - a
    # compute edge angle
    ea = atan2(v1[1]-v0[1],v1[0]-v0[0])
    h = cos(0.5*(pi-a2))
    w = sin(0.5*(pi-a2))
    vlist = [[-w,0],[0,-h],[w,0],[0,h]]

    rot = []
    for v in vlist:
        x0 = v[0]+w
        y0 = v[1]
        r = rotate([x0,y0],ea+(a2/2))
        x1 = r[0]+v0[0]
        y1 = r[1]+v0[1]
        rot.append([x1,y1])
    if flip:
        pe = [[rot[1],rot[0]],[rot[2],rot[1]],[rot[3],rot[2]],[rot[0],rot[3]]]
    else:
        pe = [[rot[0],rot[1]],[rot[1],rot[2]],[rot[2],rot[3]],[rot[3],rot[0]]]
    return [pe,e0,a]

def overlaps_polygons(pl,sp,ignore=False):
    p0 = []

    for e in pl:
        p0.append(e[0])
    poly = shapely.geometry.Polygon(p0)
    i = 0
    for p in sp:
        try:
            if poly.intersects(p):
                pintersect = poly.intersection(p)
                if (type(pintersect) is shapely.geometry.polygon.Polygon) and pintersect.area > 0.01:
                    return True
        except:
            print "overlap intersect function error"
            return True
        i += 1
    return False

# compute a new regular polygon given an edge and an order
# and add it to the end of the polygon list

def rotate(v,theta):
    xc = cos(theta) * v[0] - sin(theta) * v[1]
    yc = sin(theta) * v[0] + cos(theta) * v[1]
    return xc, yc

def get_polygon_edges(r,centre,phi,n):
    pe = []
    p0 = [centre[0]+r*sin(phi),centre[1]+ r*cos(phi)]
    for i in range(1,n):
        p1 = [centre[0]+r*sin(phi+(2*i*pi/n)), centre[1]+ r*cos(phi+(2*i*pi/n))]
        pe.append([p0,p1])
        p0 = p1
    pe.append([p1,pe[0][0]])
    return pe

def different_vertex(v0,v1):
    return (abs(v0[0] - v1.x) > 0.0001) or (abs(v0[1] - v1.y) > 0.0001)

def equal_vertex(v0,v1):
    return (abs(v0[0] - v1[0]) < 0.0001) and (abs(v0[1] - v1[1]) < 0.0001)

def different_edge(v0,v1,e):
    return (v0 != e.v0 or v1 != e.v1) and (v0 != e.v1 or v1 != e.v0)

def angle_between(v0,v1):
    return atan2(v0[0]*v1[1] - v0[1]*v1[0], v0[0]*v1[0]+v0[1]*v1[1])

# edges are equal if the vertices are equal in either order
def equal_poly_edges(e,et):
    t0 = equal_vertex(e[0],et[0]) and equal_vertex(e[1],et[1])
    t1 = equal_vertex(e[0],et[1]) and equal_vertex(e[1],et[0])
    return t0 or t1

def convert_sweep(s):
    s2 = []
    for a in s:
        if a > 0:
            s2.append("p"+str(a))
        else:
            s2.append("r"+str(-180*a/pi))
    return s2

def rh(n,m=2):
    if n > 0:
        return -m*pi/n
    else:
        n = -n
        return (-(n-m)*pi/n)

def get_plist_from_sweep(sweep,sp):
    lsp = len(sp)
    for i in range(0,lsp):
        sp0 = i
        sp1 = i
        spi = i
        plist = []
        success = True
        for a in sweep:
            if abs(sp[spi] - a) < 0.0001:
                spi += 1
                if spi >= lsp:
                    spi = 0
                sp1 = spi
            else:
                success = False
                break
        if success:
            break
    if success:
        if sp1 <= sp0:
            return sp[sp1:sp0]
        else:
            return sp[sp1:]+sp[:sp0]
    else:
        return []

def equal_sweep(s0,s1):
    if len(s0) == len(s1):
        for i in range(0,len(s0)):
            if abs(s0[i]-s1[i]) > 0.0001:
                return False
        return True
    return False

# v - raw vertex
# et - raw root edge
# poly_edges - list of raw polygon edges
def get_other_edge(v,et,poly_edges):
    for e in poly_edges:
        if not equal_poly_edges(e,et) and (equal_vertex(e[0],v) or equal_vertex(e[1],v)):
            return e
    return []

### You will want to change the code after this point

# hendecagon tiling example
# You can change n below to generate a translational unit
# for a tiling containing any n-gon
print "Computing parts for translational unit."
print
print "Decomposing 6n polygon."
n = 11
t = Tiling()
t.start_polygon(6*n)
#t.decompose6_odd(n)
t.decompose_polygon(0,6,n,position_6)
#t.decompose_zonagon(0)
t.clip_all_worms()
decomp = copy.deepcopy(t)

print
print "Computing translational unit. This may take some time."
t = Tiling()
t.universal_tiler_tu(n,decomp)

print t.pc, "polygons"

xml = t.toSVG()

# You need to change where the results are saved.
fd = open("d:/tess/output.svg","w")
fd.write(xml)
fd.close()

print "Tiling output completed."