[1093] | 1 | // Python - C extension for finite_volumes util module. |
---|
| 2 | // |
---|
| 3 | // To compile (Python2.3): |
---|
| 4 | // gcc -c util_ext.c -I/usr/include/python2.3 -o util_ext.o -Wall -O |
---|
| 5 | // gcc -shared util_ext.o -o util_ext.so |
---|
| 6 | // |
---|
| 7 | // See the module util.py |
---|
| 8 | // |
---|
| 9 | // |
---|
| 10 | // Ole Nielsen, GA 2004 |
---|
| 11 | // |
---|
| 12 | //NOTE: On 64 bit systems use long* instead of int* for Numeric arrays |
---|
| 13 | //this will also work on 32 bit systems |
---|
| 14 | |
---|
| 15 | |
---|
| 16 | #include "Python.h" |
---|
| 17 | #include "Numeric/arrayobject.h" |
---|
| 18 | #include "math.h" |
---|
| 19 | |
---|
| 20 | //Shared snippets |
---|
| 21 | #include "util_ext.h" |
---|
| 22 | |
---|
| 23 | |
---|
| 24 | |
---|
| 25 | int _point_on_line(double x, double y, |
---|
| 26 | double x0, double y0, |
---|
| 27 | double x1, double y1) { |
---|
| 28 | /*Determine whether a point is on a line segment |
---|
| 29 | |
---|
| 30 | Input: x, y, x0, x0, x1, y1: where |
---|
| 31 | point is given by x, y |
---|
| 32 | line is given by (x0, y0) and (x1, y1) |
---|
| 33 | |
---|
| 34 | */ |
---|
| 35 | |
---|
| 36 | double a0, a1, a_normal0, a_normal1, b0, b1, len_a, len_b; |
---|
| 37 | |
---|
| 38 | a0 = x - x0; |
---|
| 39 | a1 = y - y0; |
---|
| 40 | |
---|
| 41 | a_normal0 = a1; |
---|
| 42 | a_normal1 = -a0; |
---|
| 43 | |
---|
| 44 | b0 = x1 - x0; |
---|
| 45 | b1 = y1 - y0; |
---|
| 46 | |
---|
| 47 | if ( a_normal0*b0 + a_normal1*b1 == 0 ) { |
---|
| 48 | //Point is somewhere on the infinite extension of the line |
---|
| 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 && len_a <= len_b) { |
---|
| 54 | return 1; |
---|
| 55 | } else { |
---|
| 56 | return 0; |
---|
| 57 | } |
---|
| 58 | } else { |
---|
| 59 | return 0; |
---|
| 60 | } |
---|
| 61 | } |
---|
| 62 | |
---|
| 63 | |
---|
| 64 | |
---|
| 65 | int _separate_points_by_polygon(int M, // Number of points |
---|
| 66 | int N, // Number of polygon vertices |
---|
| 67 | double* points, |
---|
| 68 | double* polygon, |
---|
| 69 | long* indices, // M-Array for storage indices |
---|
| 70 | int closed, |
---|
| 71 | int verbose) { |
---|
| 72 | |
---|
| 73 | double minpx, maxpx, minpy, maxpy, x, y, px_i, py_i, px_j, py_j; |
---|
| 74 | int i, j, k, outside_index, inside_index, inside; |
---|
| 75 | |
---|
| 76 | //Find min and max of poly used for optimisation when points |
---|
| 77 | //are far away from polygon |
---|
| 78 | |
---|
| 79 | minpx = polygon[0]; maxpx = minpx; |
---|
| 80 | minpy = polygon[1]; maxpy = minpy; |
---|
| 81 | |
---|
| 82 | for (i=0; i<N; i++) { |
---|
| 83 | px_i = polygon[2*i]; |
---|
| 84 | py_i = polygon[2*i + 1]; |
---|
| 85 | |
---|
| 86 | if (px_i < minpx) minpx = px_i; |
---|
| 87 | if (px_i > maxpx) maxpx = px_i; |
---|
| 88 | if (py_i < minpy) minpy = py_i; |
---|
| 89 | if (py_i > maxpy) maxpy = py_i; |
---|
| 90 | } |
---|
| 91 | |
---|
| 92 | //Begin main loop (for each point) |
---|
| 93 | inside_index = 0; //Keep track of points inside |
---|
| 94 | outside_index = M-1; //Keep track of points outside (starting from end) |
---|
| 95 | for (k=0; k<M; k++) { |
---|
| 96 | if (verbose){ |
---|
| 97 | if (k %((M+10)/10)==0) printf("Doing %d of %d\n", k, M); |
---|
| 98 | } |
---|
| 99 | |
---|
| 100 | x = points[2*k]; |
---|
| 101 | y = points[2*k + 1]; |
---|
| 102 | |
---|
| 103 | inside = 0; |
---|
| 104 | |
---|
| 105 | //Optimisation |
---|
| 106 | if ((x > maxpx) || (x < minpx) || (y > maxpy) || (y < minpy)) { |
---|
| 107 | //Nothing |
---|
| 108 | } else { |
---|
| 109 | //Check polygon |
---|
| 110 | for (i=0; i<N; i++) { |
---|
| 111 | //printf("k,i=%d,%d\n", k, i); |
---|
| 112 | j = (i+1)%N; |
---|
| 113 | |
---|
| 114 | px_i = polygon[2*i]; |
---|
| 115 | py_i = polygon[2*i+1]; |
---|
| 116 | px_j = polygon[2*j]; |
---|
| 117 | py_j = polygon[2*j+1]; |
---|
| 118 | |
---|
| 119 | //Check for case where point is contained in line segment |
---|
| 120 | if (_point_on_line(x, y, px_i, py_i, px_j, py_j)) { |
---|
| 121 | if (closed == 1) { |
---|
| 122 | inside = 1; |
---|
| 123 | } else { |
---|
| 124 | inside = 0; |
---|
| 125 | } |
---|
| 126 | break; |
---|
| 127 | } else { |
---|
| 128 | //Check if truly inside polygon |
---|
| 129 | if ( ((py_i < y) && (py_j >= y)) || |
---|
| 130 | ((py_j < y) && (py_i >= y)) ) { |
---|
| 131 | if (px_i + (y-py_i)/(py_j-py_i)*(px_j-px_i) < x) |
---|
| 132 | inside = 1-inside; |
---|
| 133 | } |
---|
| 134 | } |
---|
| 135 | } |
---|
| 136 | } |
---|
| 137 | if (inside == 1) { |
---|
| 138 | indices[inside_index] = k; |
---|
| 139 | inside_index += 1; |
---|
| 140 | } else { |
---|
| 141 | indices[outside_index] = k; |
---|
| 142 | outside_index -= 1; |
---|
| 143 | } |
---|
| 144 | } // End k |
---|
| 145 | |
---|
| 146 | return inside_index; |
---|
| 147 | } |
---|
| 148 | |
---|
| 149 | |
---|
| 150 | |
---|
| 151 | // Gateways to Python |
---|
| 152 | PyObject *point_on_line(PyObject *self, PyObject *args) { |
---|
| 153 | // |
---|
| 154 | // point_on_line(x, y, x0, y0, x1, y1) |
---|
| 155 | // |
---|
| 156 | |
---|
| 157 | double x, y, x0, y0, x1, y1; |
---|
| 158 | int res; |
---|
| 159 | PyObject *result; |
---|
| 160 | |
---|
| 161 | // Convert Python arguments to C |
---|
| 162 | if (!PyArg_ParseTuple(args, "dddddd", &x, &y, &x0, &y0, &x1, &y1)) |
---|
| 163 | return NULL; |
---|
| 164 | |
---|
| 165 | |
---|
| 166 | // Call underlying routine |
---|
| 167 | res = _point_on_line(x, y, x0, y0, x1, y1); |
---|
| 168 | |
---|
| 169 | // Return values a and b |
---|
| 170 | result = Py_BuildValue("i", res); |
---|
| 171 | return result; |
---|
| 172 | } |
---|
| 173 | |
---|
| 174 | |
---|
| 175 | |
---|
| 176 | PyObject *separate_points_by_polygon(PyObject *self, PyObject *args) { |
---|
| 177 | //def separate_points_by_polygon(points, polygon, closed, verbose, one_point): |
---|
| 178 | // """Determine whether points are inside or outside a polygon |
---|
| 179 | // |
---|
| 180 | // Input: |
---|
| 181 | // point - Tuple of (x, y) coordinates, or list of tuples |
---|
| 182 | // polygon - list of vertices of polygon |
---|
| 183 | // closed - (optional) determine whether points on boundary should be |
---|
| 184 | // regarded as belonging to the polygon (closed = True) |
---|
| 185 | // or not (closed = False) |
---|
| 186 | |
---|
| 187 | // |
---|
| 188 | // Output: |
---|
| 189 | // indices: array of same length as points with indices of points falling |
---|
| 190 | // inside the polygon listed from the beginning and indices of points |
---|
| 191 | // falling outside listed from the end. |
---|
| 192 | // |
---|
| 193 | // count: count of points falling inside the polygon |
---|
| 194 | // |
---|
| 195 | // The indices of points inside are obtained as indices[:count] |
---|
| 196 | // The indices of points outside are obtained as indices[count:] |
---|
| 197 | // |
---|
| 198 | // Examples: |
---|
| 199 | // separate_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]] ) |
---|
| 200 | // will return the indices [0, 2, 1] as only the first and the last point |
---|
| 201 | // is inside the unit square |
---|
| 202 | // |
---|
| 203 | // Remarks: |
---|
| 204 | // The vertices may be listed clockwise or counterclockwise and |
---|
| 205 | // the first point may optionally be repeated. |
---|
| 206 | // Polygons do not need to be convex. |
---|
| 207 | // Polygons can have holes in them and points inside a hole is |
---|
| 208 | // regarded as being outside the polygon. |
---|
| 209 | // |
---|
| 210 | // |
---|
| 211 | // Algorithm is based on work by Darel Finley, |
---|
| 212 | // http://www.alienryderflex.com/polygon/ |
---|
| 213 | // |
---|
| 214 | // |
---|
| 215 | |
---|
| 216 | PyArrayObject |
---|
| 217 | *points, |
---|
| 218 | *polygon, |
---|
| 219 | *indices; |
---|
| 220 | |
---|
| 221 | int closed, verbose; //Flags |
---|
| 222 | int count, M, N; |
---|
| 223 | |
---|
| 224 | // Convert Python arguments to C |
---|
| 225 | if (!PyArg_ParseTuple(args, "OOOii", |
---|
| 226 | &points, |
---|
| 227 | &polygon, |
---|
| 228 | &indices, |
---|
| 229 | &closed, |
---|
| 230 | &verbose)) |
---|
| 231 | return NULL; |
---|
| 232 | |
---|
| 233 | M = points -> dimensions[0]; //Number of points |
---|
| 234 | N = polygon -> dimensions[0]; //Number of vertices in polygon |
---|
| 235 | |
---|
| 236 | if (verbose) printf("Got %d points and %d polygon vertices\n", M, N); |
---|
| 237 | |
---|
| 238 | //Call underlying routine |
---|
| 239 | count = _separate_points_by_polygon(M, N, |
---|
| 240 | (double*) points -> data, |
---|
| 241 | (double*) polygon -> data, |
---|
| 242 | (long*) indices -> data, |
---|
| 243 | closed, verbose); |
---|
| 244 | |
---|
| 245 | //NOTE: return number of points inside.. |
---|
| 246 | return Py_BuildValue("i", count); |
---|
| 247 | } |
---|
| 248 | |
---|
| 249 | |
---|
| 250 | |
---|
| 251 | PyObject *gradient(PyObject *self, PyObject *args) { |
---|
| 252 | // |
---|
| 253 | // a,b = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2) |
---|
| 254 | // |
---|
| 255 | |
---|
| 256 | double x0, y0, x1, y1, x2, y2, q0, q1, q2, a, b; |
---|
| 257 | PyObject *result; |
---|
| 258 | |
---|
| 259 | // Convert Python arguments to C |
---|
| 260 | if (!PyArg_ParseTuple(args, "ddddddddd", &x0, &y0, &x1, &y1, &x2, &y2, |
---|
| 261 | &q0, &q1, &q2)) |
---|
| 262 | return NULL; |
---|
| 263 | |
---|
| 264 | |
---|
| 265 | // Call underlying routine |
---|
| 266 | _gradient(x0, y0, x1, y1, x2, y2, |
---|
| 267 | q0, q1, q2, &a, &b); |
---|
| 268 | |
---|
| 269 | // Return values a and b |
---|
| 270 | result = Py_BuildValue("dd", a, b); |
---|
| 271 | return result; |
---|
| 272 | } |
---|
| 273 | |
---|
| 274 | |
---|
| 275 | // Method table for python module |
---|
| 276 | static struct PyMethodDef MethodTable[] = { |
---|
| 277 | /* The cast of the function is necessary since PyCFunction values |
---|
| 278 | * only take two PyObject* parameters, and rotate() takes |
---|
| 279 | * three. |
---|
| 280 | */ |
---|
| 281 | |
---|
| 282 | //{"rotate", (PyCFunction)rotate, METH_VARARGS | METH_KEYWORDS, "Print out"}, |
---|
| 283 | {"gradient", gradient, METH_VARARGS, "Print out"}, |
---|
| 284 | {"point_on_line", point_on_line, METH_VARARGS, "Print out"}, |
---|
| 285 | //{"inside_polygon", inside_polygon, METH_VARARGS, "Print out"}, |
---|
| 286 | {"separate_points_by_polygon", separate_points_by_polygon, |
---|
| 287 | METH_VARARGS, "Print out"}, |
---|
| 288 | {NULL, NULL, 0, NULL} /* sentinel */ |
---|
| 289 | }; |
---|
| 290 | |
---|
| 291 | |
---|
| 292 | |
---|
| 293 | // Module initialisation |
---|
| 294 | void initutil_ext(void){ |
---|
| 295 | Py_InitModule("util_ext", MethodTable); |
---|
| 296 | |
---|
| 297 | import_array(); //Necessary for handling of NumPY structures |
---|
| 298 | } |
---|
| 299 | |
---|
| 300 | |
---|
| 301 | |
---|
| 302 | |
---|
| 303 | //OBSOLETE |
---|
| 304 | |
---|
| 305 | int _inside_polygon(int M, // Number of points |
---|
| 306 | int N, // Number of polygon vertices |
---|
| 307 | double* points, |
---|
| 308 | double* polygon, |
---|
| 309 | long* indices, // M-Array for storage indices |
---|
| 310 | int closed, |
---|
| 311 | int verbose) { |
---|
| 312 | |
---|
| 313 | double minpx, maxpx, minpy, maxpy, x, y, px_i, py_i, px_j, py_j; |
---|
| 314 | int i, j, k, count, inside; |
---|
| 315 | |
---|
| 316 | //Find min and max of poly used for optimisation when points |
---|
| 317 | //are far away from polygon |
---|
| 318 | |
---|
| 319 | minpx = polygon[0]; maxpx = minpx; |
---|
| 320 | minpy = polygon[1]; maxpy = minpy; |
---|
| 321 | |
---|
| 322 | for (i=0; i<N; i++) { |
---|
| 323 | px_i = polygon[2*i]; |
---|
| 324 | py_i = polygon[2*i + 1]; |
---|
| 325 | |
---|
| 326 | if (px_i < minpx) minpx = px_i; |
---|
| 327 | if (px_i > maxpx) maxpx = px_i; |
---|
| 328 | if (py_i < minpy) minpy = py_i; |
---|
| 329 | if (py_i > maxpy) maxpy = py_i; |
---|
| 330 | } |
---|
| 331 | |
---|
| 332 | //Begin main loop (for each point) |
---|
| 333 | count = 0; |
---|
| 334 | for (k=0; k<M; k++) { |
---|
| 335 | if (verbose){ |
---|
| 336 | if (k %((M+10)/10)==0) printf("Doing %d of %d\n", k, M); |
---|
| 337 | } |
---|
| 338 | |
---|
| 339 | x = points[2*k]; |
---|
| 340 | y = points[2*k + 1]; |
---|
| 341 | |
---|
| 342 | inside = 0; |
---|
| 343 | |
---|
| 344 | //Optimisation |
---|
| 345 | if ((x > maxpx) || (x < minpx)) continue; |
---|
| 346 | if ((y > maxpy) || (y < minpy)) continue; |
---|
| 347 | |
---|
| 348 | //Check polygon |
---|
| 349 | for (i=0; i<N; i++) { |
---|
| 350 | //printf("k,i=%d,%d\n", k, i); |
---|
| 351 | j = (i+1)%N; |
---|
| 352 | |
---|
| 353 | px_i = polygon[2*i]; |
---|
| 354 | py_i = polygon[2*i+1]; |
---|
| 355 | px_j = polygon[2*j]; |
---|
| 356 | py_j = polygon[2*j+1]; |
---|
| 357 | |
---|
| 358 | //Check for case where point is contained in line segment |
---|
| 359 | if (_point_on_line(x, y, px_i, py_i, px_j, py_j)) { |
---|
| 360 | if (closed == 1) { |
---|
| 361 | inside = 1; |
---|
| 362 | } else { |
---|
| 363 | inside = 0; |
---|
| 364 | } |
---|
| 365 | break; |
---|
| 366 | } else { |
---|
| 367 | //Check if truly inside polygon |
---|
| 368 | if ( ((py_i < y) && (py_j >= y)) || |
---|
| 369 | ((py_j < y) && (py_i >= y)) ) { |
---|
| 370 | if (px_i + (y-py_i)/(py_j-py_i)*(px_j-px_i) < x) |
---|
| 371 | inside = 1-inside; |
---|
| 372 | } |
---|
| 373 | } |
---|
| 374 | } |
---|
| 375 | if (inside == 1) { |
---|
| 376 | indices[count] = k; |
---|
| 377 | count++; |
---|
| 378 | } |
---|
| 379 | } // End k |
---|
| 380 | |
---|
| 381 | return count; |
---|
| 382 | } |
---|
| 383 | |
---|
| 384 | |
---|
| 385 | |
---|
| 386 | PyObject *inside_polygon(PyObject *self, PyObject *args) { |
---|
| 387 | //def inside_polygon(point, polygon, closed, verbose, one_point): |
---|
| 388 | // """Determine whether points are inside or outside a polygon |
---|
| 389 | // |
---|
| 390 | // Input: |
---|
| 391 | // point - Tuple of (x, y) coordinates, or list of tuples |
---|
| 392 | // polygon - list of vertices of polygon |
---|
| 393 | // one_poin - If True Boolean value should be returned |
---|
| 394 | // If False, indices of points inside returned |
---|
| 395 | // closed - (optional) determine whether points on boundary should be |
---|
| 396 | // regarded as belonging to the polygon (closed = True) |
---|
| 397 | // or not (closed = False) |
---|
| 398 | |
---|
| 399 | // |
---|
| 400 | // Output: |
---|
| 401 | // If one point is considered, True or False is returned. |
---|
| 402 | // If multiple points are passed in, the function returns indices |
---|
| 403 | // of those points that fall inside the polygon |
---|
| 404 | // |
---|
| 405 | // Examples: |
---|
| 406 | // inside_polygon( [0.5, 0.5], [[0,0], [1,0], [1,1], [0,1]] ) |
---|
| 407 | // will evaluate to True as the point 0.5, 0.5 is inside the unit square |
---|
| 408 | // |
---|
| 409 | // inside_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]] ) |
---|
| 410 | // will return the indices [0, 2] as only the first and the last point |
---|
| 411 | // is inside the unit square |
---|
| 412 | // |
---|
| 413 | // Remarks: |
---|
| 414 | // The vertices may be listed clockwise or counterclockwise and |
---|
| 415 | // the first point may optionally be repeated. |
---|
| 416 | // Polygons do not need to be convex. |
---|
| 417 | // Polygons can have holes in them and points inside a hole is |
---|
| 418 | // regarded as being outside the polygon. |
---|
| 419 | // |
---|
| 420 | // |
---|
| 421 | // Algorithm is based on work by Darel Finley, |
---|
| 422 | // http://www.alienryderflex.com/polygon/ |
---|
| 423 | // |
---|
| 424 | // |
---|
| 425 | |
---|
| 426 | PyArrayObject |
---|
| 427 | *point, |
---|
| 428 | *polygon, |
---|
| 429 | *indices; |
---|
| 430 | |
---|
| 431 | int closed, verbose; //Flags |
---|
| 432 | int count, M, N; |
---|
| 433 | |
---|
| 434 | // Convert Python arguments to C |
---|
| 435 | if (!PyArg_ParseTuple(args, "OOOii", |
---|
| 436 | &point, |
---|
| 437 | &polygon, |
---|
| 438 | &indices, |
---|
| 439 | &closed, |
---|
| 440 | &verbose)) |
---|
| 441 | return NULL; |
---|
| 442 | |
---|
| 443 | M = point -> dimensions[0]; //Number of points |
---|
| 444 | N = polygon -> dimensions[0]; //Number of vertices in polygon |
---|
| 445 | |
---|
| 446 | if (verbose) printf("Got %d points and %d polygon vertices\n", M, N); |
---|
| 447 | |
---|
| 448 | //Call underlying routine |
---|
| 449 | count = _inside_polygon(M, N, |
---|
| 450 | (double*) point -> data, |
---|
| 451 | (double*) polygon -> data, |
---|
| 452 | (long*) indices -> data, |
---|
| 453 | closed, verbose); |
---|
| 454 | |
---|
| 455 | //NOTE: return number of points inside.. |
---|
| 456 | //printf("COunt=%d\n", count); |
---|
| 457 | return Py_BuildValue("i", count); |
---|
| 458 | } |
---|