import math

class DelaunayData:
    '''A class that provides fast high-level access to data returned from
    Delaunay triangulation.
    
    Object properties:
    selected_triangles: A set of triangle index numbers.
    selected_points:    A set of point index numbers.
    tri_points:         A map of triangle index -> tuple of three point ids.
    point_owners:       A map of point index -> set of triangle indexes.
    tri_neighbours:     A map of triangle index -> set of triangle indexes.
    regions:            A map of triangle index -> set of connected triangles.
    '''
    def __init__(self, selected_triangles, tri_points, tri_neighbours):
        '''Parameters:
        selected_triangles: An iterable of index numbers of triangles that
                            should be included in polygons.
        tri_points:         A sequence of 3-item sequences, containing the
                            points indexed by triangle id. All triangles should
                            be in same direction, either clockwise or
                            counter-clockwise.
        tri_neighbours:     A sequence of sequences, containing neighbour ids
                            for each triangle. If -1 occurs, it is ignored.
        '''
        self.selected_triangles = set(selected_triangles)
        self.selected_points = set()
        self.tri_points = {}
        self.point_owners = {}
        self.tri_neighbours = {}
        self.regions = {}
        self.allregions = []
        
        for triangle in selected_triangles:
            # Add the points of this triangle
            points = tuple(tri_points[triangle])
            self.selected_points.update(points)
            self.tri_points[triangle] = points
            
            # Map points back to triangle
            for point in points:
                if not self.point_owners.has_key(point):
                    self.point_owners[point] = set()
                self.point_owners[point].add(triangle)
            
            # Neighbours of this triangle
            neighbours = set()
            for neighbour in tri_neighbours[triangle]:
                if neighbour != -1 and neighbour in selected_triangles:
                    neighbours.add(neighbour)
            self.tri_neighbours[triangle] = neighbours
        
        missing_triangles = set(self.selected_triangles)
        while missing_triangles:
            triangle = missing_triangles.pop()
            region = self.collect_region(triangle)
            missing_triangles.difference_update(region)
            
            self.allregions.append(region)
            
            for triangle in region:
                self.regions[triangle] = region
            
    
    def collect_region(self, triangle):
        '''Walk the neighbours of a triangle until the whole region has
        been mapped.
        '''
        visited = set([triangle])
        remaining = set([triangle])
        
        while remaining:
            triangle = remaining.pop()
            neighbours = self.tri_neighbours[triangle]
            remaining.update(neighbours.difference(visited))
            visited.update(neighbours)
        
        return visited
    
    def is_edge_on_boundary(self, point1, point2):
        '''Returns True if the edge between point1 and point2 is on the
        boundary of a polygon, ie. not between two selected triangles.
        Returns False if points do not form an edge.
        '''
        p1_owners = self.point_owners[point1]
        p2_owners = self.point_owners[point2]
        edge_owners = p1_owners.intersection(p2_owners)
        return len(edge_owners) == 1
    
    def get_point_neighbours(self, point):
        '''Get neighbours of a point in the direction of triangles.
        Returns half of the connected edges, ie. only edges that are in
        the same direction as the triangles.
        '''
        result = set()
        for tri_index in self.point_owners[point]:
            points = self.tri_points[tri_index]
            myindex = list(points).index(point)
            neighbour = points[(myindex + 1) % 3]
            result.add(neighbour)
        return result
    
    def get_point_regions(self, point):
        '''Get the set of regions this point belongs in.'''
        regions = set()
        for triangle in self.point_owners[point]:
            regions.add(frozenset(self.regions[triangle]))
        return regions
    
    def is_point_on_boundary(self, point):
        '''Returns True if the point is on the boundary of a polygon.'''
        # get_point_neighbours returns only half of the edges, but a point
        # on boundary always has both forward- and backward-edge on the boundary.
        for point2 in self.get_point_neighbours(point):
            if self.is_edge_on_boundary(point, point2):
                return True
        return False

    def traverse_from_point(self, point, region):
        '''Start walking from a given point and stay always in the region
        specified by a set of triangle ids. Return a list of points on the
        boundary.
        '''
        result = []
        while point not in result:
            result.append(point)
            
            for nextpoint in self.get_point_neighbours(point):
                if not self.point_owners[nextpoint].intersection(region):
                    continue # Point in wrong region
            
                if self.is_edge_on_boundary(point, nextpoint):
                    point = nextpoint
                    break
            else:
                raise Exception('Lost track near point %s of polygon %s.' % (point, result))
        
        return result
     
    def traverse_all(self):
        '''Traverse all boundaries in the triangle set. Returns a
        set of tuples of point indexes. Outer boundaries are polygons
        in the direction of triangles, inner boundaries (holes) are in
        the opposite direction.
        
        Result polygons are normalized so that the lowest index is first.
        '''
        remaining_points = set(self.selected_points)
        polygons = set()
        while remaining_points:
            point = remaining_points.pop()
        
            if not self.is_point_on_boundary(point):
                continue
            
            for region in self.get_point_regions(point):
                polygon = self.traverse_from_point(point, region)
                remaining_points.difference_update(set(polygon))
                
                polygon = normalize_polygon(polygon)
                polygons.add(tuple(polygon))
        
        return polygons

class PolygonData:
    '''Post-processes a polygon.'''
    def __init__(self, polygon, points_x, points_y):
        '''Parameters:
        polygon: a list of point indexes
        points_x: point x-coordinates
        points_y: point y-coordinates
        '''
        self.points = list(polygon)
        self.loc = {}
        for point in self.points:
            self.loc[point] = (points_x[point], points_y[point])
        
        self.distancecache = {}
        self.boundarydistancecache = {}
    
    def get_distance(self, point1, point2):
        '''Get straight distance between two points.
        Results are cached.
        '''
        key = frozenset([point1, point2])
        
        if not self.distancecache.has_key(key):
            x1, y1 = self.loc[point1]
            x2, y2 = self.loc[point2]
            self.distancecache[key] = math.sqrt((x1 - x2) ** 2 + (y1 - y2) ** 2)
        
        return self.distancecache[key]
    
    def next_point(self, point, delta = 1):
        '''Get the next point in sequence.'''
        pos = self.points.index(point)
        return self.points[(pos + delta) % len(self.points)]
    
    def distance_along_boundary(self, point1, point2):
        '''Get distance along boundary between two points.
        Results are cached.
        '''
        key = frozenset([point1, point2])
        
        if not self.boundarydistancecache.has_key(key):
            dist1 = 0
            point = point1
            while point != point2:
                nextpoint = self.next_point(point)
                dist1 += self.get_distance(point, nextpoint)
                point = nextpoint
            
            dist2 = 0
            point = point2
            while point != point1:
                nextpoint = self.next_point(point)
                dist2 += self.get_distance(point, nextpoint)
                point = nextpoint
            
            self.boundarydistancecache[key] = min(dist1, dist2)
        
        return self.boundarydistancecache[key]
    
    def get_length(self):
        '''Get the total length of the polygon.'''
        result = 0
        for i in range(len(self.points)):
            point1 = self.points[i - 1]
            point2 = self.points[i]
            result += self.get_distance(point1, point2)
        return result
    
    
        
def normalize_polygon(polygon):
    '''Reorders the points in a polygon so that the lowest index is first. Retains
    relative order and direction of the points.
    '''
    minpos = polygon.index(min(polygon))
    return polygon[minpos:] + polygon[:minpos]

if __name__ == '__main__':
    print "Unit tests"
    
    assert normalize_polygon([5, 6, 0, 8]) == [0, 8, 5, 6]
    assert normalize_polygon([0, 1, 2]) == [0, 1, 2]
    assert normalize_polygon([2, 1, 0]) == [0, 2, 1]
    
    # Two squares:
    # 0-1  4-5
    # |/|  |/|
    # 3-2  7-6
    d = DelaunayData(range(4),  
        [[0, 1, 3], [1, 2, 3], [4, 5, 7], [5, 6, 7]],
        [[1], [0], [3], [2]])
    
    assert d.selected_triangles == set(range(4))
    assert d.selected_points == set(range(8))
    assert d.tri_points == {0: (0, 1, 3), 1: (1, 2, 3), 2: (4, 5, 7), 3: (5, 6, 7)}
    assert d.point_owners == {0: set([0]), 1: set([0, 1]), 2: set([1]), 3: set([0, 1]),
                              4: set([2]), 5: set([2, 3]), 6: set([3]), 7: set([2, 3])}
    assert d.tri_neighbours == {0: set([1]), 1: set([0]), 2: set([3]), 3: set([2])}
    assert d.regions == {0: set([0, 1]), 1: set([0, 1]), 2: set([2, 3]), 3: set([2, 3])}
    
    assert d.is_edge_on_boundary(0, 1) == True
    assert d.is_edge_on_boundary(1, 3) == False
    
    assert d.get_point_neighbours(0) == set([1])
    assert d.get_point_neighbours(3) == set([0, 1])
    
    assert d.get_point_regions(0) == set([frozenset([0, 1])])
    assert d.get_point_regions(5) == set([frozenset([2, 3])])
    
    assert d.traverse_from_point(0, d.regions[0]) == [0, 1, 2, 3]
    assert d.traverse_from_point(7, d.regions[2]) == [7, 4, 5, 6]
    
    for point in range(8):
        assert d.is_point_on_boundary(point)
    
    assert d.traverse_all() == set([(0, 1, 2, 3), (4, 5, 6, 7)])
    
    # Two triangles:
    # 0       1
    # | \   / |
    # |   2   |
    # | /   \ |
    # 3       4
    d = DelaunayData(range(2), [[0, 2, 3], [2, 1, 4]], [[], []])
    assert d.point_owners[2] == set([0, 1])
    assert d.regions == {0: set([0]), 1: set([1])}
    
    assert d.is_edge_on_boundary(1, 2) == True
    assert d.is_edge_on_boundary(0, 1) == False
    
    assert d.get_point_neighbours(2) == set([1, 3])
    assert d.get_point_neighbours(1) == set([4])
    
    assert d.is_point_on_boundary(2) == True
    
    assert d.get_point_regions(0) == d.get_point_regions(3) == set([frozenset([0])])
    assert d.get_point_regions(2) == set([frozenset([0]), frozenset([1])])
    assert d.get_point_regions(1) == d.get_point_regions(4) == set([frozenset([1])])
    
    assert d.traverse_from_point(0, d.regions[0]) == [0, 2, 3]
    assert d.traverse_from_point(2, d.regions[1]) == [2, 1, 4]
    
    assert d.traverse_all() == set([(0, 2, 3), (1, 4, 2)])
    
    # Four triangles:
    # 0 ----- 1
    # | \   / |
    # |   2   |
    # | /   \ |
    # 3 ----- 4
    d = DelaunayData(range(4),
        [[0, 2, 3], [0, 1, 2], [2, 1, 4], [2, 4, 3]],
        [[3, 1], [0, 2], [1, 3], [0, 2]])
    
    assert d.point_owners[2] == set([0, 1, 2, 3])
    assert d.regions[0] == d.regions[1] == d.regions[2] == d.regions[3] == set([0, 1, 2, 3])
    
    assert d.is_edge_on_boundary(0, 2) == False
    assert d.is_edge_on_boundary(1, 2) == False
    assert d.is_edge_on_boundary(4, 2) == False
    assert d.is_edge_on_boundary(3, 2) == False
    assert d.is_edge_on_boundary(3, 4) == True
    
    assert d.get_point_neighbours(2) == set([0, 1, 3, 4])
    
    assert d.is_point_on_boundary(2) == False
    assert d.is_point_on_boundary(0) == True
    
    for point in range(5):
        assert d.get_point_regions(point) == set([frozenset([0, 1, 2, 3])])
    
    assert d.traverse_from_point(4, d.regions[0]) == [4, 3, 0, 1]
    try:
        assert d.traverse_from_point(2, d.regions[0]) is None
    except:
        pass
    
    assert d.traverse_all() == set([(0, 1, 4, 3)])
    
    # 12 triangles, with a hole in the middle
    #    0--1
    #  / | /| \        0 1 2 3
    # 2--3--4--5
    # | /|  |/ |     4 5     6 7
    # 6--7--8--9
    #  \ | \| /        8 9 A B
    #    A--B
    #
    d = DelaunayData(range(12),
        [[0, 3, 2], [0, 1, 3], [1, 4, 3], [1, 5, 4],
         [2, 3, 6], [3, 7, 6], [4, 5, 8], [5, 9, 8],
         [6, 7, 10], [7, 11, 10], [7, 8, 11], [8, 9, 11]],
        [[1, 4], [0, 2], [1, 3], [2, 6],
         [0, 5], [4, 8], [3, 7], [6, 11],
         [5, 9], [8, 10], [9, 11], [9, 7]]
    )
    
    assert d.traverse_all() == set([(0, 1, 5, 9, 11, 10, 6, 2), (3, 7, 8, 4)])
    
    print "Unit tests OK"