source: inundation/ga/storm_surge/alpha_shape/alpha_shape.py @ 643

"""Alpha shape   
"""

import exceptions
from Numeric import array, Float
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"

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

        
        """ Build interpolation matrix mapping from
        function values at vertices to function values at data points

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

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

    def _set_points(self, 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)

        if len (points) <= 2:
            raise PointError, "Too few points to find an alpha shape"
        
    
    def write_boundary(self,file_name):
        """
        Write the boundary to a file
        """
        #print " this info will be in the file",boundary
        export_boundary_file(file_name, self.get_boundary(),
                             OUTPUT_FILE_TITLE, delimiter = ',')
    
    def get_boundary(self):
        """
        """
        return self.boundary
    
    def _alpha_shape_algorithm(self):
        
        # A python style guide suggests using _[function name] to
        # specify interal functions
        #- this is the first time I've used it though - DSG
        
        # 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
        
#-------------------------------------------------------------
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
        alpha_shape_via_files(point_file, boundary_file, alpha)