source: inundation/pmesh/mesh.py @ 3180

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

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