source: inundation/pmesh/mesh.py @ 3380

Last change on this file since 3380 was 3351, checked in by duncan, 19 years ago

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