Changeset 62045
- Timestamp:
- Sep 20, 2014, 7:46:53 AM (10 years ago)
- Location:
- grass/trunk
- Files:
-
- 1 edited
- 1 copied
-
include/defs/vector.h (modified) (1 diff)
-
lib/vector/Vlib/intersect2.c (copied) (copied from grass/trunk/lib/vector/Vlib/intersect.c ) (19 diffs)
Legend:
- Unmodified
- Added
- Removed
-
grass/trunk/include/defs/vector.h
r59594 r62045 440 440 struct line_pnts ***, struct line_pnts ***, int *, 441 441 int *, int); 442 int Vect_line_intersection2(struct line_pnts *, struct line_pnts *, 443 struct bound_box *, struct bound_box *, 444 struct line_pnts ***, struct line_pnts ***, int *, 445 int *, int); 442 446 int Vect_line_check_intersection(struct line_pnts *, struct line_pnts *, int); 447 int Vect_line_check_intersection2(struct line_pnts *, struct line_pnts *, int); 443 448 int Vect_line_get_intersections(struct line_pnts *, struct line_pnts *, 444 449 struct line_pnts *, int); 450 int Vect_line_get_intersections2(struct line_pnts *, struct line_pnts *, 451 struct line_pnts *, int); 445 452 char *Vect_subst_var(const char *, const struct Map_info *); 446 453 -
grass/trunk/lib/vector/Vlib/intersect2.c
r61910 r62045 57 57 coordinates, the results would be different) 58 58 59 (C) 2001-20 09by the GRASS Development Team59 (C) 2001-2014 by the GRASS Development Team 60 60 61 61 This program is free software under the GNU General Public License … … 64 64 \author Original author CERL, probably Dave Gerdes or Mike Higgins. 65 65 \author Update to GRASS 5.7 Radim Blazek. 66 \author Update to GRASS 7 Markus Metz. 66 67 */ 67 68 … … 71 72 #include <math.h> 72 73 #include <grass/vector.h> 74 #include <grass/rbtree.h> 73 75 #include <grass/glocale.h> 74 76 … … 84 86 static int ident(double x1, double y1, double x2, double y2, double thresh); 85 87 #endif 86 static int cross_seg(int id, const struct RTree_Rect *rect, int *arg); 87 static int find_cross(int id, const struct RTree_Rect *rect, int *arg); 88 89 #define D ((ax2-ax1)*(by1-by2) - (ay2-ay1)*(bx1-bx2)) 90 #define D1 ((bx1-ax1)*(by1-by2) - (by1-ay1)*(bx1-bx2)) 91 #define D2 ((ax2-ax1)*(by1-ay1) - (ay2-ay1)*(bx1-ax1)) 92 93 /*! 94 * \brief Check for intersect of 2 line segments. 95 * 96 * \param ax1,ay1,az1,ax2,ay2,az2 input line a 97 * \param bx1,by1,bz1,bx2,by2,bz2 input line b 98 * \param[out] x1,y1,z1 intersection point1 (case 2-4) 99 * \param[out] x2,y2,z2 intersection point2 (case 2-4) 100 * \param with_z use z coordinate (3D) (TODO) 101 * 102 * \return 0 - do not intersect, 103 * \return 1 - intersect at one point, 104 * \return 2 - partial overlap, 105 * \return 3 - a contains b, 106 * \return 4 - b contains a, 107 * \return 5 - identical 108 */ 109 int Vect_segment_intersection(double ax1, double ay1, double az1, double ax2, 110 double ay2, double az2, double bx1, double by1, 111 double bz1, double bx2, double by2, double bz2, 112 double *x1, double *y1, double *z1, double *x2, 113 double *y2, double *z2, int with_z) 114 { 115 static int first_3d = 1; 116 double d, d1, d2, r1, dtol, t; 117 int switched; 118 119 /* TODO: Works for points ? */ 120 121 G_debug(4, "Vect_segment_intersection()"); 122 G_debug(4, " %.15g , %.15g - %.15g , %.15g", ax1, ay1, ax2, ay2); 123 G_debug(4, " %.15g , %.15g - %.15g , %.15g", bx1, by1, bx2, by2); 124 125 /* TODO 3D */ 126 if (with_z && first_3d) { 127 G_warning(_("3D not supported by Vect_segment_intersection()")); 128 first_3d = 0; 129 } 130 131 /* Check identical segments */ 132 if ((ax1 == bx1 && ay1 == by1 && ax2 == bx2 && ay2 == by2) || 133 (ax1 == bx2 && ay1 == by2 && ax2 == bx1 && ay2 == by1)) { 134 G_debug(2, " -> identical segments"); 135 *x1 = ax1; 136 *y1 = ay1; 137 *z1 = az1; 138 *x2 = ax2; 139 *y2 = ay2; 140 *z2 = az2; 141 return 5; 142 } 143 144 /* 'Sort' lines by x, y 145 * MUST happen before D, D1, D2 are calculated */ 146 switched = 0; 147 if (bx2 < bx1) 148 switched = 1; 149 else if (bx2 == bx1) { 150 if (by2 < by1) 151 switched = 1; 152 } 153 154 if (switched) { 155 t = bx1; 156 bx1 = bx2; 157 bx2 = t; 158 t = by1; 159 by1 = by2; 160 by2 = t; 161 t = bz1; 162 bz1 = bz2; 163 bz2 = t; 164 } 165 166 switched = 0; 167 if (ax2 < ax1) 168 switched = 1; 169 else if (ax2 == ax1) { 170 if (ay2 < ay1) 171 switched = 1; 172 } 173 174 if (switched) { 175 t = ax1; 176 ax1 = ax2; 177 ax2 = t; 178 t = ay1; 179 ay1 = ay2; 180 ay2 = t; 181 t = az1; 182 az1 = az2; 183 az2 = t; 184 } 185 186 d = D; 187 d1 = D1; 188 d2 = D2; 189 190 G_debug(2, "Vect_segment_intersection(): d = %f, d1 = %f, d2 = %f", d, d1, 191 d2); 192 193 /* TODO: dtol was originaly set to 1.0e-10, which was usualy working but not always. 194 * Can it be a problem to set the tolerance to 0.0 ? */ 195 dtol = 0.0; 196 if (fabs(d) > dtol) { 197 198 G_debug(2, " -> not parallel/collinear: d1 = %f, d2 = %f", d1, d2); 199 if (d > 0) { 200 if (d1 < 0 || d1 > d || d2 < 0 || d2 > d) { 201 G_debug(2, " -> no intersection"); 202 return 0; 203 } 204 } 205 else { 206 if (d1 < d || d1 > 0 || d2 < d || d2 > 0) { 207 G_debug(2, " -> no intersection"); 208 return 0; 209 } 210 } 211 212 r1 = d1 / d; 213 214 *x1 = ax1 + r1 * (ax2 - ax1); 215 *y1 = ay1 + r1 * (ay2 - ay1); 216 *z1 = 0; 217 218 G_debug(2, " -> intersection %f, %f", *x1, *y1); 219 return 1; 220 } 221 222 /* segments are parallel or collinear */ 223 G_debug(3, " -> parallel/collinear"); 224 225 if (d1 || d2) { /* lines are parallel */ 226 G_debug(2, " -> parallel"); 227 return 0; 228 } 229 230 /* segments are colinear. check for overlap */ 231 232 /* Collinear vertical */ 233 /* original code assumed lines were not both vertical 234 * so there is a special case if they are */ 235 if (ax1 == ax2) { 236 G_debug(2, " -> collinear vertical"); 237 if (ay1 > by2 || ay2 < by1) { 238 G_debug(2, " -> no intersection"); 239 return 0; 240 } 241 242 /* end points */ 243 if (ay1 == by2) { 244 G_debug(2, " -> connected by end points"); 245 *x1 = ax1; 246 *y1 = ay1; 247 *z1 = 0; 248 249 return 1; /* endpoints only */ 250 } 251 if (ay2 == by1) { 252 G_debug(2, " -> connected by end points"); 253 *x1 = ax2; 254 *y1 = ay2; 255 *z1 = 0; 256 257 return 1; /* endpoints only */ 258 } 259 260 /* general overlap */ 261 G_debug(3, " -> vertical overlap"); 262 /* a contains b */ 263 if (ay1 <= by1 && ay2 >= by2) { 264 G_debug(2, " -> a contains b"); 265 *x1 = bx1; 266 *y1 = by1; 267 *z1 = 0; 268 *x2 = bx2; 269 *y2 = by2; 270 *z2 = 0; 271 272 return 3; 273 } 274 /* b contains a */ 275 if (ay1 >= by1 && ay2 <= by2) { 276 G_debug(2, " -> b contains a"); 277 *x1 = ax1; 278 *y1 = ay1; 279 *z1 = 0; 280 *x2 = ax2; 281 *y2 = ay2; 282 *z2 = 0; 283 284 return 4; 285 } 286 287 /* general overlap, 2 intersection points */ 288 G_debug(2, " -> partial overlap"); 289 if (by1 > ay1 && by1 < ay2) { /* b1 in a */ 290 G_debug(2, " -> b1 in a"); 291 *x1 = bx1; 292 *y1 = by1; 293 *z1 = 0; 294 *x2 = ax2; 295 *y2 = ay2; 296 *z2 = 0; 297 298 return 2; 299 } 300 if (by2 > ay1 && by2 < ay2) { /* b2 in a */ 301 G_debug(2, " -> b2 in a"); 302 *x1 = ax1; 303 *y1 = ay1; 304 *z1 = 0; 305 *x2 = bx2; 306 *y2 = by2; 307 *z2 = 0; 308 309 return 2; 310 } 311 312 /* should not be reached */ 313 G_warning(_("Vect_segment_intersection() ERROR (collinear vertical segments)")); 314 G_warning("a"); 315 G_warning("%.15g %.15g", ax1, ay1); 316 G_warning("%.15g %.15g", ax2, ay2); 317 G_warning("b"); 318 G_warning("%.15g %.15g", bx1, by1); 319 G_warning("%.15g %.15g", bx2, by2); 320 321 return 0; 322 } 323 324 /* Collinear non vertical */ 325 326 G_debug(2, " -> collinear non vertical"); 327 328 /* b is to the left or right of a */ 329 if ((bx1 > ax2) || (bx2 < ax1)) { 330 /* should not happen if segments are selected from rtree */ 331 G_debug(2, " -> no intersection"); 332 return 0; 333 } 334 335 /* there is overlap or connected end points */ 336 G_debug(2, " -> overlap/connected end points"); 337 338 /* end points */ 339 if (ax1 == bx2 && ay1 == by2) { 340 G_debug(2, " -> connected by end points"); 341 *x1 = ax1; 342 *y1 = ay1; 343 *z1 = 0; 344 345 return 1; 346 } 347 if (ax2 == bx1 && ay2 == by1) { 348 G_debug(2, " -> connected by end points"); 349 *x1 = ax2; 350 *y1 = ay2; 351 *z1 = 0; 352 353 return 1; 354 } 355 356 /* a contains b */ 357 if (ax1 <= bx1 && ax2 >= bx2) { 358 G_debug(2, " -> a contains b"); 359 *x1 = bx1; 360 *y1 = by1; 361 *z1 = 0; 362 *x2 = bx2; 363 *y2 = by2; 364 *z2 = 0; 365 366 return 3; 367 } 368 /* b contains a */ 369 if (ax1 >= bx1 && ax2 <= bx2) { 370 G_debug(2, " -> b contains a"); 371 *x1 = ax1; 372 *y1 = ay1; 373 *z1 = 0; 374 *x2 = ax2; 375 *y2 = ay2; 376 *z2 = 0; 377 378 return 4; 379 } 380 381 /* general overlap, 2 intersection points (lines are not vertical) */ 382 G_debug(2, " -> partial overlap"); 383 if (bx1 > ax1 && bx1 < ax2) { /* b1 is in a */ 384 G_debug(2, " -> b1 in a"); 385 *x1 = bx1; 386 *y1 = by1; 387 *z1 = 0; 388 *x2 = ax2; 389 *y2 = ay2; 390 *z2 = 0; 391 392 return 2; 393 } 394 if (bx2 > ax1 && bx2 < ax2) { /* b2 is in a */ 395 G_debug(2, " -> b2 in a"); 396 *x1 = ax1; 397 *y1 = ay1; 398 *z1 = 0; 399 *x2 = bx2; 400 *y2 = by2; 401 *z2 = 0; 402 403 return 2; 404 } 405 406 /* should not be reached */ 407 G_warning(_("Vect_segment_intersection() ERROR (collinear non vertical segments)")); 408 G_warning("a"); 409 G_warning("%.15g %.15g", ax1, ay1); 410 G_warning("%.15g %.15g", ax2, ay2); 411 G_warning("b"); 412 G_warning("%.15g %.15g", bx1, by1); 413 G_warning("%.15g %.15g", bx2, by2); 414 415 return 0; 416 } 88 static int cross_seg(int i, int j); 89 static int find_cross(int i, int j); 90 417 91 418 92 typedef struct … … 431 105 static CROSS *cross = NULL; 432 106 static int *use_cross = NULL; 107 108 static double rethresh = 0.000001; /* TODO */ 109 433 110 434 111 static void add_cross(int asegment, double adistance, int bsegment, … … 500 177 501 178 /* shared by Vect_line_intersection, Vect_line_check_intersection, cross_seg, find_cross */ 502 static struct line_pnts *APnts, *BPnts ;503 504 /* break segments (called by rtree search)*/505 static int cross_seg(int i d, const struct RTree_Rect *rect, int *arg)179 static struct line_pnts *APnts, *BPnts, *ABPnts[2], *IPnts; 180 181 /* break segments */ 182 static int cross_seg(int i, int j) 506 183 { 507 184 double x1, y1, z1, x2, y2, z2; 508 int i, j, ret; 509 510 /* !!! segment number for B lines is returned as +1 */ 511 i = *arg; 512 j = id - 1; 513 /* Note: -1 to make up for the +1 when data was inserted */ 185 double y1min, y1max, y2min, y2max; 186 int ret; 187 188 y1min = APnts->y[i]; 189 y1max = APnts->y[i + 1]; 190 if (APnts->y[i] > APnts->y[i + 1]) { 191 y1min = APnts->y[i + 1]; 192 y1max = APnts->y[i]; 193 } 194 195 y2min = BPnts->y[j]; 196 y2max = BPnts->y[j + 1]; 197 if (BPnts->y[j] > BPnts->y[j + 1]) { 198 y2min = BPnts->y[j + 1]; 199 y2max = BPnts->y[j]; 200 } 201 202 if (y1min > y2max || y1max < y2min) 203 return 0; 514 204 515 205 ret = Vect_segment_intersection(APnts->x[i], APnts->y[i], APnts->z[i], … … 540 230 } 541 231 232 /* event queue for Bentley-Ottmann */ 233 #define QEVT_IN 1 234 #define QEVT_OUT 2 235 #define QEVT_CRS 3 236 237 #define GET_PARENT(p, c) ((p) = (int) (((c) - 2) / 3 + 1)) 238 #define GET_CHILD(c, p) ((c) = (int) (((p) * 3) - 1)) 239 240 struct qitem 241 { 242 int l; /* line 0 - A line , 1 - B line */ 243 int s; /* segment index */ 244 int p; /* point index */ 245 int e; /* event type */ 246 }; 247 248 struct boq 249 { 250 int count; 251 int alloc; 252 struct qitem *i; 253 }; 254 255 /* compare two queue points */ 256 /* return 1 if a < b else 0 */ 257 static int cmp_q_x(struct qitem *a, struct qitem *b) 258 { 259 double x1, y1, z1, x2, y2, z2; 260 261 x1 = ABPnts[a->l]->x[a->p]; 262 y1 = ABPnts[a->l]->y[a->p]; 263 z1 = ABPnts[a->l]->z[a->p]; 264 265 x2 = ABPnts[b->l]->x[b->p]; 266 y2 = ABPnts[b->l]->y[b->p]; 267 z2 = ABPnts[b->l]->z[b->p]; 268 269 if (x1 < x2) 270 return 1; 271 if (x1 > x2) 272 return 0; 273 274 if (y1 < y2) 275 return 1; 276 if (y1 > y2) 277 return 0; 278 279 if (z1 < z2) 280 return 1; 281 if (z1 > z2) 282 return 0; 283 284 if (a->e < b->e) 285 return 1; 286 287 if (a->l < b->l) 288 return 1; 289 290 if (a->s < b->s) 291 return 1; 292 293 return 0; 294 } 295 296 /* sift up routine for min heap */ 297 int sift_up(struct boq *q, int start) 298 { 299 register int parent, child; 300 struct qitem a, *b; 301 302 child = start; 303 a = q->i[child]; 304 305 while (child > 1) { 306 GET_PARENT(parent, child); 307 308 b = &q->i[parent]; 309 /* child smaller */ 310 if (cmp_q_x(&a, b)) { 311 /* push parent point down */ 312 q->i[child] = q->i[parent]; 313 child = parent; 314 } 315 else 316 /* no more sifting up, found new slot for child */ 317 break; 318 } 319 320 /* put point in new slot */ 321 if (child < start) { 322 q->i[child] = a; 323 } 324 325 return child; 326 } 327 328 int boq_add(struct boq *q, struct qitem *i) 329 { 330 if (q->count + 2 >= q->alloc) { 331 q->alloc = q->count + 100; 332 q->i = G_realloc(q->i, q->alloc * sizeof(struct qitem)); 333 } 334 q->i[q->count + 1] = *i; 335 sift_up(q, q->count + 1); 336 337 q->count++; 338 339 return 1; 340 } 341 342 /* drop point routine for min heap */ 343 int boq_drop(struct boq *q, struct qitem *qi) 344 { 345 register int child, childr, parent; 346 register int i; 347 struct qitem *a, *b; 348 349 if (q->count == 0) 350 return 0; 351 352 *qi = q->i[1]; 353 354 if (q->count == 1) { 355 q->count = 0; 356 return 1; 357 } 358 359 /* start with root */ 360 parent = 1; 361 362 /* sift down: move hole back towards bottom of heap */ 363 364 while (GET_CHILD(child, parent) <= q->count) { 365 a = &q->i[child]; 366 i = child + 3; 367 for (childr = child + 1; childr <= q->count && childr < i; childr++) { 368 b = &q->i[childr]; 369 if (cmp_q_x(b, a)) { 370 child = childr; 371 a = &q->i[child]; 372 } 373 } 374 375 /* move hole down */ 376 q->i[parent] = q->i[child]; 377 parent = child; 378 } 379 380 /* hole is in lowest layer, move to heap end */ 381 if (parent < q->count) { 382 q->i[parent] = q->i[q->count]; 383 384 /* sift up last swapped point, only necessary if hole moved to heap end */ 385 sift_up(q, parent); 386 } 387 388 /* the actual drop */ 389 q->count--; 390 391 return 1; 392 } 393 394 /* compare two tree points */ 395 /* return -1 if a < b, 1 if a > b, 0 if a == b */ 396 static int cmp_t_y(const void *aa, const void *bb) 397 { 398 double x1, y1, z1, x2, y2, z2; 399 struct qitem *a = (struct qitem *) aa; 400 struct qitem *b = (struct qitem *) bb; 401 402 x1 = ABPnts[a->l]->x[a->p]; 403 y1 = ABPnts[a->l]->y[a->p]; 404 z1 = ABPnts[a->l]->z[a->p]; 405 406 x2 = ABPnts[b->l]->x[b->p]; 407 y2 = ABPnts[b->l]->y[b->p]; 408 z2 = ABPnts[b->l]->z[b->p]; 409 410 if (y1 < y2) 411 return -1; 412 if (y1 > y2) 413 return 1; 414 415 if (x1 < x2) 416 return -1; 417 if (x1 > x2) 418 return 1; 419 420 if (z1 < z2) 421 return -1; 422 if (z1 > z2) 423 return 1; 424 425 if (a->s < b->s) 426 return -1; 427 if (a->s > b->s) 428 return 1; 429 430 return 0; 431 } 432 433 int boq_load(struct boq *q, struct line_pnts *Pnts, 434 struct bound_box *abbox, int l, int with_z) 435 { 436 int i, loaded; 437 int vi, vo; 438 double x1, y1, z1, x2, y2, z2; 439 struct bound_box box; 440 struct qitem qi; 441 442 /* load Pnts to queue */ 443 qi.l = l; 444 loaded = 0; 445 446 for (i = 0; i < Pnts->n_points - 1; i++) { 447 x1 = Pnts->x[i]; 448 y1 = Pnts->y[i]; 449 z1 = Pnts->z[i]; 450 451 x2 = Pnts->x[i + 1]; 452 y2 = Pnts->y[i + 1]; 453 z2 = Pnts->z[i + 1]; 454 455 if (x1 == x2 && y1 == y2 && (!with_z || z1 == z2)) 456 continue; 457 458 if (x1 < x2) { 459 box.W = x1; 460 box.E = x2; 461 } 462 else { 463 box.E = x1; 464 box.W = x2; 465 } 466 if (y1 < y2) { 467 box.S = y1; 468 box.N = y2; 469 } 470 else { 471 box.N = y1; 472 box.S = y2; 473 } 474 if (z1 < z2) { 475 box.B = z1; 476 box.T = z2; 477 } 478 else { 479 box.T = z1; 480 box.B = z2; 481 } 482 box.W -= rethresh; 483 box.S -= rethresh; 484 box.B -= rethresh; 485 box.E += rethresh; 486 box.N += rethresh; 487 box.T += rethresh; 488 489 if (!Vect_box_overlap(abbox, &box)) 490 continue; 491 492 vi = i; 493 vo = i + 1; 494 495 if (x1 < x2) { 496 vi = i; 497 vo = i + 1; 498 } 499 else if (x1 > x2) { 500 vi = i + 1; 501 vo = i; 502 } 503 else { 504 if (y1 < y2) { 505 vi = i; 506 vo = i + 1; 507 } 508 else if (y1 > y2) { 509 vi = i + 1; 510 vo = i; 511 } 512 else { 513 if (z1 < z2) { 514 vi = i; 515 vo = i + 1; 516 } 517 else if (z1 > z2) { 518 vi = i + 1; 519 vo = i; 520 } 521 else { 522 G_fatal_error("Identical points"); 523 } 524 } 525 } 526 527 qi.s = i; 528 529 /* event in */ 530 qi.e = QEVT_IN; 531 qi.p = vi; 532 boq_add(q, &qi); 533 534 /* event out */ 535 qi.e = QEVT_OUT; 536 qi.p = vo; 537 boq_add(q, &qi); 538 539 loaded += 2; 540 } 541 542 return loaded; 543 } 544 542 545 /*! 543 546 * \brief Intersect 2 lines. … … 546 549 * intersection with B line. Points (Points->n_points == 1) are not 547 550 * supported. 551 * 552 * simplified Bentley–Ottmann Algorithm 548 553 * 549 554 * \param APoints first input line … … 559 564 */ 560 565 int 561 Vect_line_intersection (struct line_pnts *APoints,562 struct line_pnts *BPoints,563 struct bound_box *ABox,564 struct bound_box *BBox,565 struct line_pnts ***ALines,566 struct line_pnts ***BLines,567 int *nalines, int *nblines, int with_z)568 { 569 int i, j, k, l, last_seg, seg, last;566 Vect_line_intersection2(struct line_pnts *APoints, 567 struct line_pnts *BPoints, 568 struct bound_box *ABox, 569 struct bound_box *BBox, 570 struct line_pnts ***ALines, 571 struct line_pnts ***BLines, 572 int *nalines, int *nblines, int with_z) 573 { 574 int i, j, k, l, nl, last_seg, seg, last; 570 575 int n_alive_cross; 571 576 double dist, curdist, last_x, last_y, last_z; 572 double x, y , rethresh;577 double x, y; 573 578 struct line_pnts **XLines, *Points; 574 struct RTree *MyRTree;575 579 struct line_pnts *Points1, *Points2; /* first, second points */ 576 580 int seg1, seg2, vert1, vert2; 577 static struct RTree_Rect rect; 578 static int rect_init = 0; 579 struct bound_box box, abbox; 581 struct bound_box abbox; 582 struct boq bo_queue; 583 struct qitem qi, *found; 584 struct RB_TREE *bo_ta, *bo_tb; 585 struct RB_TRAV bo_t_trav; 586 int same = 0; 580 587 581 588 if (debug_level == -1) { … … 588 595 } 589 596 590 if (!rect_init) {591 rect.boundary = G_malloc(6 * sizeof(RectReal));592 rect_init = 6;593 }594 595 597 n_cross = 0; 596 rethresh = 0.000001; /* TODO */597 598 APnts = APoints; 598 599 BPnts = BPoints; 600 601 ABPnts[0] = APnts; 602 ABPnts[1] = BPnts; 603 604 *nalines = 0; 605 *nblines = 0; 599 606 600 607 /* RE (representation error). … … 653 660 */ 654 661 655 /* Spatial index: lines may be very long (thousands of segments) and check each segment656 * with each from second line takes a long time (n*m). Because of that, spatial index657 * is build first for the second line and segments from the first line are broken by segments658 * in bound box */659 660 662 if (!Vect_box_overlap(ABox, BBox)) { 661 663 *nalines = 0; … … 663 665 return 0; 664 666 } 665 667 668 /* overlap box of A line and B line */ 666 669 abbox = *BBox; 667 670 if (abbox.N > ABox->N) … … 684 687 abbox.T += rethresh; 685 688 abbox.B -= rethresh; 686 687 /* Create rtree for B line */ 688 MyRTree = RTreeCreateTree(-1, 0, 2); 689 RTreeSetOverflow(MyRTree, 0); 690 for (i = 0; i < BPoints->n_points - 1; i++) { 691 if (BPoints->x[i] <= BPoints->x[i + 1]) { 692 rect.boundary[0] = BPoints->x[i]; 693 rect.boundary[3] = BPoints->x[i + 1]; 694 } 695 else { 696 rect.boundary[0] = BPoints->x[i + 1]; 697 rect.boundary[3] = BPoints->x[i]; 698 } 699 700 if (BPoints->y[i] <= BPoints->y[i + 1]) { 701 rect.boundary[1] = BPoints->y[i]; 702 rect.boundary[4] = BPoints->y[i + 1]; 703 } 704 else { 705 rect.boundary[1] = BPoints->y[i + 1]; 706 rect.boundary[4] = BPoints->y[i]; 707 } 708 709 if (BPoints->z[i] <= BPoints->z[i + 1]) { 710 rect.boundary[2] = BPoints->z[i]; 711 rect.boundary[5] = BPoints->z[i + 1]; 712 } 713 else { 714 rect.boundary[2] = BPoints->z[i + 1]; 715 rect.boundary[5] = BPoints->z[i]; 716 } 717 718 box.W = rect.boundary[0] - rethresh; 719 box.S = rect.boundary[1] - rethresh; 720 box.B = rect.boundary[2] - rethresh; 721 box.E = rect.boundary[3] + rethresh; 722 box.N = rect.boundary[4] + rethresh; 723 box.T = rect.boundary[5] + rethresh; 724 725 if (Vect_box_overlap(&abbox, &box)) { 726 RTreeInsertRect(&rect, i + 1, MyRTree); /* B line segment numbers in rtree start from 1 */ 727 } 728 } 729 730 /* Break segments in A by segments in B */ 731 for (i = 0; i < APoints->n_points - 1; i++) { 732 if (APoints->x[i] <= APoints->x[i + 1]) { 733 rect.boundary[0] = APoints->x[i]; 734 rect.boundary[3] = APoints->x[i + 1]; 735 } 736 else { 737 rect.boundary[0] = APoints->x[i + 1]; 738 rect.boundary[3] = APoints->x[i]; 739 } 740 741 if (APoints->y[i] <= APoints->y[i + 1]) { 742 rect.boundary[1] = APoints->y[i]; 743 rect.boundary[4] = APoints->y[i + 1]; 744 } 745 else { 746 rect.boundary[1] = APoints->y[i + 1]; 747 rect.boundary[4] = APoints->y[i]; 748 } 749 if (APoints->z[i] <= APoints->z[i + 1]) { 750 rect.boundary[2] = APoints->z[i]; 751 rect.boundary[5] = APoints->z[i + 1]; 752 } 753 else { 754 rect.boundary[2] = APoints->z[i + 1]; 755 rect.boundary[5] = APoints->z[i]; 756 } 757 box.W = rect.boundary[0] - rethresh; 758 box.S = rect.boundary[1] - rethresh; 759 box.B = rect.boundary[2] - rethresh; 760 box.E = rect.boundary[3] + rethresh; 761 box.N = rect.boundary[4] + rethresh; 762 box.T = rect.boundary[5] + rethresh; 763 764 if (Vect_box_overlap(&abbox, &box)) { 765 j = RTreeSearch(MyRTree, &rect, (void *)cross_seg, &i); /* A segment number from 0 */ 766 } 767 } 768 769 /* Free RTree */ 770 RTreeDestroyTree(MyRTree); 689 690 if (APnts->n_points < 2 || BPnts->n_points < 2) { 691 G_fatal_error("Intersection with points is not yet supported"); 692 return 0; 693 } 694 695 if (APnts->n_points == BPnts->n_points) { 696 same = 1; 697 for (i = 0; i < APnts->n_points; i++) { 698 if (APnts->x[i] != BPnts->x[i] || 699 APnts->y[i] != BPnts->y[i] || 700 (with_z && APnts->z[i] != BPnts->z[i])) { 701 same = 0; 702 break; 703 } 704 } 705 if (same) 706 G_debug(3, "Intersecting different lines"); 707 } 708 709 /* initialize queue */ 710 bo_queue.count = 0; 711 bo_queue.alloc = 2 * (APnts->n_points + BPnts->n_points); 712 bo_queue.i = G_malloc(bo_queue.alloc * sizeof(struct qitem)); 713 714 /* load APnts to queue */ 715 boq_load(&bo_queue, APnts, &abbox, 0, with_z); 716 717 if (!same) { 718 /* load BPnts to queue */ 719 boq_load(&bo_queue, BPnts, &abbox, 1, with_z); 720 } 721 722 /* initialize search tree */ 723 bo_ta = rbtree_create(cmp_t_y, sizeof(struct qitem)); 724 if (same) 725 bo_tb = bo_ta; 726 else 727 bo_tb = rbtree_create(cmp_t_y, sizeof(struct qitem)); 728 729 /* find intersections */ 730 while (boq_drop(&bo_queue, &qi)) { 731 732 if (qi.e == QEVT_IN) { 733 /* not the original Bentley-Ottmann algorithm */ 734 if (qi.l == 0) { 735 /* test for intersection of s with all segments in T */ 736 rbtree_init_trav(&bo_t_trav, bo_tb); 737 while ((found = rbtree_traverse(&bo_t_trav))) { 738 cross_seg(qi.s, found->s); 739 } 740 741 /* insert s in T */ 742 rbtree_insert(bo_ta, &qi); 743 } 744 else { 745 /* test for intersection of s with all segments in T */ 746 rbtree_init_trav(&bo_t_trav, bo_ta); 747 while ((found = rbtree_traverse(&bo_t_trav))) { 748 cross_seg(found->s, qi.s); 749 } 750 751 /* insert s in T */ 752 rbtree_insert(bo_tb, &qi); 753 } 754 } 755 else if (qi.e == QEVT_OUT) { 756 /* remove from T */ 757 758 /* adjust p */ 759 if (qi.p == qi.s) 760 qi.p++; 761 else 762 qi.p--; 763 764 if (qi.l == 0) { 765 if (!rbtree_remove(bo_ta, &qi)) 766 G_fatal_error("RB tree error!"); 767 } 768 else { 769 if (!rbtree_remove(bo_tb, &qi)) 770 G_fatal_error("RB tree error!"); 771 } 772 } 773 } 774 G_free(bo_queue.i); 775 rbtree_destroy(bo_ta); 776 if (!same) 777 rbtree_destroy(bo_tb); 771 778 772 779 G_debug(2, "n_cross = %d", n_cross); 773 780 /* Lines do not cross each other */ 774 781 if (n_cross == 0) { 775 *nalines = 0;776 *nblines = 0;777 782 return 0; 778 783 } … … 834 839 835 840 /* l = 1 ~ line A, l = 2 ~ line B */ 836 for (l = 1; l < 3; l++) { 841 nl = 3; 842 if (same) 843 nl = 2; 844 for (l = 1; l < nl; l++) { 837 845 for (i = 0; i < n_cross; i++) 838 846 use_cross[i] = 1; … … 1089 1097 Points->y[last_seg], 1090 1098 Points->y[last_seg + 1]); 1091 last_z = (Points->z[last_seg] * sqrt(cross[i].distance[current]) + 1092 Points->z[last_seg + 1] * (sqrt(dist) - sqrt(cross[i].distance[current]))) / 1093 sqrt(dist); 1099 if (with_z) { 1100 last_z = (Points->z[last_seg] * sqrt(cross[i].distance[current]) + 1101 Points->z[last_seg + 1] * 1102 (sqrt(dist) - sqrt(cross[i].distance[current]))) / 1103 sqrt(dist); 1104 } 1094 1105 } 1095 1106 … … 1118 1129 } 1119 1130 } 1120 1121 /* clean up */1122 1131 1123 1132 return 1; 1124 1133 } 1125 1134 1126 static struct line_pnts *APnts, *BPnts, *IPnts;1127 1128 1135 static int cross_found; /* set by find_cross() */ 1129 1136 1130 /* break segments (called by rtree search)*/1131 static int find_cross(int i d, const struct RTree_Rect *rect, int *arg)1137 /* break segments */ 1138 static int find_cross(int i, int j) 1132 1139 { 1133 1140 double x1, y1, z1, x2, y2, z2; 1134 int i, j, ret; 1135 1136 /* !!! segment number for B lines is returned as +1 */ 1137 i = *arg; 1138 j = id - 1; 1139 /* Note: -1 to make up for the +1 when data was inserted */ 1141 double y1min, y1max, y2min, y2max; 1142 int ret; 1143 1144 y1min = APnts->y[i]; 1145 y1max = APnts->y[i + 1]; 1146 if (APnts->y[i] > APnts->y[i + 1]) { 1147 y1min = APnts->y[i + 1]; 1148 y1max = APnts->y[i]; 1149 } 1150 1151 y2min = BPnts->y[j]; 1152 y2max = BPnts->y[j + 1]; 1153 if (BPnts->y[j] > BPnts->y[j + 1]) { 1154 y2min = BPnts->y[j + 1]; 1155 y2max = BPnts->y[j]; 1156 } 1157 1158 if (y1min > y2max || y1max < y2min) 1159 return 0; 1140 1160 1141 1161 ret = Vect_segment_intersection(APnts->x[i], APnts->y[i], APnts->z[i], … … 1187 1207 */ 1188 1208 int 1189 Vect_line_check_intersection(struct line_pnts *APoints, 1190 struct line_pnts *BPoints, int with_z) 1191 { 1192 int i; 1193 double dist, rethresh; 1194 struct RTree *MyRTree; 1195 static struct RTree_Rect rect; 1196 static int rect_init = 0; 1197 1198 if (!rect_init) { 1199 rect.boundary = G_malloc(6 * sizeof(RectReal)); 1200 rect_init = 6; 1201 } 1202 1203 rethresh = 0.000001; /* TODO */ 1209 Vect_line_check_intersection2(struct line_pnts *APoints, 1210 struct line_pnts *BPoints, int with_z) 1211 { 1212 double dist; 1213 struct bound_box ABox, BBox, abbox; 1214 struct boq bo_queue; 1215 struct qitem qi, *found; 1216 struct RB_TREE *bo_ta, *bo_tb; 1217 struct RB_TRAV bo_t_trav; 1218 1204 1219 APnts = APoints; 1205 1220 BPnts = BPoints; … … 1273 1288 /* Take each segment from A and find if intersects any segment from B. */ 1274 1289 1275 /* Spatial index: lines may be very long (thousands of segments) and check each segment 1276 * with each from second line takes a long time (n*m). Because of that, spatial index 1277 * is build first for the second line and segments from the first line are broken by segments 1278 * in bound box */ 1279 1280 /* Create rtree for B line */ 1281 MyRTree = RTreeCreateTree(-1, 0, 2); 1282 RTreeSetOverflow(MyRTree, 0); 1283 for (i = 0; i < BPoints->n_points - 1; i++) { 1284 if (BPoints->x[i] <= BPoints->x[i + 1]) { 1285 rect.boundary[0] = BPoints->x[i]; 1286 rect.boundary[3] = BPoints->x[i + 1]; 1287 } 1288 else { 1289 rect.boundary[0] = BPoints->x[i + 1]; 1290 rect.boundary[3] = BPoints->x[i]; 1291 } 1292 1293 if (BPoints->y[i] <= BPoints->y[i + 1]) { 1294 rect.boundary[1] = BPoints->y[i]; 1295 rect.boundary[4] = BPoints->y[i + 1]; 1296 } 1297 else { 1298 rect.boundary[1] = BPoints->y[i + 1]; 1299 rect.boundary[4] = BPoints->y[i]; 1300 } 1301 1302 if (BPoints->z[i] <= BPoints->z[i + 1]) { 1303 rect.boundary[2] = BPoints->z[i]; 1304 rect.boundary[5] = BPoints->z[i + 1]; 1305 } 1306 else { 1307 rect.boundary[2] = BPoints->z[i + 1]; 1308 rect.boundary[5] = BPoints->z[i]; 1309 } 1310 1311 RTreeInsertRect(&rect, i + 1, MyRTree); /* B line segment numbers in rtree start from 1 */ 1312 } 1313 1314 /* Find intersection */ 1290 dig_line_box(APoints, &ABox); 1291 dig_line_box(BPoints, &BBox); 1292 1293 if (!Vect_box_overlap(&ABox, &BBox)) { 1294 return 0; 1295 } 1296 1297 /* overlap box */ 1298 abbox = BBox; 1299 if (abbox.N > ABox.N) 1300 abbox.N = ABox.N; 1301 if (abbox.S < ABox.S) 1302 abbox.S = ABox.S; 1303 if (abbox.E > ABox.E) 1304 abbox.E = ABox.E; 1305 if (abbox.W < ABox.W) 1306 abbox.W = ABox.W; 1307 if (abbox.T > ABox.T) 1308 abbox.T = ABox.T; 1309 if (abbox.B < ABox.B) 1310 abbox.B = ABox.B; 1311 1312 abbox.N += rethresh; 1313 abbox.S -= rethresh; 1314 abbox.E += rethresh; 1315 abbox.W -= rethresh; 1316 abbox.T += rethresh; 1317 abbox.B -= rethresh; 1318 1319 /* initialize queue */ 1320 bo_queue.count = 0; 1321 bo_queue.alloc = 2 * (APnts->n_points + BPnts->n_points); 1322 bo_queue.i = G_malloc(bo_queue.alloc * sizeof(struct qitem)); 1323 1324 /* load APnts to queue */ 1325 boq_load(&bo_queue, APnts, &abbox, 0, with_z); 1326 1327 /* load BPnts to queue */ 1328 boq_load(&bo_queue, BPnts, &abbox, 1, with_z); 1329 1330 /* initialize search tree */ 1331 bo_ta = rbtree_create(cmp_t_y, sizeof(struct qitem)); 1332 bo_tb = rbtree_create(cmp_t_y, sizeof(struct qitem)); 1333 1334 /* find intersection */ 1315 1335 cross_found = 0; 1316 1317 for (i = 0; i < APoints->n_points - 1; i++) { 1318 if (APoints->x[i] <= APoints->x[i + 1]) { 1319 rect.boundary[0] = APoints->x[i]; 1320 rect.boundary[3] = APoints->x[i + 1]; 1321 } 1322 else { 1323 rect.boundary[0] = APoints->x[i + 1]; 1324 rect.boundary[3] = APoints->x[i]; 1325 } 1326 1327 if (APoints->y[i] <= APoints->y[i + 1]) { 1328 rect.boundary[1] = APoints->y[i]; 1329 rect.boundary[4] = APoints->y[i + 1]; 1330 } 1331 else { 1332 rect.boundary[1] = APoints->y[i + 1]; 1333 rect.boundary[4] = APoints->y[i]; 1334 } 1335 if (APoints->z[i] <= APoints->z[i + 1]) { 1336 rect.boundary[2] = APoints->z[i]; 1337 rect.boundary[5] = APoints->z[i + 1]; 1338 } 1339 else { 1340 rect.boundary[2] = APoints->z[i + 1]; 1341 rect.boundary[5] = APoints->z[i]; 1342 } 1343 1344 RTreeSearch(MyRTree, &rect, (void *)find_cross, &i); /* A segment number from 0 */ 1345 1336 while (boq_drop(&bo_queue, &qi)) { 1337 1338 if (qi.e == QEVT_IN) { 1339 /* not the original Bentley-Ottmann algorithm */ 1340 if (qi.l == 0) { 1341 /* test for intersection of s with all segments in T */ 1342 rbtree_init_trav(&bo_t_trav, bo_tb); 1343 while ((found = rbtree_traverse(&bo_t_trav))) { 1344 find_cross(qi.s, found->s); 1345 1346 if (cross_found) { 1347 break; 1348 } 1349 } 1350 1351 /* insert s in T */ 1352 rbtree_insert(bo_ta, &qi); 1353 } 1354 else { 1355 /* test for intersection of s with all segments in T */ 1356 rbtree_init_trav(&bo_t_trav, bo_ta); 1357 while ((found = rbtree_traverse(&bo_t_trav))) { 1358 find_cross(found->s, qi.s); 1359 1360 if (cross_found) { 1361 break; 1362 } 1363 } 1364 1365 /* insert s in T */ 1366 rbtree_insert(bo_tb, &qi); 1367 } 1368 } 1369 else if (qi.e == QEVT_OUT) { 1370 /* remove from T */ 1371 1372 /* adjust p */ 1373 if (qi.p == qi.s) 1374 qi.p++; 1375 else 1376 qi.p--; 1377 1378 if (qi.l == 0) { 1379 if (!rbtree_remove(bo_ta, &qi)) 1380 G_fatal_error("RB tree error!"); 1381 } 1382 else { 1383 if (!rbtree_remove(bo_tb, &qi)) 1384 G_fatal_error("RB tree error!"); 1385 } 1386 } 1346 1387 if (cross_found) { 1347 1388 break; 1348 1389 } 1349 1390 } 1350 1351 /* Free RTree */1352 RTreeDestroyTree(MyRTree);1391 G_free(bo_queue.i); 1392 rbtree_destroy(bo_ta); 1393 rbtree_destroy(bo_tb); 1353 1394 1354 1395 return cross_found; … … 1369 1410 */ 1370 1411 int 1371 Vect_line_get_intersections (struct line_pnts *APoints,1372 struct line_pnts *BPoints,1373 struct line_pnts *IPoints, int with_z)1412 Vect_line_get_intersections2(struct line_pnts *APoints, 1413 struct line_pnts *BPoints, 1414 struct line_pnts *IPoints, int with_z) 1374 1415 { 1375 1416 int ret; 1376 1417 1377 1418 IPnts = IPoints; 1378 ret = Vect_line_check_intersection (APoints, BPoints, with_z);1419 ret = Vect_line_check_intersection2(APoints, BPoints, with_z); 1379 1420 1380 1421 return ret;
Note:
See TracChangeset
for help on using the changeset viewer.
