Changeset 6535 for anuga_core/source/anuga/utilities/polygon_ext.c
- Timestamp:
- Mar 17, 2009, 11:06:55 PM (15 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/utilities/polygon_ext.c
r6534 r6535 244 244 245 245 246 247 int __is_inside_triangle(double* point, 248 double* triangle, 249 int closed, 250 double rtol, 251 double atol) { 252 253 double vx, vy, v0x, v0y, v1x, v1y; 254 double a00, a10, a01, a11, b0, b1; 255 double denom, alpha, beta; 256 int i, j, res; 257 258 // v0 = C-A 259 v0x = triangle[4]-triangle[0]; 260 v0y = triangle[5]-triangle[1]; 261 262 // v1 = B-A 263 v1x = triangle[2]-triangle[0]; 264 v1y = triangle[3]-triangle[1]; 265 266 // First check if point lies wholly inside triangle 267 a00 = v0x*v0x + v0y*v0y; // innerproduct(v0, v0) 268 a01 = v0x*v1x + v0y*v1y; // innerproduct(v0, v1) 269 a10 = a01; // innerproduct(v1, v0) 270 a11 = v1x*v1x + v1y*v1y; // innerproduct(v1, v1) 271 272 denom = a11*a00 - a01*a10; 273 274 if (fabs(denom) > 0.0) { 275 // v = point-A 276 vx = point[0] - triangle[0]; 277 vy = point[1] - triangle[1]; 278 279 b0 = v0x*vx + v0y*vy; // innerproduct(v0, v) 280 b1 = v1x*vx + v1y*vy; // innerproduct(v1, v) 281 282 alpha = (b0*a11 - b1*a01)/denom; 283 beta = (b1*a00 - b0*a10)/denom; 284 285 if ((alpha > 0.0) && (beta > 0.0) && (alpha+beta < 1.0)) return 1; 286 } 287 288 if (closed) { 289 // Check if point lies on one of the edges 290 291 for (i=0; i<3; i++) { 292 j = (i+1) % 3; // Circular index into triangle array 293 res = __point_on_line(point[0], point[1], 294 triangle[2*i], triangle[2*i+1], 295 triangle[2*j], triangle[2*j+1], 296 rtol, atol); 297 if (res) return 1; 298 } 299 } 300 301 // Default return if point is outside triangle 302 return 0; 303 } 304 305 246 306 int __separate_points_by_polygon(int M, // Number of points 247 307 int N, // Number of polygon vertices … … 255 315 int i, j, k, outside_index, inside_index, inside; 256 316 257 // Find min and max of poly used for optimisation when points258 // are far away from polygon259 260 // FIXME(Ole): Pass in rtol and atol from Python317 // Find min and max of poly used for optimisation when points 318 // are far away from polygon 319 320 // FIXME(Ole): Pass in rtol and atol from Python 261 321 262 322 minpx = polygon[0]; maxpx = minpx; … … 273 333 } 274 334 275 // Begin main loop (for each point)276 inside_index = 0; // Keep track of points inside277 outside_index = M-1; // Keep track of points outside (starting from end)335 // Begin main loop (for each point) 336 inside_index = 0; // Keep track of points inside 337 outside_index = M-1; // Keep track of points outside (starting from end) 278 338 if (verbose){ 279 339 printf("Separating %d points\n", M); … … 289 349 inside = 0; 290 350 291 // Optimisation351 // Optimisation 292 352 if ((x > maxpx) || (x < minpx) || (y > maxpy) || (y < minpy)) { 293 // Nothing353 // Nothing 294 354 } else { 295 // Check polygon355 // Check polygon 296 356 for (i=0; i<N; i++) { 297 //printf("k,i=%d,%d\n", k, i);298 357 j = (i+1)%N; 299 358 … … 303 362 py_j = polygon[2*j+1]; 304 363 305 // Check for case where point is contained in line segment364 // Check for case where point is contained in line segment 306 365 if (__point_on_line(x, y, px_i, py_i, px_j, py_j, rtol, atol)) { 307 366 if (closed == 1) { … … 312 371 break; 313 372 } else { 314 // Check if truly inside polygon373 // Check if truly inside polygon 315 374 if ( ((py_i < y) && (py_j >= y)) || 316 375 ((py_j < y) && (py_i >= y)) ) { … … 356 415 res = __point_on_line(x, y, x0, y0, x1, y1, rtol, atol); 357 416 358 // Return values a and b417 // Return result 359 418 result = Py_BuildValue("i", res); 360 419 return result; … … 417 476 418 477 478 479 480 PyObject *_is_inside_triangle(PyObject *self, PyObject *args) { 481 // 482 // _is_inside_triangle(point, triangle, int(closed), rtol, atol) 483 // 484 485 486 PyArrayObject 487 *point, 488 *triangle; 489 490 double rtol, atol; 491 int closed, res; 492 493 PyObject *result; 494 495 // Convert Python arguments to C 496 if (!PyArg_ParseTuple(args, "OOidd", 497 &point, 498 &triangle, 499 &closed, 500 &rtol, 501 &atol)) { 502 503 PyErr_SetString(PyExc_RuntimeError, 504 "_is_inside_triangle could not parse input"); 505 return NULL; 506 } 507 508 // Call underlying routine 509 res = __is_inside_triangle((double*) point -> data, 510 (double*) triangle -> data, 511 closed, 512 rtol, 513 atol); 514 515 516 // Return result 517 result = Py_BuildValue("i", res); 518 return result; 519 } 520 521 522 419 523 /* 420 524 PyObject *_intersection(PyObject *self, PyObject *args) { … … 539 643 {"_interpolate_polyline", _interpolate_polyline, 540 644 METH_VARARGS, "Print out"}, 645 {"_is_inside_triangle", _is_inside_triangle, 646 METH_VARARGS, "Print out"}, 541 647 {NULL, NULL, 0, NULL} /* sentinel */ 542 648 };
Note: See TracChangeset
for help on using the changeset viewer.