[GRASS-SVN] r32656 - grass-addons/vector/v.parallel2
svn_grass at osgeo.org
svn_grass at osgeo.org
Sat Aug 9 17:15:04 EDT 2008
Author: rmatev
Date: 2008-08-09 17:15:04 -0400 (Sat, 09 Aug 2008)
New Revision: 32656
Added:
grass-addons/vector/v.parallel2/e_intersect.c
grass-addons/vector/v.parallel2/e_intersect.h
Modified:
grass-addons/vector/v.parallel2/dgraph.c
grass-addons/vector/v.parallel2/main.c
grass-addons/vector/v.parallel2/vlib_buffer.c
Log:
Everything working. Needs GMP to compile
Modified: grass-addons/vector/v.parallel2/dgraph.c
===================================================================
--- grass-addons/vector/v.parallel2/dgraph.c 2008-08-09 13:21:36 UTC (rev 32655)
+++ grass-addons/vector/v.parallel2/dgraph.c 2008-08-09 21:15:04 UTC (rev 32656)
@@ -2,11 +2,11 @@
#include <math.h>
#include <grass/Vect.h>
#include <grass/gis.h>
+#include <gmp.h>
#include "dgraph.h"
+#include "e_intersect.h"
#define LENGTH(DX, DY) (sqrt((DX*DX)+(DY*DY)))
-#define FZERO(X, TOL) (fabs(X)<TOL)
-#define FEQUAL(X, Y, TOL) (fabs(X-Y)<TOL)
#ifndef MIN
#define MIN(X,Y) ((X<Y)?X:Y)
#endif
@@ -15,268 +15,6 @@
#endif
#define PI M_PI
-/* tollerance aware version */
-/* TODO: fix all =='s left */
-int segment_intersection_2d(double ax1, double ay1, double ax2, double ay2, double bx1, double by1, double bx2, double by2,
- double *x1, double *y1, double *x2, double *y2, double tol)
-{
- double adx, ady, bdx, bdy;
- double aa, bb;
- double tola, tolb;
- double d, d1, d2, ra, rb, t;
- int switched = 0;
-
- /* TODO: Works for points ? */
- G_debug(4, "segment_intersection_2d()");
- G_debug(4, " ax1 = %.18f, ay1 = %.18f", ax1, ay1);
- G_debug(4, " ax2 = %.18f, ay2 = %.18f", ax2, ay2);
- G_debug(4, " bx1 = %.18f, by1 = %.18f", bx1, by1);
- G_debug(4, " bx2 = %.18f, by2 = %.18f", bx2, by2);
-
-
- /* Check identical lines */
- if ((FEQUAL(ax1, bx1, tol) && FEQUAL(ay1, by1, tol) && FEQUAL(ax2, bx2, tol) && FEQUAL(ay2, by2, tol)) ||
- (FEQUAL(ax1, bx2, tol) && FEQUAL(ay1, by2, tol) && FEQUAL(ax2, bx1, tol) && FEQUAL(ay2, by1, tol))) {
- G_debug (2, " -> identical segments" );
- *x1 = ax1; *y1 = ay1;
- *x2 = ax2; *y2 = ay2;
- return 5;
- }
-
- /* 'Sort' lines by x1, x2, y1, y2 */
- if ( bx1 < ax1 )
- switched = 1;
- else if ( bx1 == ax1 ) {
- if ( bx2 < ax2 ) switched = 1;
- else if ( bx2 == ax2 ) {
- if ( by1 < ay1 ) switched = 1;
- else if ( by1 == ay1 ) {
- if ( by2 < ay2 ) switched = 1; /* by2 != ay2 (would be identical */
- }
- }
- }
-
- if (switched) {
- t = ax1; ax1 = bx1; bx1 = t; t = ay1; ay1 = by1; by1 = t;
- t = ax2; ax2 = bx2; bx2 = t; t = ay2; ay2 = by2; by2 = t;
- }
-
- adx = ax2 - ax1;
- ady = ay2 - ay1;
- bdx = bx2 - bx1;
- bdy = by2 - by1;
- d = (ax2-ax1)*(by1-by2) - (ay2-ay1)*(bx1-bx2);
- d1 = (bx1-ax1)*(by1-by2) - (by1-ay1)*(bx1-bx2);
- d2 = (ax2-ax1)*(by1-ay1) - (ay2-ay1)*(bx1-ax1);
-
- G_debug(2, " d = %.18f", d);
- G_debug(2, " d1 = %.18f", d1);
- G_debug(2, " d2 = %.18f", d2);
-
- tola = tol/MAX(fabs(adx), fabs(ady));
- tolb = tol/MAX(fabs(bdx), fabs(bdy));
- G_debug(2, " tol = %.18f", tol);
- G_debug(2, " tola = %.18f", tola);
- G_debug(2, " tolb = %.18f", tolb);
- if (!FZERO(d, tol)) {
- ra = d1/d;
- rb = d2/d;
-
- G_debug(2, " not parallel/collinear: ra = %.18f", ra);
- G_debug(2, " rb = %.18f", rb);
-
- if ((ra <= -tola) || (ra >= 1+tola) || (rb <= -tolb) || (rb >= 1+tolb)) {
- G_debug(2, " no intersection");
- return 0;
- }
-
- ra = MIN(MAX(ra, 0), 1);
- *x1 = ax1 + ra*(ax2 - ax1);
- *y1 = ay1 + ra*(ay2 - ay1);
-
- G_debug(2, " intersection %.18f, %.18f", *x1, *y1);
- return 1;
- }
-
- /* segments are parallel or collinear */
- G_debug (3, " -> parallel/collinear" );
-
- if ((!FZERO(d1, tol)) || (!FZERO(d2, tol))) { /* lines are parallel */
- G_debug (2, " -> parallel" );
- return 0;
- }
-
- /* segments are colinear. check for overlap */
-
-/* aa = adx*adx + ady*ady;
- bb = bdx*bdx + bdy*bdy;
-
- t = (ax1-bx1)*bdx + (ay1-by1)*bdy;*/
-
-
- /* Collinear vertical */
- /* original code assumed lines were not both vertical
- * so there is a special case if they are */
- if (FEQUAL(ax1, ax2, tol) && FEQUAL(bx1, bx2, tol) && FEQUAL(ax1, bx1, tol)) {
- G_debug(2, " -> collinear vertical");
- if (ay1 > ay2) { t=ay1; ay1=ay2; ay2=t; } /* to be sure that ay1 < ay2 */
- if (by1 > by2) { t=by1; by1=by2; by2=t; } /* to be sure that by1 < by2 */
- if (ay1 > by2 || ay2 < by1) {
- G_debug (2, " -> no intersection");
- return 0;
- }
-
- /* end points */
- if (FEQUAL(ay1, by2, tol)) {
- *x1 = ax1; *y1 = ay1;
- G_debug(2, " -> connected by end points");
- return 1; /* endpoints only */
- }
- if (FEQUAL(ay2, by1, tol)) {
- *x1 = ax2; *y1 = ay2;
- G_debug (2, " -> connected by end points");
- return 1; /* endpoints only */
- }
-
- /* general overlap */
- G_debug(3, " -> vertical overlap");
- /* a contains b */
- if (ay1 <= by1 && ay2 >= by2) {
- G_debug (2, " -> a contains b" );
- *x1 = bx1; *y1 = by1;
- *x2 = bx2; *y2 = by2;
- if (!switched)
- return 3;
- else
- return 4;
- }
- /* b contains a */
- if ( ay1 >= by1 && ay2 <= by2 ) {
- G_debug (2, " -> b contains a" );
- *x1 = ax1; *y1 = ay1;
- *x2 = ax2; *y2 = ay2;
- if ( !switched )
- return 4;
- else
- return 3;
- }
-
- /* general overlap, 2 intersection points */
- G_debug (2, " -> partial overlap" );
- if ( by1 > ay1 && by1 < ay2 ) { /* b1 in a */
- if ( !switched ) {
- *x1 = bx1; *y1 = by1;
- *x2 = ax2; *y2 = ay2;
- } else {
- *x1 = ax2; *y1 = ay2;
- *x2 = bx1; *y2 = by1;
- }
- return 2;
- }
-
- if ( by2 > ay1 && by2 < ay2 ) { /* b2 in a */
- if ( !switched ) {
- *x1 = bx2; *y1 = by2;
- *x2 = ax1; *y2 = ay1;
- } else {
- *x1 = ax1; *y1 = ay1;
- *x2 = bx2; *y2 = by2;
- }
- return 2;
- }
-
- /* should not be reached */
- G_warning(("Vect_segment_intersection() ERROR (collinear vertical segments)"));
- G_warning("%.15g %.15g", ax1, ay1);
- G_warning("%.15g %.15g", ax2, ay2);
- G_warning("x");
- G_warning("%.15g %.15g", bx1, by1);
- G_warning("%.15g %.15g", bx2, by2);
- return 0;
- }
-
- G_debug (2, " -> collinear non vertical" );
-
- /* Collinear non vertical */
- if ( ( bx1 > ax1 && bx2 > ax1 && bx1 > ax2 && bx2 > ax2 ) ||
- ( bx1 < ax1 && bx2 < ax1 && bx1 < ax2 && bx2 < ax2 ) ) {
- G_debug (2, " -> no intersection" );
- return 0;
- }
-
- /* there is overlap or connected end points */
- G_debug (2, " -> overlap/connected end points" );
-
- /* end points */
- if ( (ax1 == bx1 && ay1 == by1) || (ax1 == bx2 && ay1 == by2) ) {
- *x1 = ax1; *y1 = ay1;
- G_debug (2, " -> connected by end points");
- return 1;
- }
- if ( (ax2 == bx1 && ay2 == by1) || (ax2 == bx2 && ay2 == by2) ) {
- *x1 = ax2; *y1 = ay2;
- G_debug (2, " -> connected by end points");
- return 1;
- }
-
- if (ax1 > ax2) { t=ax1; ax1=ax2; ax2=t; t=ay1; ay1=ay2; ay2=t; } /* to be sure that ax1 < ax2 */
- if (bx1 > bx2) { t=bx1; bx1=bx2; bx2=t; t=by1; by1=by2; by2=t; } /* to be sure that bx1 < bx2 */
-
- /* a contains b */
- if ( ax1 <= bx1 && ax2 >= bx2 ) {
- G_debug (2, " -> a contains b" );
- *x1 = bx1; *y1 = by1;
- *x2 = bx2; *y2 = by2;
- if ( !switched )
- return 3;
- else
- return 4;
- }
- /* b contains a */
- if ( ax1 >= bx1 && ax2 <= bx2 ) {
- G_debug (2, " -> b contains a" );
- *x1 = ax1; *y1 = ay1;
- *x2 = ax2; *y2 = ay2;
- if ( !switched )
- return 4;
- else
- return 3;
- }
-
- /* general overlap, 2 intersection points (lines are not vertical) */
- G_debug (2, " -> partial overlap" );
- if ( bx1 > ax1 && bx1 < ax2 ) { /* b1 is in a */
- if ( !switched ) {
- *x1 = bx1; *y1 = by1;
- *x2 = ax2; *y2 = ay2;
- } else {
- *x1 = ax2; *y1 = ay2;
- *x2 = bx1; *y2 = by1;
- }
- return 2;
- }
- if ( bx2 > ax1 && bx2 < ax2 ) { /* b2 is in a */
- if ( !switched ) {
- *x1 = bx2; *y1 = by2;
- *x2 = ax1; *y2 = ay1;
- } else {
- *x1 = ax1; *y1 = ay1;
- *x2 = bx2; *y2 = by2;
- }
- return 2;
- }
-
- /* should not be reached */
- G_warning(("segment_intersection_2d() ERROR (collinear non vertical segments)"));
- G_warning("%.15g %.15g", ax1, ay1);
- G_warning("%.15g %.15g", ax2, ay2);
- G_warning("x");
- G_warning("%.15g %.15g", bx1, by1);
- G_warning("%.15g %.15g", bx2, by2);
-
- return 0;
-}
-
struct intersection_point {
double x;
double y;
@@ -414,21 +152,43 @@
return (aa->y < bb->y)?-1:((aa->y > bb->y)?1:0);
}
+/* O(Points->n_points) time */
+double get_epsilon(struct line_pnts *Points) {
+ int i, np;
+ double min, t;
+ double *x, *y;
+
+ np = Points->n_points;
+ x = Points->x;
+ y = Points->y;
+
+ min = MAX(fabs(x[1]-x[0]), fabs(y[1]-y[0]));
+ for (i = 1; i <= np-2; i++) {
+ t = MAX(fabs(x[i+1]-x[i]), fabs(y[i+1]-y[i]));
+ if ((t > 0) && (t < min)) {
+ min = t;
+ }
+ }
+
+ /* ??? is 0.001 ok ??? */
+ return min*0.000001;
+}
+
/* currently O(n*n); future implementation O(nlogn) */
struct seg_intersections* find_all_intersections(struct line_pnts *Points) {
int i, j, np;
int group, t;
int looped;
- double EPSILON = 1e-10;
+ double EPSILON = 0.00000001;
double *x, *y;
double x1, y1, x2, y2;
- int res;
+ int res, res2;
+ double x1_, y1_, x2_, y2_;
struct seg_intersections *si;
struct seg_intersection_list *il;
struct intersection_point **sorted;
G_debug(4, "find_all_intersections()");
- G_debug(4, " EPSILON=%.18f", EPSILON);
np = Points->n_points;
x = Points->x;
@@ -436,14 +196,20 @@
si = create_si_struct(np-1);
- looped = FEQUAL(x[0], x[np-1], EPSILON) && FEQUAL(y[0], y[np-1], EPSILON);
+ looped = ((x[0] == x[np-1]) && (y[0] == y[np-1]));
G_debug(4, " looped=%d", looped);
G_debug(4, " finding intersections...");
for (i = 0; i < np-1; i++) {
for (j = i+1; j < np-1; j++) {
- G_debug(3, " checking %d-%d %d-%d", i, i+1, j, j+1);
- res = segment_intersection_2d(x[i], y[i], x[i+1], y[i+1], x[j], y[j], x[j+1], y[j+1], &x1, &y1, &x2, &y2, EPSILON);
+ G_debug(3, " checking %d-%d %d-%d", i, i+1, j, j+1);
+ res = segment_intersection_2d_e(x[i], y[i], x[i+1], y[i+1], x[j], y[j], x[j+1], y[j+1], &x1, &y1, &x2, &y2);
+/* res = segment_intersection_2d(x[i], y[i], x[i+1], y[i+1], x[j], y[j], x[j+1], y[j+1], &x1, &y1, &x2, &y2, EPSILON);*/
+/* res2 = segment_intersection_2d_e(x[i], y[i], x[i+1], y[i+1], x[j], y[j], x[j+1], y[j+1], &x1_, &y1_, &x2_, &y2_);
+ if (res != res2) {
+ G_debug(0, "res=%d res2=%d", res, res2);
+ }*/
+
G_debug(3, " intersection type = %d", res);
if (res == 1) {
add_ipoint(Points, i, j, x1, y1, si);
@@ -492,8 +258,10 @@
t = group;
for (j = i-1; j >= 0; j--) {
if (!FEQUAL(sorted[j]->x, sorted[i]->x, EPSILON))
+/* if (!almost_equal(sorted[j]->x, sorted[i]->x, 16))*/
break;
if (FEQUAL(sorted[j]->y, sorted[i]->y, EPSILON)) {
+/* if (almost_equal(sorted[j]->y, sorted[i]->y, 16)) {*/
t = sorted[j]->group;
break;
}
@@ -508,7 +276,6 @@
G_debug(4, " postprocessing...done");
/* output contents of si */
-
for (i = 0; i < si->ilcount; i++) {
G_debug(4, "%d-%d :", i, i+1);
for (j = 0; j < si->il[i].count; j++) {
Added: grass-addons/vector/v.parallel2/e_intersect.c
===================================================================
--- grass-addons/vector/v.parallel2/e_intersect.c (rev 0)
+++ grass-addons/vector/v.parallel2/e_intersect.c 2008-08-09 21:15:04 UTC (rev 32656)
@@ -0,0 +1,565 @@
+#include <stdlib.h>
+#include <math.h>
+#include <grass/gis.h>
+#include <gmp.h>
+#include "e_intersect.h"
+
+#define SWAP(a,b) {t = a; a = b; b = t;}
+
+mpf_t p11, p12, p21, p22, t1, t2;
+mpf_t dd, dda, ddb, delta;
+
+int initialized = 0;
+
+void initialize_mpf_vars() {
+ mpf_set_default_prec(2048);
+
+ mpf_init(p11);
+ mpf_init(p12);
+ mpf_init(p21);
+ mpf_init(p22);
+
+ mpf_init(t1);
+ mpf_init(t2);
+
+ mpf_init(dd);
+ mpf_init(dda);
+ mpf_init(ddb);
+ mpf_init(delta);
+
+ initialized = 1;
+}
+
+/*
+ Caclulates:
+ |a11-b11 a12-b12|
+ |a21-b21 a22-b22|
+*/
+void det22(mpf_t rop, double a11, double b11, double a12, double b12, double a21, double b21, double a22, double b22) {
+ mpf_set_d(t1, a11);
+ mpf_set_d(t2, b11);
+ mpf_sub(p11, t1, t2);
+ mpf_set_d(t1, a12);
+ mpf_set_d(t2, b12);
+ mpf_sub(p12, t1, t2);
+ mpf_set_d(t1, a21);
+ mpf_set_d(t2, b21);
+ mpf_sub(p21, t1, t2);
+ mpf_set_d(t1, a22);
+ mpf_set_d(t2, b22);
+ mpf_sub(p22, t1, t2);
+
+ mpf_mul(t1, p11, p22);
+ mpf_mul(t2, p12, p21);
+ mpf_sub(rop, t1, t2);
+
+ return;
+}
+
+void swap(double *a, double *b) {
+ double t = *a;
+ *a = *b;
+ *b = t;
+ return;
+}
+
+/* multi-precision version */
+int segment_intersection_2d_e(double ax1, double ay1, double ax2, double ay2, double bx1, double by1, double bx2, double by2,
+ double *x1, double *y1, double *x2, double *y2)
+{
+ double t;
+ double max_ax, min_ax, max_ay, min_ay;
+ double max_bx, min_bx, max_by, min_by;
+ int sgn_d, sgn_da, sgn_db;
+ int vertical;
+ int f11, f12, f21, f22;
+ mp_exp_t exp;
+ char *s;
+
+ if (!initialized)
+ initialize_mpf_vars();
+
+ /* TODO: Works for points ? */
+ G_debug(3, "segment_intersection_2d_e()");
+ G_debug(4, " ax1 = %.18f, ay1 = %.18f", ax1, ay1);
+ G_debug(4, " ax2 = %.18f, ay2 = %.18f", ax2, ay2);
+ G_debug(4, " bx1 = %.18f, by1 = %.18f", bx1, by1);
+ G_debug(4, " bx2 = %.18f, by2 = %.18f", bx2, by2);
+
+ f11 = ((ax1 == bx1) && (ay1 == by1));
+ f12 = ((ax1 == bx2) && (ay1 == by2));
+ f21 = ((ax2 == bx1) && (ay2 == by1));
+ f22 = ((ax2 == bx2) && (ay2 == by2));
+
+ /* Check for identical segments */
+ if ((f11 && f22) || (f12 && f21)) {
+ G_debug (3, " identical segments" );
+ *x1 = ax1; *y1 = ay1;
+ *x2 = ax2; *y2 = ay2;
+ return 5;
+ }
+ /* Check for identical endpoints */
+ if (f11 || f12) {
+ G_debug (3, " connected by endpoints" );
+ *x1 = ax1; *y1 = ay1;
+ return 1;
+ }
+ if (f21 || f22) {
+ G_debug (3, " connected by endpoints" );
+ *x1 = ax2; *y1 = ay2;
+ return 1;
+ }
+
+ if ((fmax(ax1, ax2) < fmin(bx1, bx2)) || (fmax(bx1, bx2) < fmin(ax1, ax2))) {
+ G_debug(3, " no intersection (disjoint bounding boxes)");
+ return 0;
+ }
+ if ((fmax(ay1, ay2) < fmin(by1, by2)) || (fmax(by1, by2) < fmin(ay1, ay2))) {
+ G_debug(3, " no intersection (disjoint bounding boxes)");
+ return 0;
+ }
+
+ det22(dd, ax2, ax1, bx1, bx2, ay2, ay1, by1, by2);
+ sgn_d = mpf_sgn(dd);
+ if (sgn_d != 0) {
+ G_debug(3, " general position");
+
+ det22(dda, bx1, ax1, bx1, bx2, by1, ay1, by1, by2);
+ sgn_da = mpf_sgn(dda);
+
+ /*mpf_div(rra, dda, dd);
+ mpf_div(rrb, ddb, dd);
+ s = mpf_get_str(NULL, &exp, 10, 40, rra);
+ G_debug(4, " ra = %sE%d", (s[0]==0)?"0":s, exp);
+ s = mpf_get_str(NULL, &exp, 10, 24, rrb);
+ G_debug(4, " rb = %sE%d", (s[0]==0)?"0":s, exp);
+ */
+
+ if (sgn_d > 0) {
+ if ((sgn_da < 0) || (mpf_cmp(dda, dd) > 0)) {
+ G_debug(3, " no intersection");
+ return 0;
+ }
+
+ det22(ddb, ax2, ax1, bx1, ax1, ay2, ay1, by1, ay1);
+ sgn_db = mpf_sgn(ddb);
+ if ((sgn_db < 0) || (mpf_cmp(ddb, dd) > 0)) {
+ G_debug(3, " no intersection");
+ return 0;
+ }
+ }
+ else { /* if sgn_d < 0 */
+ if ((sgn_da > 0) || (mpf_cmp(dda, dd) < 0)) {
+ G_debug(3, " no intersection");
+ return 0;
+ }
+
+ det22(ddb, ax2, ax1, bx1, ax1, ay2, ay1, by1, ay1);
+ sgn_db = mpf_sgn(ddb);
+ if ((sgn_db > 0) || (mpf_cmp(ddb, dd) < 0)) {
+ G_debug(3, " no intersection");
+ return 0;
+ }
+ }
+
+ G_debug(3, " ra=%.17g rb=%.17g", mpf_get_d(dda)/mpf_get_d(dd), mpf_get_d(ddb)/mpf_get_d(dd));
+ G_debug(3, " sgn_d=%d sgn_da=%d sgn_db=%d cmp(dda,dd)=%d cmp(ddb,dd)=%d", sgn_d, sgn_da, sgn_db, mpf_cmp(dda, dd), mpf_cmp(ddb, dd));
+
+ mpf_set_d(delta, ax2 - ax1);
+ mpf_mul(t1, dda, delta);
+ mpf_div(t2, t1, dd);
+ *x1 = ax1 + mpf_get_d(t2);
+
+ mpf_set_d(delta, ay2 - ay1);
+ mpf_mul(t1, dda, delta);
+ mpf_div(t2, t1, dd);
+ *y1 = ay1 + mpf_get_d(t2);
+
+ G_debug(3, " intersection %.16g, %.16g", *x1, *y1);
+ return 1;
+ }
+
+ /* segments are parallel or collinear */
+ det22(dda, bx1, ax1, bx1, bx2, by1, ay1, by1, by2);
+ sgn_da = mpf_sgn(dda);
+ if (sgn_da != 0) {
+ /* segments are parallel */
+ G_debug(3, " parallel segments");
+ return 0;
+ }
+
+ /* segments are colinear. check for overlap */
+
+ /* swap endpoints if needed */
+ /* if segments are vertical, we swap x-coords with y-coords */
+ vertical = 0;
+ if (ax1 > ax2) {
+ SWAP(ax1, ax2);
+ SWAP(ay1, ay2);
+ } else if (ax1 == ax2) {
+ vertical = 1;
+ if (ay1 > ay2)
+ SWAP(ay1, ay2);
+ SWAP(ax1, ay1);
+ SWAP(ax2, ay2);
+ }
+ if (bx1 > bx2) {
+ SWAP(bx1, bx2);
+ SWAP(by1, by2);
+ } else if (bx1 == bx2) {
+ if (by1 > by2)
+ SWAP(by1, by2);
+ SWAP(bx1, by1);
+ SWAP(bx2, by2);
+ }
+
+ G_debug(3, " collinear segments");
+
+ if ((bx2 < ax1) || (bx1 > ax2)) {
+ G_debug(3, " no intersection");
+ return 0;
+ }
+
+ /* there is overlap or connected end points */
+ G_debug(3, " overlap");
+
+ /* a contains b */
+ if ((ax1 < bx1) && (ax2 > bx2)) {
+ G_debug(3, " a contains b");
+ if (!vertical) {
+ *x1 = bx1; *y1 = by1;
+ *x2 = bx2; *y2 = by2;
+ }
+ else {
+ *x1 = by1; *y1 = bx1;
+ *x2 = by2; *y2 = bx2;
+ }
+ return 3;
+ }
+
+ /* b contains a */
+ if ((ax1 > bx1) && (ax2 < bx2)) {
+ G_debug(3, " b contains a");
+ if (!vertical) {
+ *x1 = bx1; *y1 = by1;
+ *x2 = bx2; *y2 = by2;
+ }
+ else {
+ *x1 = by1; *y1 = bx1;
+ *x2 = by2; *y2 = bx2;
+ }
+ return 4;
+ }
+
+ /* general overlap, 2 intersection points */
+ G_debug(3, " partial overlap");
+ if ((bx1 > ax1) && (bx1 < ax2)) { /* b1 is in a */
+ if (!vertical) {
+ *x1 = bx1; *y1 = by1;
+ *x2 = ax2; *y2 = ay2;
+ }
+ else {
+ *x1 = by1; *y1 = bx1;
+ *x2 = ay2; *y2 = ax2;
+ }
+ return 2;
+ }
+ if ((bx2 > ax1) && (bx2 < ax2)) { /* b2 is in a */
+ if (!vertical) {
+ *x1 = bx2; *y1 = by2;
+ *x2 = ax1; *y2 = ay1;
+ }
+ else {
+ *x1 = by2; *y1 = bx2;
+ *x2 = ay1; *y2 = ax1;
+ }
+ return 2;
+ }
+
+ /* should not be reached */
+ G_warning(("segment_intersection_2d() ERROR (should not be reached)"));
+ G_warning("%.16g %.16g", ax1, ay1);
+ G_warning("%.16g %.16g", ax2, ay2);
+ G_warning("x");
+ G_warning("%.16g %.16g", bx1, by1);
+ G_warning("%.16g %.16g", bx2, by2);
+
+ return 0;
+}
+
+/* tollerance aware version */
+/* TODO: fix all ==s left */
+int segment_intersection_2d(double ax1, double ay1, double ax2, double ay2, double bx1, double by1, double bx2, double by2,
+ double *x1, double *y1, double *x2, double *y2, double tol)
+{
+ double tola, tolb;
+ double d, d1, d2, ra, rb, t;
+ int switched = 0;
+
+ /* TODO: Works for points ? */
+ G_debug(4, "segment_intersection_2d()");
+ G_debug(4, " ax1 = %.18f, ay1 = %.18f", ax1, ay1);
+ G_debug(4, " ax2 = %.18f, ay2 = %.18f", ax2, ay2);
+ G_debug(4, " bx1 = %.18f, by1 = %.18f", bx1, by1);
+ G_debug(4, " bx2 = %.18f, by2 = %.18f", bx2, by2);
+
+
+ /* Check identical lines */
+ if ((FEQUAL(ax1, bx1, tol) && FEQUAL(ay1, by1, tol) && FEQUAL(ax2, bx2, tol) && FEQUAL(ay2, by2, tol)) ||
+ (FEQUAL(ax1, bx2, tol) && FEQUAL(ay1, by2, tol) && FEQUAL(ax2, bx1, tol) && FEQUAL(ay2, by1, tol))) {
+ G_debug (2, " -> identical segments" );
+ *x1 = ax1; *y1 = ay1;
+ *x2 = ax2; *y2 = ay2;
+ return 5;
+ }
+
+ /* 'Sort' lines by x1, x2, y1, y2 */
+ if ( bx1 < ax1 )
+ switched = 1;
+ else if ( bx1 == ax1 ) {
+ if ( bx2 < ax2 ) switched = 1;
+ else if ( bx2 == ax2 ) {
+ if ( by1 < ay1 ) switched = 1;
+ else if ( by1 == ay1 ) {
+ if ( by2 < ay2 ) switched = 1; /* by2 != ay2 (would be identical */
+ }
+ }
+ }
+
+ if (switched) {
+ t = ax1; ax1 = bx1; bx1 = t; t = ay1; ay1 = by1; by1 = t;
+ t = ax2; ax2 = bx2; bx2 = t; t = ay2; ay2 = by2; by2 = t;
+ }
+
+ d = (ax2-ax1)*(by1-by2) - (ay2-ay1)*(bx1-bx2);
+ d1 = (bx1-ax1)*(by1-by2) - (by1-ay1)*(bx1-bx2);
+ d2 = (ax2-ax1)*(by1-ay1) - (ay2-ay1)*(bx1-ax1);
+
+ G_debug(2, " d = %.18f", d);
+ G_debug(2, " d1 = %.18f", d1);
+ G_debug(2, " d2 = %.18f", d2);
+
+ tola = tol/fmax(fabs(ax2-ax1), fabs(ay2-ay1));
+ tolb = tol/fmax(fabs(bx2-bx1), fabs(by2-by1));
+ G_debug(2, " tol = %.18f", tol);
+ G_debug(2, " tola = %.18f", tola);
+ G_debug(2, " tolb = %.18f", tolb);
+ if (!FZERO(d, tol)) {
+ ra = d1/d;
+ rb = d2/d;
+
+ G_debug(2, " not parallel/collinear: ra = %.18f", ra);
+ G_debug(2, " rb = %.18f", rb);
+
+ if ((ra <= -tola) || (ra >= 1+tola) || (rb <= -tolb) || (rb >= 1+tolb)) {
+ G_debug(2, " no intersection");
+ return 0;
+ }
+
+ ra = fmin(fmax(ra, 0), 1);
+ *x1 = ax1 + ra*(ax2 - ax1);
+ *y1 = ay1 + ra*(ay2 - ay1);
+
+ G_debug(2, " intersection %.18f, %.18f", *x1, *y1);
+ return 1;
+ }
+
+ /* segments are parallel or collinear */
+ G_debug (3, " -> parallel/collinear" );
+
+ if ((!FZERO(d1, tol)) || (!FZERO(d2, tol))) { /* lines are parallel */
+ G_debug (2, " -> parallel" );
+ return 0;
+ }
+
+ /* segments are colinear. check for overlap */
+
+/* aa = adx*adx + ady*ady;
+ bb = bdx*bdx + bdy*bdy;
+
+ t = (ax1-bx1)*bdx + (ay1-by1)*bdy;*/
+
+
+ /* Collinear vertical */
+ /* original code assumed lines were not both vertical
+ * so there is a special case if they are */
+ if (FEQUAL(ax1, ax2, tol) && FEQUAL(bx1, bx2, tol) && FEQUAL(ax1, bx1, tol)) {
+ G_debug(2, " -> collinear vertical");
+ if (ay1 > ay2) { t=ay1; ay1=ay2; ay2=t; } /* to be sure that ay1 < ay2 */
+ if (by1 > by2) { t=by1; by1=by2; by2=t; } /* to be sure that by1 < by2 */
+ if (ay1 > by2 || ay2 < by1) {
+ G_debug (2, " -> no intersection");
+ return 0;
+ }
+
+ /* end points */
+ if (FEQUAL(ay1, by2, tol)) {
+ *x1 = ax1; *y1 = ay1;
+ G_debug(2, " -> connected by end points");
+ return 1; /* endpoints only */
+ }
+ if (FEQUAL(ay2, by1, tol)) {
+ *x1 = ax2; *y1 = ay2;
+ G_debug (2, " -> connected by end points");
+ return 1; /* endpoints only */
+ }
+
+ /* general overlap */
+ G_debug(3, " -> vertical overlap");
+ /* a contains b */
+ if (ay1 <= by1 && ay2 >= by2) {
+ G_debug (2, " -> a contains b" );
+ *x1 = bx1; *y1 = by1;
+ *x2 = bx2; *y2 = by2;
+ if (!switched)
+ return 3;
+ else
+ return 4;
+ }
+ /* b contains a */
+ if ( ay1 >= by1 && ay2 <= by2 ) {
+ G_debug (2, " -> b contains a" );
+ *x1 = ax1; *y1 = ay1;
+ *x2 = ax2; *y2 = ay2;
+ if ( !switched )
+ return 4;
+ else
+ return 3;
+ }
+
+ /* general overlap, 2 intersection points */
+ G_debug (2, " -> partial overlap" );
+ if ( by1 > ay1 && by1 < ay2 ) { /* b1 in a */
+ if ( !switched ) {
+ *x1 = bx1; *y1 = by1;
+ *x2 = ax2; *y2 = ay2;
+ } else {
+ *x1 = ax2; *y1 = ay2;
+ *x2 = bx1; *y2 = by1;
+ }
+ return 2;
+ }
+
+ if ( by2 > ay1 && by2 < ay2 ) { /* b2 in a */
+ if ( !switched ) {
+ *x1 = bx2; *y1 = by2;
+ *x2 = ax1; *y2 = ay1;
+ } else {
+ *x1 = ax1; *y1 = ay1;
+ *x2 = bx2; *y2 = by2;
+ }
+ return 2;
+ }
+
+ /* should not be reached */
+ G_warning(("Vect_segment_intersection() ERROR (collinear vertical segments)"));
+ G_warning("%.15g %.15g", ax1, ay1);
+ G_warning("%.15g %.15g", ax2, ay2);
+ G_warning("x");
+ G_warning("%.15g %.15g", bx1, by1);
+ G_warning("%.15g %.15g", bx2, by2);
+ return 0;
+ }
+
+ G_debug (2, " -> collinear non vertical" );
+
+ /* Collinear non vertical */
+ if ( ( bx1 > ax1 && bx2 > ax1 && bx1 > ax2 && bx2 > ax2 ) ||
+ ( bx1 < ax1 && bx2 < ax1 && bx1 < ax2 && bx2 < ax2 ) ) {
+ G_debug (2, " -> no intersection" );
+ return 0;
+ }
+
+ /* there is overlap or connected end points */
+ G_debug (2, " -> overlap/connected end points" );
+
+ /* end points */
+ if ( (ax1 == bx1 && ay1 == by1) || (ax1 == bx2 && ay1 == by2) ) {
+ *x1 = ax1; *y1 = ay1;
+ G_debug (2, " -> connected by end points");
+ return 1;
+ }
+ if ( (ax2 == bx1 && ay2 == by1) || (ax2 == bx2 && ay2 == by2) ) {
+ *x1 = ax2; *y1 = ay2;
+ G_debug (2, " -> connected by end points");
+ return 1;
+ }
+
+ if (ax1 > ax2) { t=ax1; ax1=ax2; ax2=t; t=ay1; ay1=ay2; ay2=t; } /* to be sure that ax1 < ax2 */
+ if (bx1 > bx2) { t=bx1; bx1=bx2; bx2=t; t=by1; by1=by2; by2=t; } /* to be sure that bx1 < bx2 */
+
+ /* a contains b */
+ if ( ax1 <= bx1 && ax2 >= bx2 ) {
+ G_debug (2, " -> a contains b" );
+ *x1 = bx1; *y1 = by1;
+ *x2 = bx2; *y2 = by2;
+ if ( !switched )
+ return 3;
+ else
+ return 4;
+ }
+ /* b contains a */
+ if ( ax1 >= bx1 && ax2 <= bx2 ) {
+ G_debug (2, " -> b contains a" );
+ *x1 = ax1; *y1 = ay1;
+ *x2 = ax2; *y2 = ay2;
+ if ( !switched )
+ return 4;
+ else
+ return 3;
+ }
+
+ /* general overlap, 2 intersection points (lines are not vertical) */
+ G_debug (2, " -> partial overlap" );
+ if ( bx1 > ax1 && bx1 < ax2 ) { /* b1 is in a */
+ if ( !switched ) {
+ *x1 = bx1; *y1 = by1;
+ *x2 = ax2; *y2 = ay2;
+ } else {
+ *x1 = ax2; *y1 = ay2;
+ *x2 = bx1; *y2 = by1;
+ }
+ return 2;
+ }
+ if ( bx2 > ax1 && bx2 < ax2 ) { /* b2 is in a */
+ if ( !switched ) {
+ *x1 = bx2; *y1 = by2;
+ *x2 = ax1; *y2 = ay1;
+ } else {
+ *x1 = ax1; *y1 = ay1;
+ *x2 = bx2; *y2 = by2;
+ }
+ return 2;
+ }
+
+ /* should not be reached */
+ G_warning(("segment_intersection_2d() ERROR (collinear non vertical segments)"));
+ G_warning("%.15g %.15g", ax1, ay1);
+ G_warning("%.15g %.15g", ax2, ay2);
+ G_warning("x");
+ G_warning("%.15g %.15g", bx1, by1);
+ G_warning("%.15g %.15g", bx2, by2);
+
+ return 0;
+}
+
+#define N 52 /* double's mantisa size in bits */
+/* a and b are different in at most <bits> significant digits */
+int almost_equal(double a, double b, int bits) {
+ int ea, eb, e;
+
+ if (a == b)
+ return 1;
+
+ if (a == 0 || b == 0) {
+ /* return (0 < -N+bits); */
+ return (bits > N);
+ }
+
+ frexp(a, &ea);
+ frexp(b, &eb);
+ if (ea != eb)
+ return (bits > N + abs(ea-eb));
+ frexp(a-b, &e);
+ return (e < ea-N+bits);
+}
Added: grass-addons/vector/v.parallel2/e_intersect.h
===================================================================
--- grass-addons/vector/v.parallel2/e_intersect.h (rev 0)
+++ grass-addons/vector/v.parallel2/e_intersect.h 2008-08-09 21:15:04 UTC (rev 32656)
@@ -0,0 +1,15 @@
+#ifndef E_INTERSECT_H
+#define E_INTERSECT_H
+
+#define FZERO(X, TOL) (fabs(X)<TOL)
+#define FEQUAL(X, Y, TOL) (fabs(X-Y)<TOL)
+
+int segment_intersection_2d_e(double ax1, double ay1, double ax2, double ay2, double bx1, double by1, double bx2, double by2,
+ double *x1, double *y1, double *x2, double *y2);
+
+int segment_intersection_2d(double ax1, double ay1, double ax2, double ay2, double bx1, double by1, double bx2, double by2,
+ double *x1, double *y1, double *x2, double *y2, double tol);
+
+int almost_equal(double a, double b, int bits);
+
+#endif
Modified: grass-addons/vector/v.parallel2/main.c
===================================================================
--- grass-addons/vector/v.parallel2/main.c 2008-08-09 13:21:36 UTC (rev 32655)
+++ grass-addons/vector/v.parallel2/main.c 2008-08-09 21:15:04 UTC (rev 32656)
@@ -56,7 +56,7 @@
distb_opt->key = "minordistance";
distb_opt->type = TYPE_DOUBLE;
distb_opt->required = NO;
- dista_opt->options = "0-100000000";
+ distb_opt->options = "0-100000000";
distb_opt->multiple = NO;
distb_opt->description = _("Offset along minor axis in map units");
@@ -121,7 +121,7 @@
else
side = 0;
- Vect_set_open_level (2);
+ Vect_set_open_level(2);
Vect_open_old (&In, in_opt->answer, "");
Vect_open_new (&Out, out_opt->answer, 0);
Vect_copy_head_data (&In, &Out);
@@ -155,7 +155,7 @@
}
}
else {
- parallel_line_b(Points, da, db, dalpha, round_flag->answer, 1, tolerance, &oPoints, &iPoints, &inner_count);
+ Vect_line_buffer2(Points, da, db, dalpha, round_flag->answer, 1, tolerance, &oPoints, &iPoints, &inner_count);
Vect_write_line(&Out, ltype, oPoints, Cats);
for (j = 0; j < inner_count; j++) {
Vect_write_line(&Out, ltype, iPoints[j], Cats);
Modified: grass-addons/vector/v.parallel2/vlib_buffer.c
===================================================================
--- grass-addons/vector/v.parallel2/vlib_buffer.c 2008-08-09 13:21:36 UTC (rev 32655)
+++ grass-addons/vector/v.parallel2/vlib_buffer.c 2008-08-09 21:15:04 UTC (rev 32656)
@@ -343,6 +343,7 @@
for (i = 0; i <= np-2; i++)
{
+ G_debug(4, "point %d, segment %d-%d", i, i, i+1);
/* save the old values */
if (!round) {
a0 = a1;
@@ -379,13 +380,17 @@
norm_vector(0, 0, vx, vy, &tx, &ty);
elliptic_tangent(side*tx, side*ty, da, db, dalpha, &tx, &ty);
Vect_append_point(nPoints, x[i] + wx + tx, y[i] + wy + ty, 0);
+ G_debug(4, " append point (c) x=%.16f y=%.16f", x[i] + wx + tx, y[i] + wy + ty);
Vect_append_point(nPoints, nx + tx, ny + ty, 0); /* nx == x[i] + vx, ny == y[i] + vy */
+ G_debug(4, " append point (c) x=%.16f y=%.16f", nx + tx, ny + ty);
}
if ((!turns180) && (!round) && (!inner_corner)) {
res = line_intersection(a0, b0, c0, a1, b1, c1, &rx, &ry);
- if (res == 1)
+ if (res == 1) {
Vect_append_point(nPoints, rx, ry, 0);
+ G_debug(4, " append point (o) x=%.16f y=%.16f", rx, ry);
+ }
else
G_fatal_error("unexpected result of line_intersection()");
}
@@ -408,15 +413,19 @@
nsegments = (int)(delta_phi/angular_tol) + 1;
angular_step = side*(delta_phi/nsegments);
+ phi1 += angular_step;
for (j = 1; j <= nsegments-1; j++) {
elliptic_transform(cos(phi1), sin(phi1), da, db, dalpha, &tx, &ty);
Vect_append_point(nPoints, x[i] + tx, y[i] + ty, 0);
+ G_debug(4, " append point (r) x=%.16f y=%.16f", x[i] + tx, y[i] + ty);
phi1 += angular_step;
}
}
Vect_append_point(nPoints, nx, ny, 0);
+ G_debug(4, " append point (s) x=%.16f y=%.16f", nx, ny);
Vect_append_point(nPoints, mx, my, 0);
+ G_debug(4, " append point (s) x=%.16f y=%.16f", mx, my);
}
/* close the output line */
@@ -465,7 +474,6 @@
while (1) {
Vect_append_point(nPoints, vert0->x, vert0->y, 0);
G_debug(4, "ec: v0=%d, v=%d, eside=%d, edge->v1=%d, edge->v2=%d", v0, v, eside, edge->v1, edge->v2);
- G_debug(4, "ec: edge=%X, first=%X", edge, first);
G_debug(4, "ec: append point x=%.18f y=%.18f", vert0->x, vert0->y);
/* mark current edge as visited on the appropriate side */
@@ -502,17 +510,35 @@
/* if line end is reached (no other edges at curr vertex) */
if (opt_flag) {
- if (stop_at_line_end)
+ if (stop_at_line_end) {
+ G_debug(3, " end has been reached, will stop here");
break;
+ }
else {
opt_j = 0; /* the only edge of vert is vert->edges[0] */
opt_side = -eside; /* go to the other side of the current edge */
+ G_debug(3, " end has been reached, turning around");
}
}
+ /* break condition */
if ((vert->edges[opt_j] == first) && (opt_side == side))
break;
-
+ if (opt_side == 1) {
+ if (vert->edges[opt_j]->visited_right) {
+ G_warning("next edge was visited but it is not the first one !!! breaking loop");
+ G_debug(4, "ec: v0=%d, v=%d, eside=%d, edge->v1=%d, edge->v2=%d", v, (edge->v1 == v)?(edge->v2):(edge->v1), opt_side, vert->edges[opt_j]->v1, vert->edges[opt_j]->v2);
+ break;
+ }
+ }
+ else {
+ if (vert->edges[opt_j]->visited_left) {
+ G_warning("next edge was visited but it is not the first one !!! breaking loop");
+ G_debug(4, "ec: v0=%d, v=%d, eside=%d, edge->v1=%d, edge->v2=%d", v, (edge->v1 == v)?(edge->v2):(edge->v1), opt_side, vert->edges[opt_j]->v1, vert->edges[opt_j]->v2);
+ break;
+ }
+ }
+
edge = vert->edges[opt_j];
eside = opt_side;
v0 = v;
@@ -683,6 +709,29 @@
return 0;
}
+/* returns 0 for ccw, non-zero for cw
+*/
+int get_polygon_orientation(const double *x, const double *y, int n) {
+ double x1,y1,x2,y2;
+ double area;
+
+ x2 = x[n-1];
+ y2 = y[n-1];
+
+ area = 0;
+ while (--n >= 0)
+ {
+ x1 = x2;
+ y1 = y2;
+
+ x2 = *x++;
+ y2 = *y++;
+
+ area += (y2+y1)*(x2-x1);
+ }
+ return signbit(area);
+}
+
/* internal */
void add_line_to_array(struct line_pnts *Points, struct line_pnts ***arrPoints, int *count, int *allocated, int more) {
if (*allocated == *count) {
@@ -694,29 +743,41 @@
return;
}
-void parallel_line_b(struct line_pnts *Points, double da, double db, double dalpha, int round, int caps, double tol, struct line_pnts **oPoints, struct line_pnts ***iPoints, int *inner_count) {
- struct planar_graph *pg, *pg2;
- struct line_pnts *tPoints, *sPoints, *cPoints;
+void destroy_lines_array(struct line_pnts **arr, int count) {
+ int i;
+ for (i = 0; i < count; i++)
+ Vect_destroy_line_struct(arr[i]);
+ G_free(arr);
+}
+
+/* area_outer and area_isles[i] must be closed non self-intersecting lines
+ side: 0 - auto, 1 - right, -1 left
+ */
+void buffer_lines(struct line_pnts *area_outer, struct line_pnts **area_isles, int isles_count, int side, double da, double db, double dalpha, int round, int caps, double tol, struct line_pnts **oPoints, struct line_pnts ***iPoints, int *inner_count) {
+ struct planar_graph *pg2;
+ struct line_pnts *sPoints, *cPoints;
struct line_pnts **arrPoints;
int i, count = 0;
int res, winding;
+ int auto_side;
int more = 8;
int allocated = 0;
double px, py;
- G_debug(4, "parallel_line_b()");
+ G_debug(4, "buffer_lines()");
+ auto_side = (side == 0);
+
/* initializations */
- tPoints = Vect_new_line_struct();
sPoints = Vect_new_line_struct();
cPoints = Vect_new_line_struct();
arrPoints = NULL;
- pg = pg_create(Points);
/* outer contour */
*oPoints = Vect_new_line_struct();
- extract_outer_contour(pg, 0, tPoints);
- convolution_line(tPoints, da, db, dalpha, RIGHT_SIDE, round, caps, tol, sPoints);
+ if (auto_side)
+ side = get_polygon_orientation(area_outer->x, area_outer->y, area_outer->n_points-1)?LEFT_SIDE:RIGHT_SIDE;
+ convolution_line(area_outer, da, db, dalpha, side, round, caps, tol, sPoints);
pg2 = pg_create(sPoints);
extract_outer_contour(pg2, 0, *oPoints);
res = extract_inner_contour(pg2, &winding, cPoints);
@@ -730,21 +791,22 @@
pg_destroy_struct(pg2);
/* inner contours */
- res = extract_inner_contour(pg, &winding, tPoints);
- while (res != 0) {
- convolution_line(tPoints, da, db, dalpha, RIGHT_SIDE, round, caps, tol, sPoints);
+ for (i = 0; i < isles_count; i++) {
+ if (auto_side)
+ side = get_polygon_orientation(area_isles[i]->x, area_isles[i]->y, area_isles[i]->n_points-1)?RIGHT_SIDE:LEFT_SIDE;
+ convolution_line(area_isles[i], da, db, dalpha, side, round, caps, tol, sPoints);
pg2 = pg_create(sPoints);
extract_outer_contour(pg2, 0, cPoints);
res = extract_inner_contour(pg2, &winding, cPoints);
while (res != 0) {
if (winding == -1) {
/* we need to check if the area is in the buffer.
- I simplfied convolution_line, so that it runs faster,
+ I've simplfied convolution_line(), so that it runs faster,
however that leads to ocasional problems */
- if (Vect_point_in_poly(cPoints->x[0], cPoints->y[0], tPoints)) {
+ if (Vect_point_in_poly(cPoints->x[0], cPoints->y[0], area_isles[i])) {
if (Vect_get_point_in_poly(cPoints, &px, &py) != 0)
G_fatal_error("Vect_get_point_in_poly() failed.");
- if (!point_in_buf(tPoints, px, py, da, db, dalpha)) {
+ if (!point_in_buf(area_isles[i], px, py, da, db, dalpha)) {
add_line_to_array(cPoints, &arrPoints, &count, &allocated, more);
cPoints = Vect_new_line_struct();
}
@@ -777,22 +839,160 @@
Vect_copy_xyz_to_pnts(arrPoints[count], arrPoints2[i]->x, arrPoints2[i]->y, arrPoints2[i]->z, arrPoints2[i]->n_points);
count++;
} */
-
- res = extract_inner_contour(pg, &winding, tPoints);
}
arrPoints = G_realloc(arrPoints, count*sizeof(struct line_pnts *));
*inner_count = count;
*iPoints = arrPoints;
- Vect_destroy_line_struct(tPoints);
Vect_destroy_line_struct(sPoints);
Vect_destroy_line_struct(cPoints);
+
+ return;
+}
+
+
+/*!
+ \fn void Vect_line_buffer(struct line_pnts *Points, double da, double db, double dalpha, int round, int caps, double tol, struct line_pnts **oPoints, struct line_pnts ***iPoints, int *inner_count) {
+ \brief Creates buffer around the line Points.
+ \param InPoints input line
+ \param da distance along major axis
+ \param da distance along minor axis
+ \param dalpha angle between 0x and major axis
+ \param round make corners round
+ \param caps add caps at line ends
+ \param tol maximum distance between theoretical arc and output segments
+ \param oPoints output polygon outer border (ccw order)
+ \param inner_count number of holes
+ \param iPoints array of output polygon's holes (cw order)
+*/
+void Vect_line_buffer2(struct line_pnts *Points, double da, double db, double dalpha, int round, int caps, double tol, struct line_pnts **oPoints, struct line_pnts ***iPoints, int *inner_count) {
+ struct planar_graph *pg;
+ struct line_pnts *tPoints, *outer;
+ struct line_pnts **isles;
+ int isles_count = 0;
+ int res, winding;
+ int more = 8;
+ int isles_allocated = 0;
+
+ G_debug(4, "Vect_line_buffer()");
+
+ /* initializations */
+ tPoints = Vect_new_line_struct();
+ isles = NULL;
+ pg = pg_create(Points);
+
+ /* outer contour */
+ outer = Vect_new_line_struct();
+ extract_outer_contour(pg, 0, outer);
+
+ /* inner contours */
+ res = extract_inner_contour(pg, &winding, tPoints);
+ while (res != 0) {
+ add_line_to_array(tPoints, &isles, &isles_count, &isles_allocated, more);
+ tPoints = Vect_new_line_struct();
+ res = extract_inner_contour(pg, &winding, tPoints);
+ }
+
+ buffer_lines(outer, isles, isles_count, RIGHT_SIDE, da, db, dalpha, round, caps, tol, oPoints, iPoints, inner_count);
+
+ Vect_destroy_line_struct(tPoints);
+ Vect_destroy_line_struct(outer);
+ destroy_lines_array(isles, isles_count);
pg_destroy_struct(pg);
+
+ return;
+}
+
+void Vect_area_buffer2(struct Map_info *Map, int area, double da, double db, double dalpha, int round, int caps, double tol, struct line_pnts **oPoints, struct line_pnts ***iPoints, int *inner_count) {
+ struct line_pnts *tPoints, *outer;
+ struct line_pnts **isles;
+ int isles_count = 0;
+ int i, res, winding, isle;
+ int more = 8;
+ int isles_allocated = 0;
+ double max = fmax(da, db);
+ G_debug(4, "Vect_area_buffer()");
+
+ /* initializations */
+ tPoints = Vect_new_line_struct();
+ isles_count = Vect_get_area_num_isles(Map, area);
+ isles_allocated = isles_count;
+ isles = G_malloc(isles_allocated*sizeof(struct line_pnts *));
+
+ /* outer contour */
+ outer = Vect_new_line_struct();
+ Vect_get_area_points(Map, area, outer);
+ Vect_append_point(outer, outer->x[0], outer->y[0], outer->z[0]);
+
+ /* inner contours */
+ for (i = 0; i < isles_count; i++) {
+ isle = Vect_get_area_isle(Map, area, i);
+ Vect_get_isle_points(Map, isle, tPoints);
+
+ /* Check if the isle is big enough */
+ /*
+ if (Vect_line_length(tPoints) < 2*PI*max)
+ continue;
+ */
+ Vect_append_point(tPoints, tPoints->x[0], tPoints->y[0], tPoints->z[0]);
+ add_line_to_array(tPoints, &isles, &isles_count, &isles_allocated, more);
+ tPoints = Vect_new_line_struct();
+ }
+
+ buffer_lines(outer, isles, isles_count, 0, da, db, dalpha, round, caps, tol, oPoints, iPoints, inner_count);
+
+ Vect_destroy_line_struct(tPoints);
+ Vect_destroy_line_struct(outer);
+ destroy_lines_array(isles, isles_count);
+
return;
}
+/*!
+ \fn void Vect_point_buffer(double px, double py, double da, double db, double dalpha, int round, double tol, struct line_pnts *nPoints) {
+ \brief Creates buffer around the point (px,py).
+ \param px input point x-coordinate
+ \param py input point y-coordinate
+ \param da distance along major axis
+ \param da distance along minor axis
+ \param dalpha angle between 0x and major axis
+ \param round make corners round
+ \param tol maximum distance between theoretical arc and output segments
+ \param nPoints output polygon outer border (ccw order)
+*/
+void Vect_point_buffer2(double px, double py, double da, double db, double dalpha, int round, double tol, struct line_pnts *nPoints) {
+ double tx, ty;
+ double angular_tol, angular_step, phi1;
+ int j, nsegments;
+
+ Vect_reset_line(nPoints);
+
+ dalpha *= PI/180; /* convert dalpha from degrees to radians */
+
+ if (round || (!round)) {
+ angular_tol = angular_tolerance(tol, da, db);
+
+ nsegments = (int)(2*PI/angular_tol) + 1;
+ angular_step = 2*PI/nsegments;
+
+ phi1 = 0;
+ for (j = 0; j < nsegments; j++) {
+ elliptic_transform(cos(phi1), sin(phi1), da, db, dalpha, &tx, &ty);
+ Vect_append_point(nPoints, px + tx, py + ty, 0);
+ phi1 += angular_step;
+ }
+ }
+ else {
+
+ }
+
+ /* close the output line */
+ Vect_append_point(nPoints, nPoints->x[0], nPoints->y[0], nPoints->z[0]);
+}
+
+
/*
\fn void Vect_line_parallel2 ( struct line_pnts *InPoints, double distance, double tolerance, int rm_end,
struct line_pnts *OutPoints )
@@ -815,112 +1015,3 @@
*/
return;
}
-
-/*!
- \fn void Vect_line_buffer ( struct line_pnts *InPoints, double distance, double tolerance,
- struct line_pnts *OutPoints )
- \brief Create buffer around the line line.
- Buffer is closed counter clockwise polygon.
- Warning: output line may contain loops!
- \param InPoints input line
- \param distance
- \param tolerance maximum distance between theoretical arc and polygon segments
- \param OutPoints output line
-*/
-void
-Vect_line_buffer ( struct line_pnts *InPoints, double distance, double tolerance,
- struct line_pnts *OutPoints )
-{
- double dangle;
- int side, npoints;
- static struct line_pnts *Points = NULL;
- static struct line_pnts *PPoints = NULL;
-
- distance = fabs (distance );
-
- dangle = 2 * acos( 1-tolerance/fabs(distance) ); /* angle step */
-
- if ( Points == NULL )
- Points = Vect_new_line_struct();
-
- if ( PPoints == NULL )
- PPoints = Vect_new_line_struct();
-
- /* Copy and prune input */
- Vect_reset_line ( Points );
- Vect_append_points ( Points, InPoints, GV_FORWARD );
- Vect_line_prune ( Points );
-
- Vect_reset_line ( OutPoints );
-
- npoints = Points->n_points;
- if ( npoints <= 0 ) {
- return;
- } else if ( npoints == 1 ) { /* make a circle */
- double angle, x, y;
-
- for ( angle = 0; angle < 2*PI; angle += dangle ) {
- x = Points->x[0] + distance * cos( angle );
- y = Points->y[0] + distance * sin( angle );
- Vect_append_point ( OutPoints, x, y, 0 );
- }
- /* Close polygon */
- Vect_append_point ( OutPoints, OutPoints->x[0], OutPoints->y[0], 0 );
- } else { /* 2 and more points */
- for ( side = 0; side < 2; side++ ) {
- double angle, sangle;
- double lx1, ly1, lx2, ly2;
- double x, y, nx, ny, sx, sy, ex, ey;
-
- /* Parallel on one side */
- if ( side == 0 ) {
- Vect_line_parallel ( Points, distance, tolerance, 0, PPoints );
- Vect_append_points ( OutPoints, PPoints, GV_FORWARD );
- } else {
- Vect_line_parallel ( Points, -distance, tolerance, 0, PPoints );
- Vect_append_points ( OutPoints, PPoints, GV_BACKWARD );
- }
-
- /* Arc at the end */
- /* 2 points at theend of original line */
- if ( side == 0 ) {
- lx1 = Points->x[npoints-2];
- ly1 = Points->y[npoints-2];
- lx2 = Points->x[npoints-1];
- ly2 = Points->y[npoints-1];
- } else {
- lx1 = Points->x[1];
- ly1 = Points->y[1];
- lx2 = Points->x[0];
- ly2 = Points->y[0];
- }
-
- /* normalized vector */
- norm_vector( lx1, ly1, lx2, ly2, &nx, &ny);
-
- /* starting point */
- sangle = atan2 ( -nx, ny ); /* starting angle */
- sx = lx2 + ny * distance;
- sy = ly2 - nx * distance;
-
- /* end point */
- ex = lx2 - ny * distance;
- ey = ly2 + nx * distance;
-
- Vect_append_point ( OutPoints, sx, sy, 0 );
-
- /* arc */
- for ( angle = dangle; angle < PI; angle += dangle ) {
- x = lx2 + distance * cos( sangle + angle );
- y = ly2 + distance * sin( sangle + angle );
- Vect_append_point ( OutPoints, x, y, 0 );
- }
-
- Vect_append_point ( OutPoints, ex, ey, 0 );
- }
-
- /* Close polygon */
- Vect_append_point ( OutPoints, OutPoints->x[0], OutPoints->y[0], 0 );
- }
- Vect_line_prune ( OutPoints );
-}
More information about the grass-commit
mailing list