23 #include <grass/gis.h>
24 #include <grass/Vect.h>
25 #include <grass/glocale.h>
28 static int cmp_cross(
const void *pa,
const void *pb);
29 static void add_cross(
int asegment,
double adistance,
int bsegment,
30 double bdistance,
double x,
double y);
31 static double dist2(
double x1,
double y1,
double x2,
double y2);
34 static int ident(
double x1,
double y1,
double x2,
double y2,
double thresh);
36 static int cross_seg(
int id,
int *arg);
37 static int find_cross(
int id,
int *arg);
52 #define D ((ax2-ax1)*(by1-by2) - (ay2-ay1)*(bx1-bx2))
53 #define D1 ((bx1-ax1)*(by1-by2) - (by1-ay1)*(bx1-bx2))
54 #define D2 ((ax2-ax1)*(by1-ay1) - (ay2-ay1)*(bx1-ax1))
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)
115 static int first_3d = 1;
116 double d, d1, d2, r1, r2, dtol, t;
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);
126 if (with_z && first_3d) {
127 G_warning(_(
"3D not supported by Vect_segment_intersection()"));
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");
147 else if (bx1 == ax1) {
150 else if (bx2 == ax2) {
153 else if (by1 == ay1) {
178 G_debug(2,
"Vect_segment_intersection(): d = %f, d1 = %f, d2 = %f", d, d1,
184 if (fabs(d) > dtol) {
188 G_debug(2,
" -> not parallel/collinear: r1 = %f, r2 = %f", r1, r2);
190 if (r1 < 0 || r1 > 1 || r2 < 0 || r2 > 1) {
191 G_debug(2,
" -> no intersection");
195 *x1 = ax1 + r1 * (ax2 - ax1);
196 *y1 = ay1 + r1 * (ay2 - ay1);
199 G_debug(2,
" -> intersection %f, %f", *x1, *y1);
204 G_debug(3,
" -> parallel/collinear");
216 if (ax1 == ax2 && bx1 == bx2 && ax1 == bx1) {
217 G_debug(2,
" -> collinear vertical");
228 if (ay1 > by2 || ay2 < by1) {
229 G_debug(2,
" -> no intersection");
238 G_debug(2,
" -> connected by end points");
245 G_debug(2,
" -> connected by end points");
250 G_debug(3,
" -> vertical overlap");
252 if (ay1 <= by1 && ay2 >= by2) {
253 G_debug(2,
" -> a contains b");
266 if (ay1 >= by1 && ay2 <= by2) {
267 G_debug(2,
" -> b contains a");
281 G_debug(2,
" -> partial overlap");
282 if (by1 > ay1 && by1 < ay2) {
301 if (by2 > ay1 && by2 < ay2) {
322 G_warning(_(
"Vect_segment_intersection() ERROR (collinear vertical segments)"));
332 G_debug(2,
" -> collinear non vertical");
335 if ((bx1 > ax1 && bx2 > ax1 && bx1 > ax2 && bx2 > ax2) ||
336 (bx1 < ax1 && bx2 < ax1 && bx1 < ax2 && bx2 < ax2)) {
337 G_debug(2,
" -> no intersection");
342 G_debug(2,
" -> overlap/connected end points");
345 if ((ax1 == bx1 && ay1 == by1) || (ax1 == bx2 && ay1 == by2)) {
349 G_debug(2,
" -> connected by end points");
352 if ((ax2 == bx1 && ay2 == by1) || (ax2 == bx2 && ay2 == by2)) {
356 G_debug(2,
" -> connected by end points");
378 if (ax1 <= bx1 && ax2 >= bx2) {
379 G_debug(2,
" -> a contains b");
392 if (ax1 >= bx1 && ax2 <= bx2) {
393 G_debug(2,
" -> b contains a");
407 G_debug(2,
" -> partial overlap");
408 if (bx1 > ax1 && bx1 < ax2) {
427 if (bx2 > ax1 && bx2 < ax2) {
448 G_warning(_(
"Vect_segment_intersection() ERROR (collinear non vertical segments)"));
469 static int a_cross = 0;
472 static int *use_cross =
NULL;
474 static void add_cross(
int asegment,
double adistance,
int bsegment,
475 double bdistance,
double x,
double y)
477 if (n_cross == a_cross) {
480 (
CROSS *) G_realloc((
void *)cross,
481 (a_cross + 101) *
sizeof(
CROSS));
483 (
int *)G_realloc((
void *)use_cross,
484 (a_cross + 101) *
sizeof(
int));
489 " add new cross: aseg/dist = %d/%f bseg/dist = %d/%f, x = %f y = %f",
490 asegment, adistance, bsegment, bdistance, x, y);
491 cross[n_cross].
segment[0] = asegment;
492 cross[n_cross].
distance[0] = adistance;
493 cross[n_cross].
segment[1] = bsegment;
494 cross[n_cross].
distance[1] = bdistance;
495 cross[n_cross].
x = x;
496 cross[n_cross].
y = y;
500 static int cmp_cross(
const void *pa,
const void *pb)
517 static double dist2(
double x1,
double y1,
double x2,
double y2)
523 return (dx * dx + dy * dy);
528 static int ident(
double x1,
double y1,
double x2,
double y2,
double thresh)
534 if ((dx * dx + dy * dy) <= thresh * thresh)
542 static struct line_pnts *APnts, *BPnts;
545 static int cross_seg(
int id,
int *arg)
547 double x1, y1, z1, x2, y2, z2;
556 APnts->x[i + 1], APnts->y[i + 1],
557 APnts->z[i + 1], BPnts->x[j], BPnts->y[j],
558 BPnts->z[j], BPnts->x[j + 1],
559 BPnts->y[j + 1], BPnts->z[j + 1], &x1,
560 &y1, &z1, &x2, &y2, &z2, 0);
564 G_debug(2,
" -> %d x %d: intersection type = %d", i, j, ret);
566 G_debug(3,
" in %f, %f ", x1, y1);
567 add_cross(i, 0.0, j, 0.0, x1, y1);
569 else if (ret == 2 || ret == 3 || ret == 4 || ret == 5) {
574 G_debug(3,
" in %f, %f; %f, %f", x1, y1, x2, y2);
575 add_cross(i, 0.0, j, 0.0, x1, y1);
576 add_cross(i, 0.0, j, 0.0, x2, y2);
602 struct line_pnts *BPoints,
603 struct line_pnts ***ALines,
604 struct line_pnts ***BLines,
605 int *nalines,
int *nblines,
int with_z)
607 int i, j, k, l, last_seg, seg, last;
609 double dist, curdist, last_dist, last_x, last_y, last_z;
610 double x, y, rethresh;
612 struct line_pnts **XLines, *Points;
614 struct line_pnts *Points1, *Points2;
615 int seg1, seg2, vert1, vert2;
681 for (i = 0; i < BPoints->n_points - 1; i++) {
682 if (BPoints->x[i] <= BPoints->x[i + 1]) {
684 rect.
boundary[3] = BPoints->x[i + 1];
687 rect.
boundary[0] = BPoints->x[i + 1];
691 if (BPoints->y[i] <= BPoints->y[i + 1]) {
693 rect.
boundary[4] = BPoints->y[i + 1];
696 rect.
boundary[1] = BPoints->y[i + 1];
700 if (BPoints->z[i] <= BPoints->z[i + 1]) {
702 rect.
boundary[5] = BPoints->z[i + 1];
705 rect.
boundary[2] = BPoints->z[i + 1];
713 for (i = 0; i < APoints->n_points - 1; i++) {
714 if (APoints->x[i] <= APoints->x[i + 1]) {
716 rect.
boundary[3] = APoints->x[i + 1];
719 rect.
boundary[0] = APoints->x[i + 1];
723 if (APoints->y[i] <= APoints->y[i + 1]) {
725 rect.
boundary[4] = APoints->y[i + 1];
728 rect.
boundary[1] = APoints->y[i + 1];
731 if (APoints->z[i] <= APoints->z[i + 1]) {
733 rect.
boundary[5] = APoints->z[i + 1];
736 rect.
boundary[2] = APoints->z[i + 1];
740 j =
RTreeSearch(RTree, &rect, (
void *)cross_seg, &i);
746 G_debug(2,
"n_cross = %d", n_cross);
755 for (i = 0; i < n_cross; i++) {
759 dist2(cross[i].x, cross[i].y, APoints->
x[seg], APoints->y[seg]);
765 dist2(cross[i].x, cross[i].y, APoints->
x[seg + 1],
766 APoints->y[seg + 1]);
767 if (dist < curdist) {
769 x = APoints->x[seg + 1];
770 y = APoints->y[seg + 1];
776 dist2(cross[i].x, cross[i].y, BPoints->
x[seg], BPoints->y[seg]);
777 if (dist < curdist) {
782 dist = dist2(cross[i].x, cross[i].y, BPoints->
x[seg + 1], BPoints->y[seg + 1]);
783 if (dist < curdist) {
785 x = BPoints->x[seg + 1];
786 y = BPoints->y[seg + 1];
788 if (curdist < rethresh * rethresh) {
795 for (i = 0; i < n_cross; i++) {
798 dist2(APoints->x[seg], APoints->y[seg], cross[i].
x, cross[i].
y);
801 dist2(BPoints->x[seg], BPoints->y[seg], cross[i].
x, cross[i].
y);
805 for (l = 1; l < 3; l++) {
806 for (i = 0; i < n_cross; i++)
810 XLines = G_malloc((n_cross + 1) *
sizeof(
struct line_pnts *));
813 G_debug(2,
"Clean and create array for line A");
821 G_debug(2,
"Clean and create array for line B");
830 qsort((
void *)cross,
sizeof(
char) * n_cross,
sizeof(
CROSS),
834 for (i = 0; i < n_cross; i++) {
836 " cross = %d seg1/dist1 = %d/%f seg2/dist2 = %d/%f x = %f y = %f",
837 i, cross[i].segment[current],
838 sqrt(cross[i].distance[current]),
839 cross[i].segment[second], sqrt(cross[i].distance[second]),
840 cross[i].x, cross[i].y);
844 for (i = 0; i < n_cross; i++) {
845 if (use_cross[i] == 1) {
846 j = Points1->n_points - 1;
849 if ((cross[i].segment[current] == 0 &&
850 cross[i].x == Points1->
x[0] &&
851 cross[i].
y == Points1->y[0]) ||
852 (cross[i].
segment[current] == j - 1 &&
853 cross[i].
x == Points1->x[j] &&
854 cross[i].
y == Points1->y[j])) {
856 G_debug(3,
"cross %d deleted (first/last point)", i);
872 for (i = 0; i < n_cross; i++) {
873 if (use_cross[i] == 0)
875 G_debug(3,
" is %d between colinear?", i);
877 seg1 = cross[i].
segment[current];
878 seg2 = cross[i].
segment[second];
881 if (cross[i].x == Points1->
x[seg1] &&
882 cross[i].
y == Points1->y[seg1]) {
885 else if (cross[i].x == Points1->
x[seg1 + 1] &&
886 cross[i].
y == Points1->y[seg1 + 1]) {
890 G_debug(3,
" -> is not vertex on 1. line");
897 if (cross[i].x == Points2->
x[seg2] &&
898 cross[i].
y == Points2->y[seg2]) {
901 else if (cross[i].x == Points2->
x[seg2 + 1] &&
902 cross[i].
y == Points2->y[seg2 + 1]) {
906 G_debug(3,
" -> is not vertex on 2. line");
909 G_debug(3,
" seg1/vert1 = %d/%d seg2/ver2 = %d/%d", seg1,
913 if (vert2 == 0 || vert2 == Points2->n_points - 1) {
914 G_debug(3,
" -> vertex 2 (%d) is first/last", vert2);
919 if (!((Points1->x[vert1 - 1] == Points2->x[vert2 - 1] &&
920 Points1->y[vert1 - 1] == Points2->y[vert2 - 1] &&
921 Points1->x[vert1 + 1] == Points2->x[vert2 + 1] &&
922 Points1->y[vert1 + 1] == Points2->y[vert2 + 1]) ||
923 (Points1->x[vert1 - 1] == Points2->x[vert2 + 1] &&
924 Points1->y[vert1 - 1] == Points2->y[vert2 + 1] &&
925 Points1->x[vert1 + 1] == Points2->x[vert2 - 1] &&
926 Points1->y[vert1 + 1] == Points2->y[vert2 - 1])
929 G_debug(3,
" -> previous/next are not identical");
935 G_debug(3,
" -> collinear -> remove");
954 for (i = 1; i < n_cross; i++) {
955 if (use_cross[i] == 0)
961 seg = cross[i].
segment[current];
963 G_debug(3,
" duplicate ?: cross = %d seg = %d dist = %f", i,
964 cross[i].segment[current], cross[i].distance[current]);
965 if ((cross[i].segment[current] == cross[last].segment[current] &&
966 cross[i].distance[current] == cross[last].distance[current])
967 || (cross[i].segment[current] ==
968 cross[last].segment[current] + 1 &&
969 cross[i].distance[current] == 0 &&
970 cross[i].x == cross[last].x &&
971 cross[i].y == cross[last].y)) {
972 G_debug(3,
" cross %d identical to last -> removed", i);
984 for (i = 0; i < n_cross; i++) {
985 if (use_cross[i] == 1) {
992 if (n_alive_cross > 0) {
994 use_cross[n_cross] = 1;
995 j = Points->n_points - 1;
996 cross[n_cross].
x = Points->x[j];
997 cross[n_cross].
y = Points->y[j];
998 cross[n_cross].
segment[current] = Points->n_points - 2;
1002 last_x = Points->x[0];
1003 last_y = Points->y[0];
1004 last_z = Points->z[0];
1007 for (i = 0; i <= n_cross; i++) {
1008 seg = cross[i].
segment[current];
1009 G_debug(2,
"%d seg = %d dist = %f", i, seg,
1010 cross[i].distance[current]);
1011 if (use_cross[i] == 0) {
1012 G_debug(3,
" removed -> next");
1020 G_debug(2,
" append last vert: %f %f", last_x, last_y);
1023 for (j = last_seg + 1; j <= seg; j++) {
1024 G_debug(2,
" segment j = %d", j);
1026 if ((j == last_seg + 1) && Points->x[j] == last_x &&
1027 Points->y[j] == last_y) {
1028 G_debug(2,
" -> skip (identical to last break)");
1033 G_debug(2,
" append first of seg: %f %f", Points->x[j],
1039 G_debug(2,
" append cross / last point: %f %f", cross[i].x,
1042 last_x = cross[i].
x;
1043 last_y = cross[i].
y, last_z = 0;
1047 G_debug(2,
" line is degenerate -> skipped");
1068 static struct line_pnts *APnts, *BPnts, *IPnts;
1070 static int cross_found;
1073 static int find_cross(
int id,
int *arg)
1075 double x1, y1, z1, x2, y2, z2;
1084 APnts->x[i + 1], APnts->y[i + 1],
1085 APnts->z[i + 1], BPnts->x[j], BPnts->y[j],
1086 BPnts->z[j], BPnts->x[j + 1],
1087 BPnts->y[j + 1], BPnts->z[j + 1], &x1,
1088 &y1, &z1, &x2, &y2, &z2, 0);
1099 G_warning(_(
"Error while adding point to array. Out of memory"));
1105 G_warning(_(
"Error while adding point to array. Out of memory"));
1107 G_warning(_(
"Error while adding point to array. Out of memory"));
1132 struct line_pnts *BPoints,
int with_z)
1135 double dist, rethresh;
1139 rethresh = 0.000001;
1149 if (APoints->n_points == 1 && BPoints->n_points == 1) {
1150 if (APoints->x[0] == BPoints->x[0] && APoints->y[0] == BPoints->y[0]) {
1154 &APoints->y[0],
NULL, 1))
1155 G_warning(_(
"Error while adding point to array. Out of memory"));
1159 if (APoints->z[0] == BPoints->z[0]) {
1162 &APoints->y[0], &APoints->z[0],
1164 G_warning(_(
"Error while adding point to array. Out of memory"));
1176 if (APoints->n_points == 1) {
1181 if (dist <= rethresh) {
1185 G_warning(_(
"Error while adding point to array. Out of memory"));
1193 if (BPoints->n_points == 1) {
1198 if (dist <= rethresh) {
1202 G_warning(_(
"Error while adding point to array. Out of memory"));
1218 for (i = 0; i < BPoints->n_points - 1; i++) {
1219 if (BPoints->x[i] <= BPoints->x[i + 1]) {
1221 rect.
boundary[3] = BPoints->x[i + 1];
1224 rect.
boundary[0] = BPoints->x[i + 1];
1228 if (BPoints->y[i] <= BPoints->y[i + 1]) {
1230 rect.
boundary[4] = BPoints->y[i + 1];
1233 rect.
boundary[1] = BPoints->y[i + 1];
1237 if (BPoints->z[i] <= BPoints->z[i + 1]) {
1239 rect.
boundary[5] = BPoints->z[i + 1];
1242 rect.
boundary[2] = BPoints->z[i + 1];
1252 for (i = 0; i < APoints->n_points - 1; i++) {
1253 if (APoints->x[i] <= APoints->x[i + 1]) {
1255 rect.
boundary[3] = APoints->x[i + 1];
1258 rect.
boundary[0] = APoints->x[i + 1];
1262 if (APoints->y[i] <= APoints->y[i + 1]) {
1264 rect.
boundary[4] = APoints->y[i + 1];
1267 rect.
boundary[1] = APoints->y[i + 1];
1270 if (APoints->z[i] <= APoints->z[i + 1]) {
1272 rect.
boundary[5] = APoints->z[i + 1];
1275 rect.
boundary[2] = APoints->z[i + 1];
1279 RTreeSearch(RTree, &rect, (
void *)find_cross, &i);
1307 struct line_pnts *BPoints,
1308 struct line_pnts *IPoints,
int with_z)