source: inundation/pmesh/mesh.py @ 2392

Last change on this file since 2392 was 2392, checked in by duncan, 18 years ago

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