Changeset 5403

Ignore:
Timestamp:
Jun 13, 2008, 3:48:16 PM (15 years ago)
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

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

 r5372 point_coordinates=None, verbose=False): """Interpolate linearly between values f on nodes (vertex coordinates) of a polyline to midpoints of triangles of boundary. f is the data on the polyline nodes. The mesh values representing a smooth surface are assumed to be specified in f. Inputs: f: Vector or array of data at the polyline nodes. If f is an array, interpolation will be done for each column as per underlying matrix-matrix multiplication point_coordinates: Interpolate polyline data to these positions. List of coordinate pairs [x, y] of data points or an nx2 Numeric array or a Geospatial_data object Output: Interpolated values at inputted points (z). """ # FIXME: There is an option of passing a tolerance into this if isinstance(point_coordinates, Geospatial_data): point_coordinates = point_coordinates.get_data_points( \ absolute = True) from utilities.polygon import point_on_line,point_on_line_py from utilities.polygon import point_on_line from Numeric import ones z=ones(len(point_coordinates),Float) assert point_coordinates is not None, msg msg='function value must be specified at every interpolation node' assert f.shape[0]==vertex_coordinates.shape[0],msg assert f.shape[0]==vertex_coordinates.shape[0], msg msg='Must define function value at one or more nodes' assert f.shape[0]>0,msg assert f.shape[0]>0, msg n=f.shape[0] if n==1: z=f*z msg = 'Polyline contained only one point. I need more. ', str(f) raise Exception, msg # FIXME (John): add unit test for only One vertex point. Exception should be thrown. elif n>1: for i in range(len(point_coordinates)): for j in range(n): if gauge_neighbour_id[j]>=0: if point_on_line_py(point_coordinates[i],[vertex_coordinates[j],vertex_coordinates[gauge_neighbour_id[j]]]): if point_on_line(point_coordinates[i], [vertex_coordinates[j], vertex_coordinates[gauge_neighbour_id[j]]], rtol=1.0e-6): found=True x0=vertex_coordinates[j][0]
• anuga_core/source/anuga/shallow_water/data_manager.py

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

 r5372 from anuga.geospatial_data.geospatial_data import ensure_absolute def point_on_line_py(point, line): from Numeric import fabs point = ensure_numeric(point) line = ensure_numeric(line) x=point[0];y=point[1] x0=line[0,0];y0=line[0,1] x1=line[1,0];y1=line[1,1] #from pylab import plot,show #plot(line[:,0],line[:,1]) #plot([x],[y],'o') #show() a0 = x - x0; a1 = y - y0; a_normal0 = a1; a_normal1 = -a0; b0 = x1 - x0; b1 = y1 - y0; #if ( a_normal0*b0 + a_normal1*b1 == 0.0 ): eps=200 #print 'normal',a_normal0*b0 + a_normal1*b1 if ( fabs(a_normal0*b0 + a_normal1*b1) < eps ): #for some reason (perhaps catastrophic cancellation) urs_boundary.py #example only works for eps=2 #point is somewhere on the infinite extension of the line # FIXME (Ole): Perhaps add a tolerance here instead of 0.0 len_a = sqrt(a0*a0 + a1*a1); len_b = sqrt(b0*b0 + b1*b1); if (a0*b0 + a1*b1 >= 0 and len_a <= len_b): return 1 else: return 0 else: return 0 def point_on_line(point, line): def point_on_line(point, line, rtol=0.0, atol=0.0): """Determine whether a point is on a line segment Output: """ Note: Line can be degenerate and function still works to discern coinciding points from non-coinciding. """ # FIXME(Ole): Perhaps make defaults as in allclose: rtol=1.0e-5, atol=1.0e-8 point = ensure_numeric(point) line = ensure_numeric(line) res = _point_on_line(point[0], point[1], line[0,0], line[0,1], line[1,0], line[1,1]) line[1,0], line[1,1], rtol, atol) return bool(res)
• anuga_core/source/anuga/utilities/polygon_ext.c

 r5321 int __point_on_line(double x, double y, double x0, double y0, double x1, double y1) { double x1, double y1, double rtol, double atol) { /*Determine whether a point is on a line segment double a0, a1, a_normal0, a_normal1, b0, b1, len_a, len_b; double nominator, denominator; int is_parallel; a0 = x - x0; b1 = y1 - y0; if ( a_normal0*b0 + a_normal1*b1 == 0.0 ) { //Point is somewhere on the infinite extension of the line // FIXME (Ole): Perhaps add a tolerance here instead of 0.0 nominator = fabs(a_normal0*b0 + a_normal1*b1); denominator = b0*b0 + b1*b1; // Determine if line is parallel to point vector up to a tolerance is_parallel = 0; if (denominator == 0.0) { // Use absolute tolerance if (nominator <= atol) { is_parallel = 1; } } else { // Denominator is positive - use relative tolerance if (nominator/denominator <= rtol) { is_parallel = 1; } } if (is_parallel) { // Point is somewhere on the infinite extension of the line // subject to specified absolute tolerance len_a = sqrt(a0*a0 + a1*a1); int verbose) { double minpx, maxpx, minpy, maxpy, x, y, px_i, py_i, px_j, py_j; double minpx, maxpx, minpy, maxpy, x, y, px_i, py_i, px_j, py_j, rtol=0.0, atol=0.0; int i, j, k, outside_index, inside_index, inside; //Check for case where point is contained in line segment if (__point_on_line(x, y, px_i, py_i, px_j, py_j)) { if (__point_on_line(x, y, px_i, py_i, px_j, py_j, rtol, atol)) { if (closed == 1) { inside = 1; // double x, y, x0, y0, x1, y1; double x, y, x0, y0, x1, y1, rtol, atol; int res; PyObject *result; // Convert Python arguments to C if (!PyArg_ParseTuple(args, "dddddd", &x, &y, &x0, &y0, &x1, &y1)) { if (!PyArg_ParseTuple(args, "dddddddd", &x, &y, &x0, &y0, &x1, &y1, &rtol, &atol)) { PyErr_SetString(PyExc_RuntimeError, "point_on_line could not parse input"); // Call underlying routine res = __point_on_line(x, y, x0, y0, x1, y1); res = __point_on_line(x, y, x0, y0, x1, y1, rtol, atol); // Return values a and b
• anuga_core/source/anuga/utilities/test_polygon.py

 r5288 assert not point_on_line( [40,19], [[40,20], [40,40]] ) # Degenerate line assert point_on_line( [40,19], [[40,19], [40,19]] ) assert not point_on_line( [40,19], [[40,40], [40,40]] )
Note: See TracChangeset for help on using the changeset viewer.