source: inundation/pmesh/mesh.py @ 3031

Last change on this file since 3031 was 3027, checked in by duncan, 19 years ago

directory name change to avoid name conflicts.

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