00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018 #include <stdlib.h>
00019 #include <math.h>
00020 #include "gis.h"
00021 #include "Vect.h"
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035 #define D ((ax2-ax1)*(by1-by2) - (ay2-ay1)*(bx1-bx2))
00036 #define D1 ((bx1-ax1)*(by1-by2) - (by1-ay1)*(bx1-bx2))
00037 #define D2 ((ax2-ax1)*(by1-ay1) - (ay2-ay1)*(bx1-ax1))
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00094 int Vect_segment_intersection (
00095 double ax1, double ay1, double az1, double ax2, double ay2, double az2,
00096 double bx1, double by1, double bz1, double bx2, double by2, double bz2,
00097 double *x1, double *y1, double *z1,
00098 double *x2, double *y2, double *z2,
00099 int with_z)
00100 {
00101 static int first_3d = 1;
00102 double d, d1, d2, r1, r2, dtol, t;
00103 int switched = 0;
00104
00105
00106
00107 G_debug (4, "Vect_segment_intersection()" );
00108 G_debug (4, " %.15g , %.15g - %.15g , %.15g", ax1, ay1, ax2, ay2 );
00109 G_debug (4, " %.15g , %.15g - %.15g , %.15g", bx1, by1, bx2, by2 );
00110
00111
00112 if ( with_z && first_3d ) {
00113 G_warning ( "3D not supported by Vect_segment_intersection()" );
00114 first_3d = 0;
00115 }
00116
00117
00118 if ( ( ax1 == bx1 && ay1 == by1 && ax2 == bx2 && ay2 == by2 ) ||
00119 ( ax1 == bx2 && ay1 == by2 && ax2 == bx1 && ay2 == by1 ) ) {
00120 G_debug (2, " -> identical segments" );
00121 *x1 = ax1; *y1 = ay1; *z1 = az1;
00122 *x2 = ax2; *y2 = ay2; *z2 = az2;
00123 return 5;
00124 }
00125
00126
00127 if ( bx1 < ax1 ) switched = 1;
00128 else if ( bx1 == ax1 ) {
00129 if ( bx2 < ax2 ) switched = 1;
00130 else if ( bx2 == ax2 ) {
00131 if ( by1 < ay1 ) switched = 1;
00132 else if ( by1 == ay1 ) {
00133 if ( by2 < ay2 ) switched = 1;
00134 }
00135 }
00136 }
00137 if ( switched ) {
00138 t = ax1; ax1 = bx1; bx1 = t; t = ay1; ay1 = by1; by1 = t;
00139 t = ax2; ax2 = bx2; bx2 = t; t = ay2; ay2 = by2; by2 = t;
00140 }
00141
00142 d = D;
00143 d1 = D1;
00144 d2 = D2;
00145
00146 G_debug (2, "Vect_segment_intersection(): d = %f, d1 = %f, d2 = %f", d, d1, d2 );
00147
00148
00149
00150 dtol = 0.0;
00151 if (fabs (d) > dtol) {
00152 r1 = D1/d;
00153 r2 = D2/d;
00154
00155 G_debug (2, " -> not parallel/collinear: r1 = %f, r2 = %f", r1, r2 );
00156
00157 if (r1 < 0 || r1 > 1 || r2 < 0 || r2 > 1) {
00158 G_debug (2, " -> no intersection" );
00159 return 0;
00160 }
00161
00162 *x1 = ax1 + r1 * (ax2 - ax1);
00163 *y1 = ay1 + r1 * (ay2 - ay1);
00164 *z1 = 0;
00165
00166 G_debug (2, " -> intersection %f, %f", *x1, *y1 );
00167 return 1;
00168 }
00169
00170
00171 G_debug (3, " -> parallel/collinear" );
00172
00173 if (D1 || D2) {
00174 G_debug (2, " -> parallel" );
00175 return 0;
00176 }
00177
00178
00179
00180
00181
00182
00183 if (ax1 == ax2 && bx1==bx2 && ax1==bx1) {
00184 G_debug (2, " -> collinear vertical" );
00185 if (ay1 > ay2) { t=ay1; ay1=ay2; ay2=t; }
00186 if (by1 > by2) { t=by1; by1=by2; by2=t; }
00187 if (ay1 > by2 || ay2 < by1) {
00188 G_debug (2, " -> no intersection" );
00189 return 0;
00190 }
00191
00192
00193 if (ay1 == by2) {
00194 *x1 = ax1; *y1 = ay1; *z1 = 0;
00195 G_debug (2, " -> connected by end points");
00196 return 1;
00197 }
00198 if(ay2 == by1) {
00199 *x1 = ax2; *y1 = ay2; *z1 = 0;
00200 G_debug (2, " -> connected by end points");
00201 return 1;
00202 }
00203
00204
00205 G_debug (3, " -> vertical overlap" );
00206
00207 if ( ay1 <= by1 && ay2 >= by2 ) {
00208 G_debug (2, " -> a contains b" );
00209 *x1 = bx1; *y1 = by1; *z1 = 0;
00210 *x2 = bx2; *y2 = by2; *z2 = 0;
00211 if ( !switched )
00212 return 3;
00213 else
00214 return 4;
00215 }
00216
00217 if ( ay1 >= by1 && ay2 <= by2 ) {
00218 G_debug (2, " -> b contains a" );
00219 *x1 = ax1; *y1 = ay1; *z1 = 0;
00220 *x2 = ax2; *y2 = ay2; *z2 = 0;
00221 if ( !switched )
00222 return 4;
00223 else
00224 return 3;
00225 }
00226
00227
00228 G_debug (2, " -> partial overlap" );
00229 if ( by1 > ay1 && by1 < ay2 ) {
00230 if ( !switched ) {
00231 *x1 = bx1; *y1 = by1; *z1 = 0;
00232 *x2 = ax2; *y2 = ay2; *z2 = 0;
00233 } else {
00234 *x1 = ax2; *y1 = ay2; *z1 = 0;
00235 *x2 = bx1; *y2 = by1; *z2 = 0;
00236 }
00237 return 2;
00238 }
00239 if ( by2 > ay1 && by2 < ay2 ) {
00240 if ( !switched ) {
00241 *x1 = bx2; *y1 = by2; *z1 = 0;
00242 *x2 = ax1; *y2 = ay1; *z2 = 0;
00243 } else {
00244 *x1 = ax1; *y1 = ay1; *z1 = 0;
00245 *x2 = bx2; *y2 = by2; *z2 = 0;
00246 }
00247 return 2;
00248 }
00249
00250
00251 G_warning("Vect_segment_intersection() ERROR (collinear vertical segments)");
00252 G_warning("%.15g %.15g", ax1, ay1);
00253 G_warning("%.15g %.15g", ax2, ay2);
00254 G_warning("x");
00255 G_warning("%.15g %.15g", bx1, by1);
00256 G_warning("%.15g %.15g", bx2, by2);
00257
00258 return 0;
00259 }
00260
00261 G_debug (2, " -> collinear non vertical" );
00262
00263
00264 if ( ( bx1 > ax1 && bx2 > ax1 && bx1 > ax2 && bx2 > ax2 ) ||
00265 ( bx1 < ax1 && bx2 < ax1 && bx1 < ax2 && bx2 < ax2 ) ) {
00266 G_debug (2, " -> no intersection" );
00267 return 0;
00268 }
00269
00270
00271 G_debug (2, " -> overlap/connected end points" );
00272
00273
00274 if ( (ax1 == bx1 && ay1 == by1) || (ax1 == bx2 && ay1 == by2) ) {
00275 *x1 = ax1; *y1 = ay1; *z1 = 0;
00276 G_debug (2, " -> connected by end points");
00277 return 1;
00278 }
00279 if ( (ax2 == bx1 && ay2 == by1) || (ax2 == bx2 && ay2 == by2) ) {
00280 *x1 = ax2; *y1 = ay2; *z1 = 0;
00281 G_debug (2, " -> connected by end points");
00282 return 1;
00283 }
00284
00285 if (ax1 > ax2) { t=ax1; ax1=ax2; ax2=t; t=ay1; ay1=ay2; ay2=t; }
00286 if (bx1 > bx2) { t=bx1; bx1=bx2; bx2=t; t=by1; by1=by2; by2=t; }
00287
00288
00289 if ( ax1 <= bx1 && ax2 >= bx2 ) {
00290 G_debug (2, " -> a contains b" );
00291 *x1 = bx1; *y1 = by1; *z1 = 0;
00292 *x2 = bx2; *y2 = by2; *z2 = 0;
00293 if ( !switched )
00294 return 3;
00295 else
00296 return 4;
00297 }
00298
00299 if ( ax1 >= bx1 && ax2 <= bx2 ) {
00300 G_debug (2, " -> b contains a" );
00301 *x1 = ax1; *y1 = ay1; *z1 = 0;
00302 *x2 = ax2; *y2 = ay2; *z2 = 0;
00303 if ( !switched )
00304 return 4;
00305 else
00306 return 3;
00307 }
00308
00309
00310 G_debug (2, " -> partial overlap" );
00311 if ( bx1 > ax1 && bx1 < ax2 ) {
00312 if ( !switched ) {
00313 *x1 = bx1; *y1 = by1; *z1 = 0;
00314 *x2 = ax2; *y2 = ay2; *z2 = 0;
00315 } else {
00316 *x1 = ax2; *y1 = ay2; *z1 = 0;
00317 *x2 = bx1; *y2 = by1; *z2 = 0;
00318 }
00319 return 2;
00320 }
00321 if ( bx2 > ax1 && bx2 < ax2 ) {
00322 if ( !switched ) {
00323 *x1 = bx2; *y1 = by2; *z1 = 0;
00324 *x2 = ax1; *y2 = ay1; *z2 = 0;
00325 } else {
00326 *x1 = ax1; *y1 = ay1; *z1 = 0;
00327 *x2 = bx2; *y2 = by2; *z2 = 0;
00328 }
00329 return 2;
00330 }
00331
00332
00333 G_warning("Vect_segment_intersection() ERROR (collinear non vertical segments)");
00334 G_warning("%.15g %.15g", ax1, ay1);
00335 G_warning("%.15g %.15g", ax2, ay2);
00336 G_warning("x");
00337 G_warning("%.15g %.15g", bx1, by1);
00338 G_warning("%.15g %.15g", bx2, by2);
00339
00340 return 0;
00341 }
00342
00343 typedef struct {
00344 int segment[2];
00345 double distance[2];
00346 double x,y,z;
00347 } CROSS;
00348
00349
00350 static int current;
00351 static int second;
00352
00353 static int a_cross = 0;
00354 static int n_cross;
00355 static CROSS *cross = NULL;
00356 static int *use_cross = NULL;
00357
00358 static void add_cross ( int asegment, double adistance, int bsegment, double bdistance, double x, double y)
00359 {
00360 if ( n_cross == a_cross ) {
00361
00362 cross = (CROSS *) G_realloc ( (void *) cross, (a_cross + 101) * sizeof(CROSS) );
00363 use_cross = (int *) G_realloc ( (void *) use_cross, (a_cross + 101) * sizeof(int) );
00364 a_cross += 100;
00365 }
00366
00367 G_debug(5, " add new cross: aseg/dist = %d/%f bseg/dist = %d/%f, x = %f y = %f",
00368 asegment, adistance, bsegment, bdistance, x, y );
00369 cross[n_cross].segment[0] = asegment;
00370 cross[n_cross].distance[0] = adistance;
00371 cross[n_cross].segment[1] = bsegment;
00372 cross[n_cross].distance[1] = bdistance;
00373 cross[n_cross].x = x;
00374 cross[n_cross].y = y;
00375 n_cross++;
00376 }
00377
00378 static int cmp_cross ( const void *pa, const void *pb)
00379 {
00380 CROSS *p1 = (CROSS *) pa;
00381 CROSS *p2 = (CROSS *) pb;
00382
00383 if( p1->segment[current] < p2->segment[current] ) return -1;
00384 if( p1->segment[current] > p2->segment[current] ) return 1;
00385
00386 if( p1->distance[current] < p2->distance[current] ) return -1;
00387 if( p1->distance[current] > p2->distance[current] ) return 1;
00388 return 0;
00389 }
00390
00391 static double dist2 (double x1, double y1, double x2, double y2 ) {
00392 double dx, dy;
00393 dx = x2 - x1; dy = y2 - y1;
00394 return ( dx * dx + dy * dy );
00395 }
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409
00410 static struct line_pnts *APnts, *BPnts;
00411
00412
00413 static int cross_seg(int id, int *arg)
00414 {
00415 double x1, y1, z1, x2, y2, z2;
00416 int i, j, ret;
00417
00418
00419 i = *arg;
00420 j = id - 1;
00421
00422
00423 ret = Vect_segment_intersection (
00424 APnts->x[i], APnts->y[i], APnts->z[i],
00425 APnts->x[i+1], APnts->y[i+1], APnts->z[i+1],
00426 BPnts->x[j], BPnts->y[j], BPnts->z[j],
00427 BPnts->x[j+1], BPnts->y[j+1], BPnts->z[j+1],
00428 &x1, &y1, &z1, &x2, &y2, &z2,
00429 0);
00430
00431
00432 if ( ret > 0 ) {
00433 G_debug(2, " -> %d x %d: intersection type = %d", i, j, ret );
00434 if ( ret == 1 ) {
00435 G_debug(3, " in %f, %f ", x1, y1 );
00436 add_cross ( i, 0.0, j, 0.0, x1, y1 );
00437 } else if ( ret == 2 || ret == 3 || ret == 4 || ret == 5 ) {
00438
00439
00440
00441
00442 G_debug(3, " in %f, %f; %f, %f", x1, y1, x2, y2 );
00443 add_cross ( i, 0.0, j, 0.0, x1, y1 );
00444 add_cross ( i, 0.0, j, 0.0, x2, y2 );
00445 }
00446 }
00447 return 1;
00448 }
00449
00462 int
00463 Vect_line_intersection (
00464 struct line_pnts *APoints,
00465 struct line_pnts *BPoints,
00466 struct line_pnts ***ALines,
00467 struct line_pnts ***BLines,
00468 int *nalines,
00469 int *nblines,
00470 int with_z)
00471 {
00472 int i, j, k, l, last_seg, seg, last;
00473 int n_alive_cross;
00474 double dist, curdist, last_dist, last_x, last_y, last_z;
00475 double x, y, rethresh;
00476 struct Rect rect;
00477 struct line_pnts **XLines, *Points;
00478 struct Node *RTree;
00479 struct line_pnts *Points1, *Points2;
00480 int seg1, seg2, vert1, vert2;
00481
00482 n_cross = 0;
00483 rethresh = 0.000001;
00484 APnts = APoints;
00485 BPnts = BPoints;
00486
00487
00488
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498
00499
00500
00501
00502
00503
00504
00505
00506
00507
00508
00509
00510
00511
00512
00513
00514
00515
00516
00517
00518
00519
00520
00521
00522
00523
00524
00525
00526
00527
00528
00529
00530
00531
00532
00533
00534
00535
00536
00537
00538
00539
00540
00541
00542
00543
00544
00545 RTree = RTreeNewIndex();
00546 for (i = 0; i < BPoints->n_points - 1; i++) {
00547 if ( BPoints->x[i] <= BPoints->x[i+1] ) {
00548 rect.boundary[0] = BPoints->x[i]; rect.boundary[3] = BPoints->x[i+1];
00549 } else {
00550 rect.boundary[0] = BPoints->x[i+1]; rect.boundary[3] = BPoints->x[i];
00551 }
00552
00553 if ( BPoints->y[i] <= BPoints->y[i+1] ) {
00554 rect.boundary[1] = BPoints->y[i]; rect.boundary[4] = BPoints->y[i+1];
00555 } else {
00556 rect.boundary[1] = BPoints->y[i+1]; rect.boundary[4] = BPoints->y[i];
00557 }
00558
00559 if ( BPoints->z[i] <= BPoints->z[i+1] ) {
00560 rect.boundary[2] = BPoints->z[i]; rect.boundary[5] = BPoints->z[i+1];
00561 } else {
00562 rect.boundary[2] = BPoints->z[i+1]; rect.boundary[5] = BPoints->z[i];
00563 }
00564
00565 RTreeInsertRect( &rect, i+1, &RTree, 0);
00566 }
00567
00568
00569 for (i = 0; i < APoints->n_points - 1; i++) {
00570 if ( APoints->x[i] <= APoints->x[i+1] ) {
00571 rect.boundary[0] = APoints->x[i]; rect.boundary[3] = APoints->x[i+1];
00572 } else {
00573 rect.boundary[0] = APoints->x[i+1]; rect.boundary[3] = APoints->x[i];
00574 }
00575
00576 if ( APoints->y[i] <= APoints->y[i+1] ) {
00577 rect.boundary[1] = APoints->y[i]; rect.boundary[4] = APoints->y[i+1];
00578 } else {
00579 rect.boundary[1] = APoints->y[i+1]; rect.boundary[4] = APoints->y[i];
00580 }
00581 if ( APoints->z[i] <= APoints->z[i+1] ) {
00582 rect.boundary[2] = APoints->z[i]; rect.boundary[5] = APoints->z[i+1];
00583 } else {
00584 rect.boundary[2] = APoints->z[i+1]; rect.boundary[5] = APoints->z[i];
00585 }
00586
00587 j = RTreeSearch(RTree, &rect, (void *)cross_seg, &i);
00588 }
00589
00590
00591 RTreeDestroyNode ( RTree );
00592
00593 G_debug ( 2, "n_cross = %d", n_cross );
00594
00595 if ( n_cross == 0 ) {
00596 *nalines = 0;
00597 *nblines = 0;
00598 return 0;
00599 }
00600
00601
00602 for ( i = 0; i < n_cross; i++ ) {
00603
00604 seg = cross[i].segment[0];
00605 curdist = dist2 ( cross[i].x, cross[i].y, APoints->x[seg], APoints->y[seg] );
00606 x = APoints->x[seg]; y = APoints->y[seg];
00607
00608
00609 dist = dist2 ( cross[i].x, cross[i].y, APoints->x[seg+1], APoints->y[seg+1] );
00610 if ( dist < curdist ) {
00611 curdist = dist;
00612 x = APoints->x[seg+1]; y = APoints->y[seg+1];
00613 }
00614
00615
00616 seg = cross[i].segment[1];
00617 dist = dist2 ( cross[i].x, cross[i].y, BPoints->x[seg], BPoints->y[seg] );
00618 if ( dist < curdist ) {
00619 curdist = dist;
00620 x = BPoints->x[seg]; y = BPoints->y[seg];
00621 }
00622 dist = dist2 ( cross[i].x, cross[i].y, BPoints->x[seg+1], BPoints->y[seg+1] );
00623 if ( dist < curdist ) {
00624 curdist = dist;
00625 x = BPoints->x[seg+1]; y = BPoints->y[seg+1];
00626 }
00627 if ( curdist < rethresh * rethresh ) {
00628 cross[i].x = x; cross[i].y = y;
00629 }
00630 }
00631
00632
00633 for ( i = 0; i < n_cross; i++ ) {
00634 seg = cross[i].segment[0];
00635 cross[i].distance[0] = dist2 ( APoints->x[seg], APoints->y[seg], cross[i].x, cross[i].y );
00636 seg = cross[i].segment[1];
00637 cross[i].distance[1] = dist2 ( BPoints->x[seg], BPoints->y[seg], cross[i].x, cross[i].y );
00638 }
00639
00640
00641 for ( l = 1; l < 3; l++ ) {
00642 for ( i = 0; i < n_cross; i++ ) use_cross[i] = 1;
00643
00644
00645 XLines = G_malloc ( (n_cross + 1 ) * sizeof(struct line_pnts *) );
00646
00647 if ( l == 1 ) {
00648 G_debug ( 2, "Clean and create array for line A" );
00649 Points = APoints;
00650 Points1 = APoints;
00651 Points2 = BPoints;
00652 current = 0;
00653 second = 1;
00654 } else {
00655 G_debug ( 2, "Clean and create array for line B" );
00656 Points = BPoints;
00657 Points1 = BPoints;
00658 Points2 = APoints;
00659 current = 1;
00660 second = 0;
00661 }
00662
00663
00664 qsort( (void *)cross, n_cross, sizeof(CROSS), cmp_cross );
00665
00666
00667 for ( i = 0; i < n_cross; i++ ) {
00668 G_debug ( 3, " cross = %d seg1/dist1 = %d/%f seg2/dist2 = %d/%f x = %f y = %f",
00669 i, cross[i].segment[current], sqrt(cross[i].distance[current]),
00670 cross[i].segment[second], sqrt(cross[i].distance[second]),
00671 cross[i].x, cross[i].y );
00672 }
00673
00674
00675 for ( i = 0; i < n_cross; i++ ) {
00676 if ( use_cross[i] == 1 ) {
00677 j = Points1->n_points - 1;
00678
00679
00680 if ((cross[i].segment[current] == 0 && cross[i].x == Points1->x[0] && cross[i].y == Points1->y[0]) ||
00681 (cross[i].segment[current] == j-1 && cross[i].x == Points1->x[j] && cross[i].y == Points1->y[j]))
00682 {
00683 use_cross[i] = 0;
00684 G_debug ( 3, "cross %d deleted (first/last point)", i );
00685 }
00686 }
00687 }
00688
00689
00690
00691
00692
00693
00694
00695
00696
00697
00698
00699
00700 for ( i = 0; i < n_cross; i++ ) {
00701 if ( use_cross[i] == 0 ) continue;
00702 G_debug ( 3, " is %d between colinear?", i);
00703
00704 seg1 = cross[i].segment[current];
00705 seg2 = cross[i].segment[second];
00706
00707
00708 if ( cross[i].x == Points1->x[seg1] && cross[i].y == Points1->y[seg1] ) {
00709 vert1 = seg1;
00710 } else if ( cross[i].x == Points1->x[seg1+1] && cross[i].y == Points1->y[seg1+1] ) {
00711 vert1 = seg1 + 1;
00712 } else {
00713 G_debug ( 3, " -> is not vertex on 1. line");
00714 continue;
00715 }
00716
00717
00718
00719
00720 if ( cross[i].x == Points2->x[seg2] && cross[i].y == Points2->y[seg2] ) {
00721 vert2 = seg2;
00722 } else if ( cross[i].x == Points2->x[seg2+1] && cross[i].y == Points2->y[seg2+1] ) {
00723 vert2 = seg2 + 1;
00724 } else {
00725 G_debug ( 3, " -> is not vertex on 2. line");
00726 continue;
00727 }
00728 G_debug ( 3, " seg1/vert1 = %d/%d seg2/ver2 = %d/%d", seg1, vert1, seg2, vert2);
00729
00730
00731 if ( vert2 == 0 || vert2 == Points2->n_points - 1 ) {
00732 G_debug ( 3, " -> vertex 2 (%d) is first/last", vert2);
00733 continue;
00734 }
00735
00736
00737 if ( !( ( Points1->x[vert1-1] == Points2->x[vert2-1] &&
00738 Points1->y[vert1-1] == Points2->y[vert2-1] &&
00739 Points1->x[vert1+1] == Points2->x[vert2+1] &&
00740 Points1->y[vert1+1] == Points2->y[vert2+1]) ||
00741 ( Points1->x[vert1-1] == Points2->x[vert2+1] &&
00742 Points1->y[vert1-1] == Points2->y[vert2+1] &&
00743 Points1->x[vert1+1] == Points2->x[vert2-1] &&
00744 Points1->y[vert1+1] == Points2->y[vert2-1])
00745 )
00746 ) {
00747 G_debug ( 3, " -> previous/next are not identical");
00748 continue;
00749 }
00750
00751 use_cross[i] = 0;
00752
00753 G_debug (3, " -> collinear -> remove");
00754 }
00755
00756
00757
00758
00759
00760
00761
00762
00763
00764
00765
00766
00767
00768
00769
00770
00771 last = -1;
00772 for ( i = 1; i < n_cross; i++ ) {
00773 if ( use_cross[i] == 0 ) continue;
00774 if ( last == -1 ) {
00775 last = i;
00776 continue;
00777 }
00778 seg = cross[i].segment[current];
00779
00780 G_debug (3, " duplicate ?: cross = %d seg = %d dist = %f", i, cross[i].segment[current] ,
00781 cross[i].distance[current]);
00782 if((cross[i].segment[current] == cross[last].segment[current] && cross[i].distance[current] == cross[last].distance[current]) ||
00783 (cross[i].segment[current] == cross[last].segment[current] + 1 && cross[i].distance[current] == 0 &&
00784 cross[i].x == cross[last].x && cross[i].y == cross[last].y ) )
00785 {
00786 G_debug (3, " cross %d identical to last -> removed", i );
00787 use_cross[i] = 0;
00788 } else {
00789 last = i;
00790 }
00791 }
00792
00793
00794
00795 n_alive_cross = 0;
00796 G_debug (3, " alive crosses:");
00797 for ( i = 0; i < n_cross; i++ ) {
00798 if ( use_cross[i] == 1 ) {
00799 G_debug (3, " %d", i);
00800 n_alive_cross++;
00801 }
00802 }
00803
00804 k = 0;
00805 if ( n_alive_cross > 0 ) {
00806
00807 use_cross[n_cross] = 1;
00808 j = Points->n_points - 1;
00809 cross[n_cross].x = Points->x[j];
00810 cross[n_cross].y = Points->y[j];
00811 cross[n_cross].segment[current] = Points->n_points - 2;
00812
00813 last_seg = 0;
00814 last_dist = 0;
00815 last_x = Points->x[0]; last_y = Points->y[0]; last_z = Points->z[0];
00816
00817
00818 for ( i = 0; i <= n_cross; i++ ) {
00819 seg = cross[i].segment[current];
00820 G_debug ( 2, "%d seg = %d dist = %f", i, seg, cross[i].distance[current] );
00821 if ( use_cross[i] == 0 ) {
00822 G_debug ( 3, " removed -> next" );
00823 continue;
00824 }
00825
00826 G_debug ( 2, " New line:" );
00827 XLines[k] = Vect_new_line_struct ();
00828
00829 Vect_append_point ( XLines[k], last_x, last_y, last_z);
00830 G_debug ( 2, " append last vert: %f %f", last_x, last_y );
00831
00832
00833 for ( j = last_seg + 1; j <= seg; j++ ) {
00834 G_debug ( 2, " segment j = %d", j );
00835
00836 if ( (j == last_seg + 1) && Points->x[j] == last_x && Points->y[j] == last_y ) {
00837 G_debug ( 2, " -> skip (identical to last break)");
00838 continue;
00839 }
00840 Vect_append_point ( XLines[k], Points->x[j], Points->y[j], Points->z[j]);
00841 G_debug ( 2, " append first of seg: %f %f", Points->x[j], Points->y[j] );
00842 }
00843
00844
00845 Vect_append_point ( XLines[k], cross[i].x, cross[i].y, 0.0 );
00846 G_debug ( 2, " append cross / last point: %f %f", cross[i].x, cross[i].y );
00847 last_seg = seg; last_x = cross[i].x; last_y = cross[i].y, last_z = 0;
00848
00849
00850 if ( dig_line_degenerate ( XLines[k] ) > 0 ) {
00851 G_debug ( 2, " line is degenerate -> skipped" );
00852 Vect_destroy_line_struct ( XLines[k] );
00853 } else {
00854 k++;
00855 }
00856 }
00857 }
00858 if ( l == 1 ) {
00859 *nalines = k;
00860 *ALines = XLines;
00861 } else {
00862 *nblines = k;
00863 *BLines = XLines;
00864 }
00865 }
00866
00867 return 1;
00868 }
00869
00870 static struct line_pnts *APnts, *BPnts;
00871
00872 static int cross_found;
00873
00874
00875 static int find_cross(int id, int *arg)
00876 {
00877 double x1, y1, z1, x2, y2, z2;
00878 int i, j, ret;
00879
00880
00881 i = *arg;
00882 j = id - 1;
00883
00884
00885 ret = Vect_segment_intersection (
00886 APnts->x[i], APnts->y[i], APnts->z[i],
00887 APnts->x[i+1], APnts->y[i+1], APnts->z[i+1],
00888 BPnts->x[j], BPnts->y[j], BPnts->z[j],
00889 BPnts->x[j+1], BPnts->y[j+1], BPnts->z[j+1],
00890 &x1, &y1, &z1, &x2, &y2, &z2,
00891 0);
00892
00893
00894 if ( ret > 0 ) {
00895 cross_found = 1;
00896 return 0;
00897 }
00898 return 1;
00899 }
00900
00910 int
00911 Vect_line_check_intersection (
00912 struct line_pnts *APoints,
00913 struct line_pnts *BPoints,
00914 int with_z)
00915 {
00916 int i;
00917 double dist, rethresh;
00918 struct Rect rect;
00919 struct Node *RTree;
00920
00921 rethresh = 0.000001;
00922 APnts = APoints;
00923 BPnts = BPoints;
00924
00925
00926
00927
00928 if ( APoints->n_points == 1 && BPoints->n_points == 1 ) {
00929 if ( APoints->x[0] == BPoints->x[0] && APoints->y[0] == BPoints->y[0] ) {
00930 if ( !with_z ) {
00931 return 1;
00932 } else {
00933 if ( APoints->z[0] == BPoints->z[0] )
00934 return 1;
00935 else
00936 return 0;
00937 }
00938 } else {
00939 return 0;
00940 }
00941 }
00942
00943 if ( APoints->n_points == 1 ) {
00944 Vect_line_distance ( BPoints, APoints->x[0], APoints->y[0], APoints->z[0], with_z,
00945 NULL, NULL, NULL, &dist, NULL, NULL );
00946
00947 if ( dist <= rethresh ) {
00948 return 1;
00949 } else {
00950 return 0;
00951 }
00952 }
00953
00954 if ( BPoints->n_points == 1 ) {
00955 Vect_line_distance ( APoints, BPoints->x[0], BPoints->y[0], BPoints->z[0], with_z,
00956 NULL, NULL, NULL, &dist, NULL, NULL );
00957
00958 if ( dist <= rethresh )
00959 return 1;
00960 else
00961 return 0;
00962 }
00963
00964
00965
00966
00967
00968
00969
00970
00971
00972 RTree = RTreeNewIndex();
00973 for (i = 0; i < BPoints->n_points - 1; i++) {
00974 if ( BPoints->x[i] <= BPoints->x[i+1] ) {
00975 rect.boundary[0] = BPoints->x[i]; rect.boundary[3] = BPoints->x[i+1];
00976 } else {
00977 rect.boundary[0] = BPoints->x[i+1]; rect.boundary[3] = BPoints->x[i];
00978 }
00979
00980 if ( BPoints->y[i] <= BPoints->y[i+1] ) {
00981 rect.boundary[1] = BPoints->y[i]; rect.boundary[4] = BPoints->y[i+1];
00982 } else {
00983 rect.boundary[1] = BPoints->y[i+1]; rect.boundary[4] = BPoints->y[i];
00984 }
00985
00986 if ( BPoints->z[i] <= BPoints->z[i+1] ) {
00987 rect.boundary[2] = BPoints->z[i]; rect.boundary[5] = BPoints->z[i+1];
00988 } else {
00989 rect.boundary[2] = BPoints->z[i+1]; rect.boundary[5] = BPoints->z[i];
00990 }
00991
00992 RTreeInsertRect( &rect, i+1, &RTree, 0);
00993 }
00994
00995
00996 cross_found = 0;
00997
00998 for (i = 0; i < APoints->n_points - 1; i++) {
00999 if ( APoints->x[i] <= APoints->x[i+1] ) {
01000 rect.boundary[0] = APoints->x[i]; rect.boundary[3] = APoints->x[i+1];
01001 } else {
01002 rect.boundary[0] = APoints->x[i+1]; rect.boundary[3] = APoints->x[i];
01003 }
01004
01005 if ( APoints->y[i] <= APoints->y[i+1] ) {
01006 rect.boundary[1] = APoints->y[i]; rect.boundary[4] = APoints->y[i+1];
01007 } else {
01008 rect.boundary[1] = APoints->y[i+1]; rect.boundary[4] = APoints->y[i];
01009 }
01010 if ( APoints->z[i] <= APoints->z[i+1] ) {
01011 rect.boundary[2] = APoints->z[i]; rect.boundary[5] = APoints->z[i+1];
01012 } else {
01013 rect.boundary[2] = APoints->z[i+1]; rect.boundary[5] = APoints->z[i];
01014 }
01015
01016 RTreeSearch(RTree, &rect, (void *)find_cross, &i);
01017
01018 if ( cross_found ) {
01019 break;
01020 }
01021 }
01022
01023
01024 RTreeDestroyNode ( RTree );
01025
01026 return cross_found;
01027 }
01028