# source:inundation/pyvolution/util.py@1905

Last change on this file since 1905 was 1905, checked in by ole, 18 years ago

Better docstring in file_function

File size: 25.8 KB
Line
1"""This module contains various auxiliary function used by pyvolution.
2
3It is also a clearing house for functions that may later earn a module
4of their own.
5"""
6
7#FIXME: Move polygon stuff out
8
9def angle(v):
10    """Compute angle between e1 (the unit vector in the x-direction)
11    and the specified vector
12    """
13
14    from math import acos, pi, sqrt
15    from Numeric import sum, array
16
17    l = sqrt( sum (array(v)**2))
18    v1 = v/l
19    v2 = v/l
20
21
22    theta = acos(v1)
23
24    #try:
25    #   theta = acos(v1)
26    #except ValueError, e:
27    #    print 'WARNING (util.py): Angle acos(%s) failed: %s' %(str(v1), e)
28    #
29    #    #FIXME, hack to avoid acos(1.0) Value error
30    #    # why is it happening?
31    #    # is it catching something we should avoid?
32    #    #FIXME (Ole): When did this happen? We need a unit test to
33    #    #reveal this condition or otherwise remove the hack
34    #
35    #    s = 1e-6
36    #    if (v1+s >  1.0) and (v1-s < 1.0) :
37    #        theta = 0.0
38    #    elif (v1+s >  -1.0) and (v1-s < -1.0):
39    #        theta = 3.1415926535897931
40    #    print 'WARNING (util.py): angle v1 is %f, setting acos to %f '\
41    #          %(v1, theta)
42
43    if v2 < 0:
45        theta = 2*pi-theta
46
47    return theta
48
49
50def anglediff(v0, v1):
51    """Compute difference between angle of vector x0, y0 and x1, y1.
52    This is used for determining the ordering of vertices,
53    e.g. for checking if they are counter clockwise.
54
55    Always return a positive value
56    """
57
58    from math import pi
59
60    a0 = angle(v0)
61    a1 = angle(v1)
62
63    #Ensure that difference will be positive
64    if a0 < a1:
65        a0 += 2*pi
66
67    return a0-a1
68
69
70def mean(x):
71    from Numeric import sum
72    return sum(x)/len(x)
73
74
75def point_on_line(x, y, x0, y0, x1, y1):
76    """Determine whether a point is on a line segment
77
78    Input: x, y, x0, x0, x1, y1: where
79        point is given by x, y
80        line is given by (x0, y0) and (x1, y1)
81
82    """
83
84    from Numeric import array, dot, allclose
85    from math import sqrt
86
87    a = array([x - x0, y - y0])
88    a_normal = array([a, -a])
89
90    b = array([x1 - x0, y1 - y0])
91
92    if dot(a_normal, b) == 0:
93        #Point is somewhere on the infinite extension of the line
94
95        len_a = sqrt(sum(a**2))
96        len_b = sqrt(sum(b**2))
97        if dot(a, b) >= 0 and len_a <= len_b:
98           return True
99        else:
100           return False
101    else:
102      return False
103
104def ensure_numeric(A, typecode = None):
105    """Ensure that sequence is a Numeric array.
106    Inputs:
107        A: Sequance. If A is already a Numeric array it will be returned
108                     unaltered
109                     If not, an attempt is made to convert it to a Numeric
110                     array
111        typecode: Numeric type. If specified, use this in the conversion.
112                                If not, let Numeric decide
113
114    This function is necessary as array(A) can cause memory overflow.
115    """
116
117    from Numeric import ArrayType, array
118
119    if typecode is None:
120        if type(A) == ArrayType:
121            return A
122        else:
123            return array(A)
124    else:
125        if type(A) == ArrayType:
126            if A.typecode == typecode:
127                return array(A)  #FIXME: Shouldn't this be just return A?
128            else:
129                return A.astype(typecode)
130        else:
131            return array(A).astype(typecode)
132
133
134
135def file_function(filename,
136                  domain = None,
137                  quantities = None,
138                  interpolation_points = None, verbose = False):
139    """Read time history of spatial data from NetCDF file and return
140    a callable object.
141
142    Input variables:
143
144    filename - Name of sww or tms file
145
146       If the file has extension 'sww' then it is assumed to be spatio-temporal
147       or temporal and the callable object will have the form f(t,x,y) or f(t)
148       depending on whether the file contains spatial data
149
150       If the file has extension 'tms' then it is assumed to be temporal only
151       and the callable object will have the form f(t)
152
153       Either form will return interpolated values based on the input file
154       using the underlying interpolation_function.
155
156    domain - Associated domain object
157       If domain is specified, model time (domain.starttime)
158       will be checked and possibly modified.
159
160       All times are assumed to be in UTC
161
162       All spatial information is assumed to be in absolute UTM coordinates.
163
164    quantities - the name of the quantity to be interpolated or a
165                 list of quantity names. The resulting function will return
166                 a tuple of values - one for each quantity
167
168    interpolation_points - list of absolute UTM coordinates for points at
169    which values are sought
170
171
172
173    See Interpolation function for further documentation
174    """
175
176
177    #FIXME (OLE): Should check origin of domain against that of file
178    #In fact, this is where origin should be converted to that of domain
179    #Also, check that file covers domain fully.
180    #If we use the suggested Point_set class for interpolation points
181    #here it would be easier
182
183    #Take into account:
184    #- domain's georef
185    #- sww file's georef
186    #- interpolation points as absolute UTM coordinates
187
188
189    assert type(filename) == type(''),\
190               'First argument to File_function must be a string'
191
192    try:
193        fid = open(filename)
194    except Exception, e:
195        msg = 'File "%s" could not be opened: Error="%s"'\
196                  %(filename, e)
197        raise msg
198
200    fid.close()
201
202    if quantities is None:
203        if domain is not None:
204            quantities = domain.conserved_quantities
205
206
207
208    if line[:3] == 'CDF':
209        return get_netcdf_file_function(filename, domain, quantities,
210                                        interpolation_points, verbose = verbose)
211    else:
212        raise 'Must be a NetCDF File'
213
214
215
216def get_netcdf_file_function(filename, domain=None, quantity_names=None,
217                             interpolation_points=None, verbose = False):
218    """Read time history of spatial data from NetCDF sww file and
219    return a callable object f(t,x,y)
220    which will return interpolated values based on the input file.
221
222    If domain is specified, model time (domain.starttime)
223    will be checked and possibly modified
224
225    All times are assumed to be in UTC
226
227    See Interpolation function for further documetation
228
229
230    """
231
232    #verbose = True
233
234    #FIXME: Check that model origin is the same as file's origin
235    #(both in UTM coordinates)
236    #If not - modify those from file to match domain
237    #Take this code from e.g. dem2pts in data_manager.py
238    #FIXME: Use geo_reference to read and write xllcorner...
239
240
241    import time, calendar, types
242    from config import time_format
243    from Scientific.IO.NetCDF import NetCDFFile
244    from Numeric import array, zeros, Float, alltrue, concatenate, reshape
245    from util import ensure_numeric
246
247    #Open NetCDF file
248    if verbose: print 'Reading', filename
249    fid = NetCDFFile(filename, 'r')
250
251    if type(quantity_names) == types.StringType:
252        quantity_names = [quantity_names]
253
254    if quantity_names is None or len(quantity_names) < 1:
255        #If no quantities are specified get quantities from file
256        #x, y, time are assumed as the independent variables so
257        #they are excluded as potentiol quantities
258        quantity_names = []
259        for name in fid.variables.keys():
260            if name not in ['x', 'y', 'time']:
261                quantity_names.append(name)
262
263    if len(quantity_names) < 1:
264        msg = 'ERROR: At least one independent value must be specified'
265        raise msg
266
267
268    if interpolation_points is not None:
269        interpolation_points = ensure_numeric(interpolation_points)
270
271
272
273    #Now assert that requested quantitites (and the independent ones)
274    #are present in file
275    missing = []
276    for quantity in ['time'] + quantity_names:
277        if not fid.variables.has_key(quantity):
278            missing.append(quantity)
279
280    if len(missing) > 0:
281        msg = 'Quantities %s could not be found in file %s'\
282              %(str(missing), filename)
283        fid.close()
284        raise msg
285
286    #Decide whether this data has a spatial dimension
287    spatial = True
288    for quantity in ['x', 'y']:
289        if not fid.variables.has_key(quantity):
290            spatial = False
291
292    if filename[-3:] == 'tms' and spatial is True:
293        msg = 'Files of type tms must not contain spatial information'
294        raise msg
295
296    if filename[-3:] == 'sww' and spatial is False:
297        msg = 'Files of type sww must contain spatial information'
298        raise msg
299
300    #Get first timestep
301    try:
302        starttime = fid.starttime
303    except ValueError:
304        msg = 'Could not read starttime from file %s' %filename
305        raise msg
306
307    #Get variables
308    if verbose: print 'Get variables'
309    time = fid.variables['time'][:]
310
311    if spatial:
312        #Get origin
313        xllcorner = fid.xllcorner
314        yllcorner = fid.yllcorner
315        zone = fid.zone
316
317        x = fid.variables['x'][:]
318        y = fid.variables['y'][:]
319        triangles = fid.variables['volumes'][:]
320
321        x = reshape(x, (len(x),1))
322        y = reshape(y, (len(y),1))
323        vertex_coordinates = concatenate((x,y), axis=1) #m x 2 array
324
325        if interpolation_points is not None:
327            interpolation_points[:,0] -= xllcorner
328            interpolation_points[:,1] -= yllcorner
329
330
331
332
333    if domain is not None:
334        #Update domain.startime if it is *earlier* than starttime
335        if starttime > domain.starttime:
336            msg = 'WARNING: Start time as specified in domain (%f)'\
337                  %domain.starttime
338            msg += ' is earlier than the starttime of file %s (%f).'\
339                     %(filename, starttime)
340            msg += ' Modifying domain starttime accordingly.'
341
342            if verbose: print msg
343            domain.starttime = starttime #Modifying model time
344            if verbose: print 'Domain starttime is now set to %f'\
345               %domain.starttime
346
347
348        #If domain.startime is *later* than starttime,
349        #move time back - relative to domain's time
350        if domain.starttime > starttime:
351            time = time - domain.starttime + starttime
352
353        #FIXME Use method in geo to reconcile
354        #if spatial:
355        #assert domain.geo_reference.xllcorner == xllcorner
356        #assert domain.geo_reference.yllcorner == yllcorner
357        #assert domain.geo_reference.zone == zone
358
359    if verbose:
360        print 'File_function data obtained from: %s' %filename
361        print '  References:'
362        #print '    Datum ....' #FIXME
363        if spatial:
364            print '    Lower left corner: [%f, %f]'\
365                  %(xllcorner, yllcorner)
366        print '    Start time:   %f' %starttime
367
368
369    #Produce values for desired data points at
370    #each timestep
371
372    quantities = {}
373    for i, name in enumerate(quantity_names):
374        quantities[name] = fid.variables[name][:]
375    fid.close()
376
377
378    from least_squares import Interpolation_function
379
380    if not spatial:
381        vertex_coordinates = triangles = interpolation_points = None
382
383    return Interpolation_function(time,
384                                  quantities,
385                                  quantity_names,
386                                  vertex_coordinates,
387                                  triangles,
388                                  interpolation_points,
389                                  verbose = verbose)
390
391
392
393
394
396    """Read simple xya file, possibly with a header in the first line, just
397    x y [attributes]
398    separated by whitespace
399
400    Format can be either ASCII or NetCDF
401
402    Return list of points, list of attributes and
403    attribute names if present otherwise None
404    """
406    #FIXME: Probably obsoleted by similar function in load_ASCII
408
409
411
412    points_dict = import_points_file(filename)
413    return points_dict['pointlist'], points_dict['attributelist']
414
415
416#####################################
417#POLYGON STUFF
418#
419#FIXME: All these should be put into new module polygon.py
420
421
422
423#These have been redefined to use separate_by_polygon which will
424#put all inside indices in the first part of the indices array and
425#outside indices in the last
426
427def inside_polygon(points, polygon, closed = True, verbose = False):
428    """See separate_points_by_polygon for documentation
429    """
430
431    from Numeric import array, Float, reshape
432
433    if verbose: print 'Checking input to inside_polygon'
434    try:
435        points = ensure_numeric(points, Float)
436    except:
437        msg = 'Points could not be converted to Numeric array'
438        raise msg
439
440    try:
441        polygon = ensure_numeric(polygon, Float)
442    except:
443        msg = 'Polygon could not be converted to Numeric array'
444        raise msg
445
446
447
448    if len(points.shape) == 1:
449        one_point = True
450        points = reshape(points, (1,2))
451    else:
452        one_point = False
453
454    indices, count = separate_points_by_polygon(points, polygon,
455                                                closed, verbose)
456
457    if one_point:
458        return count == 1
459    else:
460        return indices[:count]
461
462def outside_polygon(points, polygon, closed = True, verbose = False):
463    """See separate_points_by_polygon for documentation
464    """
465
466    from Numeric import array, Float, reshape
467
468    if verbose: print 'Checking input to outside_polygon'
469    try:
470        points = ensure_numeric(points, Float)
471    except:
472        msg = 'Points could not be converted to Numeric array'
473        raise msg
474
475    try:
476        polygon = ensure_numeric(polygon, Float)
477    except:
478        msg = 'Polygon could not be converted to Numeric array'
479        raise msg
480
481
482
483    if len(points.shape) == 1:
484        one_point = True
485        points = reshape(points, (1,2))
486    else:
487        one_point = False
488
489    indices, count = separate_points_by_polygon(points, polygon,
490                                                closed, verbose)
491
492    if one_point:
493        return count != 1
494    else:
495        return indices[count:][::-1]  #return reversed
496
497
498def separate_points_by_polygon(points, polygon,
499                               closed = True, verbose = False):
500    """Determine whether points are inside or outside a polygon
501
502    Input:
503       points - Tuple of (x, y) coordinates, or list of tuples
504       polygon - list of vertices of polygon
505       closed - (optional) determine whether points on boundary should be
506       regarded as belonging to the polygon (closed = True)
507       or not (closed = False)
508
509    Outputs:
510       indices: array of same length as points with indices of points falling
511       inside the polygon listed from the beginning and indices of points
512       falling outside listed from the end.
513
514       count: count of points falling inside the polygon
515
516       The indices of points inside are obtained as indices[:count]
517       The indices of points outside are obtained as indices[count:]
518
519
520    Examples:
521       U = [[0,0], [1,0], [1,1], [0,1]] #Unit square
522
523       separate_points_by_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U)
524       will return the indices [0, 2, 1] and count == 2 as only the first
525       and the last point are inside the unit square
526
527    Remarks:
528       The vertices may be listed clockwise or counterclockwise and
529       the first point may optionally be repeated.
530       Polygons do not need to be convex.
531       Polygons can have holes in them and points inside a hole is
532       regarded as being outside the polygon.
533
534    Algorithm is based on work by Darel Finley,
535    http://www.alienryderflex.com/polygon/
536
537    """
538
539    from Numeric import array, Float, reshape, Int, zeros
540
541
542    #Input checks
543    try:
544        points = ensure_numeric(points, Float)
545    except:
546        msg = 'Points could not be converted to Numeric array'
547        raise msg
548
549    try:
550        polygon = ensure_numeric(polygon, Float)
551    except:
552        msg = 'Polygon could not be converted to Numeric array'
553        raise msg
554
555    assert len(polygon.shape) == 2,\
556       'Polygon array must be a 2d array of vertices'
557
558    assert polygon.shape == 2,\
559       'Polygon array must have two columns'
560
561    assert len(points.shape) == 2,\
562       'Points array must be a 2d array'
563
564    assert points.shape == 2,\
565       'Points array must have two columns'
566
567    N = polygon.shape #Number of vertices in polygon
568    M = points.shape  #Number of points
569
570    px = polygon[:,0]
571    py = polygon[:,1]
572
573    #Used for an optimisation when points are far away from polygon
574    minpx = min(px); maxpx = max(px)
575    minpy = min(py); maxpy = max(py)
576
577
578    #Begin main loop
579    indices = zeros(M, Int)
580
581    inside_index = 0    #Keep track of points inside
582    outside_index = M-1 #Keep track of points outside (starting from end)
583
584    for k in range(M):
585
586        if verbose:
587            if k %((M+10)/10)==0: print 'Doing %d of %d' %(k, M)
588
589        #for each point
590        x = points[k, 0]
591        y = points[k, 1]
592
593        inside = False
594
595        if not x > maxpx or x < minpx or y > maxpy or y < minpy:
596            #Check polygon
597            for i in range(N):
598                j = (i+1)%N
599
600                #Check for case where point is contained in line segment
601                if point_on_line(x, y, px[i], py[i], px[j], py[j]):
602                    if closed:
603                        inside = True
604                    else:
605                        inside = False
606                    break
607                else:
608                    #Check if truly inside polygon
609                    if py[i] < y and py[j] >= y or\
610                       py[j] < y and py[i] >= y:
611                        if px[i] + (y-py[i])/(py[j]-py[i])*(px[j]-px[i]) < x:
612                            inside = not inside
613
614        if inside:
615            indices[inside_index] = k
616            inside_index += 1
617        else:
618            indices[outside_index] = k
619            outside_index -= 1
620
621    return indices, inside_index
622
623
624def separate_points_by_polygon_c(points, polygon,
625                                 closed = True, verbose = False):
626    """Determine whether points are inside or outside a polygon
627
628    C-wrapper
629    """
630
631    from Numeric import array, Float, reshape, zeros, Int
632
633
634    if verbose: print 'Checking input to separate_points_by_polygon'
635    #Input checks
636    try:
637        points = ensure_numeric(points, Float)
638    except:
639        msg = 'Points could not be converted to Numeric array'
640        raise msg
641
642    #if verbose: print 'Checking input to separate_points_by_polygon 2'
643    try:
644        polygon = ensure_numeric(polygon, Float)
645    except:
646        msg = 'Polygon could not be converted to Numeric array'
647        raise msg
648
649    if verbose: print 'check'
650
651    assert len(polygon.shape) == 2,\
652       'Polygon array must be a 2d array of vertices'
653
654    assert polygon.shape == 2,\
655       'Polygon array must have two columns'
656
657    assert len(points.shape) == 2,\
658       'Points array must be a 2d array'
659
660    assert points.shape == 2,\
661       'Points array must have two columns'
662
663    N = polygon.shape #Number of vertices in polygon
664    M = points.shape  #Number of points
665
666    from util_ext import separate_points_by_polygon
667
668    if verbose: print 'Allocating array for indices'
669
670    indices = zeros( M, Int )
671
672    if verbose: print 'Calling C-version of inside poly'
673    count = separate_points_by_polygon(points, polygon, indices,
674                                       int(closed), int(verbose))
675
676    return indices, count
677
678
679
680class Polygon_function:
681    """Create callable object f: x,y -> z, where a,y,z are vectors and
682    where f will return different values depending on whether x,y belongs
683    to specified polygons.
684
685    To instantiate:
686
687       Polygon_function(polygons)
688
689    where polygons is a list of tuples of the form
690
691      [ (P0, v0), (P1, v1), ...]
692
693      with Pi being lists of vertices defining polygons and vi either
694      constants or functions of x,y to be applied to points with the polygon.
695
696    The function takes an optional argument, default which is the value
697    (or function) to used for points not belonging to any polygon.
698    For example:
699
700       Polygon_function(polygons, default = 0.03)
701
702    If omitted the default value will be 0.0
703
704    Note: If two polygons overlap, the one last in the list takes precedence
705
706    FIXME? : Currently, coordinates specified here are assumed to be relative to the origin (georeference) used by domain. This function is more general than domain so perhaps it'd be an idea to allow the optional argument georeference???
707
708    """
709
710    def __init__(self, regions, default = 0.0):
711
712        try:
713            len(regions)
714        except:
715            msg = 'Polygon_function takes a list of pairs (polygon, value). Got %s' %polygons
716            raise msg
717
718
719        T = regions
720        try:
721            a = len(T)
722        except:
723            msg = 'Polygon_function takes a list of pairs (polygon, value). Got %s' %polygons
724            raise msg
725
726        assert a == 2, 'Must have two component each: %s' %T
727
728        self.regions = regions
729        self.default = default
730
731
732    def __call__(self, x, y):
733        from util import inside_polygon
734        from Numeric import ones, Float, concatenate, array, reshape, choose
735
736        x = array(x).astype(Float)
737        y = array(y).astype(Float)
738
739        N = len(x)
740        assert len(y) == N
741
742        points = concatenate( (reshape(x, (N, 1)),
743                               reshape(y, (N, 1))), axis=1 )
744
745        if callable(self.default):
746            z = self.default(x,y)
747        else:
748            z = ones(N, Float) * self.default
749
750        for polygon, value in self.regions:
751            indices = inside_polygon(points, polygon)
752
753            #FIXME: This needs to be vectorised
754            if callable(value):
755                for i in indices:
756                    xx = array([x[i]])
757                    yy = array([y[i]])
758                    z[i] = value(xx, yy)
759            else:
760                for i in indices:
761                    z[i] = value
762
763        return z
764
766    """Read points assumed to form a polygon
767       There must be exactly two numbers in each line
768    """
769
770    #Get polygon
771    fid = open(filename)
773    fid.close()
774    polygon = []
775    for line in lines:
776        fields = line.split(',')
777        polygon.append( [float(fields), float(fields)] )
778
779    return polygon
780
781def populate_polygon(polygon, number_of_points, seed = None, exclude = None):
782    """Populate given polygon with uniformly distributed points.
783
784    Input:
785       polygon - list of vertices of polygon
786       number_of_points - (optional) number of points
787       seed - seed for random number generator (default=None)
788       exclude - list of polygons (inside main polygon) from where points should be excluded
789
790    Output:
791       points - list of points inside polygon
792
793    Examples:
794       populate_polygon( [[0,0], [1,0], [1,1], [0,1]], 5 )
795       will return five randomly selected points inside the unit square
796    """
797
798    from random import uniform, seed
799
800    seed(seed)
801
802    points = []
803
804    #Find outer extent of polygon
805    max_x = min_x = polygon
806    max_y = min_y = polygon
807    for point in polygon[1:]:
808        x = point
809        if x > max_x: max_x = x
810        if x < min_x: min_x = x
811        y = point
812        if y > max_y: max_y = y
813        if y < min_y: min_y = y
814
815
816    while len(points) < number_of_points:
817        x = uniform(min_x, max_x)
818        y = uniform(min_y, max_y)
819
820        append = False
821        if inside_polygon( [x,y], polygon ):
822
823            append = True
824
825            #Check exclusions
826            if exclude is not None:
827                for ex_poly in exclude:
828                    if inside_polygon( [x,y], ex_poly ):
829                        append = False
830
831
832        if append is True:
833            points.append([x,y])
834
835    return points
836
837
838
839def tsh2sww(filename, verbose=True):
840    """
841    to check if a tsh/msh file 'looks' good.
842    """
843
844    #FIXME Move to data_manager
845    from shallow_water import Domain
846    from pmesh2domain import pmesh_to_domain_instance
847    import time, os
848    from data_manager import get_dataobject
849    from os import sep, path
850    #from util import mean
851
852    if verbose == True:print 'Creating domain from', filename
853    domain = pmesh_to_domain_instance(filename, Domain)
854    if verbose == True:print "Number of triangles = ", len(domain)
855
856    domain.smooth = True
857    domain.format = 'sww'   #Native netcdf visualisation format
858    file_path, filename = path.split(filename)
859    filename, ext = path.splitext(filename)
860    domain.filename = filename
861    domain.reduction = mean
862    if verbose == True:print "file_path",file_path
863    if file_path == "":file_path = "."
865
866    if verbose == True:
867        print "Output written to " + domain.get_datadir() + sep + \
868              domain.filename + "." + domain.format
869    sww = get_dataobject(domain)
870    sww.store_connectivity()
871    sww.store_timestep('stage')
872
873####################################################################
874#Python versions of function that are also implemented in util_gateway.c
875#
876
877def gradient_python(x0, y0, x1, y1, x2, y2, q0, q1, q2):
878    """
879    """
880
881    det = (y2-y0)*(x1-x0) - (y1-y0)*(x2-x0)
882    a = (y2-y0)*(q1-q0) - (y1-y0)*(q2-q0)
883    a /= det
884
885    b = (x1-x0)*(q2-q0) - (x2-x0)*(q1-q0)
886    b /= det
887
888    return a, b
889
890
891def gradient2_python(x0, y0, x1, y1, q0, q1):
893    in the direction orthogonal to (x1-x0), (y1-y0)
894    """
895
896    #Old code
897    #det = x0*y1 - x1*y0
898    #if det != 0.0:
899    #    a = (y1*q0 - y0*q1)/det
900    #    b = (x0*q1 - x1*q0)/det
901
902    #Correct code (ON)
903    det = (x1-x0)**2 + (y1-y0)**2
904    if det != 0.0:
905        a = (x1-x0)*(q1-q0)/det
906        b = (y1-y0)*(q1-q0)/det
907
908    return a, b
909
910
911##############################################
912#Initialise module
913
914import compile
915if compile.can_use_C_extension('util_ext.c'):