"""Alpha shape   
"""

import exceptions
from Numeric import array, Float, divide_safe, sqrt
from load_mesh.loadASCII import load_xya_file, export_boundary_file

class PointError(exceptions.Exception): pass

OUTPUT_FILE_TITLE = "# The alpha shape boundary defined by point index pairs of edges"
INF = pow(10,9)

def alpha_shape_via_files(point_file, boundary_file, alpha= None):
    
    from load_mesh.loadASCII import load_xya_file
    
    point_dict = load_xya_file(point_file)
    points = point_dict['pointlist']
    #title_string = point_dict['title']
    
    alpha = Alpha_Shape(points)
    alpha.write_boundary(boundary_file)

class Alpha_Shape:

    def __init__(self, points, alpha = None):

        """ 
        Inputs:
        
          points: List of coordinate pairs [[x1, y1],[x2, y2]..] 

          alpha: alpha shape parameter
          
        """
        self._set_points(points)
        self._alpha_shape_algorithm(alpha)
            

    def _set_points(self, points):
        """
        """
        # print "setting points"
        if len (points) <= 2:
            raise PointError, "Too few points to find an alpha shape"
        
        #Convert input to Numeric arrays
        self.points = array(points).astype(Float)

    
    def write_boundary(self,file_name):
        """
        Write the boundary to a file
        """
        #print " this info will be in the file" 
        export_boundary_file(file_name, self.get_boundary(),
                             OUTPUT_FILE_TITLE, delimiter = ',')
    
    def get_boundary(self):
        """
        """
        #print "getting boundary" 
        return self.boundary

    def get_delaunay(self):
        """
        """
        return self.deltri
    
            
    def _alpha_shape_algorithm(self,alpha):
        
        # A python style guide suggests using _[function name] to
        # specify internal functions
        # - this is the first time I've used it though - DSG

        #print "starting alpha shape algorithm" 

        # Build Delaunay triangulation
        import triang
        points = []
        seglist = []
        holelist = []
        regionlist = []
        pointattlist = []
        segattlist = []
        trilist = []

        points = [(pt[0], pt[1]) for pt in self.points]
        pointattlist = [ [] for i in range(len(points)) ] 
        mode = "Qzcn"
        print "computing delaunay triangulation ... \n" 
        tridata = triang.genMesh(points,seglist,holelist,regionlist,pointattlist,segattlist,trilist,mode)
        #print tridata
        #print "point attribute list: ", tridata['generatedpointattributelist'],"\n"
        #print "hull segments: ", tridata['generatedsegmentlist'], "\n"
        self.deltri = tridata['generatedtrianglelist']
        self.deltrinbr = tridata['generatedtriangleneighborlist']
        self.hulledges = tridata['generatedsegmentlist']                              

        print "Number of delaunay triangles: ", len(self.deltri), "\n"
        #print "deltrinbrs: ", self.deltrinbr, "\n"

        # Build Alpha table
        print "Building alpha table ... \n" 
        self._tri_circumradius()
        print "Largest circumradius ", max(self.triradius)
        self._edge_intervals()
        self._vertex_intervals()

        if alpha==None:
            # Find optimum alpha
            # Ken Clarkson's hull program uses smallest alpha so that
            # every vertex is non-singular so... 
            self.optimum_alpha = max([ interval[0] for interval in self.vertexinterval])
            print "optimum alpha ", self.optimum_alpha

            #alpha_tri = self.get_alpha_triangles(self.optimum_alpha)
            #print "alpha shape triangles ", alpha_tri

            reg_edge = self.get_regular_edges(self.optimum_alpha)
            #print "alpha boundary edges", reg_edge
        else:
            reg_edge = self.get_regular_edges(alpha)
    
        self.boundary = [self.edge[k] for k in reg_edge]
        #print "alpha boundary edges ", self.boundary
            
        # this produces a baaad boundary
        #boundary = []
        #for point_index in range(len(self.points)-1):
        #    boundary.append([point_index,point_index +1])
        #boundary.append([len(self.points)-1,0])
        #self.boundary = boundary
        return

    def _tri_circumradius(self):

        # compute circumradii of the delaunay triangles
        x = self.points[:,0]
        y = self.points[:,1]
        ind1 = [self.deltri[j][0] for j in range(len(self.deltri))]
        ind2 = [self.deltri[j][1] for j in range(len(self.deltri))]
        ind3 = [self.deltri[j][2] for j in range(len(self.deltri))]

        x1 = array([ x[j] for j in ind1 ])
        y1 = array([ y[j] for j in ind1 ])
        x2 = array([ x[j] for j in ind2 ])
        y2 = array([ y[j] for j in ind2 ])
        x3 = array([ x[j] for j in ind3 ])
        y3 = array([ y[j] for j in ind3 ])

        x21 = x2-x1
        x31 = x3-x1
        y21 = y2-y1
        y31 = y3-y1

        dist21 = x21*x21 + y21*y21
        dist31 = x31*x31 + y31*y31

        denom = x21*y31 - x31*y21
        #print "denom = ", denom
                
        # dx/2, dy/2 give circumcenter relative to x1,y1. 
        #dx = (y31*dist21 - y21*dist31)/denom
        #dy = (x21*dist31 - x31*dist21)/denom
        # from Numeric import divide_safe
        try:
            dx = divide_safe(y31*dist21 - y21*dist31,denom)
            dy = divide_safe(x21*dist31 - x31*dist21,denom)
        except ZeroDivisionError:
            print "Found some degenerate triangles."
            raise
            # really need to take care of cases when denom = 0.
            # eg: something like:
            # print "Handling degenerate triangles."
            # ... add some random displacments to trouble points? 
            
        self.triradius = 0.5*sqrt(dx*dx + dy*dy)
        #print "triangle radii", self.triradius
        return

    def _edge_intervals(self):
        # for each edge, find triples
        # (length/2, min_adj_triradius, max_adj_triradius) if unattached
        # (min_adj_triradius, min_adj_triradius, max_adj_triradius) if attached.

        # It should be possible to rewrite this routine in an array-friendly form
        # like _tri_circumradius()  if we need to speed things up. 

        edges = []
        edgenbrs = []
        edgeinterval = []
        for t in range(len(self.deltri)):
            tri = self.deltri[t]
            trinbr = self.deltrinbr[t]
            dx = array([self.points[tri[(i+1)%3],0] - self.points[tri[(i+2)%3],0] for i in [0,1,2]])
            dy = array([self.points[tri[(i+1)%3],1] - self.points[tri[(i+2)%3],1] for i in [0,1,2]])
            elen = sqrt(dx*dx+dy*dy)
            #angle = array([(-dx[(i+1)%3]*dx[(i+2)%3]-dy[(i+1)%3]*dy[(i+2)%3])/(elen[(i+1)%3]*elen[(i+2)%3]) for i in [0,1,2]])
            #angle = arccos(angle)
            # really only need sign - not angle value:
            anglesign = array([(-dx[(i+1)%3]*dx[(i+2)%3]-dy[(i+1)%3]*dy[(i+2)%3]) for i in [0,1,2]])
            
            #print "dx ",dx,"\n"
            #print "dy ",dy,"\n"
            #print "edge lengths of triangle ",t,"\t",elen,"\n"
            #print "angles ",angle,"\n"
                          
            for i in [0,1,2]:
                j = (i+1)%3
                k = (i+2)%3
                if trinbr[i]==-1:
                    edges.append((tri[j], tri[k]))
                    edgenbrs.append((t, -1))
                    edgeinterval.append([0.5*elen[i], self.triradius[t], INF])
                elif (tri[j]<tri[k]):
                    edges.append((tri[j], tri[k]))
                    edgenbrs.append((t, trinbr[i]))
                    edgeinterval.append([0.5*elen[i],\
                                         min(self.triradius[t],self.triradius[trinbr[i]]),\
                                             max(self.triradius[t],self.triradius[trinbr[i]]) ])
                else:
                    continue
                #if angle[i]>pi/2:
                if anglesign[i] < 0: 
                    edgeinterval[-1][0] = edgeinterval[-1][1]
                    
        self.edge = edges
        self.edgenbr = edgenbrs
        self.edgeinterval = edgeinterval
        #print "edges: ",edges, "\n"
        #print "edge nbrs:", edgenbrs ,"\n"
        #print "edge intervals: ",edgeinterval , "\n" 

    def _vertex_intervals(self):
        # for each vertex find pairs
        # (min_adj_triradius, max_adj_triradius)
        nv = len(self.points)
        vertexnbrs = [ [] for i in range(nv)]
        vertexinterval = [ [] for i in range(nv)] 
        for t in range(len(self.deltri)):
            for j in self.deltri[t]:
                vertexnbrs[j].append(t)
        for h in range(len(self.hulledges)):
            for j in self.hulledges[h]:
                vertexnbrs[j].append(-1)
        
        for i in range(nv):
            radii = [ self.triradius[t] for t in vertexnbrs[i] if t>=0 ]
            try:
                vertexinterval[i] = [min(radii), max(radii)]
                if vertexnbrs[i][-1]==-1:
                    vertexinterval[i][1]=INF
            except ValueError:
                print "Vertex %i is isolated?"%i
                print "coords: ",self.points[i]
                print "Vertex nbrs: ", vertexnbrs[i]
                print "nbr radii: ",radii
                vertexinterval[i] = [0,INF]
                pass

        self.vertexnbr = vertexnbrs
        self.vertexinterval = vertexinterval
        #print "vertex nbrs ", vertexnbrs, "\n"
        #print "vertex intervals ",vertexinterval, "\n"
        
    def get_alpha_triangles(self,alpha):
        """
        Given an alpha value,
        return indices of triangles in the alpha-shape
        """
        def tri_rad_lta(k):
            return self.triradius[k]<=alpha

        return filter(tri_rad_lta, range(len(self.triradius)))

    def get_regular_edges(self,alpha):
        """
        Given and alpha value,
        return the edges on the boundary of the alpha-shape
        """
        def reg_edge(k):
            return self.edgeinterval[k][1]<=alpha and self.edgeinterval[k][2]>alpha

        return filter(reg_edge, range(len(self.edgeinterval)))

        
   
#-------------------------------------------------------------
if __name__ == "__main__":
    """
    Load in a data point file.
    Determine the alpha shape boundary
    Save the boundary to a file.

    usage: alpha_shape.py point_file.xya boundary_file.bnd [alpha]
    
    The alpha value is optional.
    """
    
    import os, sys
    usage = "usage: %s point_file.xya boundary_file.bnd [alpha]"%os.path.basename(sys.argv[0])
    # I made up the .bnd affix. Other ideas welcome. -DSG 
    if len(sys.argv) < 3:
        print usage
    else:
        point_file = sys.argv[1]
        boundary_file = sys.argv[2]
        if len(sys.argv) > 4:
            alpha = sys.argv[3]
        else:
            alpha = None

        print "about to call alpha shape routine \n"
        alpha_shape_via_files(point_file, boundary_file, alpha)