[GRASS-SVN] r32228 - grass-addons/vector/v.parallel2
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed Jul 23 11:15:36 EDT 2008
Author: rmatev
Date: 2008-07-23 11:15:36 -0400 (Wed, 23 Jul 2008)
New Revision: 32228
Added:
grass-addons/vector/v.parallel2/dgraph.c
grass-addons/vector/v.parallel2/dgraph.h
Modified:
grass-addons/vector/v.parallel2/main.c
grass-addons/vector/v.parallel2/vlib_buffer.c
Log:
Implemented contour extraction algorithm
Added: grass-addons/vector/v.parallel2/dgraph.c
===================================================================
--- grass-addons/vector/v.parallel2/dgraph.c (rev 0)
+++ grass-addons/vector/v.parallel2/dgraph.c 2008-07-23 15:15:36 UTC (rev 32228)
@@ -0,0 +1,691 @@
+#include <stdlib.h>
+#include <math.h>
+#include <grass/Vect.h>
+#include <grass/gis.h>
+#include "dgraph.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
+#ifndef MAX
+ #define MAX(X,Y) ((X>Y)?X:Y)
+#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(("Vect_segment_intersection() 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;
+ int group; /* IPs with very similar dist will be in the same group */
+};
+
+struct seg_intersection {
+ int with; /* second segment */
+ int ip; /* index of the IP */
+ double dist; /* distance from first point of first segment to intersection point (IP) */
+};
+
+struct seg_intersection_list {
+ int count;
+ int allocated;
+ struct seg_intersection *a;
+};
+
+struct seg_intersections {
+ int ipcount;
+ int ipallocated;
+ struct intersection_point *ip;
+ int ilcount;
+ struct seg_intersection_list *il;
+ int group_count;
+};
+
+struct seg_intersections* create_si_struct(int segments_count) {
+ struct seg_intersections *si;
+ int i;
+
+ si = G_malloc(sizeof(struct seg_intersections));
+ si->ipcount = 0;
+ si->ipallocated = segments_count + 16;
+ si->ip = G_malloc((si->ipallocated)*sizeof(struct intersection_point));
+ si->ilcount = segments_count;
+ si->il = G_malloc(segments_count*sizeof(struct seg_intersection_list));
+ for (i = 0; i < segments_count; i++) {
+ si->il[i].count = 0;
+ si->il[i].allocated = 0;
+ si->il[i].a = NULL;
+ }
+
+ return si;
+}
+
+void destroy_si_struct(struct seg_intersections *si) {
+ int i;
+
+ for (i = 0; i < si->ilcount; i++)
+ G_free(si->il[i].a);
+ G_free(si->il);
+ G_free(si->ip);
+ G_free(si);
+
+ return;
+}
+
+/* internal use */
+void add_ipoint1(struct seg_intersection_list *il, int with, double dist, int ip) {
+ struct seg_intersection *s;
+
+ if (il->count == il->allocated) {
+ il->allocated += 4;
+ il->a = G_realloc(il->a, (il->allocated)*sizeof(struct seg_intersection));
+ }
+ s = &(il->a[il->count]);
+ s->with = with;
+ s->ip = ip;
+ s->dist = dist;
+ il->count++;
+
+ return;
+}
+
+/* adds intersection point to the structure */
+void add_ipoint(struct line_pnts *Points, int first_seg, int second_seg, double x, double y, struct seg_intersections *si) {
+ struct intersection_point *t;
+ int ip;
+
+ G_debug(4, "add_ipoint()");
+
+ if (si->ipcount == si->ipallocated) {
+ si->ipallocated += 16;
+ si->ip = G_realloc(si->ip, (si->ipallocated)*sizeof(struct intersection_point));
+ }
+ ip = si->ipcount;
+ t = &(si->ip[ip]);
+ t->x = x;
+ t->y = y;
+ t->group = -1;
+ si->ipcount++;
+
+ add_ipoint1(&(si->il[first_seg]), second_seg, LENGTH((Points->x[first_seg] - x), (Points->y[first_seg] - y)), ip);
+ if (second_seg >= 0)
+ add_ipoint1(&(si->il[second_seg]), first_seg, LENGTH((Points->x[second_seg] - x), (Points->y[second_seg] - y)), ip);
+}
+
+void sort_intersection_list(struct seg_intersection_list *il) {
+ int n, i, j, min;
+ struct seg_intersection t;
+
+ G_debug(4, "sort_intersection_list()");
+
+ n = il->count;
+ G_debug(4, " n=%d", n);
+ for (i = 0; i < n-1; i++) {
+ min = i;
+ for (j = i+1; j < n; j++) {
+ if (il->a[j].dist < il->a[min].dist) {
+ min = j;
+ }
+ }
+ if (min != i) {
+ t = il->a[i];
+ il->a[i] = il->a[min];
+ il->a[min] = t;
+ }
+ }
+
+ return;
+}
+
+static int compare(const void *a, const void *b)
+{
+ struct intersection_point *aa, *bb;
+ aa = *((struct intersection_point **)a);
+ bb = *((struct intersection_point **)b);
+
+ if (aa->x < bb->x)
+ return -1;
+ else if (aa->x > bb->x)
+ return 1;
+ else
+ return (aa->y < bb->y)?-1:((aa->y > bb->y)?1:0);
+}
+
+/* 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 *x, *y;
+ double x1, y1, x2, y2;
+ int res;
+ 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;
+ y = Points->y;
+
+ si = create_si_struct(np-1);
+
+ looped = FEQUAL(x[0], x[np-1], EPSILON) && FEQUAL(y[0], y[np-1], EPSILON);
+ 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, " intersection type = %d", res);
+ if (res == 1) {
+ add_ipoint(Points, i, j, x1, y1, si);
+ } else if ((res >= 2) && (res <= 5)) {
+ add_ipoint(Points, i, j, x1, y1, si);
+ add_ipoint(Points, i, j, x2, y2, si);
+ }
+ }
+ }
+ if (!looped) {
+ /* these are not really intersection points */
+ add_ipoint(Points, 0, -1, Points->x[0], Points->y[0], si);
+ add_ipoint(Points, np-2, -1, Points->x[np-1], Points->y[np-1], si);
+ }
+ G_debug(4, " finding intersections...done");
+
+ G_debug(4, " postprocessing...");
+ if (si->ipallocated > si->ipcount) {
+ si->ipallocated = si->ipcount;
+ si->ip = G_realloc(si->ip, (si->ipcount)*sizeof(struct intersection_point));
+ }
+ for (i = 0; i < si->ilcount; i++) {
+ il = &(si->il[i]);
+ if (il->allocated > il->count) {
+ il->allocated = il->count;
+ il->a = G_realloc(il->a, (il->count)*sizeof(struct seg_intersection));
+ }
+
+ if (il->count > 0) {
+ sort_intersection_list(il);
+ /* is it ok to use qsort here ? */
+ }
+ }
+
+ /* si->ip will not be reallocated any more so we can use pointers */
+ sorted = G_malloc((si->ipcount)*sizeof(struct intersection_point *));
+ for (i = 0; i < si->ipcount; i++)
+ sorted[i] = &(si->ip[i]);
+
+ qsort(sorted, si->ipcount, sizeof(struct intersection_point *), compare);
+
+ /* assign groups */
+ group = 0; /* next available group number */
+ for (i = 0; i < si->ipcount; i++) {
+
+ t = group;
+ for (j = i-1; j >= 0; j--) {
+ if (!FEQUAL(sorted[j]->x, sorted[i]->x, EPSILON))
+ break;
+ if (FEQUAL(sorted[j]->y, sorted[i]->y, EPSILON)) {
+ t = sorted[j]->group;
+ break;
+ }
+ }
+ G_debug(4, " group=%d, ip=%d", t, (int)(sorted[i] - &(si->ip[0])));
+ sorted[i]->group = t;
+ if (t == group)
+ group++;
+ }
+ si->group_count = group;
+
+ 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++) {
+ G_debug(4, " %d-%d, group=%d", si->il[i].a[j].with, si->il[i].a[j].with+1, si->ip[si->il[i].a[j].ip].group);
+ G_debug(4, " dist=%.18f", si->il[i].a[j].dist);
+ G_debug(4, " x=%.18f, y=%.18f", si->ip[si->il[i].a[j].ip].x, si->ip[si->il[i].a[j].ip].y);
+ }
+ }
+
+ return si;
+}
+
+/* create's graph with n vertices and allocates memory for e edges */
+/* trying to add more than e edges, produces fatal error */
+struct planar_graph* pg_create_struct(int n, int e) {
+ struct planar_graph *pg;
+ int i;
+
+ pg = G_malloc(sizeof(struct planar_graph));
+ pg->vcount = n;
+ pg->v = G_malloc(n*sizeof(struct pg_vertex));
+ memset(pg->v, 0, n*sizeof(struct pg_vertex));
+ pg->ecount = 0;
+ pg->eallocated = MAX(e, 0);
+ pg->e = NULL;
+ pg->e = G_malloc(e*sizeof(struct pg_edge));
+
+ return pg;
+}
+
+void pg_destroy_struct(struct planar_graph *pg) {
+ int i;
+
+ for (i = 0; i < pg->vcount; i++) {
+ G_free(pg->v[i].edges);
+ G_free(pg->v[i].angles);
+ }
+ G_free(pg->v);
+ G_free(pg->e);
+ G_free(pg);
+}
+
+/* v1 and v2 must be valid and v1 must be smaller than v2 (v1 < v2) */
+int pg_existsedge(struct planar_graph *pg, int v1, int v2) {
+ struct pg_vertex *v;
+ struct pg_edge *e;
+ int i, ecount;
+
+ if (pg->v[v1].ecount <= pg->v[v2].ecount)
+ v = &(pg->v[v1]);
+ else
+ v = &(pg->v[v2]);
+
+ ecount = v->ecount;
+ for (i = 0; i < ecount; i++) {
+ e = v->edges[i];
+ if ((e->v1 == v1) && (e->v2 == v2))
+ return 1;
+ }
+
+ return 0;
+}
+
+/* for internal use */
+void pg_addedge1(struct pg_vertex *v, struct pg_edge *e) {
+ if (v->ecount == v->eallocated) {
+ v->eallocated += 4;
+ v->edges = G_realloc(v->edges, (v->eallocated)*sizeof(struct pg_edge *));
+ }
+ v->edges[v->ecount] = e;
+ v->ecount++;
+}
+
+void pg_addedge(struct planar_graph *pg, int v1, int v2) {
+ struct pg_edge *e;
+ int t;
+
+ G_debug(4, "pg_addedge(), v1=%d, v2=%d", v1, v2);
+
+ if ((v1 == v2) || (v1 < 0) || (v1 >= pg->vcount) || (v2 < 0) || (v2 >= pg->vcount)) {
+ G_warning(" v1=%d, v2=%d, pg->vcount=%d", v1, v2, pg->vcount);
+ G_fatal_error(" pg_addedge(), v1 and/or v2 is invalid");
+ return;
+ }
+
+ if (v2 < v1) {
+ t = v1;
+ v1 = v2;
+ v2 = t;
+ }
+
+ if (pg_existsedge(pg, v1, v2))
+ return;
+
+ if (pg->ecount == pg->eallocated) {
+ G_fatal_error("Trying to add more edges to the planar_graph than the initial allocation size allows");
+ }
+ e = &(pg->e[pg->ecount]);
+ e->v1 = v1;
+ e->v2 = v2;
+ e->visited_left = 0;
+ e->visited_right = 0;
+ pg->ecount++;
+ pg_addedge1(&(pg->v[v1]), e);
+ pg_addedge1(&(pg->v[v2]), e);
+
+ return;
+}
+
+struct planar_graph* pg_create(struct line_pnts *Points) {
+ struct seg_intersections *si;
+ struct planar_graph *pg;
+ struct intersection_point *ip;
+ struct pg_vertex *vert;
+ struct pg_edge *edge;
+ int i, j, t, v;
+
+ G_debug(4, "pg_create()");
+
+ si = find_all_intersections(Points);
+ pg = pg_create_struct(si->group_count, 2*(si->ipcount));
+
+ for (i = 0; i < si->ipcount; i++) {
+ ip = &(si->ip[i]);
+ t = ip->group;
+ pg->v[t].x = ip->x;
+ pg->v[t].y = ip->y;
+ }
+
+ for (i = 0; i < si->ilcount; i++) {
+ v = si->ip[si->il[i].a[0].ip].group;
+ for (j = 1; j < si->il[i].count; j++) {
+ t = si->ip[si->il[i].a[j].ip].group;
+ if (t != v) {
+ pg_addedge(pg, t, v);
+ v = t;
+ }
+ }
+ }
+
+ for (i = 0; i < pg->vcount; i++) {
+ vert = &(pg->v[i]);
+ vert->angles = G_malloc((vert->ecount)*sizeof(double));
+ for (j = 0; j < vert->ecount; j++) {
+ edge = vert->edges[j];
+ t = (edge->v1 != i)?(edge->v1):(edge->v2);
+ vert->angles[j] = atan2(pg->v[t].y - vert->y, pg->v[t].x - vert->x);
+ }
+ }
+
+ destroy_si_struct(si);
+ /*
+ I'm not sure if shrinking of the allocated memory always preserves it's physical place.
+ That's why I don't want to do this:
+ if (pg->ecount < pg->eallocated) {
+ pg->eallocated = pg->ecount;
+ pg->e = G_realloc(pg->e, (pg->ecount)*sizeof(struct pg_edge));
+ }
+ */
+
+ /* very time consuming */
+ /*
+ for (i = 0; i < pg->vcount; i++) {
+ if (pg->v[i].ecount < pg->v[i].eallocated) {
+ pg->v[i].eallocated = pg->v[i].ecount;
+ pg->v[i].edges = G_realloc(pg->v[i].edges, (pg->v[i].ecount)*sizeof(struct pg_edges));
+ }
+ }
+ */
+
+ /* output pg */
+ for (i = 0; i < pg->vcount; i++) {
+ G_debug(4, " vertex %d (%g, %g)", i, pg->v[i].x, pg->v[i].y);
+ for (j = 0; j < pg->v[i].ecount; j++) {
+ G_debug(4, " edge %d-%d", pg->v[i].edges[j]->v1, pg->v[i].edges[j]->v2);
+ }
+ }
+
+ return pg;
+}
Added: grass-addons/vector/v.parallel2/dgraph.h
===================================================================
--- grass-addons/vector/v.parallel2/dgraph.h (rev 0)
+++ grass-addons/vector/v.parallel2/dgraph.h 2008-07-23 15:15:36 UTC (rev 32228)
@@ -0,0 +1,35 @@
+#ifndef DGRAPH_H
+#define DGRAPH_H
+
+/* pg comes from "planar graph" */
+struct pg_edge {
+ int v1; /* v1 should be smaller than v2 */
+ int v2;
+ char visited_left;
+ char visited_right;
+};
+
+struct pg_vertex {
+ double x; /* coordinates */
+ double y;
+ int ecount; /* number of neighbours */
+ int eallocated; /* size of the array bellow */
+ struct pg_edge **edges; /* array of pointers */
+ double *angles; /* precalculated angles with Ox */
+};
+
+struct planar_graph {
+ int vcount; /* number of vertices */
+ struct pg_vertex *v;
+ int ecount;
+ int eallocated;
+ struct pg_edge *e;
+};
+
+struct planar_graph* pg_create_struct(int n, int e);
+void pg_destroy_struct(struct planar_graph *pg);
+int pg_existsedge(struct planar_graph *pg, int v1, int v2);
+void pg_addedge(struct planar_graph *pg, int v1, int v2);
+struct planar_graph* pg_create(struct line_pnts *Points);
+
+#endif
Modified: grass-addons/vector/v.parallel2/main.c
===================================================================
--- grass-addons/vector/v.parallel2/main.c 2008-07-23 14:55:18 UTC (rev 32227)
+++ grass-addons/vector/v.parallel2/main.c 2008-07-23 15:15:36 UTC (rev 32228)
@@ -6,15 +6,13 @@
*
* PURPOSE: Create parallel lines
*
- * COPYRIGHT: (C) 2001 by the GRASS Development Team
+ * COPYRIGHT: (C) 2008 by the GRASS Development Team
*
* This program is free software under the
* GNU General Public License (>=v2).
* Read the file COPYING that comes with GRASS
* for details.
*
- * TODO/BUG: The vector segments are not properly connected
- * and should be snapped.
**************************************************************/
#include <stdlib.h>
#include <math.h>
@@ -139,7 +137,7 @@
for ( line = 1; line <= nlines; line++ ) {
int ltype;
- G_percent ( line, nlines, 1 );
+ G_percent(line, nlines, 1);
ltype = Vect_read_line ( &In, Points, Cats, line);
@@ -168,8 +166,14 @@
}
}
+/* Vect_build_partial(&Out, GV_BUILD_BASE, stderr);
+ Vect_snap_lines(&Out, GV_BOUNDARY, 1e-7, NULL, stderr);
+ Vect_break_lines(&Out, GV_BOUNDARY, NULL, stderr);
+ Vect_remove_duplicates(&Out, GV_BOUNDARY, NULL, stderr);
+ Vect_build_partial (&Out, GV_BUILD_NONE, NULL); */
+ Vect_build (&Out, stderr);
+
Vect_close (&In);
- Vect_build (&Out, stderr);
Vect_close (&Out);
exit(EXIT_SUCCESS) ;
Modified: grass-addons/vector/v.parallel2/vlib_buffer.c
===================================================================
--- grass-addons/vector/v.parallel2/vlib_buffer.c 2008-07-23 14:55:18 UTC (rev 32227)
+++ grass-addons/vector/v.parallel2/vlib_buffer.c 2008-07-23 15:15:36 UTC (rev 32228)
@@ -1,6 +1,6 @@
/* Functions: ...
**
-** Author: Radim Blazek, Rosen Matev; June 2008
+** Author: Radim Blazek, Rosen Matev; July 2008
**
**
*/
@@ -8,6 +8,7 @@
#include <math.h>
#include <grass/Vect.h>
#include <grass/gis.h>
+#include "dgraph.h"
#define LENGTH(DX, DY) (sqrt((DX*DX)+(DY*DY)))
#ifndef MIN
@@ -157,7 +158,7 @@
*/
void parallel_line(struct line_pnts *Points, double da, double db, double dalpha, int side, int round, int caps, int looped, double tol, struct line_pnts *nPoints)
{
- int i, i_minus_1, j, res, np;
+ int i, j, res, np;
double *x, *y;
double tx, ty, vx, vy, wx, wy, nx, ny, mx, my, rx, ry;
double vx1, vy1, wx1, wy1;
@@ -165,16 +166,24 @@
double phi1, phi2, delta_phi;
double nsegments, angular_tol, angular_step;
double cosa, sina, r;
- int inner_corner, turns180, loop_flag;
+ int inner_corner, turns180;
G_debug(4, "parallel_line()");
- Vect_reset_line(nPoints);
+ if (looped && 0) {
+ /* start point != end point */
+ return;
+ }
+ Vect_reset_line(nPoints);
+
+ if (looped) {
+ Vect_append_point(Points, Points->x[1], Points->y[1], Points->z[1]);
+ }
np = Points->n_points;
x = Points->x;
y = Points->y;
-
+
if ((np == 0) || (np == 1))
return;
@@ -186,22 +195,17 @@
side = (side >= 0)?(1):(-1); /* normalize variable */
dalpha *= PI/180; /* convert dalpha from degrees to radians */
angular_tol = angular_tolerance(tol, da, db);
- loop_flag = 1;
- for (i = 0; (i < np)&&(loop_flag); i++)
+ for (i = 0; i < np-1; i++)
{
- i_minus_1 = i - 1;
- if (i == np-1) {
- if (!looped) {
- break;
- }
- else {
- i = 0;
- i_minus_1 = np - 2;
- loop_flag = 0;
- }
- }
+ /* save the old values */
+ a0 = a1;
+ b0 = b1;
+ c0 = c1;
+ wx = vx;
+ wy = vy;
+
norm_vector(x[i], y[i], x[i+1], y[i+1], &tx, &ty);
elliptic_tangent(side*tx, side*ty, da, db, dalpha, &vx, &vy);
@@ -212,286 +216,185 @@
my = y[i+1] + vy;
line_coefficients(nx, ny, mx, my, &a1, &b1, &c1);
-
- if ((i == 0) && loop_flag) {
- if (!looped) {
+
+ if (i == 0) {
+ if (!looped)
Vect_append_point(nPoints, nx, ny, 0);
- G_debug(4, "append point x=%.18f y=%.18f i=%d", nx, ny, i);
+ continue;
+ }
+
+ delta_phi = atan2(ty, tx) - atan2(y[i]-y[i-1], x[i]-x[i-1]);
+ if (delta_phi > PI)
+ delta_phi -= 2*PI;
+ else if (delta_phi <= -PI)
+ delta_phi += 2*PI;
+ /* now delta_phi is in [-pi;pi] */
+ turns180 = (fabs(fabs(delta_phi) - PI) < 1e-15);
+ inner_corner = (side*delta_phi <= 0) && (!turns180);
+
+ if ((turns180) && (!(caps && round))) {
+ if (caps) {
+ norm_vector(0, 0, vx, vy, &tx, &ty);
+ elliptic_tangent(side*tx, side*ty, da, db, dalpha, &tx, &ty);
}
+ else {
+ tx = 0;
+ ty = 0;
+ }
+ Vect_append_point(nPoints, x[i] + wx + tx, y[i] + wy + ty, 0);
+ Vect_append_point(nPoints, nx + tx, ny + ty, 0); /* nx == x[i] + vx, ny == y[i] + vy */
}
+ else if ((!round) || inner_corner) {
+ res = line_intersection(a0, b0, c0, a1, b1, c1, &rx, &ry);
+/* if (res == 0) {
+ G_debug(4, "a0=%.18f, b0=%.18f, c0=%.18f, a1=%.18f, b1=%.18f, c1=%.18f", a0, b0, c0, a1, b1, c1);
+ G_fatal_error("Two consequtive line segments are parallel, but not on one straight line! This should never happen.");
+ return;
+ } */
+ if (res == 1) {
+ Vect_append_point(nPoints, rx, ry, 0);
+ }
+ }
else {
- delta_phi = atan2(ty, tx) - atan2(y[i_minus_1+1]-y[i_minus_1], x[i_minus_1+1]-x[i_minus_1]);
- if (delta_phi > PI)
- delta_phi -= 2*PI;
- else if (delta_phi <= -PI)
+ /* we should draw elliptical arc for outside corner */
+
+ /* inverse transforms */
+ elliptic_transform(wx, wy, 1/da, 1/db, dalpha, &wx1, &wy1);
+ elliptic_transform(vx, vy, 1/da, 1/db, dalpha, &vx1, &vy1);
+
+ phi1 = atan2(wy1, wx1);
+ phi2 = atan2(vy1, vx1);
+ delta_phi = side*(phi2 - phi1);
+
+ /* make delta_phi in [0, 2pi] */
+ if (delta_phi < 0)
delta_phi += 2*PI;
- /* now delta_phi is in [-pi;pi] */
- turns180 = (fabs(fabs(delta_phi) - PI) < 1e-15);
- inner_corner = (side*delta_phi <= 0) && (!turns180);
- if ((turns180) && (!(caps && round))) {
- if (caps) {
- norm_vector(0, 0, vx, vy, &tx, &ty);
- elliptic_tangent(side*tx, side*ty, da, db, dalpha, &tx, &ty);
- }
- else {
- tx = 0;
- ty = 0;
- }
- Vect_append_point(nPoints, x[i] + wx + tx, y[i] + wy + ty, 0);
- Vect_append_point(nPoints, nx + tx, ny + ty, 0); /* nx == x[i] + vx, ny == y[i] + vy */
+ nsegments = (int)(delta_phi/angular_tol) + 1;
+ angular_step = side*(delta_phi/nsegments);
+
+ for (j = 0; j <= nsegments; j++) {
+ elliptic_transform(cos(phi1), sin(phi1), da, db, dalpha, &tx, &ty);
+ Vect_append_point(nPoints, x[i] + tx, y[i] + ty, 0);
+ phi1 += angular_step;
}
- else if ((!round) || inner_corner) {
- res = line_intersection(a0, b0, c0, a1, b1, c1, &rx, &ry);
-/* if (res == 0) {
- G_debug(4, "a0=%.18f, b0=%.18f, c0=%.18f, a1=%.18f, b1=%.18f, c1=%.18f", a0, b0, c0, a1, b1, c1);
- G_fatal_error("Two consequtive line segments are parallel, but not on one straight line! This should never happen.");
- return;
- } */
- if (res == 1) {
- Vect_append_point(nPoints, rx, ry, 0);
- }
- }
- else {
- /* we should draw elliptical arc for outside corner */
-
- /* inverse transforms */
- elliptic_transform(wx, wy, 1/da, 1/db, dalpha, &wx1, &wy1);
- elliptic_transform(vx, vy, 1/da, 1/db, dalpha, &vx1, &vy1);
-
- phi1 = atan2(wy1, wx1);
- phi2 = atan2(vy1, vx1);
- delta_phi = side*(phi2 - phi1);
-
- /* make delta_phi in [0, 2pi] */
- if (delta_phi < 0)
- delta_phi += 2*PI;
-
- nsegments = (int)(delta_phi/angular_tol) + 1;
- angular_step = side*(delta_phi/nsegments);
-
- for (j = 0; j <= nsegments; j++) {
- elliptic_transform(cos(phi1), sin(phi1), da, db, dalpha, &tx, &ty);
- Vect_append_point(nPoints, x[i] + tx, y[i] + ty, 0);
- phi1 += angular_step;
- }
- }
}
- if ((i == np-2) && (!looped)) {
+ if ((!looped) && (i == np-2)) {
Vect_append_point(nPoints, mx, my, 0);
}
-
- /* save the old values */
- a0 = a1;
- b0 = b1;
- c0 = c1;
- wx = vx;
- wy = vy;
}
if (looped) {
- i = nPoints->n_points - 1;
- Vect_line_insert_point(nPoints, 0, nPoints->x[i], nPoints->y[i], 0);
+ Vect_append_point(nPoints, nPoints->x[0], nPoints->y[0], nPoints->z[0]);
}
+
Vect_line_prune ( nPoints );
-}
-
-static void split_at_intersections(struct line_pnts *Points, struct line_pnts *nPoints)
-{
- #define MAX_SEGMENT_INTERSECTIONS 10 /* we should make dynamic array here? */
- #define EPSILON 1e-10 /* */
- int i, j, k, res, np;
- double *x, *y; /* input coordinates */
- double x1, y1, z1, x2, y2, z2; /* intersection points */
- double len; /* current i-th line length */
- double tdist, min;
- double px[MAX_SEGMENT_INTERSECTIONS], py[MAX_SEGMENT_INTERSECTIONS], pdist[MAX_SEGMENT_INTERSECTIONS];
- int pcount, min_i;
-
- G_debug(1, "split_at_intersections()");
- Vect_reset_line(nPoints);
-
- Vect_line_prune(Points);
- np = Points->n_points;
- x = Points->x;
- y = Points->y;
-
- if ((np == 0) || (np == 1))
- return;
-
- for (i = 0; i < np-1; i++) {
- Vect_append_point(nPoints, x[i], y[i], 0);
- G_debug(1, "appended point x=%f y=%f", x[i], y[i]);
-
-
- /* find all intersections */
- pcount = 0;
- len = LENGTH((x[i]-x[i+1]), (y[i]-y[i+1]));
- for (j = 0; j < np-1; j++) {
- /* exclude previous (i-1)th, current i-th, and next (i+1)th line */
- if ((j < i-1) || (j > i+1)) {
- res = Vect_segment_intersection(x[i], y[i], 0, x[i+1], y[i+1], 0,
- x[j], y[j], 0, x[j+1], y[j+1], 0,
- &x1, &y1, &z1, &x2, &y2, &z2, 0);
- G_debug(1, "intersection=%d", res);
- if (res == 1) {
- /* generall intersection */
- tdist = LENGTH((x[i]-x1), (y[i]-y1));
- /* we only want intersections on the inside */
- G_debug(1, "x1=%f y1=%f tdist=%f len=%f, eps=%f", x1, y1, tdist, len, EPSILON);
- if ((tdist > EPSILON) && (tdist < len - EPSILON)) {
- px[pcount] = x1;
- py[pcount] = y1;
- pdist[pcount] = tdist;
- pcount++;
- }
- }
- /* TODO: implement cases of overlappings */
- }
- }
-
- /* now we shoud output interserction points ordered by distance to (x[i],y[i]) */
- /* we do it in O(pcount^2) time, since pcount is usually small it's okey */
- G_debug(1, "pcount=%d", pcount);
- for (j = 0; j < pcount; j++) {
- min = 2*len;
- for (k = 0; k < pcount; k++) {
- if (pdist[k] < min) {
- min = pdist[k];
- min_i = k;
- }
- }
- pdist[min_i] = 4*len;
- Vect_append_point(nPoints, px[min_i], py[min_i], 0);
- G_debug(1, "appended point x=%f y=%f", px[min_i], py[min_i]);
-
- }
+ if (looped) {
+ Vect_line_delete_point(Points, Points->n_points-1);
}
- Vect_append_point(nPoints, x[np-1], y[np-1], 0);
- G_debug(1, "appended point x=%f y=%f", x[i], y[i]);
- Vect_line_prune(nPoints);
}
/*
-* IMPORTANT: split_at_intersections() must be applied to input line before calling extract_contour().
-* When visited != NULL, extract_contour() marks the sides of lines as visited along the contour.
-* visited->a must have at least Points->n_points elements
-*
-* side: side >= 0 - right contour, side < 0 - left contour
-* returns: -1 when contour is a loop; 0 or Points->n_points-1 when contour finishes at line end (depending on which one)
+* side: side >= 0 - extracts contour on right side of edge, side < 0 - extracts contour on left side of edge
+* returns: the side of the output line, where we should draw parallel line
*/
-static int extract_contour(struct line_pnts *Points, int first_point, int first_step, int side, int stop_at_line_end, struct line_pnts *nPoints, struct visited_segments *visited) {
- int i, is, j, k, step, np;
- int ret;
- double *x, *y;
- double cx, cy, tx, ty;
- double u, v, len, new_len, new_step;
- int intersection_flag;
- double x1, y1, z1, x2, y2, z2;
+static int extract_contour(struct planar_graph *pg, struct pg_edge *first, int side, int stop_at_line_end, struct line_pnts *nPoints) {
+ int i, j, ret;
+ int v; /* current vertex number */
+ int v0;
+ int eside; /* side of the current edge */
+ double eangle; /* current edge angle with Ox (according to the current direction) */
+ struct pg_vertex *vert; /* current vertex */
+ struct pg_vertex *vert0; /* last vertex */
+ struct pg_edge *edge; /* current edge; must be edge of vert */
+ int cs; /* on which side are we turning along the contour */
double opt_angle, tangle;
- int opt_j, opt_step, opt_flag;
+ int opt_j, opt_side, opt_flag;
- G_debug(4, "extract_contour(): first_point=%d, first_step=%d, side=%d, stop_at_line_end=%d", first_point, first_step, side, stop_at_line_end);
+ G_debug(4, "extract_contour(): v1=%d, v2=%d, side=%d, stop_at_line_end=%d", first->v1, first->v2, side, stop_at_line_end);
Vect_reset_line(nPoints);
-
- np = Points->n_points;
- x = Points->x;
- y = Points->y;
- G_debug(4, "ec: Points->n_points=%d", np);
- if ((np == 0) || (np == 1))
- return;
-
- /* normalize parameter side */
- if (side >= 0)
- side = 1;
- else if (side < 0)
- side = -1;
-
- i = first_point;
- step = first_step;
+ edge = first;
+ eside = (side >= 0)?(1):(-1);
+ v = edge->v1; /* we might select v1 or v2 if we always want ccw output line */
+ v0 = edge->v2;
+ vert = &(pg->v[v]);
+ vert0 = &(pg->v[v0]);
+ /*for (i = 0; i < pg->v[v0].ecount; i++)
+ if (pg->v[v0].edges[i] == edge) {
+ eangle = pg->v[v0].angles[i];
+ break;
+ }
+ */
+ eangle = atan2(vert->y - vert0->y, vert->x - vert0->x);
+ cs = -eside; /* because we chose to go v2-->v1 */
while (1) {
- /* precessing segment AB: A=(x[i],y[i]) B=(x[i+step],y[i+step]) */
- Vect_append_point(nPoints, x[i], y[i], 0);
- G_debug(4, "ec: append point x=%.18f y=%.18f", x[i], y[i]);
- is = i + step;
+ 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: first->v1=%d, first->v2=%d", first->v1, first->v2);
+ G_debug(4, "ec: append point x=%.18f y=%.18f", vert0->x, vert0->y);
+
+ /* mark current edge as visited on the appropriate side */
+ if (eside == 1)
+ edge->visited_right = 1;
+ else
+ edge->visited_left = 1;
- if (visited != NULL) {
- j = (3-step*side)/2; /* j is 1 when (step*side == 1), and j is 2 when (step*side == -1) */
- k = MIN(i, is);
- if ((visited->a[k] & j) == 0)
- visited->n++;
- visited->a[k] |= j;
- }
-
- cx = x[i] - x[is];
- cy = y[i] - y[is];
opt_flag = 1;
- for (j = 0; j < np-1; j++) {
- /* exclude current segment AB */
- if (j != i - (1-step)/2) {
- if ((fabs(x[j]-x[is]) < EPSILON) && (fabs(y[j]-y[is]) < EPSILON)) {
- tx = x[j+1] - x[j];
- ty = y[j+1] - y[j];
- tangle = atan2(ty, tx) - atan2(cy, cx);
- if (tangle < 0)
- tangle += 2*PI;
- /* now tangle is in [0, 2PI) */
-
- if (opt_flag || (side*tangle < side*opt_angle)) {
- opt_j = j;
- opt_step = 1;
- opt_angle = tangle;
- opt_flag = 0;
- }
+ for (j = 0; j < vert->ecount; j++) {
+ /* exclude current edge */
+ if (vert->edges[j] != edge) {
+ tangle = vert->angles[j] - eangle;
+ if (tangle < -PI)
+ tangle += 2*PI;
+ else if (tangle > PI)
+ tangle -= 2*PI;
+ /* now tangle is in (-PI, PI) */
+
+ if (opt_flag || (cs*tangle < cs*opt_angle)) {
+ opt_j = j;
+ opt_side = (vert->edges[j]->v1 == v)?(cs):(-cs);
+ opt_angle = tangle;
+ opt_flag = 0;
}
- else if ((fabs(x[j+1]-x[is]) < EPSILON) && (fabs(y[j+1]-y[is]) < EPSILON)) {
- tx = x[j] - x[j+1];
- ty = y[j] - y[j+1];
- tangle = atan2(ty, tx) - atan2(cy, cx);
- if (tangle < 0)
- tangle += 2*PI;
- /* now tangle is in [0, 2PI) */
-
- if (opt_flag || (side*tangle < side*opt_angle)) {
- opt_j = j+1;
- opt_step = -1;
- opt_angle = tangle;
- opt_flag = 0;
- }
-
- }
}
}
- /* if line end is reached */
+// G_debug(4, "ec: opt: side=%d opt_flag=%d opt_angle=%.18f opt_j=%d opt_step=%d", side, opt_flag, opt_angle, opt_j, opt_step);
+
+ /* if line end is reached (no other edges at curr vertex) */
if (opt_flag) {
if (stop_at_line_end) {
- Vect_append_point(nPoints, x[is], y[is], 0);
- ret = is;
+ ret = 2;
break;
}
else {
- opt_j = is;
- opt_step = -step;
+ opt_j = 0; /* the only edge of vert is vert->edges[0] */
+ opt_side = -eside; /* go to the other side of the curr edge */
}
}
- if ((opt_j == first_point) && (opt_step == first_step)) {
- Vect_append_point(nPoints, x[is], y[is], 0);
- ret = -1;
+ if ((vert->edges[opt_j] == first) && (opt_side == side)) {
+ ret = 1;
break;
}
-
- i = opt_j;
- step = opt_step;
- }
- Vect_line_prune(nPoints);
- return ret;
+ edge = vert->edges[opt_j];
+ eside = opt_side;
+ v0 = v;
+ v = (edge->v1 == v)?(edge->v2):(edge->v1);
+ vert0 = vert;
+ vert = &(pg->v[v]);
+ eangle = vert0->angles[opt_j];
+ }
+ Vect_append_point(nPoints, vert->x, vert->y, 0);
+ G_debug(4, "ec: append point x=%.18f y=%.18f", vert->x, vert->y);
+
+ return cs;
}
/*
@@ -502,128 +405,91 @@
* side: side > 0 - right contour, side < 0 - left contour, side = 0 - outer contour
* if side != 0 and there's only one contour, the function returns it
*
-* IMPORTANT: split_at_intersections() must be applied to input line before calling extract_outer_contour
* TODO: Implement side != 0 feature;
-* returns: returns on which side of output line is the exterior of the input line
+* retunrs: the side of the output line, where we should draw parallel line
*/
-static int extract_outer_contour(struct line_pnts *Points, int side, struct line_pnts *nPoints, struct visited_segments *visited) {
- int i, np, res, ret;
- double *x, *y;
+int extract_outer_contour(struct planar_graph *pg, int side, struct line_pnts *nPoints) {
+ int i;
+ int flag;
+ int v;
+ struct pg_vertex *vert;
+ struct pg_edge *edge;
double min_x, min_angle, ta;
- int t1, t2, j1=-1, j2=-1;
G_debug(4, "extract_outer_contour()");
- np = Points->n_points;
- x = Points->x;
- y = Points->y;
-
if (side != 0) {
- res = extract_contour(Points, 0, 1, side, 1, nPoints, NULL); /* should save visited and later restore it if needed */
- /* if we have finished at the other end */
- if (res == Points->n_points-1) {
- /* !!! WE ARE STILL NOT SURE WHETHER THIS IS PROPER OUTER CONTOUR.
- * !!! CONSIDER THE CASE WHERE BOTH LINE ENDS ARE IN THE SAME LOOP.
- */
- return side;
- }
+ G_fatal_error(" side != 0 feature not implemented");
+ return;
}
/* find a line segment which is on the outer contour */
- min_x = 1e300;
- for (i = 0; i < np-1; i++) {
- if (x[i] <= x[i+1]) {
- t1 = i;
- t2 = i+1;
+ flag = 1;
+ for (i = 0; i < pg->vcount; i++) {
+ if (flag || (pg->v[i].x < min_x)) {
+ v = i;
+ min_x = pg->v[i].x;
+ flag = 0;
}
- else {
- t1 = i+1;
- t2 = i;
+ }
+ vert = &(pg->v[v]);
+
+ flag = 1;
+ for (i = 0; i < vert->ecount; i++) {
+ if (flag || (vert->angles[i] < min_angle)) {
+ edge = vert->edges[i];
+ min_angle = vert->angles[i];
+ flag = 0;
}
- G_debug(4, "t1=%d (%f,%f) t2=%d (%f,%f)", t1, x[t1], y[t1], t2, x[t2], y[t2]);
- if (x[t1] < min_x) {
- min_x = x[t1];
- min_angle = atan2(y[t2]-y[t1], x[t2]-x[t1]);
- j1 = t1;
- j2 = t2;
- G_debug(4, "j1=%d j2=%d", j1, j2);
- }
- else if ((x[t1] == min_x) && (y[t1] == y[j1])) {
- ta = atan2(y[t2]-y[t1], x[t2]-x[t1]);
- if (ta < min_angle) {
- min_angle = ta;
- j1 = t1;
- j2 = t2;
- G_debug(4, "j1=%d j2=%d", j1, j2);
- }
- }
}
- G_debug(4, "j1=%d, j2-j1=%d, side=%d", j1, j2-j1, RIGHT_SIDE);
- res = extract_contour(Points, j1, j2-j1, RIGHT_SIDE, 0, nPoints, visited);
-
- return RIGHT_SIDE;
+ return extract_contour(pg, edge, (edge->v1 == v)?RIGHT_SIDE:LEFT_SIDE, 0, nPoints);
}
/*
-* Extracts contours which are not visited. Generates counterclockwise closed lines.
-* IMPORTANT: split_at_intersections() must be applied to input line before calling extract_inner_contour
-* IMPORTANT: the outer contour must be visited, so that extract_inner_contour() doesn't return it
-* returns: 0 when there are no more inner contours; otherwise, returns on which side of output line is the interior of the loop
+* Extracts contours which are not visited.
+* IMPORTANT: the outer contour must be visited (you should call extract_outer_contour() to do that),
+* so that extract_inner_contour() doesn't return it
+* returns: 0 when there are no more inner contours; otherwise, the side of the output line, where we should draw parallel line
*/
-int extract_inner_contour(struct line_pnts *Points, struct line_pnts *nPoints, struct visited_segments *visited) {
- int i, j, np, res;
- double *x, *y;
+int extract_inner_contour(struct planar_graph *pg, struct line_pnts *nPoints) {
+ int i;
G_debug(4, "extract_inner_contour()");
- np = Points->n_points;
- x = Points->x;
- y = Points->y;
-
- for (i = 0; i < np-1; i++) {
- /* if right side is not visited */
- if ((visited->a[i] & 1) == 0) {
- res = extract_contour(Points, i+1, -1, LEFT_SIDE, 0, nPoints, visited);
- return LEFT_SIDE;
- }
- /* if left side is not visited */
- if ((visited->a[i] & 2) == 0) {
- res = extract_contour(Points, i, 1, LEFT_SIDE, 0, nPoints, visited);
- return LEFT_SIDE;
- }
+ for (i = 0; i < pg->ecount; i++) {
+ if (!(pg->e[i].visited_right))
+ return extract_contour(pg, &(pg->e[i]), RIGHT_SIDE, 0, nPoints);
+
+ if (!(pg->e[i].visited_left))
+ return extract_contour(pg, &(pg->e[i]), LEFT_SIDE, 0, nPoints);
}
return 0;
}
void extract_all_inner_contours(struct line_pnts *Points, struct line_pnts ***iPoints, int *inner_count) {
- struct line_pnts *Points2, *tPoints;
+ struct planar_graph *pg;
+ struct line_pnts *tPoints;
struct line_pnts **arrPoints;
- int i, side, count;
+ int i, res, count;
int more = 8;
int allocated = more;
- struct visited_segments visited;
G_debug(4, "extract_all_inner_contours()");
- /* initializations and spliting of the input line */
+ /* initializations */
tPoints = Vect_new_line_struct();
- Points2 = Vect_new_line_struct();
- split_at_intersections(Points, Points2);
-
- visited.a = (char *)G_malloc(sizeof(char)*Points2->n_points);
- memset(visited.a, 0, sizeof(char)*Points2->n_points);
-
arrPoints = G_malloc(allocated*sizeof(struct line_pnts *));
+ pg = pg_create(Points);
+
+ /* extract outer contour so that we have visited edges appropriately */
+ extract_outer_contour(pg, 0, tPoints);
- /* extract outer contour so that we have filled visited appropriately */
- side = extract_outer_contour(Points2, 0, tPoints, &visited);
-
/* inner contours */
count = 0;
- side = extract_inner_contour(Points2, tPoints, &visited);
- while (side != 0) {
+ res = extract_inner_contour(pg, tPoints);
+ while (res != 0) {
if (allocated < count+1) {
allocated += more;
arrPoints = G_realloc(arrPoints, allocated*sizeof(struct line_pnts *));
@@ -632,26 +498,19 @@
Vect_copy_xyz_to_pnts(arrPoints[count], tPoints->x, tPoints->y, tPoints->z, tPoints->n_points);
count++;
- side = extract_inner_contour(Points2, tPoints, &visited);
+ res = extract_inner_contour(pg, tPoints);
}
arrPoints = G_realloc(arrPoints, count*sizeof(struct line_pnts *));
*inner_count = count;
*iPoints = arrPoints;
- G_free(visited.a);
+ pg_destroy_struct(pg);
Vect_destroy_line_struct(tPoints);
- Vect_destroy_line_struct(Points2);
return;
}
-double perp_dist(double x1, double y1, double x2, double y2, double px, double py) {
- double a,b,c;
- line_coefficients(x1, y1, x2, y2, &a, &b, &c);
- return fabs((a*px+b*py+c)/sqrt(a*a+b*b));
-}
-
/* point_in_buf - test if point px,py is in d buffer of Points
** returns: 1 in buffer
** 0 not in buffer
@@ -699,8 +558,8 @@
d = dig_distance2_point_to_line(nx + tx, ny + ty, 0, vx, vy, 0, wx, wy, 0,
0, NULL, NULL, NULL, NULL, NULL);
-/* G_debug(4, "pdist = %g", perp_dist(vx, vy ,wx, wy, px, py));
- G_debug(4, "sqrt(d)*da = %g, len' = %g, olen = %g", sqrt(d)*da, da*LENGTH(tx,ty), LENGTH((px-nx),(py-ny)));*/
+
+/* G_debug(4, "sqrt(d)*da = %g, len' = %g, olen = %g", sqrt(d)*da, da*LENGTH(tx,ty), LENGTH((px-nx),(py-ny)));*/
if (d <= 1)
return 1;
}
@@ -717,48 +576,47 @@
}
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 line_pnts *Points2, *tPoints, *sPoints;
+ struct planar_graph *pg, *pg2;
+ struct line_pnts *tPoints, *sPoints;
struct line_pnts **arrPoints;
struct line_pnts **arrPoints2;
- int i, side, count, count2;
- int res;
+ int i, count, count2;
+ int side, res;
int more = 8;
- int allocated = more;
+ int allocated = 0;
double px, py;
- struct visited_segments visited;
G_debug(4, "parallel_line_b()");
- /* initializations and spliting of the input line */
+ /* initializations */
tPoints = Vect_new_line_struct();
sPoints = Vect_new_line_struct();
- Points2 = Vect_new_line_struct();
- split_at_intersections(Points, Points2);
+ arrPoints = NULL;
+ pg = pg_create(Points);
- visited.a = (char *)G_malloc(sizeof(char)*Points2->n_points);
- memset(visited.a, 0, sizeof(char)*Points2->n_points);
-
- arrPoints = G_malloc(allocated*sizeof(struct line_pnts *));
-
/* outer contour */
- side = extract_outer_contour(Points2, 0, tPoints, &visited);
*oPoints = Vect_new_line_struct();
+ side = extract_outer_contour(pg, 0, tPoints);
+/* parallel_line(tPoints, da, db, dalpha, side, round, caps, LOOPED_LINE, tol, *oPoints);*/
parallel_line(tPoints, da, db, dalpha, side, round, caps, LOOPED_LINE, tol, sPoints);
- split_at_intersections(sPoints, tPoints);
- extract_outer_contour(tPoints, 0, *oPoints, NULL);
+ pg2 = pg_create(sPoints);
+ extract_outer_contour(pg2, 0, *oPoints);
+ pg_destroy_struct(pg2);
/* inner contours */
count = 0;
- side = extract_inner_contour(Points2, tPoints, &visited);
+ side = extract_inner_contour(pg, tPoints);
while (side != 0) {
parallel_line(tPoints, da, db, dalpha, side, round, caps, LOOPED_LINE, tol, sPoints);
-
extract_all_inner_contours(sPoints, &arrPoints2, &count2);
-
for (i = 0; i < count2; i++) {
res = Vect_line_check_intersection(tPoints, arrPoints2[i], 0);
- if (res != 0)
+ if (res != 0)
continue;
+
+ res = Vect_point_in_poly(arrPoints2[i]->x[0], arrPoints2[i]->y[0], tPoints);
+ if (res == 0)
+ continue;
res = Vect_get_point_in_poly(arrPoints2[i], &px, &py);
if (res != 0)
@@ -766,8 +624,7 @@
if (point_in_buf(tPoints, px, py, da, db, dalpha))
continue;
- /* passed all tests, add new island */
- if (allocated < count+1) {
+ if (allocated == count) {
allocated += more;
arrPoints = G_realloc(arrPoints, allocated*sizeof(struct line_pnts *));
}
@@ -776,17 +633,16 @@
count++;
}
- side = extract_inner_contour(Points2, tPoints, &visited);
+ side = extract_inner_contour(pg, tPoints);
}
-
+
arrPoints = G_realloc(arrPoints, count*sizeof(struct line_pnts *));
*inner_count = count;
*iPoints = arrPoints;
- G_free(visited.a);
Vect_destroy_line_struct(tPoints);
Vect_destroy_line_struct(sPoints);
- Vect_destroy_line_struct(Points2);
+ pg_destroy_struct(pg);
return;
}
@@ -922,6 +778,3 @@
}
Vect_line_prune ( OutPoints );
}
-
-
-
More information about the grass-commit
mailing list