source: branches/numpy/anuga/fit_interpolate/fit.py @ 6304

Last change on this file since 6304 was 6304, checked in by rwilson, 14 years ago

Initial commit of numpy changes. Still a long way to go.

File size: 26.3 KB
Line 
1"""Least squares fitting.
2
3   Implements a penalised least-squares fit.
4   putting point data onto the mesh.
5
6   The penalty term (or smoothing term) is controlled by the smoothing
7   parameter alpha.
8   With a value of alpha=0, the fit function will attempt
9   to interpolate as closely as possible in the least-squares sense.
10   With values alpha > 0, a certain amount of smoothing will be applied.
11   A positive alpha is essential in cases where there are too few
12   data points.
13   A negative alpha is not allowed.
14   A typical value of alpha is 1.0e-6
15
16
17   Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
18   Geoscience Australia, 2004.
19
20   TO DO
21   * test geo_ref, geo_spatial
22
23   IDEAS
24   * (DSG-) Change the interface of fit, so a domain object can
25      be passed in. (I don't know if this is feasible). If could
26      save time/memory.
27"""
28import types
29
30from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
31from anuga.caching import cache           
32from anuga.geospatial_data.geospatial_data import Geospatial_data, \
33     ensure_absolute
34from anuga.fit_interpolate.general_fit_interpolate import FitInterpolate
35from anuga.utilities.sparse import Sparse, Sparse_CSR
36from anuga.utilities.polygon import inside_polygon, is_inside_polygon
37from anuga.fit_interpolate.search_functions import search_tree_of_vertices
38
39from anuga.utilities.cg_solve import conjugate_gradient
40from anuga.utilities.numerical_tools import ensure_numeric, gradient
41from anuga.config import default_smoothing_parameter as DEFAULT_ALPHA
42
43import exceptions
44class TooFewPointsError(exceptions.Exception): pass
45class VertsWithNoTrianglesError(exceptions.Exception): pass
46
47import numpy as num
48
49
50class Fit(FitInterpolate):
51   
52    def __init__(self,
53                 vertex_coordinates=None,
54                 triangles=None,
55                 mesh=None,
56                 mesh_origin=None,
57                 alpha = None,
58                 verbose=False,
59                 max_vertices_per_cell=None):
60
61
62        """
63        Fit data at points to the vertices of a mesh.
64
65        Inputs:
66
67          vertex_coordinates: List of coordinate pairs [xi, eta] of
68              points constituting a mesh (or an m x 2 numeric array or
69              a geospatial object)
70              Points may appear multiple times
71              (e.g. if vertices have discontinuities)
72
73          triangles: List of 3-tuples (or a numeric array) of
74              integers representing indices of all vertices in the mesh.
75
76          mesh_origin: A geo_reference object or 3-tuples consisting of
77              UTM zone, easting and northing.
78              If specified vertex coordinates are assumed to be
79              relative to their respective origins.
80
81          max_vertices_per_cell: Number of vertices in a quad tree cell
82          at which the cell is split into 4.
83
84          Note: Don't supply a vertex coords as a geospatial object and
85              a mesh origin, since geospatial has its own mesh origin.
86
87
88        Usage,
89        To use this in a blocking way, call  build_fit_subset, with z info,
90        and then fit, with no point coord, z info.
91       
92        """
93        # Initialise variabels
94        if alpha is None:
95            self.alpha = DEFAULT_ALPHA
96        else:   
97            self.alpha = alpha
98           
99        FitInterpolate.__init__(self,
100                 vertex_coordinates,
101                 triangles,
102                 mesh,
103                 mesh_origin,
104                 verbose,
105                 max_vertices_per_cell)
106       
107        m = self.mesh.number_of_nodes # Nbr of basis functions (vertices)
108       
109        self.AtA = None
110        self.Atz = None
111
112        self.point_count = 0
113        if self.alpha <> 0:
114            if verbose: print 'Building smoothing matrix'
115            self._build_smoothing_matrix_D()
116           
117        bd_poly = self.mesh.get_boundary_polygon()   
118        self.mesh_boundary_polygon = ensure_numeric(bd_poly)
119           
120    def _build_coefficient_matrix_B(self,
121                                  verbose = False):
122        """
123        Build final coefficient matrix
124
125        Precon
126        If alpha is not zero, matrix D has been built
127        Matrix Ata has been built
128        """
129
130        if self.alpha <> 0:
131            #if verbose: print 'Building smoothing matrix'
132            #self._build_smoothing_matrix_D()
133            self.B = self.AtA + self.alpha*self.D
134        else:
135            self.B = self.AtA
136
137        # Convert self.B matrix to CSR format for faster matrix vector
138        self.B = Sparse_CSR(self.B)
139
140    def _build_smoothing_matrix_D(self):
141        """Build m x m smoothing matrix, where
142        m is the number of basis functions phi_k (one per vertex)
143
144        The smoothing matrix is defined as
145
146        D = D1 + D2
147
148        where
149
150        [D1]_{k,l} = \int_\Omega
151           \frac{\partial \phi_k}{\partial x}
152           \frac{\partial \phi_l}{\partial x}\,
153           dx dy
154
155        [D2]_{k,l} = \int_\Omega
156           \frac{\partial \phi_k}{\partial y}
157           \frac{\partial \phi_l}{\partial y}\,
158           dx dy
159
160
161        The derivatives \frac{\partial \phi_k}{\partial x},
162        \frac{\partial \phi_k}{\partial x} for a particular triangle
163        are obtained by computing the gradient a_k, b_k for basis function k
164        """
165       
166        # FIXME: algorithm might be optimised by computing local 9x9
167        # "element stiffness matrices:
168
169        m = self.mesh.number_of_nodes # Nbr of basis functions (1/vertex)
170
171        self.D = Sparse(m,m)
172
173        # For each triangle compute contributions to D = D1+D2
174        for i in range(len(self.mesh)):
175
176            # Get area
177            area = self.mesh.areas[i]
178
179            # Get global vertex indices
180            v0 = self.mesh.triangles[i,0]
181            v1 = self.mesh.triangles[i,1]
182            v2 = self.mesh.triangles[i,2]
183
184            # Get the three vertex_points
185            xi0 = self.mesh.get_vertex_coordinate(i, 0)
186            xi1 = self.mesh.get_vertex_coordinate(i, 1)
187            xi2 = self.mesh.get_vertex_coordinate(i, 2)
188
189            # Compute gradients for each vertex
190            a0, b0 = gradient(xi0[0], xi0[1], xi1[0], xi1[1], xi2[0], xi2[1],
191                              1, 0, 0)
192
193            a1, b1 = gradient(xi0[0], xi0[1], xi1[0], xi1[1], xi2[0], xi2[1],
194                              0, 1, 0)
195
196            a2, b2 = gradient(xi0[0], xi0[1], xi1[0], xi1[1], xi2[0], xi2[1],
197                              0, 0, 1)
198
199            # Compute diagonal contributions
200            self.D[v0,v0] += (a0*a0 + b0*b0)*area
201            self.D[v1,v1] += (a1*a1 + b1*b1)*area
202            self.D[v2,v2] += (a2*a2 + b2*b2)*area
203
204            # Compute contributions for basis functions sharing edges
205            e01 = (a0*a1 + b0*b1)*area
206            self.D[v0,v1] += e01
207            self.D[v1,v0] += e01
208
209            e12 = (a1*a2 + b1*b2)*area
210            self.D[v1,v2] += e12
211            self.D[v2,v1] += e12
212
213            e20 = (a2*a0 + b2*b0)*area
214            self.D[v2,v0] += e20
215            self.D[v0,v2] += e20
216
217    def get_D(self):
218        return self.D.todense()
219
220
221
222    def _build_matrix_AtA_Atz(self,
223                              point_coordinates,
224                              z,
225                              verbose = False):
226        """Build:
227        AtA  m x m  interpolation matrix, and,
228        Atz  m x a  interpolation matrix where,
229        m is the number of basis functions phi_k (one per vertex)
230        a is the number of data attributes
231
232        This algorithm uses a quad tree data structure for fast binning of
233        data points.
234
235        If Ata is None, the matrices AtA and Atz are created.
236
237        This function can be called again and again, with sub-sets of
238        the point coordinates.  Call fit to get the results.
239       
240        Preconditions
241        z and points are numeric
242        Point_coordindates and mesh vertices have the same origin.
243
244        The number of attributes of the data points does not change
245        """
246       
247        # Build n x m interpolation matrix
248        if self.AtA == None:
249            # AtA and Atz need to be initialised.
250            m = self.mesh.number_of_nodes
251            if len(z.shape) > 1:
252                att_num = z.shape[1]
253                self.Atz = num.zeros((m,att_num), num.float)
254            else:
255                att_num = 1
256                self.Atz = num.zeros((m,), num.float)
257            assert z.shape[0] == point_coordinates.shape[0] 
258
259            AtA = Sparse(m,m)
260            # The memory damage has been done by now.
261        else:
262             AtA = self.AtA #Did this for speed, did ~nothing
263        self.point_count += point_coordinates.shape[0]
264
265
266        inside_indices = inside_polygon(point_coordinates,
267                                        self.mesh_boundary_polygon,
268                                        closed=True,
269                                        verbose=False) # Too much output if True
270
271       
272        n = len(inside_indices)
273
274        # Compute matrix elements for points inside the mesh
275        triangles = self.mesh.triangles # Shorthand
276        for d, i in enumerate(inside_indices):
277            # For each data_coordinate point
278            # if verbose and d%((n+10)/10)==0: print 'Doing %d of %d' %(d, n)
279            x = point_coordinates[i]
280           
281            element_found, sigma0, sigma1, sigma2, k = \
282                           search_tree_of_vertices(self.root, self.mesh, x)
283           
284            if element_found is True:
285                j0 = triangles[k,0] #Global vertex id for sigma0
286                j1 = triangles[k,1] #Global vertex id for sigma1
287                j2 = triangles[k,2] #Global vertex id for sigma2
288
289                sigmas = {j0:sigma0, j1:sigma1, j2:sigma2}
290                js     = [j0,j1,j2]
291
292                for j in js:
293                    self.Atz[j] +=  sigmas[j]*z[i]
294                    #print "self.Atz building", self.Atz
295                    #print "self.Atz[j]", self.Atz[j]
296                    #print " sigmas[j]", sigmas[j]
297                    #print "z[i]",z[i]
298                    #print "result", sigmas[j]*z[i]
299                   
300                    for k in js:
301                        AtA[j,k] += sigmas[j]*sigmas[k]
302            else:
303                # FIXME(Ole): This is the message referred to in ticket:314
304               
305                flag = is_inside_polygon(x,
306                                         self.mesh_boundary_polygon,
307                                         closed=True,
308                                         verbose=False) # Too much output if True               
309                msg = 'Point (%f, %f) is not inside mesh boundary' % tuple(x)
310                assert flag is True, msg               
311               
312                minx = min(self.mesh_boundary_polygon[:,0])
313                maxx = max(self.mesh_boundary_polygon[:,0])               
314                miny = min(self.mesh_boundary_polygon[:,1])
315                maxy = max(self.mesh_boundary_polygon[:,1])
316                msg = 'Could not find triangle for point %s. ' % str(x) 
317                msg += 'Mesh boundary extent is (%.f, %.f), (%.f, %.f)'\
318                    % (minx, maxx, miny, maxy)
319               
320
321                raise Exception(msg)
322            self.AtA = AtA
323
324       
325    def fit(self, point_coordinates_or_filename=None, z=None,
326            verbose=False,
327            point_origin=None,
328            attribute_name=None,
329            max_read_lines=500):
330        """Fit a smooth surface to given 1d array of data points z.
331
332        The smooth surface is computed at each vertex in the underlying
333        mesh using the formula given in the module doc string.
334
335        Inputs:
336        point_coordinates: The co-ordinates of the data points.
337              List of coordinate pairs [x, y] of
338              data points or an nx2 numeric array or a Geospatial_data object
339              or points file filename
340          z: Single 1d vector or array of data at the point_coordinates.
341         
342        """
343       
344        # Use blocking to load in the point info
345        if type(point_coordinates_or_filename) == types.StringType:
346            msg = "Don't set a point origin when reading from a file"
347            assert point_origin is None, msg
348            filename = point_coordinates_or_filename
349
350            G_data = Geospatial_data(filename,
351                                     max_read_lines=max_read_lines,
352                                     load_file_now=False,
353                                     verbose=verbose)
354
355            for i, geo_block in enumerate(G_data):
356                if verbose is True and 0 == i%200: 
357                    # The time this will take
358                    # is dependant on the # of Triangles
359                       
360                    print 'Processing Block %d' %i
361                    # FIXME (Ole): It would be good to say how many blocks
362                    # there are here. But this is no longer necessary
363                    # for pts files as they are reported in geospatial_data
364                    # I suggest deleting this verbose output and make
365                    # Geospatial_data more informative for txt files.
366                    #
367                    # I still think so (12/12/7, Ole).
368           
369
370                   
371                # Build the array
372
373                points = geo_block.get_data_points(absolute=True)
374                z = geo_block.get_attributes(attribute_name=attribute_name)
375                self.build_fit_subset(points, z, verbose=verbose)
376
377               
378            point_coordinates = None
379        else:
380            point_coordinates =  point_coordinates_or_filename
381           
382        if point_coordinates is None:
383            if verbose: print 'Warning: no data points in fit'
384            #assert self.AtA != None, 'no interpolation matrix'
385            #assert self.Atz != None
386            assert not self.AtA is None, 'no interpolation matrix'
387            assert not self.Atz is None
388           
389            # FIXME (DSG) - do  a message
390        else:
391            point_coordinates = ensure_absolute(point_coordinates,
392                                                geo_reference=point_origin)
393            # if isinstance(point_coordinates,Geospatial_data) and z is None:
394            # z will come from the geo-ref
395            self.build_fit_subset(point_coordinates, z, verbose)
396
397        # Check sanity
398        m = self.mesh.number_of_nodes # Nbr of basis functions (1/vertex)
399        n = self.point_count
400        if n<m and self.alpha == 0.0:
401            msg = 'ERROR (least_squares): Too few data points\n'
402            msg += 'There are only %d data points and alpha == 0. ' %n
403            msg += 'Need at least %d\n' %m
404            msg += 'Alternatively, set smoothing parameter alpha to a small '
405            msg += 'positive value,\ne.g. 1.0e-3.'
406            raise TooFewPointsError(msg)
407
408        self._build_coefficient_matrix_B(verbose)
409        loners = self.mesh.get_lone_vertices()
410        # FIXME  - make this as error message.
411        # test with
412        # Not_yet_test_smooth_att_to_mesh_with_excess_verts.
413        if len(loners)>0:
414            msg = 'WARNING: (least_squares): \nVertices with no triangles\n'
415            msg += 'All vertices should be part of a triangle.\n'
416            msg += 'In the future this will be inforced.\n'
417            msg += 'The following vertices are not part of a triangle;\n'
418            msg += str(loners)
419            print msg
420            #raise VertsWithNoTrianglesError(msg)
421       
422       
423        return conjugate_gradient(self.B, self.Atz, self.Atz,
424                                  imax=2*len(self.Atz) )
425
426       
427    def build_fit_subset(self, point_coordinates, z=None, attribute_name=None,
428                              verbose=False):
429        """Fit a smooth surface to given 1d array of data points z.
430
431        The smooth surface is computed at each vertex in the underlying
432        mesh using the formula given in the module doc string.
433
434        Inputs:
435        point_coordinates: The co-ordinates of the data points.
436              List of coordinate pairs [x, y] of
437              data points or an nx2 numeric array or a Geospatial_data object
438        z: Single 1d vector or array of data at the point_coordinates.
439        attribute_name: Used to get the z values from the
440              geospatial object if no attribute_name is specified,
441              it's a bit of a lucky dip as to what attributes you get.
442              If there is only one attribute it will be that one.
443
444        """
445
446        # FIXME(DSG-DSG): Check that the vert and point coords
447        # have the same zone.
448        if isinstance(point_coordinates,Geospatial_data):
449            point_coordinates = point_coordinates.get_data_points( \
450                absolute = True)
451       
452        # Convert input to numeric arrays
453        if z is not None:
454            z = ensure_numeric(z, num.float)
455        else:
456            msg = 'z not specified'
457            assert isinstance(point_coordinates,Geospatial_data), msg
458            z = point_coordinates.get_attributes(attribute_name)
459
460        point_coordinates = ensure_numeric(point_coordinates, num.float)
461        self._build_matrix_AtA_Atz(point_coordinates, z, verbose)
462
463
464############################################################################
465
466def fit_to_mesh(point_coordinates, # this can also be a points file name
467                vertex_coordinates=None,
468                triangles=None,
469                mesh=None,
470                point_attributes=None,
471                alpha=DEFAULT_ALPHA,
472                verbose=False,
473                acceptable_overshoot=1.01, 
474                # FIXME: Move to config - this value is assumed in caching test
475                # FIXME(Ole): Just realised that this was never implemented (29 Jan 2009). I suggest removing it altogether.
476                mesh_origin=None,
477                data_origin=None,
478                max_read_lines=None,
479                attribute_name=None,
480                use_cache=False):
481    """Wrapper around internal function _fit_to_mesh for use with caching.
482   
483    """
484   
485    args = (point_coordinates, )
486    kwargs = {'vertex_coordinates': vertex_coordinates,
487              'triangles': triangles,
488              'mesh': mesh,
489              'point_attributes': point_attributes,
490              'alpha': alpha,
491              'verbose': verbose,
492              'acceptable_overshoot': acceptable_overshoot,
493              'mesh_origin': mesh_origin,
494              'data_origin': data_origin,
495              'max_read_lines': max_read_lines,
496              'attribute_name': attribute_name
497              }
498
499    if use_cache is True:
500        if isinstance(point_coordinates, basestring):
501            # We assume that point_coordinates is the name of a .csv/.txt
502            # file which must be passed onto caching as a dependency
503            # (in case it has changed on disk)
504            dep = [point_coordinates]
505        else:
506            dep = None
507
508           
509        #from caching import myhash
510        #import copy
511        #print args
512        #print kwargs
513        #print 'hashing:'
514        #print 'args', myhash( (args, kwargs) )
515        #print 'again', myhash( copy.deepcopy( (args, kwargs)) )       
516       
517        #print 'mesh hash', myhash( kwargs['mesh'] )       
518       
519        #print '-------------------------'
520        #print 'vertices hash', myhash( kwargs['mesh'].nodes )
521        #print 'triangles hash', myhash( kwargs['mesh'].triangles )
522        #print '-------------------------'       
523       
524        #for key in mesh.__dict__:
525        #    print key, myhash(mesh.__dict__[key])
526       
527        #for key in mesh.quantities.keys():
528        #    print key, myhash(mesh.quantities[key])
529       
530        #import sys; sys.exit()
531           
532        return cache(_fit_to_mesh,
533                     args, kwargs,
534                     verbose=verbose,
535                     compression=False,
536                     dependencies=dep)
537    else:
538        return apply(_fit_to_mesh,
539                     args, kwargs)
540
541def _fit_to_mesh(point_coordinates, # this can also be a points file name
542                 vertex_coordinates=None,
543                 triangles=None,
544                 mesh=None,
545                 point_attributes=None,
546                 alpha=DEFAULT_ALPHA,
547                 verbose=False,
548                 acceptable_overshoot=1.01,
549                 mesh_origin=None,
550                 data_origin=None,
551                 max_read_lines=None,
552                 attribute_name=None):
553    """
554    Fit a smooth surface to a triangulation,
555    given data points with attributes.
556
557
558        Inputs:
559        vertex_coordinates: List of coordinate pairs [xi, eta] of
560              points constituting a mesh (or an m x 2 numeric array or
561              a geospatial object)
562              Points may appear multiple times
563              (e.g. if vertices have discontinuities)
564
565          triangles: List of 3-tuples (or a numeric array) of
566          integers representing indices of all vertices in the mesh.
567
568          point_coordinates: List of coordinate pairs [x, y] of data points
569          (or an nx2 numeric array). This can also be a .csv/.txt/.pts
570          file name.
571
572          alpha: Smoothing parameter.
573
574          acceptable overshoot: NOT IMPLEMENTED
575          controls the allowed factor by which
576          fitted values
577          may exceed the value of input data. The lower limit is defined
578          as min(z) - acceptable_overshoot*delta z and upper limit
579          as max(z) + acceptable_overshoot*delta z
580         
581
582          mesh_origin: A geo_reference object or 3-tuples consisting of
583              UTM zone, easting and northing.
584              If specified vertex coordinates are assumed to be
585              relative to their respective origins.
586         
587
588          point_attributes: Vector or array of data at the
589                            point_coordinates.
590
591    """
592
593    if mesh is None:
594        # FIXME(DSG): Throw errors if triangles or vertex_coordinates
595        # are None
596           
597        #Convert input to numeric arrays
598        triangles = ensure_numeric(triangles, num.int)
599        vertex_coordinates = ensure_absolute(vertex_coordinates,
600                                             geo_reference = mesh_origin)
601
602        if verbose: print 'FitInterpolate: Building mesh'       
603        mesh = Mesh(vertex_coordinates, triangles)
604        mesh.check_integrity()
605   
606   
607    interp = Fit(mesh=mesh,
608                 verbose=verbose,
609                 alpha=alpha)
610
611    vertex_attributes = interp.fit(point_coordinates,
612                                   point_attributes,
613                                   point_origin=data_origin,
614                                   max_read_lines=max_read_lines,
615                                   attribute_name=attribute_name,
616                                   verbose=verbose)
617
618       
619    # Add the value checking stuff that's in least squares.
620    # Maybe this stuff should get pushed down into Fit.
621    # at least be a method of Fit.
622    # Or intigrate it into the fit method, saving teh max and min's
623    # as att's.
624   
625    return vertex_attributes
626
627
628#def _fit(*args, **kwargs):
629#    """Private function for use with caching. Reason is that classes
630#    may change their byte code between runs which is annoying.
631#    """
632#   
633#    return Fit(*args, **kwargs)
634
635
636def fit_to_mesh_file(mesh_file, point_file, mesh_output_file,
637                     alpha=DEFAULT_ALPHA, verbose= False,
638                     expand_search = False,
639                     precrop = False,
640                     display_errors = True):
641    """
642    Given a mesh file (tsh) and a point attribute file, fit
643    point attributes to the mesh and write a mesh file with the
644    results.
645
646    Note: the points file needs titles.  If you want anuga to use the tsh file,
647    make sure the title is elevation.
648
649    NOTE: Throws IOErrors, for a variety of file problems.
650   
651    """
652
653    from load_mesh.loadASCII import import_mesh_file, \
654         export_mesh_file, concatinate_attributelist
655
656
657    try:
658        mesh_dict = import_mesh_file(mesh_file)
659    except IOError,e:
660        if display_errors:
661            print "Could not load bad file. ", e
662        raise IOError  #Could not load bad mesh file.
663   
664    vertex_coordinates = mesh_dict['vertices']
665    triangles = mesh_dict['triangles']
666    if isinstance(mesh_dict['vertex_attributes'], num.ndarray):
667        old_point_attributes = mesh_dict['vertex_attributes'].tolist()
668    else:
669        old_point_attributes = mesh_dict['vertex_attributes']
670
671    if isinstance(mesh_dict['vertex_attribute_titles'], num.ndarray):
672        old_title_list = mesh_dict['vertex_attribute_titles'].tolist()
673    else:
674        old_title_list = mesh_dict['vertex_attribute_titles']
675
676    if verbose: print 'tsh file %s loaded' %mesh_file
677
678    # load in the points file
679    try:
680        geo = Geospatial_data(point_file, verbose=verbose)
681    except IOError,e:
682        if display_errors:
683            print "Could not load bad file. ", e
684        raise IOError  #Re-raise exception 
685
686    point_coordinates = geo.get_data_points(absolute=True)
687    title_list,point_attributes = concatinate_attributelist( \
688        geo.get_all_attributes())
689
690    if mesh_dict.has_key('geo_reference') and \
691           not mesh_dict['geo_reference'] is None:
692        mesh_origin = mesh_dict['geo_reference'].get_origin()
693    else:
694        mesh_origin = None
695
696    if verbose: print "points file loaded"
697    if verbose: print "fitting to mesh"
698    f = fit_to_mesh(point_coordinates,
699                    vertex_coordinates,
700                    triangles,
701                    None,
702                    point_attributes,
703                    alpha = alpha,
704                    verbose = verbose,
705                    data_origin = None,
706                    mesh_origin = mesh_origin)
707    if verbose: print "finished fitting to mesh"
708
709    # convert array to list of lists
710    new_point_attributes = f.tolist()
711    #FIXME have this overwrite attributes with the same title - DSG
712    #Put the newer attributes last
713    if old_title_list <> []:
714        old_title_list.extend(title_list)
715        #FIXME can this be done a faster way? - DSG
716        for i in range(len(old_point_attributes)):
717            old_point_attributes[i].extend(new_point_attributes[i])
718        mesh_dict['vertex_attributes'] = old_point_attributes
719        mesh_dict['vertex_attribute_titles'] = old_title_list
720    else:
721        mesh_dict['vertex_attributes'] = new_point_attributes
722        mesh_dict['vertex_attribute_titles'] = title_list
723
724    if verbose: print "exporting to file ", mesh_output_file
725
726    try:
727        export_mesh_file(mesh_output_file, mesh_dict)
728    except IOError,e:
729        if display_errors:
730            print "Could not write file. ", e
731        raise IOError
Note: See TracBrowser for help on using the repository browser.