Changeset 5403


Ignore:
Timestamp:
Jun 13, 2008, 3:48:16 PM (14 years ago)
Author:
ole
Message:

Added tolerances for point_on_line and changed interpolate_polyline to use a general relative tolerance.

Location:
anuga_core/source/anuga
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/fit_interpolate/interpolate.py

    r5372 r5403  
    9696                             point_coordinates=None,
    9797                             verbose=False):
    98  
     98        """Interpolate linearly between values f on nodes (vertex coordinates) of a polyline to midpoints of triangles
     99        of boundary.
     100
     101        f is the data on the polyline nodes.
     102
     103        The mesh values representing a smooth surface are
     104        assumed to be specified in f.
     105
     106        Inputs:
     107          f: Vector or array of data at the polyline nodes.
     108              If f is an array, interpolation will be done for each column as
     109              per underlying matrix-matrix multiplication
     110          point_coordinates: Interpolate polyline data to these positions.
     111              List of coordinate pairs [x, y] of
     112              data points or an nx2 Numeric array or a Geospatial_data object
     113             
     114        Output:
     115          Interpolated values at inputted points (z).
     116        """
     117       
     118
     119        # FIXME: There is an option of passing a tolerance into this
     120       
    99121        if isinstance(point_coordinates, Geospatial_data):
    100122            point_coordinates = point_coordinates.get_data_points( \
    101123                absolute = True)
    102124 
    103         from utilities.polygon import point_on_line,point_on_line_py
     125        from utilities.polygon import point_on_line
    104126        from Numeric import ones
    105127        z=ones(len(point_coordinates),Float)
     
    108130        assert point_coordinates is not None, msg
    109131        msg='function value must be specified at every interpolation node'
    110         assert f.shape[0]==vertex_coordinates.shape[0],msg
     132        assert f.shape[0]==vertex_coordinates.shape[0], msg
    111133        msg='Must define function value at one or more nodes'
    112         assert f.shape[0]>0,msg
     134        assert f.shape[0]>0, msg
    113135
    114136        n=f.shape[0]
    115137        if n==1:
    116138            z=f*z
     139            msg = 'Polyline contained only one point. I need more. ', str(f)
     140            raise Exception, msg
     141           
     142        # FIXME (John): add unit test for only One vertex point. Exception should be thrown.
     143
     144       
    117145        elif n>1:
    118146            for i in range(len(point_coordinates)):
     
    120148                for j in range(n):
    121149                    if gauge_neighbour_id[j]>=0:
    122                         if point_on_line_py(point_coordinates[i],[vertex_coordinates[j],vertex_coordinates[gauge_neighbour_id[j]]]):
     150                        if point_on_line(point_coordinates[i],
     151                                         [vertex_coordinates[j], vertex_coordinates[gauge_neighbour_id[j]]],
     152                                         rtol=1.0e-6):
    123153                            found=True
    124154                            x0=vertex_coordinates[j][0]
  • anuga_core/source/anuga/shallow_water/data_manager.py

    r5391 r5403  
    54135413    #Title
    54145414    fid.write('location_index'+d+'lat'+d+ 'long' +d+ 'Easting' +d+ \
    5415                   'Northing' +d+ 'depth m' + "\n")
     5415              'Northing' +d+ 'depth m' + "\n")
    54165416    i = 0
    54175417    for depth, point_utm, lat, long in map(None, depths,
  • anuga_core/source/anuga/utilities/polygon.py

    r5372 r5403  
    1717from anuga.geospatial_data.geospatial_data import ensure_absolute
    1818
    19 def point_on_line_py(point, line):
    20     from Numeric import fabs
    21     point = ensure_numeric(point)
    22     line = ensure_numeric(line)
    23 
    24 
    25     x=point[0];y=point[1]
    26     x0=line[0,0];y0=line[0,1]
    27     x1=line[1,0];y1=line[1,1]
    28     #from pylab import plot,show
    29     #plot(line[:,0],line[:,1])
    30     #plot([x],[y],'o')
    31     #show()
    32     a0 = x - x0;
    33     a1 = y - y0;
    34 
    35     a_normal0 = a1;
    36     a_normal1 = -a0;
    37 
    38     b0 = x1 - x0;
    39     b1 = y1 - y0;
    40 
    41     #if ( a_normal0*b0 + a_normal1*b1 == 0.0 ):
    42     eps=200
    43     #print 'normal',a_normal0*b0 + a_normal1*b1
    44     if ( fabs(a_normal0*b0 + a_normal1*b1) < eps ):
    45         #for some reason (perhaps catastrophic cancellation) urs_boundary.py
    46         #example only works for eps=2
    47         #point is somewhere on the infinite extension of the line
    48         # FIXME (Ole): Perhaps add a tolerance here instead of 0.0
    49 
    50         len_a = sqrt(a0*a0 + a1*a1);
    51         len_b = sqrt(b0*b0 + b1*b1);
    52 
    53         if (a0*b0 + a1*b1 >= 0 and len_a <= len_b):
    54             return 1
    55         else:
    56             return 0
    57     else:
    58         return 0
    59 
    60 
    61 def point_on_line(point, line):
     19
     20def point_on_line(point, line, rtol=0.0, atol=0.0):
    6221    """Determine whether a point is on a line segment
    6322
     
    6827
    6928    Output:
    70        
    71     """
     29
     30    Note: Line can be degenerate and function still works to discern coinciding points from non-coinciding.
     31    """
     32
     33    # FIXME(Ole): Perhaps make defaults as in allclose: rtol=1.0e-5, atol=1.0e-8
    7234
    7335    point = ensure_numeric(point)
    7436    line = ensure_numeric(line)
    7537
    76 
    7738    res = _point_on_line(point[0], point[1],
    7839                         line[0,0], line[0,1],
    79                          line[1,0], line[1,1])
     40                         line[1,0], line[1,1],
     41                         rtol, atol)
    8042   
    8143    return bool(res)
  • anuga_core/source/anuga/utilities/polygon_ext.c

    r5321 r5403  
    2121int __point_on_line(double x, double y,
    2222                    double x0, double y0,
    23                     double x1, double y1) {
     23                    double x1, double y1,
     24                    double rtol,
     25                    double atol) {
    2426  /*Determine whether a point is on a line segment
    2527
     
    3133
    3234  double a0, a1, a_normal0, a_normal1, b0, b1, len_a, len_b;
     35  double nominator, denominator;
     36  int is_parallel;
    3337
    3438  a0 = x - x0;
     
    4145  b1 = y1 - y0;
    4246
    43   if ( a_normal0*b0 + a_normal1*b1 == 0.0 ) {
    44     //Point is somewhere on the infinite extension of the line
    45     // FIXME (Ole): Perhaps add a tolerance here instead of 0.0
     47  nominator = fabs(a_normal0*b0 + a_normal1*b1);
     48  denominator = b0*b0 + b1*b1;
     49 
     50  // Determine if line is parallel to point vector up to a tolerance
     51  is_parallel = 0;
     52  if (denominator == 0.0) {
     53    // Use absolute tolerance
     54    if (nominator <= atol) {
     55      is_parallel = 1;
     56    }
     57  } else {
     58    // Denominator is positive - use relative tolerance
     59    if (nominator/denominator <= rtol) {
     60      is_parallel = 1;
     61    }   
     62  }
     63   
     64  if (is_parallel) {
     65    // Point is somewhere on the infinite extension of the line
     66    // subject to specified absolute tolerance
    4667
    4768    len_a = sqrt(a0*a0 + a1*a1);
     
    173194                                int verbose) {
    174195
    175   double minpx, maxpx, minpy, maxpy, x, y, px_i, py_i, px_j, py_j;
     196  double minpx, maxpx, minpy, maxpy, x, y, px_i, py_i, px_j, py_j, rtol=0.0, atol=0.0;
    176197  int i, j, k, outside_index, inside_index, inside;
    177198
     
    223244
    224245        //Check for case where point is contained in line segment
    225         if (__point_on_line(x, y, px_i, py_i, px_j, py_j)) {
     246        if (__point_on_line(x, y, px_i, py_i, px_j, py_j, rtol, atol)) {
    226247          if (closed == 1) {
    227248            inside = 1;
     
    260281  //
    261282
    262   double x, y, x0, y0, x1, y1;
     283  double x, y, x0, y0, x1, y1, rtol, atol;
    263284  int res;
    264285  PyObject *result;
    265286
    266287  // Convert Python arguments to C
    267   if (!PyArg_ParseTuple(args, "dddddd", &x, &y, &x0, &y0, &x1, &y1)) {
     288  if (!PyArg_ParseTuple(args, "dddddddd", &x, &y, &x0, &y0, &x1, &y1, &rtol, &atol)) {
    268289    PyErr_SetString(PyExc_RuntimeError,
    269290                    "point_on_line could not parse input");   
     
    273294
    274295  // Call underlying routine
    275   res = __point_on_line(x, y, x0, y0, x1, y1);
     296  res = __point_on_line(x, y, x0, y0, x1, y1, rtol, atol);
    276297
    277298  // Return values a and b
  • anuga_core/source/anuga/utilities/test_polygon.py

    r5288 r5403  
    165165        assert not point_on_line( [40,19], [[40,20], [40,40]] )
    166166
     167        # Degenerate line
     168        assert point_on_line( [40,19], [[40,19], [40,19]] )
     169        assert not point_on_line( [40,19], [[40,40], [40,40]] )       
    167170
    168171
Note: See TracChangeset for help on using the changeset viewer.