source: inundation/pmesh/mesh.py @ 3315

Last change on this file since 3315 was 3193, checked in by duncan, 19 years ago

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