source: inundation/pmesh/mesh.py @ 2299

Last change on this file since 2299 was 2297, checked in by duncan, 19 years ago

new pmesh method, add_hole_from_polygon

File size: 125.0 KB
Line 
1#!/usr/bin/env python
2#
3"""General 2D triangular classes for triangular mesh generation.
4
5   Note: A .index attribute is added to objects such as vertices and
6   segments, often when reading and writing to files.  This information
7   should not be used as percistant information.  It is not the 'index' of
8   an element in a list.
9
10   
11   Copyright 2003/2004
12   Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
13   Geoscience Australia
14"""
15import load_mesh.loadASCII
16import alpha_shape.alpha_shape
17
18import sys
19import math
20import triangle.triang as triang
21import re
22import os
23import pickle
24
25import types
26import exceptions
27
28import load_mesh
29from coordinate_transforms.geo_reference import Geo_reference,DEFAULT_ZONE
30from utilities.polygon import point_in_polygon
31
32SET_COLOUR='red'
33
34# FIXME (DSG-ALL)When the dependancies between pyvolution and utils
35#are sorted out this can be changed   
36def gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2):
37    """
38    """
39   
40    det = (y2-y0)*(x1-x0) - (y1-y0)*(x2-x0)           
41    a = (y2-y0)*(q1-q0) - (y1-y0)*(q2-q0)
42    a /= det
43
44    b = (x1-x0)*(q2-q0) - (x2-x0)*(q1-q0)
45    b /= det           
46
47    return a, b
48
49##############################################
50#Initialise module
51
52#FIXME (DSG) Don't have this dependent on pyvolution!
53#Changed by Ole: Let util.py make the decision
54#from pyvolution.util import gradient
55#import compile
56#if compile.can_use_C_extension('util_ext.c'):
57#    from util_ext import gradient, point_on_line
58#else:
59#    gradient = gradient_python
60
61
62# 1st and third values must be the same
63# FIXME: maybe make this a switch that the user can change? - DSG
64initialconversions = ['', 'exterior','']
65
66#from os import sep
67#sys.path.append('..'+sep+'pmesh')
68#print "sys.path",sys.path
69
70class MeshObject:
71    """
72    An abstract superclass for the basic mesh objects, eg vertex, segment,
73    triangle.
74    """
75    def __init__(self):
76        pass
77   
78class Point(MeshObject): 
79    """
80    Define a point in a 2D space.
81    """
82    def __init__(self,X,Y):
83        __slots__ = ['x','y']
84        self.x=X
85        self.y=Y
86       
87    def DistanceToPoint(self, OtherPoint):
88        """
89        Returns the distance from this point to another
90        """
91        SumOfSquares = ((self.x - OtherPoint.x)**2) + ((self.y - OtherPoint.y)**2)
92        return math.sqrt(SumOfSquares)
93       
94    def IsInsideCircle(self, Center, Radius):
95        """
96        Return 1 if this point is inside the circle,
97        0 otherwise
98        """
99       
100        if (self.DistanceToPoint(Center)<Radius):
101            return 1
102        else:
103            return 0
104       
105    def __repr__(self):
106        return "(%f,%f)" % (self.x,self.y) 
107
108    def cmp_xy(self, point):
109        if self.x < point.x:
110            return -1
111        elif self.x > point.x:
112            return 1
113        else:           
114            if self.y < point.y:
115                return -1
116            elif self.y > point.y:
117                return 1
118            else:
119                return 0
120       
121    def same_x_y(self, point):
122        if self.x == point.x and self.y == point.y:
123            return True
124        else:
125            return False
126       
127           
128
129class Vertex(Point):
130    """
131    A point on the mesh.
132    Object attributes based on the Triangle program
133    """
134    def __init__(self,X,Y, attributes = None):
135        __slots__ = ['x','y','attributes']
136       
137        assert (type(X) == types.FloatType or type(X) == types.IntType)
138        assert (type(Y) == types.FloatType or type(Y) == types.IntType)
139        self.x=X
140        self.y=Y       
141        self.attributes=[] 
142       
143        if attributes is None:
144            self.attributes=[]
145        else:
146            self.attributes=attributes
147   
148
149    def setAttributes(self,attributes):
150        """
151        attributes is a list.
152        """
153        self.attributes = attributes
154       
155    VERTEXSQUARESIDELENGTH = 6
156    def draw(self, canvas, tags, colour = 'black',scale = 1, xoffset = 0, yoffset =0, ):
157        x =  scale*(self.x + xoffset)
158        y = -1*scale*(self.y + yoffset)  # - since for a canvas - is up
159        #print "draw x:", x
160        #print "draw y:", y
161        cornerOffset= self.VERTEXSQUARESIDELENGTH/2
162
163        # A hack to see the vert tags
164        # note: there will be many tags, since tags will not be removed
165        #when zooming
166        #canvas.create_text(x+ 2*cornerOffset,
167        #                   y+ 2*cornerOffset,
168        #                        text=tags)
169       
170        return canvas.create_rectangle(x-cornerOffset,
171                                       y-cornerOffset,
172                                       x+cornerOffset,
173                                       y+cornerOffset,
174                                       tags = tags,
175                                       outline=colour,
176                                       fill = 'white')
177   
178        #return tags
179     
180    def __repr__(self):
181        return "[(%f,%f),%r]" % (self.x,self.y,self.attributes)
182   
183class Hole(Point):
184    """
185    A region of the mesh were no triangles are generated.
186    Defined by a point in the hole enclosed by segments.
187    """
188
189    HOLECORNERLENGTH = 6
190   
191    def draw(self, canvas, tags, colour = 'purple',scale = 1, xoffset = 0, yoffset =0, ):
192        x =  scale*(self.x + xoffset)
193        y = -1*scale*(self.y + yoffset)  # - since for a canvas - is up
194        #print "draw x:", x
195        #print "draw y:", y
196        cornerOffset= self.HOLECORNERLENGTH/2
197        return canvas.create_oval(x-cornerOffset,
198                                       y-cornerOffset,
199                                       x+cornerOffset,
200                                       y+cornerOffset,
201                                       tags = tags,
202                                       outline=colour,
203                                       fill = 'white')
204   
205class Region(Point):
206    """
207    A region of the mesh, defined by a point in the region
208    enclosed by segments. Used to tag areas.
209    """
210    CROSSLENGTH = 6
211    TAG = 0
212    MAXAREA = 1
213   
214    def __init__(self,X,Y, tag = None, maxArea = None):
215        """Precondition: tag is a string and maxArea is a real
216        """
217        # This didn't work. 
218        #super(Region,self)._init_(self,X,Y)
219        self.x=X
220        self.y=Y   
221        self.attributes =[] # index 0 is the tag string
222                            #optoinal index 1 is the max triangle area
223                            #NOTE the size of this attribute is assumed
224                            # to be 1 or 2 in regionstrings2int
225        if tag is None:
226            self.attributes.append("")
227        else:
228            self.attributes.append(tag) #this is a string
229           
230        if maxArea is not None:
231            self.setMaxArea(maxArea) # maxArea is a number
232           
233    def getTag(self,):
234        return self.attributes[self.TAG]
235   
236    def setTag(self,tag):
237        self.attributes[self.TAG] = tag
238       
239    def getMaxArea(self):
240        """ Returns the Max Area of a Triangle or
241        None, if the Max Area has not been set.
242        """
243        if self.isMaxArea():
244            return self.attributes[1]
245        else:
246            return None
247   
248    def setMaxArea(self,MaxArea):
249        if self.isMaxArea(): 
250            self.attributes[self.MAXAREA] = float(MaxArea)
251        else:
252            self.attributes.append( float(MaxArea) )
253   
254    def deleteMaxArea(self):
255        if self.isMaxArea():
256            self.attributes.pop(self.MAXAREA)
257           
258    def isMaxArea(self):
259        return len(self.attributes)> 1
260   
261    def draw(self, canvas, tags, scale=1, xoffset = 0, yoffset =0, colour = "black"):
262        """
263        Draw a black cross, returning the objectID
264        """
265        x =  scale*(self.x + xoffset)
266        y = -1*scale*(self.y + yoffset) 
267        cornerOffset= self.CROSSLENGTH/2
268        return canvas.create_polygon(x,
269                                     y-cornerOffset,
270                                     x,
271                                     y,
272                                     x+cornerOffset,
273                                     y,
274                                     x,
275                                     y,
276                                     x,
277                                     y+cornerOffset,
278                                     x,
279                                     y,
280                                     x-cornerOffset,
281                                     y,
282                                     x,
283                                     y,
284                                     tags = tags,
285                                     outline = colour,fill = '')
286   
287    def __repr__(self):
288        if self.isMaxArea():
289            area = self.getMaxArea() 
290            return "(%f,%f,%s,%f)" % (self.x,self.y,
291                                      self.getTag(), self.getMaxArea())
292        else:
293            return "(%f,%f,%s)" % (self.x,self.y,
294                                   self.getTag())
295       
296class Triangle(MeshObject):
297    """
298    A triangle element, defined by 3 vertices.
299    Attributes based on the Triangle program.
300    """
301
302    def __init__(self, vertex1, vertex2, vertex3, attribute = None, neighbors = None ):
303        """
304        Vertices, the initial arguments, are listed in counterclockwise order.
305        """
306        self.vertices= [vertex1,vertex2, vertex3 ]
307       
308        if attribute is None:
309            self.attribute =""
310        else:
311            self.attribute = attribute #this is a string
312           
313        if neighbors is None:
314            self.neighbors=[]
315        else:
316            self.neighbors=neighbors
317
318    def replace(self,new_triangle):
319        self = new_triangle
320
321    def longestSideID(self):
322        ax = self.vertices[0].x
323        ay = self.vertices[0].y
324       
325        bx = self.vertices[1].x
326        by = self.vertices[1].y
327       
328        cx = self.vertices[2].x
329        cy = self.vertices[2].y
330
331        lenA = ((cx-bx)**2+(cy-by)**2)**0.5
332        lenB = ((ax-cx)**2+(ay-cy)**2)**0.5
333        lenC = ((bx-ax)**2+(by-ay)**2)**0.5
334 
335        len = [lenA,lenB,lenC]
336        return len.index(max(len))
337
338    def rotate(self,offset):
339        """
340        permute the order of the sides of the triangle
341        offset must be 0,1 or 2
342        """
343
344        if offset == 0:
345            pass
346        else:
347            if offset == 1:
348                self.vertices = [self.vertices[1],self.vertices[2],self.vertices[0]]
349                self.neighbors = [self.neighbors[1],self.neighbors[2],self.neighbors[0]]
350            if offset == 2:
351                self.vertices = [self.vertices[2],self.vertices[0],self.vertices[1]]
352                self.neighbors = [self.neighbors[2],self.neighbors[0],self.neighbors[1]]
353
354    def rotate_longest_side(self):
355        self.rotate(self.longestSideID())
356
357    def getVertices(self):
358        return self.vertices
359   
360    def get_vertices(self):
361        """
362        Return a list of the vertices.  The x and y values will be absolute
363        Easting and Northings for the zone of the current geo_ref.
364        """
365        return self.vertices
366   
367    def calcArea(self):
368        ax = self.vertices[0].x
369        ay = self.vertices[0].y
370       
371        bx = self.vertices[1].x
372        by = self.vertices[1].y
373       
374        cx = self.vertices[2].x
375        cy = self.vertices[2].y
376       
377        return abs((bx*ay-ax*by)+(cx*by-bx*cy)+(ax*cy-cx*ay))/2
378   
379    def calcP(self):
380        #calculate the perimeter
381        ax = self.vertices[0].x
382        ay = self.vertices[0].y
383       
384        bx = self.vertices[1].x
385        by = self.vertices[1].y
386       
387        cx = self.vertices[2].x
388        cy = self.vertices[2].y
389
390        a =  ((cx-bx)**2+(cy-by)**2)**0.5 
391        b =  ((ax-cx)**2+(ay-cy)**2)**0.5 
392        c =  ((bx-ax)**2+(by-ay)**2)**0.5 
393
394        return a+b+c
395           
396    def setNeighbors(self,neighbor1 = None, neighbor2 = None, neighbor3 = None):
397        """
398        neighbor1 is the triangle opposite vertex1 and so on.
399        Null represents no neighbor
400        """
401        self.neighbors = [neighbor1, neighbor2, neighbor3]
402
403    def setAttribute(self,attribute):
404        """
405        neighbor1 is the triangle opposite vertex1 and so on.
406        Null represents no neighbor
407        """
408        self.attribute = attribute #this is a string
409       
410    def __repr__(self):
411        return "[%s,%s]" % (self.vertices,self.attribute)
412       
413
414    def draw(self, canvas, tags, scale=1, xoffset = 0, yoffset =0, colour = "green"):
415        """
416        Draw a triangle, returning the objectID
417        """
418        return canvas.create_polygon(scale*(self.vertices[1].x + xoffset),
419                                     scale*-1*(self.vertices[1].y + yoffset),
420                                     scale*(self.vertices[0].x + xoffset),
421                                     scale*-1*(self.vertices[0].y + yoffset),
422                                     scale*(self.vertices[2].x + xoffset),
423                                     scale*-1*(self.vertices[2].y + yoffset),
424                                     tags = tags,
425                                     outline = colour,fill = '')
426
427class Segment(MeshObject):
428    """
429    Segments are edges whose presence in the triangulation is enforced.
430   
431    """
432    def __init__(self, vertex1, vertex2, tag = None ):
433        """
434        Each segment is specified by listing the vertices of its endpoints
435        The vertices are Vertex objects
436        """
437        assert(vertex1 != vertex2)
438        self.vertices = [vertex1,vertex2 ]
439       
440        if tag is None:
441            self.tag = self.__class__.default
442        else:
443            self.tag = tag #this is a string
444       
445    def __repr__(self):
446        return "[%s,%s]" % (self.vertices,self.tag)
447           
448       
449    def draw(self, canvas, tags,scale=1 , xoffset=0 , yoffset=0,colour='blue' ):
450        x=[]
451        y=[]
452        for end in self.vertices:
453            #end.draw(canvas,scale, xoffset, yoffset ) # draw the vertices
454            x.append(scale*(end.x + xoffset))
455            y.append(-1*scale*(end.y + yoffset)) # - since for a canvas - is up
456       
457        return canvas.create_line(x[0], y[0], x[1], y[1],
458                                  tags = tags,fill=colour)
459    def set_tag(self,tag):
460        self.tag = tag
461       
462    # Class methods
463    def set_default_tag(cls, default):
464        cls.default = default
465   
466    def get_default_tag(cls):
467        return cls.default
468   
469    set_default_tag = classmethod(set_default_tag) 
470    get_default_tag = classmethod(get_default_tag)
471
472Segment.set_default_tag("")       
473
474class Mesh:
475    """
476    Representation of a 2D triangular mesh.
477    User attributes describe the mesh region/segments/vertices/attributes
478
479    mesh attributes describe the mesh that is produced eg triangles and vertices.
480   
481    The Mesh holds user information to define the
482    """
483
484    def __repr__(self):
485        return """
486        mesh Triangles: %s 
487        mesh Attribute Titles: %s 
488        mesh Segments: %s 
489        mesh Vertices: %s 
490        user Segments: %s 
491        user Vertices: %s 
492        holes: %s 
493        regions: %s""" % (self.meshTriangles,
494                                self.attributeTitles,
495                                self.meshSegments,
496                                self.meshVertices,
497                                self.getUserSegments(),
498                                self.userVertices,
499                                self.holes,
500                                self.regions) 
501   
502    def __init__(self,
503                 userSegments=None,
504                 userVertices=None,
505                 holes=None,
506                 regions=None,
507                 geo_reference=None):
508        self.meshTriangles=[] 
509        self.attributeTitles=[] 
510        self.meshSegments=[]
511        self.meshVertices=[]
512
513        self.setID={}
514        #a dictionary of names.
515        #multiple sets are allowed, but the gui does not yet
516        #support this
517       
518        self.setID['None']=0
519        #contains the names of the sets pointing to the indexes
520        #in the list.
521       
522        self.sets=[[]]
523        #Contains the lists of triangles (triangle sets)
524
525       
526        self.visualise_graph = True
527
528        if userSegments is None:
529            self.userSegments=[]
530        else:
531            self.userSegments=userSegments
532        self.alphaUserSegments=[]
533           
534        if userVertices is None:
535            self.userVertices=[]
536        else:
537            self.userVertices=userVertices
538           
539        if holes is None:
540            self.holes=[]
541        else:
542            self.holes=holes
543           
544        if regions is None:
545            self.regions=[]
546        else:
547            self.regions=regions
548
549        if geo_reference is None:
550            self.geo_reference = Geo_reference(DEFAULT_ZONE,0,0) 
551        else:
552            self.geo_reference = geo_reference
553           
554    def __cmp__(self,other):
555       
556        # A dic for the initial m
557        dic = self.Mesh2triangList()
558        dic_mesh = self.Mesh2MeshList()
559        for element in dic_mesh.keys():
560            dic[element] = dic_mesh[element]
561        for element in dic.keys():
562            dic[element].sort()
563           
564        # A dic for the exported/imported m
565        dic_other = other.Mesh2triangList()
566        dic_mesh = other.Mesh2MeshList()
567        for element in dic_mesh.keys():
568            dic_other[element] = dic_mesh[element]
569        for element in dic.keys():
570            dic_other[element].sort()
571
572        #print "dsg************************8"
573        #print "dic ",dic
574        #print "*******8"
575        #print "mesh",dic_other
576        #print "dic.__cmp__(dic_o)",dic.__cmp__(dic_other)
577        #print "dsg************************8"
578       
579        return (dic.__cmp__(dic_other))
580   
581    def addUserPoint(self, pointType, x,y):
582        if pointType == Vertex:
583            point = self.addUserVertex(x,y)
584        if pointType == Hole:
585            point = self._addHole(x,y)
586        if pointType == Region:
587            point = self._addRegion(x,y)
588        return point
589   
590    def addUserVertex(self, x,y):
591        v=Vertex(x, y)
592        self.userVertices.append(v)
593        return v
594
595    def _addHole(self, x,y):
596        h=Hole(x, y)
597        self.holes.append(h)
598        return h
599   
600    def add_hole(self, x,y, geo_reference=None):
601        """
602        adds a point, which represents a hole.
603
604        The point data can have it's own geo_refernece.
605        If geo_reference is None the data is asumed to be absolute
606        """
607        [[x,y]] = self.geo_reference.change_points_geo_ref([x,y],
608                                                 points_geo_ref=geo_reference)
609        return self._addHole(x, y)
610
611    def _addRegion(self, x,y):
612        h=Region(x, y)
613        self.regions.append(h)
614        return h
615   
616    def add_region(self, x,y, geo_reference=None):
617        """
618        adds a point, which represents a region.
619
620        The point data can have it's own geo_refernece.
621        If geo_reference is None the data is asumed to be absolute
622        """
623        [[x,y]] = self.geo_reference.change_points_geo_ref([x,y],
624                                                 points_geo_ref=geo_reference)
625        return self._addRegion(x, y)
626
627    # Depreciated
628    def addRegionEN(self, x,y):
629        print "depreciated, use add_region"
630        return self.add_region(x,y)
631
632   
633   
634    def add_hole_from_polygon(self, polygon, tags=None, geo_reference=None):
635        """
636        Add a polygon with tags to the current mesh, as a region.
637        The maxArea of the region can be specified.
638
639        If a geo_reference of the polygon points is given, this is used.
640        If not;
641        The x,y info is assumed to be Easting and Northing, absolute,
642        for the meshes zone.
643
644        polygon a list of points, in meters that describe the polygon
645             (e.g. [[x1,y1],[x2,y2],...]
646        tags (e.g.{'wall':[0,1,3],'ocean':[2]})
647
648        This returns the region instance, so if the user whats to modify
649        it they can.       
650        """
651        return self._add_area_from_polygon(polygon, tags=tags,
652                                           hole=True,
653                                           geo_reference=geo_reference
654                                           )
655
656       
657    def add_region_from_polygon(self, polygon, tags=None,
658                                max_triangle_area=None, geo_reference=None):
659        """
660        Add a polygon with tags to the current mesh, as a region.
661        The maxArea of the region can be specified.
662
663        If a geo_reference of the polygon points is given, this is used.
664        If not;
665        The x,y info is assumed to be Easting and Northing, absolute,
666        for the meshes zone.
667
668        polygon a list of points, in meters that describe the polygon
669             (e.g. [[x1,y1],[x2,y2],...]
670        tags (e.g.{'wall':[0,1,3],'ocean':[2]})
671
672        This returns the region instance (if a max_triangle_area is given),
673        so if the user whats to modify it they can.     
674        """
675        return self._add_area_from_polygon(polygon, tags=tags,
676                                max_triangle_area=max_triangle_area,
677                                      geo_reference=geo_reference)
678       
679    def _add_area_from_polygon(self, polygon, tags=None,
680                                max_triangle_area=None,
681                               geo_reference=None,
682                               hole=False):
683        """
684        Add a polygon with tags to the current mesh, as a region.
685        The maxArea of the region can be specified.
686
687        If a geo_reference of the polygon points is given, this is used.
688        If not;
689        The x,y info is assumed to be Easting and Northing, absolute,
690        for the meshes zone.
691
692        polygon a list of points, in meters that describe the polygon
693             (e.g. [[x1,y1],[x2,y2],...]
694        tags (e.g.{'wall':[0,1,3],'ocean':[2]})
695
696        This returns the region instance, so if the user whats to modify
697        it they can.
698       
699        """
700        #get absolute values
701        if geo_reference is not None:
702            polygon = geo_reference.get_absolute(polygon)
703        # polygon is now absolute
704        #print "polygon  should be absolute",polygon
705       
706        #create points, segs and tags
707        region_dict = {}
708        region_dict['points'] = polygon
709       
710        #Create segments
711        #E.g. [[0,1], [1,2], [2,3], [3,0]]
712        #from polygon
713        #[0,1,2,3]
714        segments = []
715        N = len(polygon)
716        for i in range(N):
717            lo = i
718            hi = (lo + 1) % N
719            segments.append( [lo, hi] ) 
720        region_dict['segments'] = segments
721
722
723        #Create tags
724        #E.g. ['wall', 'wall', 'ocean', 'wall']
725        # from a dic
726        #{'wall':[0,1,3],'ocean':[2]}
727        segment_tags = ['']*N
728        if tags is not None:
729            for key in tags:
730                indices = tags[key]
731                for i in indices:
732                    segment_tags[i] = key
733        region_dict['segment_tags'] = segment_tags
734   
735        self.addVertsSegs(region_dict) #this is assuming absolute geos
736
737        inner = None
738       
739        if max_triangle_area is not None and hole is False:
740            #get inner point - absolute values
741            inner_point = point_in_polygon(polygon)
742            inner = self.add_region(inner_point[0], inner_point[1],
743                                    geo_reference=None) 
744            inner.setMaxArea(max_triangle_area)
745
746        if hole is True:
747            #get inner point - absolute values
748            inner_point = point_in_polygon(polygon)
749            inner = self.add_hole(inner_point[0], inner_point[1],
750                                    geo_reference=None) 
751           
752        return inner
753       
754    def getUserVertices(self):
755        return self.userVertices
756   
757    def getUserSegments(self):
758        allSegments = self.userSegments + self.alphaUserSegments
759        #print "self.userSegments",self.userSegments
760        #print "self.alphaUserSegments",self.alphaUserSegments
761        #print "allSegments",allSegments
762        return allSegments
763   
764    def deleteUserSegments(self,seg):
765        if self.userSegments.count(seg) == 0:
766            self.alphaUserSegments.remove(seg)
767            pass
768        else:
769            self.userSegments.remove(seg)
770           
771    def clearUserSegments(self):
772        self.userSegments = []
773        self.alphaUserSegments = []
774               
775    def getTriangulation(self):
776        return self.meshTriangles
777   
778    def getMeshVertices(self):
779        return self.meshVertices
780 
781    def getMeshSegments(self):
782        return self.meshSegments
783   
784    def getHoles(self):
785        return self.holes
786   
787    def getRegions(self):
788        return self.regions
789   
790    def isTriangulation(self):
791        if self.meshVertices == []:
792            return False 
793        else:
794            return True
795   
796    def addUserSegment(self, v1,v2):
797        """
798        PRECON: A segment between the two vertices is not already present.
799        Check by calling isUserSegmentNew before calling this function.
800       
801        """
802        s=Segment( v1,v2)
803        self.userSegments.append(s)
804        return s
805       
806    def generate_mesh(self,
807                      maximum_triangle_area=None,
808                      minimum_triangle_angle=28.0,
809                      verbose=False):
810        if verbose is True:
811            silent = ''
812        else:
813            silent = 'Q'
814        self.generateMesh(mode = silent +"pzq"+str(minimum_triangle_angle)
815                                  +"a"+str(maximum_triangle_area)
816                                  +"a")
817        #The last a is so areas for regions will be used
818       
819    def generateMesh(self, mode = None, maxArea = None, minAngle= None,
820                     isRegionalMaxAreas = True):
821        """
822        Based on the current user vaules, holes and regions
823        generate a new mesh
824        mode is a string that sets conditions on the mesh generations
825        see triangle_instructions.txt for a definition of the commands
826        PreCondition: maxArea is a double
827        """
828        #print "mode ",mode
829        if mode == None:
830            self.mode = ""
831        else:
832            self.mode = mode
833       
834        if not self.mode.find('p'):
835            self.mode += 'p' #p - read a planar straight line graph.
836            #there must be segments to use this switch
837            # TODO throw an aception if there aren't seg's
838            # it's more comlex than this.  eg holes
839        if not self.mode.find('z'):
840            self.mode += 'z' # z - Number all items starting from zero
841                             # (rather than one)
842        if self.mode.find('n'):
843            self.mode += 'n' # n - output a list of neighboring triangles
844        if not self.mode.find('A'):
845            self.mode += 'A' # A - output region attribute list for triangles
846
847        if not self.mode.find('V') and not self.mode.find('Q'):
848            self.mode += 'V' # V - output info about what Triangle is doing
849       
850        if not self.mode.find('q') and minAngle is not None:
851            #   print "**********8minAngle******** ",minAngle
852            min_angle = 'q' + str(minAngle)
853            self.mode += min_angle # z - Number all items starting from zero
854                             # (rather than one)
855        if maxArea != None:
856            self.mode += 'a' + str(maxArea)
857        #FIXME (DSG-DSG) This isn't explained.  I don't think it works
858        # well with  maxArea = None as well.
859        if isRegionalMaxAreas:
860            self.mode += 'a'
861           
862        meshDict = self.Mesh2triangList()
863        #print "!@!@ This is going to triangle   !@!@"
864        #print meshDict
865        #print "!@!@ This is going to triangle   !@!@"
866
867        #print "meshDict['segmenttaglist']", meshDict['segmenttaglist']
868        [meshDict['segmenttaglist'],
869         segconverter] =  segment_strings2ints(meshDict['segmenttaglist'],
870                                             initialconversions)
871        #print "regionlist",meshDict['regionlist']
872        [meshDict['regionlist'],
873         regionconverter] =  region_strings2ints(meshDict['regionlist'])
874        #print "regionlist",meshDict['regionlist']
875        #print "meshDict['segmenttaglist']", meshDict['segmenttaglist'
876        generatedMesh = triang.genMesh(
877                              meshDict['pointlist'],
878                              meshDict['segmentlist'],
879                              meshDict['holelist'],
880                              meshDict['regionlist'],
881                              meshDict['pointattributelist'],
882                              meshDict['segmenttaglist'],
883                              [],  # since the trianglelist isn't used
884                              self.mode)
885        #print "generated",generatedMesh
886        generatedMesh['generatedsegmentmarkerlist'] = \
887             segment_ints2strings(generatedMesh['generatedsegmentmarkerlist'],
888                                  segconverter)
889        #print "processed gen",generatedMesh['generatedsegmentmarkerlist']
890        generatedMesh['generatedtriangleattributelist'] = \
891         region_ints2strings(generatedMesh['generatedtriangleattributelist'],
892                                  regionconverter)
893
894
895        if len(generatedMesh['generatedpointattributelist'][0])==0:
896            self.attributeTitles = []
897        generatedMesh['generatedpointattributetitlelist']= \
898                                            self.attributeTitles
899
900        self.setTriangulation(generatedMesh)
901   
902    def clearTriangulation(self):
903
904        #Clear the current generated mesh values
905        self.meshTriangles=[] 
906        self.meshSegments=[]
907        self.meshVertices=[]
908
909    def removeDuplicatedUserVertices(self):
910        """Pre-condition: There are no user segments
911        This function will keep the first duplicate
912        """
913        assert self.getUserSegments() == []
914        self.userVertices, counter =  self.removeDuplicatedVertices(
915            self.userVertices)
916        return counter
917   
918    def removeDuplicatedVertices(self, Vertices):
919        """
920        This function will keep the first duplicate, remove all others
921        Precondition: Each vertex has a dupindex, which is the list
922        index.
923
924        Note: this removes vertices that have the same x,y values,
925        not duplicate instances in the Vertices list.
926        """
927        remove = []
928        index = 0
929        for v in Vertices:
930            v.dupindex = index
931            index += 1
932        t = list(Vertices)
933        t.sort(Point.cmp_xy)
934   
935        length = len(t)
936        behind = 0
937        ahead  = 1
938        counter = 0
939        while ahead < length:
940            b = t[behind]
941            ah = t[ahead]
942            if (b.y == ah.y and b.x == ah.x):
943                remove.append(ah.dupindex) 
944            behind += 1
945            ahead += 1
946
947        # remove the duplicate vertices
948        remove.sort()
949        remove.reverse()
950        for i in remove:
951            Vertices.pop(i)
952            pass
953
954        #Remove the attribute that this function added
955        for v in Vertices:
956            del v.dupindex
957        return Vertices,counter
958   
959    def thinoutVertices(self, delta):
960        """Pre-condition: There are no user segments
961        This function will keep the first duplicate
962        """
963        assert self.getUserSegments() == []
964        #t = self.userVertices
965        #self.userVertices =[]
966        boxedVertices = {}
967        thinnedUserVertices =[]
968        delta = round(delta,1)
969       
970        for v in self.userVertices :
971            # tag is the center of the boxes
972            tag = (round(v.x/delta,0)*delta,round(v.y/delta,0)*delta)
973            #this creates a dict of lists of faces, indexed by tag
974            boxedVertices.setdefault(tag,[]).append(v)
975
976        for [tag,verts] in boxedVertices.items():
977            min = delta
978            bestVert = None
979            tagVert = Vertex(tag[0],tag[1])
980            for v in verts:
981                dist = v.DistanceToPoint(tagVert)
982                if (dist<min):
983                    min = dist
984                    bestVert = v
985            thinnedUserVertices.append(bestVert)
986        self.userVertices =thinnedUserVertices
987       
988           
989    def isUserSegmentNew(self, v1,v2):
990        identicalSegs= [x for x in self.getUserSegments() \
991                        if (x.vertices[0] == v1 and x.vertices[1] == v2)
992        or (x.vertices[0] == v2 and x.vertices[1] == v1) ]
993       
994        return len(identicalSegs) == 0
995
996       
997    def deleteSegsOfVertex(self, delVertex):
998        """
999        Delete this vertex and any segments that connect to it.
1000        """
1001        #Find segments that connect to delVertex
1002        deletedSegments = []
1003        for seg in self.getUserSegments():
1004            if (delVertex in seg.vertices):
1005                deletedSegments.append(seg)
1006        # Delete segments that connect to delVertex
1007        for seg in deletedSegments:
1008            self.deleteUserSegments(seg)
1009        return deletedSegments
1010
1011   
1012    def deleteMeshObject(self, MeshObject):
1013        """
1014        Returns a list of all objects that were removed
1015        """
1016        deletedObs = []
1017        if isinstance(MeshObject, Vertex ):
1018            deletedObs = self.deleteSegsOfVertex(MeshObject)
1019            deletedObs.append(MeshObject)
1020            self.userVertices.remove(MeshObject)
1021        elif isinstance(MeshObject, Segment):
1022            deletedObs.append(MeshObject)
1023            self.deleteUserSegments(MeshObject)
1024        elif isinstance(MeshObject, Hole):
1025            deletedObs.append(MeshObject)
1026            self.holes.remove(MeshObject)
1027        elif isinstance(MeshObject, Region):
1028            deletedObs.append(MeshObject)
1029            self.regions.remove(MeshObject)         
1030        return deletedObs
1031                                                 
1032    def Mesh2triangList(self, userVertices=None,
1033                        userSegments=None,
1034                        holes=None,
1035                        regions=None):
1036        """
1037        Convert the Mesh to a dictionary of the lists needed for the
1038        triang module
1039        points list: [(x1,y1),(x2,y2),...] (Tuples of doubles)
1040        pointattributelist: [(a11,a12,...),(a21,a22),...] (Tuples of doubles)
1041        segment list: [(point1,point2),(p3,p4),...] (Tuples of integers)
1042        hole list: [(x1,y1),...](Tuples of doubles, one inside each hole)
1043        regionlist: [ (x1,y1,tag, max area),...] (Tuple of 3-4 doubles)
1044       
1045        Note, this adds an index attribute to the user Vertex objects.
1046
1047        Used to produce output to triangle
1048        """
1049        if userVertices is None:
1050            userVertices = self.getUserVertices()
1051        if userSegments is None:
1052            userSegments = self.getUserSegments()
1053        if holes is None:
1054            holes = self.getHoles()
1055        if regions is None:
1056            regions = self.getRegions()
1057           
1058        meshDict = {}
1059       
1060        pointlist=[]
1061        pointattributelist=[]
1062        index = 0
1063        for vertex in userVertices:
1064            vertex.index = index
1065            pointlist.append((vertex.x,vertex.y))
1066            pointattributelist.append((vertex.attributes))
1067           
1068            index += 1
1069        meshDict['pointlist'] = pointlist
1070        meshDict['pointattributelist'] = pointattributelist
1071
1072        segmentlist=[]
1073        segmenttaglist=[]
1074        for seg in userSegments:
1075            segmentlist.append((seg.vertices[0].index,seg.vertices[1].index))
1076            segmenttaglist.append(seg.tag)
1077        meshDict['segmentlist'] =segmentlist
1078        meshDict['segmenttaglist'] =segmenttaglist
1079       
1080        holelist=[]
1081        for hole in holes:
1082            holelist.append((hole.x,hole.y)) 
1083        meshDict['holelist'] = holelist
1084       
1085        regionlist=[]
1086        for region in regions:
1087            if (region.getMaxArea() != None): 
1088                regionlist.append((region.x,region.y,region.getTag(),
1089                               region.getMaxArea()))
1090            else:
1091                regionlist.append((region.x,region.y,region.getTag()))
1092        meshDict['regionlist'] = regionlist
1093        #print "*(*("
1094        #print meshDict
1095        #print meshDict['regionlist']
1096        #print "*(*("
1097        return meshDict
1098                                               
1099    def Mesh2MeshList(self):
1100        """
1101        Convert the Mesh to a dictionary of lists describing the
1102        triangulation variables;
1103        generated point list: [(x1,y1),(x2,y2),...] (Tuples of doubles)
1104        generated point attribute list: [(a11,a12,...),(a21,a22),...]
1105            (Tuples of doubles)
1106        generated point attribute title list:[A1Title, A2Title ...]
1107            (list of strings)
1108        generated segment list: [(point1,point2),(p3,p4),...]
1109            (Tuples of integers)
1110        generated segment tag list: [tag,tag,...] list of strings
1111
1112        generated triangle list: [(p1,p2,p3), (p4,p5,p6),....] tuple of points
1113
1114        generated triangle attribute list: [s1,s2,...] list of strings
1115
1116        generated triangle neighbor list: [(t1,t2,t3), (t4,t5,t6),....]
1117            tuple of triangles
1118       
1119        Used to produce .tsh file
1120        """
1121       
1122        meshDict = {}       
1123        pointlist=[]
1124        pointattributelist=[]
1125
1126
1127        self.maxVertexIndex=0
1128        for vertex in self.meshVertices:
1129            vertex.index = self.maxVertexIndex
1130            pointlist.append((vertex.x,vertex.y))
1131            pointattributelist.append((vertex.attributes))           
1132            self.maxVertexIndex += 1
1133
1134        meshDict['generatedpointlist'] = pointlist
1135        meshDict['generatedpointattributelist'] = pointattributelist
1136        meshDict['generatedpointattributetitlelist'] = self.attributeTitles
1137        #segments
1138        segmentlist=[]
1139        segmenttaglist=[]
1140        for seg in self.meshSegments:
1141            segmentlist.append((seg.vertices[0].index,seg.vertices[1].index))
1142            segmenttaglist.append(seg.tag)
1143        meshDict['generatedsegmentlist'] =segmentlist
1144        meshDict['generatedsegmenttaglist'] =segmenttaglist
1145
1146        # Make sure that the indexation is correct
1147        index = 0
1148        for tri in self.meshTriangles:
1149            tri.index = index           
1150            index += 1
1151
1152        trianglelist = []
1153        triangleattributelist = []
1154        triangleneighborlist = []
1155        for tri in self.meshTriangles: 
1156            trianglelist.append((tri.vertices[0].index,tri.vertices[1].index,
1157                                 tri.vertices[2].index)) 
1158            triangleattributelist.append([tri.attribute])
1159            neighborlist = [-1,-1,-1]
1160            for neighbor,index in map(None,tri.neighbors,
1161                                      range(len(tri.neighbors))):
1162                if neighbor:   
1163                    neighborlist[index] = neighbor.index
1164            triangleneighborlist.append(neighborlist)
1165       
1166        meshDict['generatedtrianglelist'] = trianglelist
1167        meshDict['generatedtriangleattributelist'] = triangleattributelist
1168        meshDict['generatedtriangleneighborlist'] = triangleneighborlist
1169       
1170        #print "mesh.Mesh2MeshList*)*)"
1171        #print meshDict
1172        #print "mesh.Mesh2MeshList*)*)"
1173
1174        return meshDict
1175
1176                               
1177    def Mesh2MeshDic(self):
1178        """
1179        Convert the user and generated info of a mesh to a dictionary
1180        structure
1181        """
1182        dic = self.Mesh2triangList()
1183        dic_mesh = self.Mesh2MeshList()
1184        for element in dic_mesh.keys():
1185            dic[element] = dic_mesh[element]
1186        return dic
1187   
1188    def setTriangulation(self, genDict):
1189        """
1190        Set the mesh attributes given a dictionary of the lists
1191        returned from the triang module       
1192        generated point list: [(x1,y1),(x2,y2),...] (Tuples of doubles) 
1193        generated point attribute list:[(P1att1,P1attt2, ...),
1194            (P2att1,P2attt2,...),...]
1195        generated point attribute title list:[A1Title, A2Title ...]
1196            (list of strings)
1197        generated segment list: [(point1,point2),(p3,p4),...]
1198            (Tuples of integers)
1199        generated segment marker list: [S1Tag, S2Tag, ...] (list of ints)
1200        triangle list:  [(point1,point2, point3),(p5,p4, p1),...]
1201            (Tuples of integers)
1202        triangle neighbor list: [(triangle1,triangle2, triangle3),
1203            (t5,t4, t1),...] (Tuples of integers)
1204            -1 means there's no triangle neighbor
1205        triangle attribute list: [(T1att), (T2att), ...]
1206            (list of a list of strings)
1207        """
1208        #Clear the current generated mesh values
1209        self.meshTriangles=[]
1210        self.attributeTitles=[]
1211        self.meshSegments=[]
1212        self.meshVertices=[]
1213
1214        #print "mesh.setTriangulation@#@#@#"
1215        #print genDict
1216        #print "@#@#@#"
1217
1218        self.maxVertexIndex = 0
1219        for point in genDict['generatedpointlist']:
1220            v=Vertex(point[0], point[1])
1221            v.index =  self.maxVertexIndex
1222            self.maxVertexIndex +=1
1223            self.meshVertices.append(v)
1224
1225        self.attributeTitles = genDict['generatedpointattributetitlelist']
1226
1227        index = 0
1228        for seg,marker in map(None,genDict['generatedsegmentlist'],
1229                              genDict['generatedsegmentmarkerlist']):
1230            segObject = Segment( self.meshVertices[seg[0]],
1231                           self.meshVertices[seg[1]], tag = marker )
1232            segObject.index = index
1233            index +=1
1234            self.meshSegments.append(segObject)
1235
1236        index = 0
1237        for triangle in genDict['generatedtrianglelist']:
1238            tObject =Triangle( self.meshVertices[triangle[0]],
1239                        self.meshVertices[triangle[1]],
1240                        self.meshVertices[triangle[2]] )
1241            tObject.index = index
1242            index +=1
1243            self.meshTriangles.append(tObject)
1244
1245        index = 0
1246        for att in genDict['generatedtriangleattributelist']:
1247            if att == []:
1248                self.meshTriangles[index].setAttribute("")
1249            else:
1250                self.meshTriangles[index].setAttribute(att[0])
1251            index += 1
1252           
1253        index = 0
1254        for att in genDict['generatedpointattributelist']:
1255            if att == None:
1256                self.meshVertices[index].setAttributes([])
1257            else:
1258                self.meshVertices[index].setAttributes(att)
1259            index += 1
1260   
1261        index = 0
1262        for triangle in genDict['generatedtriangleneighborlist']:
1263            # Build a list of triangle object neighbors
1264            ObjectNeighbor = []
1265            for neighbor in triangle:
1266                if ( neighbor != -1):
1267                    ObjectNeighbor.append(self.meshTriangles[neighbor])
1268                else:
1269                    ObjectNeighbor.append(None)
1270            self.meshTriangles[index].setNeighbors(ObjectNeighbor[0],
1271                                                   ObjectNeighbor[1],
1272                                                   ObjectNeighbor[2])
1273            index += 1
1274
1275
1276    def setMesh(self, genDict):
1277        """
1278        Set the user Mesh attributes given a dictionary of the lists
1279        point list: [(x1,y1),(x2,y2),...] (Tuples of doubles) 
1280        point attribute list:[(P1att1,P1attt2, ...),(P2att1,P2attt2,...),...]
1281        segment list: [(point1,point2),(p3,p4),...] (Tuples of integers)
1282        segment tag list: [S1Tag, S2Tag, ...] (list of ints)
1283        region list: [(x1,y1),(x2,y2),...] (Tuples of doubles)
1284        region attribute list: ["","reservoir",""] list of strings
1285        region max area list:[real, None, Real,...] list of None and reals
1286       
1287        mesh is an instance of a mesh object
1288        """
1289        #Clear the current user mesh values
1290        self.clearUserSegments()
1291        self.userVertices=[]
1292        self.Holes=[]
1293        self.Regions=[]
1294
1295        #print "mesh.setMesh@#@#@#"
1296        #print genDict
1297        #print "@#@#@#"
1298       
1299        #index = 0
1300        for point in genDict['pointlist']:
1301            v=Vertex(point[0], point[1])
1302            #v.index = index
1303            #index +=1
1304            self.userVertices.append(v)
1305
1306        #index = 0
1307        for seg,tag in map(None,genDict['segmentlist'],
1308                           genDict['segmenttaglist']):
1309            segObject = Segment( self.userVertices[seg[0]],
1310                           self.userVertices[seg[1]], tag = tag )
1311            #segObject.index = index
1312            #index +=1
1313            self.userSegments.append(segObject)
1314
1315# Remove the loading of attribute info.
1316# Have attribute info added using least_squares in pyvolution
1317#         index = 0
1318#         for att in genDict['pointattributelist']:
1319#             if att == None:
1320#                 self.userVertices[index].setAttributes([])
1321#             else:
1322#                 self.userVertices[index].setAttributes(att)
1323#            index += 1
1324       
1325        #index = 0
1326        for point in genDict['holelist']:
1327            h=Hole(point[0], point[1])
1328            #h.index = index
1329            #index +=1
1330            self.holes.append(h)
1331
1332        #index = 0
1333        for reg,att,maxArea in map(None,
1334                                   genDict['regionlist'],
1335                                   genDict['regionattributelist'],
1336                                   genDict['regionmaxarealist']):
1337            Object = Region( reg[0],
1338                             reg[1],
1339                             tag = att,
1340                             maxArea = maxArea)
1341            #Object.index = index
1342            #index +=1
1343            self.regions.append(Object)
1344       
1345    def addVertsSegs(self, outlineDict):
1346        """
1347        Add out-line (user Mesh) attributes given a dictionary of the lists
1348        points: [(x1,y1),(x2,y2),...] (Tuples of doubles) 
1349        segments: [(point1,point2),(p3,p4),...] (Tuples of integers)
1350        segment_tags: [S1Tag, S2Tag, ...] (list of strings)
1351
1352        Assume the values are in Eastings and Northings, with no reference
1353        point. eg absolute
1354        """
1355        if not outlineDict.has_key('segment_tags'):
1356            outlineDict['segment_tags'] = []
1357            for i in range(len(outlineDict['segments'])):
1358                outlineDict['segment_tags'].append('')
1359        #print "outlineDict['segment_tags']",outlineDict['segment_tags']
1360        #print "outlineDict['points']",outlineDict['points']
1361        #print "outlineDict['segments']",outlineDict['segments']
1362       
1363        localUserVertices = []
1364        #index = 0
1365        for point in outlineDict['points']:
1366            v=Vertex(point[0]-self.geo_reference.xllcorner,
1367                     point[1]-self.geo_reference.yllcorner)
1368            #v.index = index
1369            #index +=1
1370            self.userVertices.append(v)
1371            localUserVertices.append(v)
1372           
1373        #index = 0
1374        for seg,seg_tag in map(None,outlineDict['segments'],
1375                       outlineDict['segment_tags']):
1376            segObject = Segment( localUserVertices[seg[0]],
1377                           localUserVertices[seg[1]] )
1378            if not seg_tag == '':
1379                segObject.set_tag(seg_tag)
1380            #segObject.index = index
1381            #index +=1
1382            self.userSegments.append(segObject)
1383            #DSG!!!
1384           
1385    def TestautoSegment(self):
1386        newsegs = []
1387        s1 = Segment(self.userVertices[0],
1388                               self.userVertices[1])
1389        s2 = Segment(self.userVertices[0],
1390                               self.userVertices[2])
1391        s3 = Segment(self.userVertices[2],
1392                               self.userVertices[1])
1393        if self.isUserSegmentNew(s1.vertices[0],s1.vertices[1]):
1394            newsegs.append(s1)
1395        if self.isUserSegmentNew(s2.vertices[0],s2.vertices[1]):
1396            newsegs.append(s2)
1397        if self.isUserSegmentNew(s3.vertices[0],s3.vertices[1]):
1398            newsegs.append(s3)
1399        #DSG!!!
1400        self.userSegments.extend(newsegs)
1401        return newsegs
1402
1403   
1404    def savePickle(self, currentName):
1405        fd = open(currentName, 'w')
1406        pickle.dump(self,fd)
1407        fd.close()
1408
1409    def autoSegmentHull(self):
1410        """
1411        initially work by running an executable
1412        Later compile the c code with a python wrapper.
1413
1414        Precon: There must be 3 or more vertices in the userVertices structure
1415        """
1416        newsegs = []
1417        inputfile = 'hull_in.txt'
1418        outputfile = inputfile + '-alf'
1419        #write vertices to file
1420        fd = open(inputfile,'w')
1421        for v in self.userVertices:
1422            fd.write(str(v.x))
1423            fd.write(' ')
1424            fd.write(str(v.y))
1425            fd.write('\n')
1426        fd.close()
1427       
1428        #run hull executable
1429        #warning need to compile hull for the current operating system
1430        command = 'hull.exe -A -i ' + inputfile
1431        os.system(command)
1432       
1433        #read results into this object
1434        fd = open(outputfile)
1435        lines = fd.readlines()
1436        fd.close()
1437        #print "(*(*(*("
1438        #print lines
1439        #print "(*(*(*("
1440        lines.pop(0) #remove the first (title) line
1441        for line in lines:
1442            vertindexs = line.split()
1443            #print 'int(vertindexs[0])', int(vertindexs[0])
1444            #print 'int(vertindexs[1])', int(vertindexs[1])
1445            #print 'self.userVertices[int(vertindexs[0])]' ,self.userVertices[int(vertindexs[0])]
1446            #print 'self.userVertices[int(vertindexs[1])]' ,self.userVertices[int(vertindexs[1])]
1447            v1 = self.userVertices[int(vertindexs[0])]
1448            v2 = self.userVertices[int(vertindexs[1])]
1449           
1450            if self.isUserSegmentNew(v1,v2):
1451                newseg = Segment(v1, v2)
1452                newsegs.append(newseg)
1453            #DSG!!!
1454        self.userSegments.extend(newsegs)
1455        return newsegs
1456    def autoSegmentFilter(self,raw_boundary=True,
1457                    remove_holes=False,
1458                    smooth_indents=False,
1459                    expand_pinch=False):
1460        """
1461        Precon: There is a self.shape
1462        """
1463        #FIXME remove the precon.  Internally check
1464        return self._boundary2mesh(raw_boundary=raw_boundary,
1465                    remove_holes=remove_holes,
1466                    smooth_indents=smooth_indents,
1467                    expand_pinch=expand_pinch)
1468   
1469       
1470   
1471    def autoSegment(self, alpha = None,
1472                    raw_boundary=True,
1473                    remove_holes=False,
1474                    smooth_indents=False,
1475                    expand_pinch=False): 
1476        """
1477        Precon: There must be 3 or more vertices in the userVertices structure
1478        """
1479        self._createBoundary(alpha=alpha)
1480        return self._boundary2mesh(raw_boundary=raw_boundary,
1481                    remove_holes=remove_holes,
1482                    smooth_indents=smooth_indents,
1483                    expand_pinch=expand_pinch)
1484
1485    def _createBoundary(self,alpha=None):
1486        """
1487        """
1488        points=[]
1489        for vertex in self.getUserVertices():
1490            points.append((vertex.x,vertex.y))
1491        self.shape = alpha_shape.alpha_shape.Alpha_Shape(points, alpha = alpha)
1492
1493
1494    def _boundary2mesh(self, raw_boundary=True,
1495                    remove_holes=False,
1496                    smooth_indents=False,
1497                    expand_pinch=False):
1498        """
1499        Precon there must be a shape object.
1500        """
1501        self.shape.set_boundary_type(raw_boundary=raw_boundary,
1502                                 remove_holes=remove_holes,
1503                                 smooth_indents=smooth_indents,
1504                                 expand_pinch=expand_pinch)
1505        boundary_segs = self.shape.get_boundary()
1506        #print "boundary_segs",boundary_segs
1507        segs2delete = self.alphaUserSegments
1508        #FIXME(DSG-DSG) this algorithm needs comments
1509        new_segs = {}
1510        #alpha_segs = []
1511        #user_segs = []
1512        for seg in boundary_segs:
1513            v1 = self.userVertices[int(seg[0])]
1514            v2 = self.userVertices[int(seg[1])]
1515            boundary_seg = Segment(v1, v2)
1516            new_segs[(v1,v2)] = boundary_seg
1517
1518        for user_seg in self.userSegments:
1519            if new_segs.has_key((user_seg.vertices[0],
1520                                user_seg.vertices[1])):
1521                del new_segs[user_seg.vertices[0],
1522                                user_seg.vertices[1]]
1523            elif new_segs.has_key((user_seg.vertices[1],
1524                                user_seg.vertices[0])):
1525                del new_segs[user_seg.vertices[1],
1526                                user_seg.vertices[0]]
1527               
1528        optimum_alpha = self.shape.get_alpha()
1529        alpha_segs_no_user_segs  = new_segs.values()
1530        self.alphaUserSegments = alpha_segs_no_user_segs
1531        return alpha_segs_no_user_segs, segs2delete, optimum_alpha
1532   
1533    def _boundary2mesh_old(self, raw_boundary=True,
1534                    remove_holes=False,
1535                    smooth_indents=False,
1536                    expand_pinch=False):
1537        """
1538        Precon there must be a shape object.
1539        """
1540        self.shape.set_boundary_type(raw_boundary=raw_boundary,
1541                                 remove_holes=remove_holes,
1542                                 smooth_indents=smooth_indents,
1543                                 expand_pinch=expand_pinch)
1544        boundary_segs = self.shape.get_boundary()
1545        #print "boundary_segs",boundary_segs
1546        segs2delete = self.alphaUserSegments
1547
1548        #FIXME(DSG-DSG) this algorithm needs comments
1549        #FIXME(DSG-DSG) can it be sped up?  It's slow
1550        new_segs = []
1551        alpha_segs = []
1552        user_segs = []
1553        for seg in boundary_segs:
1554            v1 = self.userVertices[int(seg[0])]
1555            v2 = self.userVertices[int(seg[1])]
1556            alpha_seg = self.representedAlphaUserSegment(v1, v2)
1557            user_seg = self.representedUserSegment(v1, v2)
1558            #DSG!!!
1559            assert not(not (alpha_seg == None) and not (user_seg == None))
1560            if not alpha_seg == None:
1561                alpha_segs.append(alpha_seg)
1562            elif not user_seg  == None:
1563                user_segs.append(user_seg)
1564            else:
1565                unique_seg = Segment(v1, v2)
1566                new_segs.append(unique_seg)
1567               
1568            for seg in alpha_segs:
1569                try:
1570                    segs2delete.remove(seg)
1571                except:
1572                    pass
1573       
1574        self.alphaUserSegments = []
1575        self.alphaUserSegments.extend(new_segs)
1576        self.alphaUserSegments.extend(alpha_segs)
1577
1578        optimum_alpha = self.shape.get_alpha()
1579        # need to draw newsegs
1580        return new_segs, segs2delete, optimum_alpha
1581   
1582    def representedAlphaUserSegment(self, v1,v2):
1583        identicalSegs= [x for x in self.alphaUserSegments \
1584                        if (x.vertices[0] == v1 and x.vertices[1] == v2)
1585        or (x.vertices[0] == v2 and x.vertices[1] == v1) ]
1586
1587        if identicalSegs == []:
1588            return None
1589        else:
1590            # Only return the first one.
1591            return identicalSegs[0]
1592   
1593    def representedUserSegment(self, v1,v2):
1594        identicalSegs= [x for x in self.userSegments \
1595                        if (x.vertices[0] == v1 and x.vertices[1] == v2)
1596        or (x.vertices[0] == v2 and x.vertices[1] == v1) ]
1597
1598        if identicalSegs == []:
1599            return None
1600        else:
1601            # Only return the first one.
1602            return identicalSegs[0]
1603       
1604    def joinVertices(self):
1605        """
1606        Return list of segments connecting all userVertices
1607        in the order they were given
1608       
1609        Precon: There must be 3 or more vertices in the userVertices structure
1610        """
1611
1612        newsegs = []
1613       
1614        v1 = self.userVertices[0]
1615        for v2 in self.userVertices[1:]:
1616            if self.isUserSegmentNew(v1,v2):           
1617                newseg = Segment(v1, v2)
1618                newsegs.append(newseg)
1619            v1 = v2
1620
1621        #Connect last point to the first
1622        v2 = self.userVertices[0]       
1623        if self.isUserSegmentNew(v1,v2):           
1624                newseg = Segment(v1, v2)
1625                newsegs.append(newseg)
1626           
1627
1628        #Update list of user segments       
1629        #DSG!!!
1630        self.userSegments.extend(newsegs)               
1631        return newsegs
1632   
1633    def normaliseMesh(self,scale, offset, height_scale):
1634        [xmin, ymin, xmax, ymax] = self.boxsize()
1635        [attmin0, attmax0] = self.maxMinVertAtt(0)
1636        #print "[attmin0, attmax0]" ,[attmin0, attmax0]
1637        [attmin1, attmax1] = self.maxMinVertAtt(1)
1638        #print [xmin, ymin, xmax, ymax]
1639        xrange = xmax - xmin
1640        yrange = ymax - ymin
1641        if xrange > yrange:
1642            min,max = xmin, xmax
1643        else:
1644            min,max = ymin, ymax
1645           
1646        for obj in self.getUserVertices():
1647            obj.x = (obj.x - xmin)/(max- min)*scale + offset
1648            obj.y = (obj.y - ymin)/(max- min)*scale + offset
1649            if len(obj.attributes)  > 0 and attmin0 != attmax0:
1650                obj.attributes[0] = (obj.attributes[0]-attmin0)/ \
1651                                    (attmax0-attmin0)*height_scale
1652            if len(obj.attributes) > 1 and attmin1 != attmax1:
1653                obj.attributes[1] = (obj.attributes[1]-attmin1)/ \
1654                                    (attmax1-attmin1)*height_scale
1655           
1656        for obj in self.getMeshVertices():
1657            obj.x = (obj.x - xmin)/(max- min)*scale + offset
1658            obj.y = (obj.y - ymin)/(max- min)*scale + offset
1659            if len(obj.attributes)  > 0 and attmin0 != attmax0:
1660                obj.attributes[0] = (obj.attributes[0]-attmin0)/ \
1661                                    (attmax0-attmin0)*height_scale
1662            if len(obj.attributes) > 1 and attmin1 != attmax1:
1663                obj.attributes[1] = (obj.attributes[1]-attmin1)/ \
1664                                    (attmax1-attmin1)*height_scale
1665               
1666        for obj in self.getHoles():
1667            obj.x = (obj.x - xmin)/(max- min)*scale + offset
1668            obj.y = (obj.y - ymin)/(max- min)*scale + offset
1669        for obj in self.getRegions():
1670            obj.x = (obj.x - xmin)/(max- min)*scale + offset
1671            obj.y = (obj.y - ymin)/(max- min)*scale + offset
1672        [xmin, ymin, xmax, ymax] = self.boxsize()
1673        #print [xmin, ymin, xmax, ymax]
1674     
1675    def boxsizeVerts(self):
1676        """
1677        Returns a list of verts denoting a box or triangle that contains
1678        verts on the xmin, ymin, xmax and ymax axis.
1679        Structure: list of verts
1680        """
1681        # FIXME dsg!!! large is a hack
1682        #You want the kinds package, part of Numeric:
1683        #In [2]: import kinds
1684       
1685        #In [3]: kinds.default_float_kind.M
1686        #kinds.default_float_kind.MAX         kinds.default_float_kind.MIN
1687    #kinds.default_float_kind.MAX_10_EXP  kinds.default_float_kind.MIN_10_EXP
1688        #kinds.default_float_kind.MAX_EXP     kinds.default_float_kind.MIN_EXP
1689
1690        #In [3]: kinds.default_float_kind.MIN
1691        #Out[3]: 2.2250738585072014e-308
1692
1693        large = 1e100
1694        xmin= large
1695        xmax=-large
1696        ymin= large
1697        ymax=-large
1698        for vertex in self.userVertices:
1699            if vertex.x < xmin:
1700                xmin = vertex.x
1701                xminVert = vertex
1702            if vertex.x > xmax:
1703                xmax = vertex.x
1704                xmaxVert = vertex
1705               
1706            if vertex.y < ymin:
1707                ymin = vertex.y
1708                yminVert = vertex
1709            if vertex.y > ymax:
1710                ymax = vertex.y
1711                ymaxVert = vertex
1712        verts, count = self.removeDuplicatedVertices([xminVert,
1713                                                      xmaxVert,
1714                                                      yminVert,
1715                                                      ymaxVert])
1716         
1717        return verts
1718   
1719    def boxsize(self):
1720        """
1721        Returns a list denoting a box that contains the entire
1722        structure of vertices
1723        Structure: [xmin, ymin, xmax, ymax]
1724        """
1725        # FIXME dsg!!! large is a hack
1726        #You want the kinds package, part of Numeric:
1727        #In [2]: import kinds
1728       
1729        #In [3]: kinds.default_float_kind.M
1730        #kinds.default_float_kind.MAX         kinds.default_float_kind.MIN
1731    #kinds.default_float_kind.MAX_10_EXP  kinds.default_fltesting oat_kind.MIN_10_EXP
1732        #kinds.default_float_kind.MAX_EXP     kinds.default_float_kind.MIN_EXP
1733
1734        #In [3]: kinds.default_float_kind.MIN
1735        #Out[3]: 2.2250738585072014e-308
1736
1737        large = 1e100
1738        xmin= large
1739        xmax=-large
1740        ymin= large
1741        ymax=-large
1742        for vertex in self.userVertices:
1743            if vertex.x < xmin:
1744                xmin = vertex.x
1745            if vertex.x > xmax:
1746                xmax = vertex.x
1747               
1748            if vertex.y < ymin:
1749                ymin = vertex.y
1750            if vertex.y > ymax:
1751                ymax = vertex.y
1752        return [xmin, ymin, xmax, ymax]
1753 
1754    def maxMinVertAtt(self, iatt):
1755        """
1756        Returns a list denoting a box that contains the entire structure
1757        of vertices
1758        Structure: [xmin, ymin, xmax, ymax]
1759        """
1760        # FIXME dsg!!! large is a hacktesting
1761        #You want the kinds package, part of Numeric:
1762        #In [2]: import kinds
1763       
1764        #In [3]: kinds.default_float_kind.M
1765        #kinds.default_float_kind.MAX         kinds.default_float_kind.MIN
1766    #kinds.default_float_kind.MAX_10_EXP  kinds.default_float_kind.MIN_10_EXP
1767        #kinds.default_float_kind.MAX_EXP     kinds.default_float_kind.MIN_EXP
1768
1769        #In [3]: kinds.default_float_kind.MIN
1770        #Out[3]: 2.2250738585072014e-308
1771
1772        large = 1e100
1773        min= large
1774        max=-large
1775        for vertex in self.userVertices:
1776            if len(vertex.attributes) > iatt:
1777                if vertex.attributes[iatt] < min:
1778                    min = vertex.attributes[iatt]
1779                if vertex.attributes[iatt] > max:
1780                    max = vertex.attributes[iatt] 
1781        for vertex in self.meshVertices:
1782            if len(vertex.attributes) > iatt:
1783                if vertex.attributes[iatt] < min:
1784                    min = vertex.attributes[iatt]
1785                if vertex.attributes[iatt] > max:
1786                    max = vertex.attributes[iatt] 
1787        return [min, max]
1788   
1789    def scaleoffset(self, WIDTH, HEIGHT):
1790        """
1791        Returns a list denoting the scale and offset terms that need to be
1792        applied when converting  mesh co-ordinates onto grid co-ordinates
1793        Structure: [scale, xoffset, yoffset]
1794        """   
1795        OFFSET = 0.05*min([WIDTH, HEIGHT])
1796        [xmin, ymin, xmax, ymax] = self.boxsize()
1797        SCALE = min([0.9*WIDTH, 0.9*HEIGHT])/max([xmax-xmin, ymax-ymin])
1798       
1799        if SCALE*xmin < OFFSET:
1800            xoffset = abs(SCALE*xmin) + OFFSET
1801        if SCALE*xmax > WIDTH - OFFSET:
1802            xoffset= -(SCALE*xmax - WIDTH + OFFSET)
1803        if SCALE*ymin < OFFSET:
1804            b = abs(SCALE*ymin)+OFFSET
1805        if SCALE*ymax > HEIGHT-OFFSET:
1806            b = -(SCALE*ymax - HEIGHT + OFFSET)
1807        yoffset = HEIGHT - b
1808        return [SCALE, xoffset, yoffset]
1809
1810         
1811    def exportASCIIobj(self,ofile):
1812        """
1813        export a file, ofile, with the format
1814         lines:  v <x> <y> <first attribute>
1815        f <vertex #>  <vertex #> <vertex #> (of the triangles)
1816        """
1817        fd = open(ofile,'w')
1818        self.writeASCIIobj(fd)   
1819        fd.close()
1820
1821
1822    def writeASCIIobj(self,fd):
1823        fd.write(" # Triangulation as an obj file\n")
1824        numVert = str(len(self.meshVertices))
1825       
1826        index1 = 1 
1827        for vert in self.meshVertices:
1828            vert.index1 = index1
1829            index1 += 1
1830           
1831            fd.write("v "
1832                     + str(vert.x) + " "
1833                     + str(vert.y) + " "
1834                     + str(vert.attributes[0]) + "\n")
1835           
1836        for tri in self.meshTriangles:
1837            fd.write("f "
1838                     + str(tri.vertices[0].index1) + " " 
1839                     + str(tri.vertices[1].index1) + " " 
1840                     + str(tri.vertices[2].index1) + "\n")
1841           
1842    def exportASCIIsegmentoutlinefile(self,ofile):
1843        """Write the boundary user mesh info, eg
1844         vertices that are connected to segments,
1845         segments
1846        """
1847       
1848        verts = {}
1849        for seg in self.getUserSegments():
1850            verts[seg.vertices[0]] = seg.vertices[0]
1851            verts[seg.vertices[1]] = seg.vertices[1]
1852        meshDict = self.Mesh2IOOutlineDict(userVertices=verts.values())
1853        load_mesh.loadASCII.export_mesh_file(ofile,meshDict)
1854       
1855        # exportASCIImeshfile   - this function is used
1856    def export_mesh_file(self,ofile):
1857        """
1858        export a file, ofile, with the format
1859        """
1860       
1861        dict = self.Mesh2IODict()
1862        load_mesh.loadASCII.export_mesh_file(ofile,dict)
1863
1864    def exportPointsFile(self,ofile):
1865        """
1866        export a points (.xya or .pts)  file, ofile.
1867       
1868        """
1869       
1870        mesh_dict = self.Mesh2IODict()
1871        point_dict = {}
1872        point_dict['attributelist'] = {} #this will need to be expanded..
1873                                         # if attributes are brought back in.
1874        point_dict['geo_reference'] = self.geo_reference
1875        if mesh_dict['vertices'] == []:
1876            point_dict['pointlist'] = mesh_dict['points']
1877        else:
1878            point_dict['pointlist'] = mesh_dict['vertices']
1879
1880        load_mesh.loadASCII.export_points_file(ofile,point_dict)
1881
1882       
1883########### IO CONVERTERS ##################
1884        """
1885        The dict fromat for IO with .tsh files is;
1886        (the triangulation)
1887        vertices: [[x1,y1],[x2,y2],...] (lists of doubles)
1888        vertex_attributes: [[a11,a12,...],[a21,a22],...] (lists of doubles)
1889        vertex_attribute_titles:[A1Title, A2Title ...] (A list of strings)
1890        segments: [[v1,v2],[v3,v4],...] (lists of integers)
1891        segment_tags : [tag,tag,...] list of strings
1892        triangles : [(v1,v2,v3), (v4,v5,v6),....] lists of points
1893        triangle tags: [s1,s2,...] A list of strings
1894        triangle neighbors: [[t1,t2,t3], [t4,t5,t6],..] lists of triangles
1895       
1896        (the outline)   
1897        points: [[x1,y1],[x2,y2],...] (lists of doubles)
1898        point_attributes: [[a11,a12,...],[a21,a22],...] (lists of doubles)
1899        outline_segments: [[point1,point2],[p3,p4],...] (lists of integers)
1900        outline_segment_tags : [tag1,tag2,...] list of strings
1901        holes : [[x1,y1],...](List of doubles, one inside each hole region)
1902        regions : [ [x1,y1],...] (List of 4 doubles)
1903        region_tags : [tag1,tag2,...] (list of strings)
1904        region_max_areas: [ma1,ma2,...] (A list of doubles)
1905        {Convension: A -ve max area means no max area}
1906       
1907        """
1908     
1909
1910                               
1911    def Mesh2IODict(self):
1912        """
1913        Convert the triangulation and outline info of a mesh to a dictionary
1914        structure
1915        """
1916        dict = self.Mesh2IOTriangulationDict()
1917        dict_mesh = self.Mesh2IOOutlineDict()
1918        for element in dict_mesh.keys():
1919            dict[element] = dict_mesh[element]
1920
1921        # add the geo reference
1922        dict['geo_reference'] = self.geo_reference
1923        return dict
1924   
1925    def Mesh2IOTriangulationDict(self):
1926        """
1927        Convert the Mesh to a dictionary of lists describing the
1928        triangulation variables;
1929       
1930        Used to produce .tsh file
1931        """
1932       
1933        meshDict = {}       
1934        vertices=[]
1935        vertex_attributes=[]
1936
1937        self.maxVertexIndex=0
1938        for vertex in self.meshVertices:
1939            vertex.index = self.maxVertexIndex
1940            vertices.append([vertex.x,vertex.y])
1941            vertex_attributes.append(vertex.attributes)           
1942            self.maxVertexIndex += 1
1943
1944        meshDict['vertices'] = vertices
1945        meshDict['vertex_attributes'] = vertex_attributes
1946        meshDict['vertex_attribute_titles'] = self.attributeTitles
1947        #segments
1948        segments=[]
1949        segment_tags=[]
1950        for seg in self.meshSegments:
1951            segments.append([seg.vertices[0].index,seg.vertices[1].index])
1952            segment_tags.append(seg.tag)
1953        meshDict['segments'] =segments
1954        meshDict['segment_tags'] =segment_tags
1955
1956        # Make sure that the indexation is correct
1957        index = 0
1958        for tri in self.meshTriangles:
1959            tri.index = index           
1960            index += 1
1961
1962        triangles = []
1963        triangle_tags = []
1964        triangle_neighbors = []
1965        for tri in self.meshTriangles: 
1966            triangles.append([tri.vertices[0].index,
1967                              tri.vertices[1].index,
1968                              tri.vertices[2].index]) 
1969            triangle_tags.append(tri.attribute)
1970            neighborlist = [-1,-1,-1]
1971            for neighbor,index in map(None,tri.neighbors,
1972                                      range(len(tri.neighbors))):
1973                if neighbor:   
1974                    neighborlist[index] = neighbor.index
1975            triangle_neighbors.append(neighborlist)
1976       
1977        meshDict['triangles'] = triangles
1978        meshDict['triangle_tags'] = triangle_tags
1979        meshDict['triangle_neighbors'] = triangle_neighbors
1980       
1981        #print "mesh.Mesh2IOTriangulationDict*)*)"
1982        #print meshDict
1983        #print "mesh.Mesh2IOTriangulationDict*)*)"
1984
1985        return meshDict
1986
1987                                                     
1988    def Mesh2IOOutlineDict(self, userVertices=None,
1989                        userSegments=None,
1990                        holes=None,
1991                        regions=None):
1992        """
1993        Convert the mesh outline to a dictionary of the lists needed for the
1994        triang module;
1995       
1996        Note, this adds an index attribute to the user Vertex objects.
1997
1998        Used to produce .tsh file and output to triangle
1999        """
2000        if userVertices is None:
2001            userVertices = self.getUserVertices()
2002        if userSegments is None:
2003            userSegments = self.getUserSegments()
2004        if holes is None:
2005            holes = self.getHoles()
2006        if regions is None:
2007            regions = self.getRegions()
2008           
2009        meshDict = {}
2010        #print "userVertices",userVertices
2011        #print "userSegments",userSegments
2012        pointlist=[]
2013        pointattributelist=[]
2014        index = 0
2015        for vertex in userVertices:
2016            vertex.index = index
2017            pointlist.append([vertex.x,vertex.y])
2018            pointattributelist.append(vertex.attributes)
2019           
2020            index += 1
2021        meshDict['points'] = pointlist
2022        meshDict['point_attributes'] = pointattributelist
2023
2024        segmentlist=[]
2025        segmenttaglist=[]
2026        for seg in userSegments:
2027            segmentlist.append([seg.vertices[0].index,seg.vertices[1].index])
2028            segmenttaglist.append(seg.tag)
2029        meshDict['outline_segments'] =segmentlist
2030        meshDict['outline_segment_tags'] =segmenttaglist
2031       
2032        holelist=[]
2033        for hole in holes:
2034            holelist.append([hole.x,hole.y]) 
2035        meshDict['holes'] = holelist
2036       
2037        regionlist=[]
2038        regiontaglist = []
2039        regionmaxarealist = []
2040        for region in regions:
2041            regionlist.append([region.x,region.y])
2042            regiontaglist.append(region.getTag())
2043           
2044            if (region.getMaxArea() != None): 
2045                regionmaxarealist.append(region.getMaxArea())
2046            else:
2047                regionmaxarealist.append( load_mesh.loadASCII.NOMAXAREA)
2048        meshDict['regions'] = regionlist
2049        meshDict['region_tags'] = regiontaglist
2050        meshDict['region_max_areas'] = regionmaxarealist
2051        #print "*(*("
2052        #print meshDict
2053        #print meshDict['regionlist']
2054        #print "*(*("
2055        return meshDict
2056
2057    def IOTriangulation2Mesh(self, genDict):
2058        """
2059        Set the mesh attributes given an tsh IO dictionary
2060        """
2061        #Clear the current generated mesh values
2062        self.meshTriangles=[]
2063        self.attributeTitles=[]
2064        self.meshSegments=[]
2065        self.meshVertices=[]
2066
2067        #print "mesh.setTriangulation@#@#@#"
2068        #print genDict
2069        #print "@#@#@#"
2070
2071        self.maxVertexIndex = 0
2072        for point in genDict['vertices']:
2073            v=Vertex(point[0], point[1])
2074            v.index =  self.maxVertexIndex
2075            self.maxVertexIndex +=1
2076            self.meshVertices.append(v)
2077
2078        self.attributeTitles = genDict['vertex_attribute_titles']
2079
2080        index = 0
2081        for seg,tag in map(None,genDict['segments'],genDict['segment_tags']):
2082            segObject = Segment( self.meshVertices[seg[0]],
2083                           self.meshVertices[seg[1]], tag = tag )
2084            segObject.index = index
2085            index +=1
2086            self.meshSegments.append(segObject)
2087
2088        index = 0
2089        for triangle in genDict['triangles']:
2090            tObject =Triangle( self.meshVertices[triangle[0]],
2091                        self.meshVertices[triangle[1]],
2092                        self.meshVertices[triangle[2]] )
2093            tObject.index = index
2094            index +=1
2095            self.meshTriangles.append(tObject)
2096
2097        index = 0
2098        for att in genDict['triangle_tags']:
2099            if att == []:
2100                self.meshTriangles[index].setAttribute("")
2101            else:
2102                self.meshTriangles[index].setAttribute(att)
2103            index += 1
2104           
2105        index = 0
2106        for att in genDict['vertex_attributes']:
2107            if att == None:
2108                self.meshVertices[index].setAttributes([])
2109            else:
2110                self.meshVertices[index].setAttributes(att)
2111            index += 1
2112   
2113        index = 0
2114        for triangle in genDict['triangle_neighbors']:
2115            # Build a list of triangle object neighbors
2116            ObjectNeighbor = []
2117            for neighbor in triangle:
2118                if ( neighbor != -1):
2119                    ObjectNeighbor.append(self.meshTriangles[neighbor])
2120                else:
2121                    ObjectNeighbor.append(None)
2122            self.meshTriangles[index].setNeighbors(ObjectNeighbor[0],
2123                                                   ObjectNeighbor[1],
2124                                                   ObjectNeighbor[2])
2125            index += 1
2126
2127
2128    def IOOutline2Mesh(self, genDict):
2129        """
2130        Set the outline (user Mesh attributes) given a IO tsh dictionary
2131       
2132        mesh is an instance of a mesh object
2133        """
2134        #Clear the current user mesh values
2135        self.clearUserSegments()
2136        self.userVertices=[]
2137        self.Holes=[]
2138        self.Regions=[]
2139
2140        #print "mesh.IOOutline2Mesh@#@#@#"
2141        #print "genDict",genDict
2142        #print "@#@#@#"
2143       
2144        #index = 0
2145        for point in genDict['points']:
2146            v=Vertex(point[0], point[1])
2147            #v.index = index
2148            #index +=1
2149            self.userVertices.append(v)
2150
2151        #index = 0
2152        for seg,tag in map(None,genDict['outline_segments'],
2153                           genDict['outline_segment_tags']):
2154            segObject = Segment( self.userVertices[seg[0]],
2155                           self.userVertices[seg[1]], tag = tag )
2156            #segObject.index = index
2157            #index +=1
2158            self.userSegments.append(segObject)
2159
2160# Remove the loading of attribute info.
2161# Have attribute info added using least_squares in pyvolution
2162#         index = 0
2163#         for att in genDict['point_attributes']:
2164#             if att == None:
2165#                 self.userVertices[index].setAttributes([])
2166#             else:
2167#                 self.userVertices[index].setAttributes(att)
2168#            index += 1
2169       
2170        #index = 0
2171        for point in genDict['holes']:
2172            h=Hole(point[0], point[1])
2173            #h.index = index
2174            #index +=1
2175            self.holes.append(h)
2176
2177        #index = 0
2178        for reg,att,maxArea in map(None,
2179                                   genDict['regions'],
2180                                   genDict['region_tags'],
2181                                   genDict['region_max_areas']):
2182            if maxArea > 0:  # maybe I should ref NOMAXAREA? Prob' not though
2183                Object = Region( reg[0],
2184                                 reg[1],
2185                                 tag = att,
2186                                 maxArea = maxArea)
2187            else:
2188                Object = Region( reg[0],
2189                                 reg[1],
2190                                 tag = att)
2191               
2192            #Object.index = index
2193            #index +=1
2194            self.regions.append(Object)
2195 
2196############################################
2197
2198       
2199    def refineSet(self,setName):
2200        Triangles = self.sets[self.setID[setName]]
2201        Refine(self,Triangles)
2202
2203    def selectAllTriangles(self):
2204        A=[]
2205        A.extend(self.meshTriangles)
2206        if not('All' in self.setID.keys()):
2207            self.setID['All']=len(self.sets)
2208            self.sets.append(A)
2209        else:
2210            self.sets[self.setID['All']]=A
2211        return 'All'
2212        # and objectIDs
2213
2214
2215    def clearSelection(self):
2216        A = []
2217        if not('None' in self.setID.keys()):
2218            self.setID['None']=len(self.sets)
2219            self.sets.append(A)
2220        return 'None'
2221
2222    def drawSet(self,canvas,setName,SCALE,colour=SET_COLOUR):
2223    #FIXME Draws over previous triangles - may bloat canvas
2224        Triangles = self.sets[self.setID[setName]]
2225        for triangle in Triangles:
2226            triangle.draw(canvas,1,
2227                          scale = SCALE,
2228                          colour = colour)
2229           
2230    def undrawSet(self,canvas,setName,SCALE,colour='green'):
2231    #FIXME Draws over previous lines - may bloat canvas
2232        Triangles = self.sets[self.setID[setName]]
2233        for triangle in Triangles:
2234            triangle.draw(canvas,1,
2235                          scale = SCALE,
2236                          colour = colour)
2237
2238    def weed(self,Vertices,Segments):
2239        #Depreciated
2240        #weed out existing duplicates
2241        print 'len(self.getUserSegments())'
2242        print len(self.getUserSegments())
2243        print 'len(self.getUserVertices())'
2244        print len(self.getUserVertices())
2245
2246        point_keys = {}
2247        for vertex in Vertices:
2248            point = (vertex.x,vertex.y)
2249            point_keys[point]=vertex
2250        #inlined would looks very ugly
2251
2252        line_keys = {}
2253        for segment in Segments:
2254            vertex1 = segment.vertices[0]
2255            vertex2 = segment.vertices[1]
2256            point1 = (vertex1.x,vertex1.y)
2257            point2 = (vertex2.x,vertex2.y)
2258            segment.vertices[0]=point_keys[point1]
2259            segment.vertices[1]=point_keys[point2]
2260            vertex1 = segment.vertices[0]
2261            vertex2 = segment.vertices[1]
2262            point1 = (vertex1.x,vertex1.y)
2263            point2 = (vertex2.x,vertex2.y)
2264            line1 = (point1,point2)
2265            line2 = (point2,point1)
2266            if not (line_keys.has_key(line1) \
2267                 or line_keys.has_key(line2)):
2268                 line_keys[line1]=segment
2269        Vertices=point_keys.values()
2270        Segments=line_keys.values()
2271        return Vertices,Segments
2272
2273    def segs_to_dict(self,segments):
2274        dict={}
2275        for segment in segments:
2276            vertex1 = segment.vertices[0]
2277            vertex2 = segment.vertices[1]
2278            point1 = (vertex1.x,vertex1.y)
2279            point2 = (vertex2.x,vertex2.y)
2280            line = (point1,point2)
2281            dict[line]=segment
2282        return dict
2283
2284    def seg2line(self,s):
2285        return ((s.vertices[0].x,s.vertices[0].y,)\
2286                (s.vertices[1].x,s.vertices[1].y))
2287
2288    def line2seg(self,line,tag=None):
2289        point0 = self.point2ver(line[0])
2290        point1 = self.point2ver(line[1])
2291        return Segment(point0,point1,tag=tag)
2292
2293    def ver2point(self,vertex):
2294        return (vertex.x,vertex.y)
2295
2296    def point2ver(self,point):
2297        return Vertex(point[0],point[1])
2298
2299    def smooth_polySet(self,min_radius=0.05):
2300        #for all pairs of connecting segments:
2301        #    propose a new segment that replaces the 2
2302
2303        #    If the difference between the new segment
2304        #    and the old lines is small: replace the
2305        #    old lines.
2306
2307        seg2line = self.seg2line
2308        ver2point= self.ver2point
2309        line2seg = self.line2seg
2310        point2ver= self.point2ver
2311
2312        #create dictionaries of lines -> segments
2313        userSegments = self.segs_to_dict(self.userSegments)
2314        alphaSegments = self.segs_to_dict(self.alphaUserSegments)
2315
2316        #lump user and alpha segments
2317        for key in alphaSegments.keys():
2318            userSegments[key]=alphaSegments[key]
2319
2320        #point_keys = tuple -> vertex
2321        #userVertices = vertex -> [line,line] - lines from that node
2322        point_keys = {}
2323        userVertices={}
2324        for vertex in self.getUserVertices():
2325            point = ver2point(vertex)
2326            if not point_keys.has_key(point):
2327                point_keys[point]=vertex
2328                userVertices[vertex]=[]
2329        for key in userSegments.keys():
2330            line = key
2331            point_0 = key[0]
2332            point_1 = key[1]
2333            userVertices[point_keys[point_0]].append(line)
2334            userVertices[point_keys[point_1]].append(line)
2335
2336        for point in point_keys.keys():
2337            try:
2338            #removed keys can cause keyerrors
2339                vertex = point_keys[point]
2340                lines = userVertices[vertex]
2341   
2342                #if there are 2 lines on the node
2343                if len(lines)==2:
2344                    line_0 = lines[0]
2345                    line_1 = lines[1]
2346   
2347                    #if the tags are the the same on the 2 lines
2348                    if userSegments[line_0].tag == userSegments[line_1].tag:
2349                        tag = userSegments[line_0].tag
2350   
2351                        #point_a is one of the next nodes, point_b is the other
2352                        if point==line_0[0]:
2353                            point_a = line_0[1]
2354                        if point==line_0[1]:
2355                            point_a = line_0[0]
2356                        if point==line_1[0]:
2357                            point_b = line_1[1]
2358                        if point==line_1[1]:
2359                            point_b = line_1[0]
2360   
2361   
2362                        #line_2 is proposed
2363                        line_2 = (point_a,point_b)
2364
2365                        #calculate the area of the triangle between
2366                        #the two existing segments and the proposed
2367                        #new segment
2368                        ax = point_a[0]
2369                        ay = point_a[1]
2370                        bx = point_b[0]
2371                        by = point_b[1]
2372                        cx = point[0]
2373                        cy = point[1]
2374                        area=abs((bx*ay-ax*by)+(cx*by-bx*cy)+(ax*cy-cx*ay))/2
2375
2376                        #calculate the perimeter
2377                        len_a =  ((cx-bx)**2+(cy-by)**2)**0.5 
2378                        len_b =  ((ax-cx)**2+(ay-cy)**2)**0.5 
2379                        len_c =  ((bx-ax)**2+(by-ay)**2)**0.5 
2380                        perimeter = len_a+len_b+len_c
2381
2382                        #calculate the radius
2383                        r = area/(2*perimeter)
2384
2385                        #if the radius is small: then replace the existing
2386                        #segments with the new one
2387                        if r < min_radius:
2388                            if len_c < min_radius: append = False
2389                            else: append = True
2390                            #if the new seg is also time, don't add it
2391                            if append:
2392                                segment = self.line2seg(line_2,tag=tag)
2393
2394                            list_a=userVertices[point_keys[point_a]]
2395                            list_b=userVertices[point_keys[point_b]]
2396
2397                            if line_0 in list_a:
2398                                list_a.remove(line_0)
2399                            else:
2400                                list_a.remove(line_1)
2401
2402                            if line_0 in list_b:
2403                                list_b.remove(line_0)
2404                            else:
2405                                list_b.remove(line_1)
2406
2407                            if append:
2408                                list_a.append(line_2)
2409                                list_b.append(line_2)
2410                            else:
2411                                if len(list_a)==0:
2412                                    userVertices.pop(point_keys[point_a])
2413                                    point_keys.pop(point_a)
2414                                if len(list_b)==0:
2415                                    userVertices.pop(point_keys[point_b])
2416                                    point_keys.pop(point_b)
2417
2418                            userVertices.pop(point_keys[point])
2419                            point_keys.pop(point)
2420                            userSegments.pop(line_0)
2421                            userSegments.pop(line_1)
2422
2423                            if append:
2424                                userSegments[line_2]=segment
2425            except:
2426                pass
2427
2428        #self.userVerticies = userVertices.keys()
2429        #self.userSegments = []
2430        #for key in userSegments.keys():
2431        #    self.userSegments.append(userSegments[key])
2432        #self.alphaUserSegments = []
2433
2434        self.userVerticies = []
2435        self.userSegments = []
2436        self.alphaUserSegments = []
2437
2438        return userVertices,userSegments,alphaSegments
2439
2440    def triangles_to_polySet(self,setName):
2441        #self.smooth_polySet()
2442
2443        seg2line = self.seg2line
2444        ver2point= self.ver2point
2445        line2seg = self.line2seg
2446        point2ver= self.point2ver
2447
2448        from Numeric import array,allclose
2449        #turn the triangles into a set
2450        Triangles = self.sets[self.setID[setName]]
2451        Triangles_dict = {}
2452        for triangle in Triangles:
2453            Triangles_dict[triangle]=None
2454 
2455
2456        #create a dict of points to vertexes (tuple -> object)
2457        #also create a set of vertexes (object -> True)
2458        point_keys = {}
2459        userVertices={}
2460        for vertex in self.getUserVertices():
2461            point = ver2point(vertex)
2462            if not point_keys.has_key(point):
2463                point_keys[point]=vertex
2464                userVertices[vertex]=True
2465
2466        #create a dict of lines to segments (tuple -> object)
2467        userSegments = self.segs_to_dict(self.userSegments)
2468        #append the userlines in an affine linespace
2469        affine_lines = Affine_Linespace()
2470        for line in userSegments.keys():
2471            affine_lines.append(line)
2472        alphaSegments = self.segs_to_dict(self.alphaUserSegments)
2473        for line in alphaSegments.keys():
2474            affine_lines.append(line)
2475       
2476        for triangle in Triangles:
2477            for i in (0,1,2):
2478                #for every triangles neighbour:
2479                if not Triangles_dict.has_key(triangle.neighbors[i]):
2480                #if the neighbour is not in the set:
2481                    a = triangle.vertices[i-1]
2482                    b = triangle.vertices[i-2]
2483                    #Get possible matches:
2484                    point_a = ver2point(a)
2485                    point_b = ver2point(b)
2486                    midpoint = ((a.x+b.x)/2,(a.y+b.y)/2)
2487                    line = (point_a,point_b)
2488                    tag = None
2489
2490
2491                    #this bit checks for matching lines
2492                    possible_lines = affine_lines[line] 
2493                    possible_lines = unique(possible_lines)
2494                    found = 0                           
2495                    for user_line in possible_lines:
2496                        if self.point_on_line(midpoint,user_line):
2497                            found+=1
2498                            assert found<2
2499                            if userSegments.has_key(user_line):
2500                                parent_segment = userSegments.pop(user_line)
2501                            if alphaSegments.has_key(user_line):
2502                                parent_segment = alphaSegments.pop(user_line)
2503                            tag = parent_segment.tag
2504                            offspring = [line]
2505                            offspring.extend(self.subtract_line(user_line,line))
2506                            affine_lines.remove(user_line)
2507                            for newline in offspring:
2508                                line_vertices = []
2509                                for point in newline:
2510                                    if point_keys.has_key(point):
2511                                        vert = point_keys[point]
2512                                    else:
2513                                        vert = Vertex(point[0],point[1])
2514                                        userVertices[vert]=True
2515                                        point_keys[point]=vert
2516                                    line_vertices.append(vert)
2517                                segment = Segment(line_vertices[0],line_vertices[1],tag)
2518                                userSegments[newline]=segment
2519                                affine_lines.append(newline)
2520                            #break
2521                    assert found<2
2522
2523
2524
2525                    #if no matching lines
2526                    if not found:
2527                        line_vertices = []
2528                        for point in line:
2529                            if point_keys.has_key(point):
2530                                vert = point_keys[point]
2531                            else:
2532                                vert = Vertex(point[0],point[1])
2533                                userVertices[vert]=True
2534                                point_keys[point]=vert
2535                            line_vertices.append(vert)
2536                        segment = Segment(line_vertices[0],line_vertices[1],tag)
2537                        userSegments[line]=segment
2538                        affine_lines.append(line)
2539       
2540        self.userVerticies = []
2541        self.userSegments = []
2542        self.alphaUserSegments = []
2543
2544        return userVertices,userSegments,alphaSegments
2545
2546    def subtract_line(self,parent,child):
2547        #Subtracts child from parent
2548        #Requires that the child is a
2549        #subline of parent to work.
2550
2551        from Numeric import allclose,dot,array
2552        A= parent[0]
2553        B= parent[1]
2554        a = child[0]
2555        b = child[1]
2556
2557        A_array = array(parent[0])
2558        B_array = array(parent[1])
2559        a_array   = array(child[0])
2560        b_array   = array(child[1])
2561
2562        assert not A == B
2563        assert not a == b
2564
2565        answer = []
2566
2567        #if the new line does not share a
2568        #vertex with the old one
2569        if not (allclose(A_array,a_array)\
2570             or allclose(B_array,b_array)\
2571             or allclose(A_array,b_array)\
2572             or allclose(a_array,B_array)):
2573            if dot(A_array-a_array,A_array-a_array) \
2574            < dot(A_array-b_array,A_array-b_array):
2575                sibling1 = (A,a)
2576                sibling2 = (B,b)
2577                return [sibling1,sibling2]
2578            else:
2579                sibling1 = (A,b)
2580                sibling2 = (B,a)
2581                return [sibling1,sibling2]
2582
2583        elif allclose(A_array,a_array):
2584            if allclose(B_array,b_array):
2585                return []
2586            else:
2587                sibling = (b,B)
2588                return [sibling]
2589        elif allclose(B_array,b_array):
2590            sibling = (a,A)
2591            return [sibling]
2592
2593        elif allclose(A_array,b_array):
2594            if allclose(B,a):
2595                return []
2596            else:
2597                sibling = (a,B)
2598                return [sibling]
2599        elif allclose(a_array,B_array):
2600            sibling = (b,A)
2601            return [sibling]
2602
2603    def point_on_line(self,point,line):
2604        #returns true within a tolerance of 3 degrees
2605        x=point[0]
2606        y=point[1]
2607        x0=line[0][0]
2608        x1=line[1][0]
2609        y0=line[0][1]
2610        y1=line[1][1]
2611        from Numeric import array, dot, allclose
2612        from math import sqrt
2613        tol = 3. #DEGREES
2614        tol = tol*3.1415/180
2615
2616        a = array([x - x0, y - y0]) 
2617        a_normal = array([a[1], -a[0]])     
2618        len_a_normal = sqrt(sum(a_normal**2)) 
2619
2620        b = array([x1 - x0, y1 - y0])         
2621        len_b = sqrt(sum(b**2)) 
2622   
2623        if abs(dot(a_normal, b)/(len_b*len_a_normal))< tol:
2624            #Point is somewhere on the infinite extension of the line
2625
2626            len_a = sqrt(sum(a**2))                                     
2627            if dot(a, b) >= 0 and len_a <= len_b:
2628               return True
2629            else:   
2630               return False
2631        else:
2632          return False
2633
2634    def line_length(self,line):
2635        x0=line[0][0]
2636        x1=line[1][0]
2637        y0=line[0][1]
2638        y1=line[1][1]
2639        return ((x1-x0)**2-(y1-y0)**2)**0.5     
2640
2641    def threshold(self,setName,min=None,max=None,attribute_name = 'elevation'):
2642        """
2643        threshold using  d
2644        """
2645        triangles = self.sets[self.setID[setName]]
2646        A = []
2647
2648        if attribute_name in self.attributeTitles:
2649            i = self.attributeTitles.index(attribute_name)
2650        else: i = -1#no attribute
2651        if not max == None:
2652            for t in triangles:
2653                if (min<self.av_att(t,i)<max):
2654                    A.append(t)
2655        else:
2656            for t in triangles:
2657                if (min<self.av_att(t,i)):
2658                    A.append(t)
2659        self.sets[self.setID[setName]] = A
2660
2661    def general_threshold(self,setName,min=None,max=None\
2662              ,attribute_name = 'elevation',function=None):
2663        """
2664        Thresholds the triangles
2665        """
2666        from visual.graph import arange,ghistogram,color as colour
2667        triangles = self.sets[self.setID[setName]]
2668        A = []
2669        data=[]
2670        #data is for the graph
2671
2672        if attribute_name in self.attributeTitles:
2673            i = self.attributeTitles.index(attribute_name)
2674        else: i = -1
2675        if not max == None:
2676            for t in triangles:
2677                value=function(t,i)
2678                if (min<value<max):
2679                    A.append(t)
2680                data.append(value)
2681        else:
2682            for t in triangles:
2683                value=function(t,i)
2684                if (min<value):
2685                    A.append(t)
2686                data.append(value)
2687        self.sets[self.setID[setName]] = A
2688
2689        if self.visualise_graph:
2690            if len(data)>0:
2691                max=data[0]
2692                min=data[0]
2693                for value in data:
2694                    if value > max:
2695                        max = value
2696                    if value < min:
2697                        min = value
2698
2699                inc = (max-min)/100
2700
2701                histogram = ghistogram(bins=arange(min,max,inc),\
2702                             color = colour.red)
2703                histogram.plot(data=data)
2704       
2705    def av_att(self,triangle,i):
2706        if i==-1: return 1
2707        else:
2708            #evaluates the average attribute of the vertices of a triangle.
2709            V = triangle.getVertices()
2710            a0 = (V[0].attributes[i])
2711            a1 = (V[1].attributes[i])
2712            a2 = (V[2].attributes[i])
2713            return (a0+a1+a2)/3
2714
2715    def Courant_ratio(self,triangle,index):
2716        """
2717        Uses the courant threshold
2718        """
2719        e = self.av_att(triangle,index)
2720        A = triangle.calcArea()
2721        P = triangle.calcP()
2722        r = A/(2*P)
2723        e = max(0.1,abs(e))
2724        return r/e**0.5
2725
2726    def Gradient(self,triangle,index):
2727        V = triangle.vertices
2728        x0, y0, x1, y1, x2, y2, q0, q1, q2 = V[0].x,V[0].y,V[1].x,V[1].y,V[2].x,V[2].y,V[0].attributes[index],V[1].attributes[index],V[2].attributes[index]
2729        grad_x,grad_y = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
2730        if ((grad_x**2)+(grad_y**2))**(0.5)<0:
2731            print ((grad_x**2)+(grad_y**2))**(0.5)
2732        return ((grad_x**2)+(grad_y**2))**(0.5)
2733   
2734
2735    def append_triangle(self,triangle):
2736        self.meshTriangles.append(triangle)
2737
2738    def replace_triangle(self,triangle,replacement):
2739        i = self.meshTriangles.index(triangle)
2740        self.meshTriangles[i]=replacement
2741        assert replacement in self.meshTriangles
2742
2743def importUngenerateFile(ofile):
2744    """
2745    import a file, ofile, with the format
2746    [poly]
2747    poly format:
2748    First line:  <# of vertices> <x centroid> <y centroid>
2749    Following lines: <x> <y>
2750    last line:  "END"
2751
2752    Note: These are clockwise.
2753    """
2754    fd = open(ofile,'r')
2755    Dict = readUngenerateFile(fd)
2756    fd.close()
2757    return Dict
2758
2759def readUngenerateFile(fd):
2760    """
2761    import a file, ofile, with the format
2762    [poly]
2763    poly format:
2764    First line:  <# of polynomial> <x centroid> <y centroid>
2765    Following lines: <x> <y>
2766    last line:  "END"
2767    """
2768    END_DELIMITER = 'END\n'
2769   
2770    points = []
2771    segments = []
2772   
2773    isEnd = False
2774    line = fd.readline() #not used <# of polynomial> <x> <y>
2775    while not isEnd:
2776        line = fd.readline()
2777        fragments = line.split()
2778        vert = [float(fragments.pop(0)),float(fragments.pop(0))]
2779        points.append(vert)
2780        PreviousVertIndex = len(points)-1
2781        firstVertIndex = PreviousVertIndex
2782       
2783        line = fd.readline() #Read the next line
2784        while line <> END_DELIMITER: 
2785            #print "line >" + line + "<"
2786            fragments = line.split()
2787            vert = [float(fragments.pop(0)),float(fragments.pop(0))]
2788            points.append(vert)
2789            thisVertIndex = len(points)-1
2790            segment = [PreviousVertIndex,thisVertIndex]
2791            segments.append(segment)
2792            PreviousVertIndex = thisVertIndex
2793            line = fd.readline() #Read the next line
2794            i =+ 1
2795        # If the last and first segments are the same,
2796        # Remove the last segment and the last vertex
2797        # then add a segment from the second last vert to the 1st vert
2798        thisVertIndex = len(points)-1
2799        firstVert = points[firstVertIndex]
2800        thisVert = points[thisVertIndex]
2801        #print "firstVert",firstVert
2802        #print "thisVert",thisVert
2803        if (firstVert[0] == thisVert[0] and firstVert[1] == thisVert[1]):
2804            points.pop()
2805            segments.pop()
2806            thisVertIndex = len(points)-1
2807            segments.append([thisVertIndex, firstVertIndex])
2808       
2809        line = fd.readline() # read <# of polynomial> <x> <y> OR END
2810        #print "line >>" + line + "<<"
2811        if line == END_DELIMITER:
2812            isEnd = True
2813   
2814    #print "points", points       
2815    #print "segments", segments
2816    ungenerated_dict = {}
2817    ungenerated_dict['points'] = points
2818    ungenerated_dict['segments'] = segments
2819    return ungenerated_dict
2820
2821def importMeshFromFile(ofile):
2822    """returns a mesh object, made from a .xya/.pts or .tsh/.msh file
2823    Often raises IOError,RuntimeError
2824    """
2825    newmesh = None
2826    if (ofile[-4:]== ".xya" or ofile[-4:]== ".pts"):
2827        dict = load_mesh.loadASCII.import_points_file(ofile)
2828        dict['points'] = dict['pointlist']
2829        dict['outline_segments'] = []
2830        dict['outline_segment_tags'] = []
2831        dict['regions'] = []
2832        dict['region_tags'] = []
2833        dict['region_max_areas'] = []
2834        dict['holes'] = [] 
2835        newmesh= Mesh(geo_reference = dict['geo_reference'])
2836        newmesh.IOOutline2Mesh(dict)
2837        counter = newmesh.removeDuplicatedUserVertices()
2838        if (counter >0):
2839            print "%i duplicate vertices removed from dataset" % (counter)
2840    elif (ofile[-4:]== ".tsh" or ofile[-4:]== ".msh"):
2841        dict = load_mesh.loadASCII.import_mesh_file(ofile)
2842        #print "********"
2843        #print "zq mesh.dict",dict
2844        #print "********"
2845        newmesh= Mesh()
2846        newmesh.IOOutline2Mesh(dict)
2847        newmesh.IOTriangulation2Mesh(dict)
2848    else:
2849        raise RuntimeError
2850   
2851    if dict.has_key('geo_reference') and not dict['geo_reference'] == None:
2852        newmesh.geo_reference = dict['geo_reference']
2853    return newmesh
2854
2855def loadPickle(currentName):
2856    fd = open(currentName)
2857    mesh = pickle.load(fd)
2858    fd.close()
2859    return mesh
2860   
2861def square_outline(side_length = 1,up = "top", left = "left", right = "right",
2862                   down = "bottom", regions = False):
2863   
2864        a = Vertex (0,0)
2865        b = Vertex (0,side_length)
2866        c = Vertex (side_length,0)
2867        d = Vertex (side_length,side_length)
2868     
2869        s2 = Segment(b,d, tag = up)
2870        s3 = Segment(b,a, tag = left)
2871        s4 = Segment(d,c, tag = right)
2872        s5 = Segment(a,c, tag = down)
2873
2874        if regions:
2875            e =  Vertex (side_length/2,side_length/2)
2876            s6 = Segment(a,e, tag = down + left)
2877            s7 = Segment(b,e, tag = up + left)
2878            s8 = Segment(c,e, tag = down + right)
2879            s9 = Segment(d,e, tag = up + right)
2880            r1 = Region(side_length/2,3.*side_length/4, tag = up)
2881            r2 = Region(1.*side_length/4,side_length/2, tag = left)
2882            r3 = Region(3.*side_length/4,side_length/2, tag = right)
2883            r4 = Region(side_length/2,1.*side_length/4, tag = down)
2884            mesh = Mesh(userVertices=[a,b,c,d,e],
2885                        userSegments=[s2,s3,s4,s5,s6,s7,s8,s9],
2886                        regions = [r1,r2,r3,r4])
2887        else:
2888            mesh = Mesh(userVertices=[a,b,c,d],
2889                 userSegments=[s2,s3,s4,s5])
2890     
2891        return mesh
2892
2893   
2894
2895def region_strings2ints(region_list):
2896    """Given a list of (x_int,y_int,tag_string) lists it returns a list of
2897    (x_int,y_int,tag_int) and a list to convert the tag_int's back to
2898    the tag_strings 
2899    """
2900    # Make sure "" has an index of 0
2901    region_list.reverse()
2902    region_list.append((1.0,2.0,""))
2903    region_list.reverse()
2904    convertint2string = []
2905    for i in xrange(len(region_list)):
2906        convertint2string.append(region_list[i][2])
2907        if len(region_list[i]) == 4: # there's an area value
2908            region_list[i] = (region_list[i][0],
2909                              region_list[i][1],i,region_list[i][3])
2910        elif len(region_list[i]) == 3: # no area value
2911            region_list[i] = (region_list[i][0],region_list[i][1],i)
2912        else:
2913            print "The region list has a bad size"
2914            # raise an error ..
2915            raise Error
2916
2917    #remove "" from the region_list
2918    region_list.pop(0)
2919   
2920    return [region_list, convertint2string]
2921
2922def region_ints2strings(region_list,convertint2string):
2923    """Reverses the transformation of region_strings2ints
2924    """
2925    if region_list[0] != []:
2926        for i in xrange(len(region_list)):
2927            region_list[i] = [convertint2string[int(region_list[i][0])]]
2928    return region_list
2929
2930def segment_ints2strings(intlist, convertint2string):
2931    """Reverses the transformation of segment_strings2ints """
2932    stringlist = []
2933    for x in intlist:
2934        stringlist.append(convertint2string[x])
2935    return stringlist
2936
2937def segment_strings2ints(stringlist, preset):
2938    """Given a list of strings return a list of 0 to n ints which represent
2939    the strings and a converting list of the strings, indexed by 0 to n ints.
2940    Also, input an initial converting list of the strings
2941    Note, the converting list starts off with
2942    ["internal boundary", "external boundary", "internal boundary"]
2943    example input and output
2944    input ["a","b","a","c"],["c"]
2945    output [[2, 1, 2, 0], ['c', 'b', 'a']]
2946
2947    the first element in the converting list will be
2948    overwritten with "".
2949    ?This will become the third item in the converting list?
2950   
2951    # note, order the initial converting list is important,
2952     since the index = the int tag
2953    """
2954    nodups = unique(stringlist)
2955    # note, order is important, the index = the int tag
2956    #preset = ["internal boundary", "external boundary"]
2957    #Remove the preset tags from the list with no duplicates
2958    nodups = [x for x in nodups if x not in preset]
2959
2960    try:
2961        nodups.remove("") # this has to go to zero
2962    except ValueError:
2963        pass
2964   
2965    # Add the preset list at the beginning of no duplicates
2966    preset.reverse()
2967    nodups.extend(preset)
2968    nodups.reverse()
2969
2970       
2971    convertstring2int = {}
2972    convertint2string = []
2973    index = 0
2974    for x in nodups:
2975        convertstring2int[x] = index
2976        convertint2string.append(x)
2977        index += 1
2978    convertstring2int[""] = 0
2979
2980    intlist = []
2981    for x in stringlist:
2982        intlist.append(convertstring2int[x])
2983    return [intlist, convertint2string]
2984
2985
2986
2987
2988def unique(s):
2989    """Return a list of the elements in s, but without duplicates.
2990
2991    For example, unique([1,2,3,1,2,3]) is some permutation of [1,2,3],
2992    unique("abcabc") some permutation of ["a", "b", "c"], and
2993    unique(([1, 2], [2, 3], [1, 2])) some permutation of
2994    [[2, 3], [1, 2]].
2995
2996    For best speed, all sequence elements should be hashable.  Then
2997    unique() will usually work in linear time.
2998
2999    If not possible, the sequence elements should enjoy a total
3000    ordering, and if list(s).sort() doesn't raise TypeError it's
3001    assumed that they do enjoy a total ordering.  Then unique() will
3002    usually work in O(N*log2(N)) time.
3003
3004    If that's not possible either, the sequence elements must support
3005    equality-testing.  Then unique() will usually work in quadratic
3006    time.
3007    """
3008
3009    n = len(s)
3010    if n == 0:
3011        return []
3012
3013    # Try using a dict first, as that's the fastest and will usually
3014    # work.  If it doesn't work, it will usually fail quickly, so it
3015    # usually doesn't cost much to *try* it.  It requires that all the
3016    # sequence elements be hashable, and support equality comparison.
3017    u = {}
3018    try:
3019        for x in s:
3020            u[x] = 1
3021    except TypeError:
3022        del u  # move on to the next method
3023    else:
3024        return u.keys()
3025
3026    # We can't hash all the elements.  Second fastest is to sort,
3027    # which brings the equal elements together; then duplicates are
3028    # easy to weed out in a single pass.
3029    # NOTE:  Python's list.sort() was designed to be efficient in the
3030    # presence of many duplicate elements.  This isn't true of all
3031    # sort functions in all languages or libraries, so this approach
3032    # is more effective in Python than it may be elsewhere.
3033    try:
3034        t = list(s)
3035        t.sort()
3036    except TypeError:
3037        del t  # move on to the next method
3038    else:
3039        assert n > 0
3040        last = t[0]
3041        lasti = i = 1
3042        while i < n:
3043            if t[i] != last:
3044                t[lasti] = last = t[i]
3045                lasti += 1
3046            i += 1
3047        return t[:lasti]
3048
3049    # Brute force is all that's left.
3050    u = []
3051    for x in s:
3052        if x not in u:
3053            u.append(x)
3054    return u
3055
3056"""Refines triangles
3057
3058   Implements the #triangular bisection?# algorithm.
3059 
3060   
3061"""
3062
3063def Refine(mesh, triangles):
3064    """
3065    Given a general_mesh, and a triangle number, split
3066    that triangle in the mesh in half. Then to prevent
3067    vertices and edges from meeting, keep refining
3068    neighbouring triangles until the mesh is clean.
3069    """
3070    state = BisectionState(mesh)
3071    for triangle in triangles:
3072        if not state.refined_triangles.has_key(triangle):
3073            triangle.rotate_longest_side()
3074            state.start(triangle)
3075            Refine_mesh(mesh, state)
3076
3077def Refine_mesh(mesh, state):
3078    """
3079    """
3080    state.getState(mesh)
3081    refine_triangle(mesh,state)
3082    state.evolve()
3083    if not state.end:
3084        Refine_mesh(mesh,state)
3085
3086def refine_triangle(mesh,state):
3087    split(mesh,state.current_triangle,state.new_point)
3088    if state.case == 'one':
3089        state.r[3]=state.current_triangle#triangle 2
3090
3091        new_triangle_id = len(mesh.meshTriangles)-1
3092        new_triangle = mesh.meshTriangles[new_triangle_id]
3093
3094        split(mesh,new_triangle,state.old_point)
3095        state.r[2]=new_triangle#triangle 1.2
3096        state.r[4]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 1.1
3097        r = state.r
3098        state.repairCaseOne()
3099
3100    if state.case == 'two':
3101        state.r[2]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 1
3102
3103        new_triangle = state.current_triangle
3104
3105        split(mesh,new_triangle,state.old_point)
3106
3107        state.r[3]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 2.1
3108        state.r[4]=new_triangle#triangle 2.2
3109        r = state.r
3110
3111        state.repairCaseTwo()
3112
3113    if state.case == 'vertex':
3114        state.r[2]=state.current_triangle#triangle 2
3115        state.r[3]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 1
3116        r = state.r
3117        state.repairCaseVertex()
3118       
3119    if state.case == 'start':
3120        state.r[2]=mesh.meshTriangles[len(mesh.meshTriangles)-1]#triangle 1
3121        state.r[3]=state.current_triangle#triangle 2
3122
3123    if state.next_case == 'boundary':
3124        state.repairCaseBoundary()
3125
3126
3127def split(mesh, triangle, new_point):
3128        """
3129        Given a mesh, triangle_id and a new point,
3130        split the corrosponding triangle into two
3131        new triangles and update the mesh.
3132        """
3133
3134        new_triangle1 = Triangle(new_point,triangle.vertices[0],triangle.vertices[1],attribute = triangle.attribute, neighbors = None)
3135        new_triangle2 = Triangle(new_point,triangle.vertices[2],triangle.vertices[0],attribute = triangle.attribute, neighbors = None)
3136
3137        new_triangle1.setNeighbors(triangle.neighbors[2],None,new_triangle2)
3138        new_triangle2.setNeighbors(triangle.neighbors[1],new_triangle1,None)
3139
3140        mesh.meshTriangles.append(new_triangle1)
3141
3142        triangle.vertices = new_triangle2.vertices
3143        triangle.neighbors = new_triangle2.neighbors
3144
3145
3146class State:
3147
3148    def __init__(self):
3149        pass
3150
3151class BisectionState(State):
3152
3153
3154    def __init__(self,mesh):
3155        self.len = len(mesh.meshTriangles)
3156        self.refined_triangles = {}
3157        self.mesh = mesh
3158        self.current_triangle = None
3159        self.case = 'start'
3160        self.end = False
3161        self.r = [None,None,None,None,None]
3162
3163    def start(self, triangle):
3164        self.current_triangle = triangle
3165        self.case = 'start'
3166        self.end = False
3167        self.r = [None,None,None,None,None]
3168
3169    def getState(self,mesh):
3170        if not self.case == 'vertex':
3171            self.new_point=self.getNewVertex(mesh, self.current_triangle)
3172            #self.neighbour=self.getNeighbour(mesh, self.current_triangle)
3173            self.neighbour = self.current_triangle.neighbors[0]
3174            if not self.neighbour is None:
3175                self.neighbour.rotate_longest_side()
3176            self.next_case = self.get_next_case(mesh,self.neighbour,self.current_triangle)
3177        if self.case == 'vertex':
3178            self.new_point=self.old_point
3179
3180
3181    def evolve(self):
3182        if self.case == 'vertex':
3183            self.end = True
3184
3185        self.last_case = self.case
3186        self.case = self.next_case
3187
3188        self.old_point = self.new_point
3189        self.current_triangle = self.neighbour
3190
3191        if self.case == 'boundary':
3192            self.end = True
3193        self.refined_triangles[self.r[2]]=1
3194        self.refined_triangles[self.r[3]]=1
3195        if not self.r[4] is None:
3196            self.refined_triangles[self.r[4]]=1
3197        self.r[0]=self.r[2]
3198        self.r[1]=self.r[3]
3199
3200
3201    def getNewVertex(self,mesh,triangle):
3202        coordinate1 = triangle.vertices[1]
3203        coordinate2 = triangle.vertices[2]
3204        a = ([coordinate1.x*1.,coordinate1.y*1.])
3205        b = ([coordinate2.x*1.,coordinate2.y*1.])
3206        attributes = []
3207        for i in range(len(coordinate1.attributes)):
3208            att = (coordinate1.attributes[i]+coordinate2.attributes[i])/2
3209            attributes.append(att)
3210        new_coordinate = [((a[0]-b[0])/2+b[0]),((a[1]-b[1])/2+b[1])]
3211        newVertex = Vertex(new_coordinate[0],new_coordinate[1], attributes = attributes)
3212        mesh.maxVertexIndex+=1
3213        newVertex.index = mesh.maxVertexIndex
3214        mesh.meshVertices.append(newVertex)
3215        return newVertex
3216
3217    def get_next_case(self, mesh,neighbour,triangle):
3218        """
3219        Given the locations of two neighbouring triangles,
3220        examine the interior indices of their vertices (i.e.
3221        0,1 or 2) to determine what how the neighbour needs
3222        to be refined.
3223        """
3224        if (neighbour is None):
3225                next_case = 'boundary'
3226        else:
3227                if triangle.vertices[1].x==neighbour.vertices[2].x:
3228                    if triangle.vertices[1].y==neighbour.vertices[2].y:
3229                        next_case = 'vertex'
3230                if triangle.vertices[1].x==neighbour.vertices[0].x:
3231                    if triangle.vertices[1].y==neighbour.vertices[0].y:
3232                        next_case = 'two'
3233                if triangle.vertices[1].x==neighbour.vertices[1].x:
3234                    if triangle.vertices[1].y==neighbour.vertices[1].y:
3235                        next_case = 'one'
3236        return next_case
3237
3238
3239
3240    def repairCaseVertex(self):
3241
3242        r = self.r
3243
3244
3245        self.link(r[0],r[2])
3246        self.repair(r[0])
3247
3248        self.link(r[1],r[3])
3249        self.repair(r[1])
3250
3251        self.repair(r[2])
3252
3253        self.repair(r[3])
3254
3255
3256    def repairCaseOne(self):
3257        r = self.rkey
3258
3259
3260        self.link(r[0],r[2])
3261        self.repair(r[0])
3262
3263        self.link(r[1],r[4])
3264        self.repair(r[1])
3265
3266        self.repair(r[4])
3267
3268    def repairCaseTwo(self):
3269        r = self.r
3270
3271        self.link(r[0],r[4])
3272        self.repair(r[0])
3273
3274        self.link(r[1],r[3])
3275        self.repair(r[1])
3276
3277        self.repair(r[4])
3278
3279    def repairCaseBoundary(self):
3280        r = self.r
3281        self.repair(r[2])
3282        self.repair(r[3])
3283
3284
3285
3286    def repair(self,triangle):
3287        """
3288        Given a triangle that knows its neighbours, this will
3289        force the neighbours to comply.
3290
3291        However, it needs to compare the vertices of triangles
3292        for this implementation
3293
3294        But it doesn't work for invalid neighbour structures
3295        """
3296        n=triangle.neighbors
3297        for i in (0,1,2):
3298            if not n[i] is None:
3299                for j in (0,1,2):#to find which side of the list is broken
3300                    if not (n[i].vertices[j] in triangle.vertices):
3301                    #ie if j is the side of n that needs fixing
3302                        k = j
3303                n[i].neighbors[k]=triangle
3304
3305
3306
3307    def link(self,triangle1,triangle2):
3308        """
3309        make triangle1 neighbors point to t
3310                #count = 0riangle2
3311        """
3312        count = 0
3313        for i in (0,1,2):#to find which side of the list is broken
3314            if not (triangle1.vertices[i] in triangle2.vertices):
3315                j = i
3316                count+=1
3317        assert count == 1
3318        triangle1.neighbors[j]=triangle2
3319
3320class Discretised_Tuple_Set:
3321    """
3322    if a={(0.0):[(0.01),(0.02)],(0.2):[(0.17)]}
3323    a[(0.01)]=a[(0.0)]=[(0.01),(0.02)]
3324    a[(10000)]=[] #NOT KEYERROR
3325
3326    a.append[(0.01)]
3327    => {0.0:[(0.01),(0.02),(0.01)],0.2:[(0.17)]}
3328
3329    #NOT IMPLIMENTED
3330    a.remove[(0.01)]
3331    => {(0.0):[(0.02),(0.01)],0.2:[(0.17)]}
3332
3333    a.remove[(0.17)]
3334    => {(0.0):[(0.02),(0.01)],0.2:[]}
3335    #NOT IMPLIMENTED
3336    at a.precision = 2:
3337        a.round_up_rel[0.0]=
3338        a.round_flat[0.0]=
3339        a.round_down_rel[0.0]=
3340
3341        a.up((0.1,2.04))=
3342
3343    If t_rel = 0, nothing gets rounded into
3344    two bins. If t_rel = 0.5, everything does.
3345
3346    Ideally, precision can be set high, so that multiple
3347    entries are rarely in the same bin. And t_rel should
3348    be low (<0.1 for 1 dimension!,<(0.1/n) for small n!!)
3349    so that it is rare to put items in mutiple bins.
3350
3351    Ex bins per entry = product(a,b,c...,n)
3352    a = 1 or 2 s.t. Ex(a) = 1+2*t_rel
3353    b = 1 or 2 ...
3354
3355    BUT!!! to avoid missing the right bin:
3356    (-10)**(precision+1)*t_rel must be greater than the
3357    greatest possible variation that an identical element
3358    can display.
3359
3360
3361    Note that if tol = 0.5 (the max allowed) 0.6 will round to .7 and .5
3362    but not .6 - this looks wrong, but note that *everything* will round,
3363    so .6 wont be missed as everything close to it will check in .7 and .5.
3364    """
3365    def __init__(self,p_rel = 6,t_rel = 0.01):
3366        self.__p_rel__ = p_rel
3367        self.__t_rel__ = t_rel
3368
3369        self.__p_abs__ = p_rel+1
3370        self.__t_abs__ = t_rel
3371
3372        assert t_rel <= 0.5
3373        self.__items__ = {}
3374        from math import frexp
3375        self.frexp = frexp
3376        roundings = [self.round_up_rel,\
3377        self.round_down_rel,self.round_flat_rel,\
3378        self.round_down_abs,self.round_up_abs,\
3379        self.round_flat_abs]#
3380
3381        self.roundings = roundings
3382
3383    def __repr__(self):
3384        return '%s'%self.__items__
3385
3386    def rounded_keys(self,key):
3387        key = tuple(key)
3388        keys = [key]
3389        keys = self.__rounded_keys__(key)
3390        return (keys)
3391
3392    def __rounded_keys__(self,key):
3393        keys = []
3394        rounded_key=list(key)
3395        rounded_values=list(key)
3396
3397        roundings = list(self.roundings)
3398
3399        #initialise rounded_values
3400        round = roundings.pop(0)
3401        for i in range(len(rounded_values)):
3402            rounded_key[i]=round(key[i])
3403            rounded_values[i]={}
3404            rounded_values[i][rounded_key[i]]=None
3405        keys.append(tuple(rounded_key))
3406       
3407        for round in roundings:
3408            for i in range(len(rounded_key)):
3409                rounded_value=round(key[i])
3410                if not rounded_values[i].has_key(rounded_value):
3411                    #ie unless round_up_rel = round_down_rel
3412                    #so the keys stay unique
3413                    for j in range(len(keys)):
3414                        rounded_key = list(keys[j])
3415                        rounded_key[i]=rounded_value
3416                        keys.append(tuple(rounded_key))
3417        return keys
3418
3419    def append(self,item):
3420        keys = self.rounded_keys(item)
3421        for key in keys:
3422            if self.__items__.has_key(key):
3423                self.__items__[key].append(item)
3424            else:
3425                self.__items__[key]=[item]
3426
3427    def __getitem__(self,key):
3428        answer = []
3429        keys = self.rounded_keys(key)
3430        for key in keys:
3431            if self.__items__.has_key(key):
3432                answer.extend(self.__items__[key])
3433        #if len(answer)==0:
3434        #    raise KeyError#FIXME or return KeyError
3435        #                  #FIXME or just return []?
3436        else:
3437            return answer #FIXME or unique(answer)?
3438
3439    def __delete__(self,item):
3440        keys = self.rounded_keys(item)
3441        answer = False
3442        #if any of the possible keys contains
3443        #a list, return true
3444        for key in keys:       
3445            if self.__items__.has_key(key):
3446                if item in self.__items__[key]:
3447                    self.__items__[key].remove(item)
3448
3449    def remove(self,item):
3450        self.__delete__(item)
3451
3452    def __contains__(self,item):
3453
3454        keys = self.rounded_keys(item)
3455        answer = False
3456        #if any of the possible keys contains
3457        #a list, return true
3458        for key in keys:       
3459            if self.__items__.has_key(key):
3460                if item in self.__items__[key]:
3461                    answer = True
3462        return answer
3463
3464
3465    def has_item(self,item):
3466        return self.__contains__(item)
3467
3468    def round_up_rel2(self,value):
3469         t_rel=self.__t_rel__
3470         #Rounding up the value
3471         m,e = self.frexp(value)
3472         m = m/2
3473         e = e + 1
3474         #m is the mantissa, e the exponent
3475         # 0.5 < |m| < 1.0
3476         m = m+t_rel*(10**-(self.__p_rel__))
3477         #bump m up
3478         m = round(m,self.__p_rel__)
3479         return m*(2.**e)
3480
3481    def round_down_rel2(self,value):
3482         t_rel=self.__t_rel__
3483         #Rounding down the value
3484         m,e = self.frexp(value)
3485         m = m/2
3486         e = e + 1
3487         #m is the mantissa, e the exponent
3488         # 0.5 < m < 1.0
3489         m = m-t_rel*(10**-(self.__p_rel__))
3490         #bump the |m| down, by 5% or whatever
3491         #self.p_rel dictates
3492         m = round(m,self.__p_rel__)
3493         return m*(2.**e)
3494
3495    def round_flat_rel2(self,value):
3496    #redundant
3497         m,e = self.frexp(value)
3498         m = m/2
3499         e = e + 1
3500         m = round(m,self.__p_rel__)
3501         return m*(2.**e)
3502
3503    def round_up_rel(self,value):
3504         t_rel=self.__t_rel__
3505         #Rounding up the value
3506         m,e = self.frexp(value)
3507         #m is the mantissa, e the exponent
3508         # 0.5 < |m| < 1.0
3509         m = m+t_rel*(10**-(self.__p_rel__))
3510         #bump m up
3511         m = round(m,self.__p_rel__)
3512         return m*(2.**e)
3513
3514    def round_down_rel(self,value):
3515         t_rel=self.__t_rel__
3516         #Rounding down the value
3517         m,e = self.frexp(value)
3518         #m is the mantissa, e the exponent
3519         # 0.5 < m < 1.0
3520         m = m-t_rel*(10**-(self.__p_rel__))
3521         #bump the |m| down, by 5% or whatever
3522         #self.p_rel dictates
3523         m = round(m,self.__p_rel__)
3524         return m*(2.**e)
3525
3526    def round_flat_rel(self,value):
3527    #redundant
3528         m,e = self.frexp(value)
3529         m = round(m,self.__p_rel__)
3530         return m*(2.**e)
3531
3532    def round_up_abs(self,value):
3533         t_abs=self.__t_abs__
3534         #Rounding up the value
3535         m = value+t_abs*(10**-(self.__p_abs__))
3536         #bump m up
3537         m = round(m,self.__p_abs__)
3538         return m
3539
3540    def round_down_abs(self,value):
3541         t_abs=self.__t_abs__
3542         #Rounding down the value
3543         m = value-t_abs*(10**-(self.__p_abs__))
3544         #bump the |m| down, by 5% or whatever
3545         #self.p_rel dictates
3546         m = round(m,self.__p_abs__)
3547         return m
3548
3549    def round_flat_abs(self,value):
3550    #redundant?
3551         m = round(value,self.__p_abs__)
3552         return m
3553
3554    def keys(self):
3555        return self.__items__.keys()
3556
3557
3558class Mapped_Discretised_Tuple_Set(Discretised_Tuple_Set):
3559    """
3560    This is a discretised tuple set, but
3561    based on a mapping. The mapping MUST
3562    return a sequence.
3563
3564    example:
3565    def weight(animal):
3566        return [animal.weight]
3567   
3568    a = Mapped_Discretised_Tuple_Set(weight)
3569    a.append[cow]
3570    a.append[fox]
3571    a.append[horse]
3572
3573    a[horse]    -> [cow,horse]
3574    a[dog]      -> [fox]
3575    a[elephant] -> []
3576    """
3577    def __init__(self,mapping,p_rel = 6, t_rel=0.01):
3578        Discretised_Tuple_Set.__init__\
3579         (self,p_rel,t_rel = t_rel)
3580        self.mapping = mapping
3581
3582    def rounded_keys(self,key):
3583        mapped_key = tuple(self.mapping(key))
3584        keys = self.__rounded_keys__(mapped_key)
3585        return keys
3586
3587class Affine_Linespace(Mapped_Discretised_Tuple_Set):
3588    """
3589    The affine linespace creates a way to record and compare lines.
3590    Precision is a bit of a hack, but it creates a way to avoid
3591    misses caused by round offs (between two lines of different
3592    lenghts, the short one gets rounded off more).
3593    I am starting to think that a quadratic search would be faster.
3594    Nearly.
3595    """
3596    def __init__(self,p_rel=4,t_rel=0.2):
3597        Mapped_Discretised_Tuple_Set.__init__\
3598            (self,self.affine_line,\
3599            p_rel=p_rel,t_rel=t_rel)
3600
3601        roundings = \
3602        [self.round_down_rel,self.round_up_rel,self.round_flat_rel]
3603        self.roundings = roundings
3604        #roundings = \
3605        #[self.round_down_abs,self.round_up_abs,self.round_flat_abs]
3606        #self.roundings = roundings
3607
3608    def affine_line(self,line):
3609        point_1 = line[0]
3610        point_2 = line[1]
3611        #returns the equation of a line
3612        #between two points, in the from
3613        #(a,b,-c), as in ax+by-c=0
3614        #or line *dot* (x,y,1) = (0,0,0)
3615
3616        #Note that it normalises the line
3617        #(a,b,-c) so Line*Line = 1.
3618        #This does not change the mathematical
3619        #properties, but it makes comparism
3620        #easier.
3621
3622        #There are probably better algorithms.
3623        x1 = point_1[0]
3624        x2 = point_2[0]
3625        y1 = point_1[1]
3626        y2 = point_2[1]
3627        dif_x = x1-x2
3628        dif_y = y1-y2
3629
3630        if dif_x == dif_y == 0:
3631            msg = 'points are the same'
3632            raise msg
3633        elif abs(dif_x)>=abs(dif_y):
3634            alpha = (-dif_y)/dif_x
3635            #a = alpha * b
3636            b = -1.
3637            c = (x1*alpha+x2*alpha+y1+y2)/2.
3638            a = alpha*b
3639        else:
3640            beta = dif_x/(-dif_y)
3641            #b = beta * a
3642            a = 1.
3643            c = (x1+x2+y1*beta+y2*beta)/2.
3644            b = beta*a
3645        mag = abs(a)+abs(b)
3646        #This does not change the mathematical
3647        #properties, but it makes comparism possible.
3648
3649        #note that the gradient is b/a, or (a/b)**-1.
3650        #so
3651
3652        #if a == 0:
3653        #    sign_a = 1.
3654        #else:
3655        #    sign_a = a/((a**2)**0.5)
3656        #if b == 0:
3657        #    sign_b = 1.
3658        #else:
3659        #    sign_b = b/((b**2)**0.5)
3660        #if c == 0:
3661        #    sign_c = 1.
3662        #else:
3663        #    sign_c = c/((c**2)**0.5)
3664        #a = a/mag*sign_a
3665        #b = b/mag*sign_b
3666        #c = c/mag*sign_c
3667        a = a/mag
3668        b = b/mag
3669        c = c/mag
3670        return a,b,c
3671
3672
3673# FIXME (DSG-DSG)
3674# To do-
3675# Create a clear interface. eg
3676# have the interface methods more at the top of this file and add comments
3677# for the interface functions/methods, use function_name (not functionName),
3678
3679#Currently
3680#function_name methods assume absolute values.  Geo-refs can be passed in.
3681#
3682
3683# instead of functionName
3684if __name__ == "__main__":
3685    #from mesh import *
3686    # THIS CAN BE DELETED
3687    m = Mesh()
3688    dict = importUngenerateFile("ungen_test.txt")
3689    m.addVertsSegs(dict)
3690    print m3
Note: See TracBrowser for help on using the repository browser.