[GRASS-SVN] r34511 - in grass/branches/develbranch_6/vector: . v.buffer2

svn_grass at osgeo.org svn_grass at osgeo.org
Wed Nov 26 16:23:44 EST 2008


Author: neteler
Date: 2008-11-26 16:23:44 -0500 (Wed, 26 Nov 2008)
New Revision: 34511

Added:
   grass/branches/develbranch_6/vector/v.buffer2/
   grass/branches/develbranch_6/vector/v.buffer2/Makefile
   grass/branches/develbranch_6/vector/v.buffer2/description.html
   grass/branches/develbranch_6/vector/v.buffer2/dgraph.c
   grass/branches/develbranch_6/vector/v.buffer2/dgraph.h
   grass/branches/develbranch_6/vector/v.buffer2/e_intersect.c
   grass/branches/develbranch_6/vector/v.buffer2/e_intersect.h
   grass/branches/develbranch_6/vector/v.buffer2/local_proto.h
   grass/branches/develbranch_6/vector/v.buffer2/main.c
   grass/branches/develbranch_6/vector/v.buffer2/vlib_buffer.c
Log:
moved here from GRASS Addons

Added: grass/branches/develbranch_6/vector/v.buffer2/Makefile
===================================================================
--- grass/branches/develbranch_6/vector/v.buffer2/Makefile	                        (rev 0)
+++ grass/branches/develbranch_6/vector/v.buffer2/Makefile	2008-11-26 21:23:44 UTC (rev 34511)
@@ -0,0 +1,13 @@
+MODULE_TOPDIR = ../..
+
+PGM = v.buffer
+
+LIBES = $(VECTLIB) $(DBMILIB) $(GISLIB)
+DEPENDENCIES= $(VECTDEP) $(DBMIDEP) $(GISDEP)
+EXTRA_INC = $(VECT_INC)
+EXTRA_CFLAGS = $(VECT_CFLAGS)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
+


Property changes on: grass/branches/develbranch_6/vector/v.buffer2/Makefile
___________________________________________________________________
Name: svn:eol-style
   + native

Added: grass/branches/develbranch_6/vector/v.buffer2/description.html
===================================================================
--- grass/branches/develbranch_6/vector/v.buffer2/description.html	                        (rev 0)
+++ grass/branches/develbranch_6/vector/v.buffer2/description.html	2008-11-26 21:23:44 UTC (rev 34511)
@@ -0,0 +1,58 @@
+<h2>DESCRIPTION</h2>
+
+<em>v.buffer</em> creates a buffer around features of given <b>type</b>, which
+have a category in the given <b>layer</b>. The <b>tolerance</b> controls
+the number of vector segments being generated (the smaller the value, the more
+vector segments are generated).
+
+<h2>NOTES</h2>
+
+<!-- TRUE?? -->
+Attributes are not transferred due to potential buffer overlap, which
+cannot be resolved automatically.
+
+<h2>EXAMPLES</h2>
+
+<h3>Buffer around input lines</h3>
+
+<div class="code"><pre>
+v.buffer input=map output=buffer type=line buffer=100
+</pre></div>
+
+<h3>Circles around input points</h3>
+
+<div class="code"><pre>
+v.buffer input=pointsmap output=circles type=point buffer=1000 
+</pre></div>
+
+<h3>Non-overlapping circles around input points with attribute transfer</h3>
+
+<div class="code"><pre>
+v.buffer input=archsites output=circles type=point buffer=200 
+# change original points to centroids: 
+v.type in=archsites out=archcentroids type=point,centroid 
+# patch circles and centroids: 
+v.patch in=archcentroids,circles out=circles_db 
+# attach attributes, either use 
+# db.copy ... 
+# or link to the original table: 
+v.db.connect map=circles_db table=archsites field=1 key=cat driver=dbf \
+database='$GISDBASE/$LOCATION_NAME/$MAPSET/dbf'
+</pre></div>
+
+<h2>SEE ALSO</h2>
+
+<em>
+<a HREF="r.buffer.html">r.buffer</a>,
+<a HREF="v.extract.html">v.extract</a>,
+<a HREF="v.type.html">v.type</a>,
+<a HREF="v.patch.html">v.patch</a>,
+<a HREF="v.db.connect.html">v.db.connect</a>
+</em>
+
+<h2>AUTHORS</h2>
+
+Radim Blazek<br>
+Rewritten by Rosen Matev (with support through the Google Summer of Code program 2008)
+<p>
+<i>Last changed: $Date: 2008-03-16 16:29:33 +0100 (Sun, 16 Mar 2008) $</i>

Added: grass/branches/develbranch_6/vector/v.buffer2/dgraph.c
===================================================================
--- grass/branches/develbranch_6/vector/v.buffer2/dgraph.c	                        (rev 0)
+++ grass/branches/develbranch_6/vector/v.buffer2/dgraph.c	2008-11-26 21:23:44 UTC (rev 34511)
@@ -0,0 +1,454 @@
+#include <stdlib.h>
+#include <grass/Vect.h>
+#include <grass/gis.h>
+#include "dgraph.h"
+#include "e_intersect.h"
+
+#define LENGTH(DX, DY) (sqrt((DX*DX)+(DY*DY)))
+#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
+
+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);
+}        
+
+/* 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 = 0.00000001;
+    double *x, *y;
+    double x1, y1, x2, y2;
+    int res;
+    /*int res2
+    double x1_, y1_, x2_, y2_, z1_, z2_;*/
+    struct seg_intersections *si;
+    struct seg_intersection_list *il;
+    struct intersection_point **sorted;
+    
+    G_debug(3, "find_all_intersections()");
+    
+    np = Points->n_points;
+    x = Points->x;
+    y = Points->y;
+    
+    si = create_si_struct(np-1);
+    
+    looped = ((x[0] == x[np-1]) && (y[0] == y[np-1]));
+    G_debug(3, "    looped=%d", looped);
+    
+    G_debug(3, "    finding intersections...");
+    for (i = 0; i < np-1; i++) {
+        for (j = i+1; j < np-1; j++) {
+            G_debug(4, "        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);
+/*            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) || ((res != 0) && (x1!=x1_ || y1!=y1_)) ) {
+                G_debug(0, "exact=%d orig=%d", res, res2);
+                segment_intersection_2d_test(x[i], y[i], x[i+1], y[i+1], x[j], y[j], x[j+1], y[j+1], &x1, &y1, &x2, &y2);
+            }
+*/            
+            G_debug(4, "        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(3, "    finding intersections...done");
+    
+    G_debug(3, "    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))
+/*            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;
+            }
+        }
+        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(3, "    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;
+    
+    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 */
+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)) || ((e->v1 == v2) && (e->v2 == v1)))
+            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;
+    
+    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_fatal_error("    pg_addedge(), v1 and/or v2 is invalid");
+        return;
+    }
+    
+    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->winding_left = 0; /* winding is undefined if the corresponding side is not visited */
+    e->winding_right = 0;
+    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(3, "pg_create()");
+    
+    si = find_all_intersections(Points);
+    pg = pg_create_struct(si->group_count, 2*(si->ipcount));
+
+    /* set vertices info */
+    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;
+    }
+    
+    /* add all edges */
+    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, v, t); /* edge direction is v ---> t */
+                v = t;
+            }
+        }
+    }
+    
+    /* precalculate angles with 0x */
+    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;
+}


Property changes on: grass/branches/develbranch_6/vector/v.buffer2/dgraph.c
___________________________________________________________________
Name: svn:eol-style
   + native

Added: grass/branches/develbranch_6/vector/v.buffer2/dgraph.h
===================================================================
--- grass/branches/develbranch_6/vector/v.buffer2/dgraph.h	                        (rev 0)
+++ grass/branches/develbranch_6/vector/v.buffer2/dgraph.h	2008-11-26 21:23:44 UTC (rev 34511)
@@ -0,0 +1,38 @@
+#ifndef DGRAPH_H
+#define DGRAPH_H
+
+/* pg comes from "planar graph" */
+/* every edge is directed. Nevertheless, we can visit it on both sides */
+struct pg_edge {
+    int v1; /* first vertex */
+    int v2; /* second vertex */
+    char visited_left;
+    char visited_right;
+    char winding_left; /* winding numbers */
+    char winding_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


Property changes on: grass/branches/develbranch_6/vector/v.buffer2/dgraph.h
___________________________________________________________________
Name: svn:eol-style
   + native

Added: grass/branches/develbranch_6/vector/v.buffer2/e_intersect.c
===================================================================
--- grass/branches/develbranch_6/vector/v.buffer2/e_intersect.c	                        (rev 0)
+++ grass/branches/develbranch_6/vector/v.buffer2/e_intersect.c	2008-11-26 21:23:44 UTC (rev 34511)
@@ -0,0 +1,918 @@
+#include <stdlib.h>
+#include <math.h>
+#include <grass/gis.h>
+#include "e_intersect.h"
+
+#define SWAP(a,b) {t = a; a = b; b = t;}
+#ifndef MIN
+    #define MIN(X,Y) ((X<Y)?X:Y)
+#endif
+#ifndef MAX
+    #define MAX(X,Y) ((X>Y)?X:Y)
+#endif    
+#define D (ax2-ax1)*(by1-by2) - (ay2-ay1)*(bx1-bx2)
+#define DA (bx1-ax1)*(by1-by2) - (by1-ay1)*(bx1-bx2)
+#define DB (ax2-ax1)*(by1-ay1) - (ay2-ay1)*(bx1-ax1)
+
+
+#ifdef ASDASDASFDSAFFDAS
+mpf_t p11, p12, p21, p22, t1, t2;
+mpf_t dd, dda, ddb, delta;
+mpf_t rra, rrb;
+
+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);
+    
+    mpf_init(rra);
+    mpf_init(rrb);
+    
+    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 ((MAX(ax1, ax2) < MIN(bx1, bx2)) || (MAX(bx1, bx2) < MIN(ax1, ax2))) {
+        G_debug(3, "    no intersection (disjoint bounding boxes)");
+        return 0;
+    }
+    if ((MAX(ay1, ay2) < MIN(by1, by2)) || (MAX(by1, by2) < MIN(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;
+}
+#endif
+
+/* OLD */
+/* tollerance aware version */
+/* TODO: fix all ==s left */
+int segment_intersection_2d_tol(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  = %.18g", d);
+    G_debug(2, "    d1 = %.18g", d1);
+    G_debug(2, "    d2 = %.18g", d2);
+    
+    tola = tol/MAX(fabs(ax2-ax1), fabs(ay2-ay1));
+    tolb = tol/MAX(fabs(bx2-bx1), fabs(by2-by1));
+    G_debug(2, "    tol  = %.18g", tol);
+    G_debug(2, "    tola = %.18g", tola);
+    G_debug(2, "    tolb = %.18g", tolb);
+    if (!FZERO(d, tol)) {
+        ra = d1/d;
+        rb = d2/d;
+        
+        G_debug(2, "    not parallel/collinear: ra = %.18g", ra);
+        G_debug(2, "                            rb = %.18g", 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;
+}
+
+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)
+{
+    const int DLEVEL = 4;
+    double t;
+    int vertical;
+    int f11, f12, f21, f22;
+    double d, da, db;
+    
+    /* TODO: Works for points ? */
+    G_debug(DLEVEL, "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);
+
+    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(DLEVEL, "    identical segments" );
+        *x1 = ax1; *y1 = ay1;
+        *x2 = ax2; *y2 = ay2;
+        return 5;
+    }
+    /* Check for identical endpoints */
+    if (f11 || f12) {
+        G_debug(DLEVEL, "    connected by endpoints" );
+        *x1 = ax1; *y1 = ay1;
+        return 1;
+    }
+    if (f21 || f22) {
+        G_debug (DLEVEL, "    connected by endpoints" );
+        *x1 = ax2; *y1 = ay2;
+        return 1;
+    }
+
+    if ((MAX(ax1, ax2) < MIN(bx1, bx2)) || (MAX(bx1, bx2) < MIN(ax1, ax2))) {
+        G_debug(DLEVEL, "    no intersection (disjoint bounding boxes)");
+        return 0;
+    }
+    if ((MAX(ay1, ay2) < MIN(by1, by2)) || (MAX(by1, by2) < MIN(ay1, ay2))) {
+        G_debug(DLEVEL, "    no intersection (disjoint bounding boxes)");
+        return 0;
+    }
+        
+    d  = D;
+    if (d != 0) {
+        G_debug(DLEVEL, "    general position");
+
+        da = DA;
+        
+        /*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 (d > 0) {
+            if ((da < 0) || (da > d)) {
+                G_debug(DLEVEL, "        no intersection");
+                return 0;
+            }
+            
+            db = DB;
+            if ((db < 0) || (db > d)) {
+                G_debug(DLEVEL, "        no intersection");
+                return 0;
+            }
+        }
+        else { /* if d < 0 */
+            if ((da > 0) || (da < d)) {
+                G_debug(DLEVEL, "        no intersection");
+                return 0;
+            }
+            
+            db = DB;
+            if ((db > 0) || (db < d)) {
+                G_debug(DLEVEL, "        no intersection");
+                return 0;
+            }
+        }
+
+        /*G_debug(DLEVEL, "        ra=%.17g rb=%.17g", mpf_get_d(dda)/mpf_get_d(dd), mpf_get_d(ddb)/mpf_get_d(dd));*/
+        /*G_debug(DLEVEL, "        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));*/
+    
+        *x1 = ax1 + (ax2 - ax1)*da/d;
+        *y1 = ay1 + (ay2 - ay1)*da/d;
+        
+        G_debug(DLEVEL, "        intersection %.16g, %.16g", *x1, *y1);
+        return 1;
+    }
+    
+    /* segments are parallel or collinear */
+    da = DA;
+    db = DB;
+    if ((da != 0) || (db != 0)) {
+        /* segments are parallel */
+        G_debug(DLEVEL, "    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(DLEVEL, "    collinear segments");
+
+    if ((bx2 < ax1) || (bx1 > ax2)) {
+        G_debug(DLEVEL, "        no intersection");
+        return 0;
+    }
+
+    /* there is overlap or connected end points */
+    G_debug(DLEVEL, "        overlap");
+    
+    /* a contains b */
+    if ((ax1 < bx1) && (ax2 > bx2)) {
+        G_debug(DLEVEL, "            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(DLEVEL, "            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(DLEVEL, "        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;
+}
+
+#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);
+}
+
+#ifdef ASDASDFASFEAS
+int segment_intersection_2d_test(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;
+    double d, da, db, ra, rb;
+    
+    if (!initialized)
+        initialize_mpf_vars();
+    
+    /* TODO: Works for points ? */
+    G_debug(0, "segment_intersection_2d_test()"); 
+    G_debug(0, "    ax1  = %.18e, ay1  = %.18e", ax1, ay1);
+    G_debug(0, "    ax2  = %.18e, ay2  = %.18e", ax2, ay2);
+    G_debug(0, "    bx1  = %.18e, by1  = %.18e", bx1, by1);
+    G_debug(0, "    bx2  = %.18e, by2  = %.18e", 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 (0, "    identical segments" );
+        *x1 = ax1; *y1 = ay1;
+        *x2 = ax2; *y2 = ay2;
+        return 5;
+    }
+    /* Check for identical endpoints */
+    if (f11 || f12) {
+        G_debug (0, "    connected by endpoints" );
+        *x1 = ax1; *y1 = ay1;
+        return 1;
+    }
+    if (f21 || f22) {
+        G_debug (0, "    connected by endpoints" );
+        *x1 = ax2; *y1 = ay2;
+        return 1;
+    }
+
+    if ((MAX(ax1, ax2) < MIN(bx1, bx2)) || (MAX(bx1, bx2) < MIN(ax1, ax2))) {
+        G_debug(0, "    no intersection (disjoint bounding boxes)");
+        return 0;
+    }
+    if ((MAX(ay1, ay2) < MIN(by1, by2)) || (MAX(by1, by2) < MIN(ay1, ay2))) {
+        G_debug(0, "    no intersection (disjoint bounding boxes)");
+        return 0;
+    }
+    
+    d  = (ax2-ax1)*(by1-by2) - (ay2-ay1)*(bx1-bx2);
+    da = (bx1-ax1)*(by1-by2) - (by1-ay1)*(bx1-bx2);
+    db = (ax2-ax1)*(by1-ay1) - (ay2-ay1)*(bx1-ax1);
+    
+    det22(dd,  ax2, ax1, bx1, bx2, ay2, ay1, by1, by2);
+    sgn_d = mpf_sgn(dd);
+    s = mpf_get_str(NULL, &exp, 10, 40, dd);
+    G_debug(0, "    dd = %sE%d", (s[0]==0)?"0":s, exp);
+    G_debug(0, "    d = %.18E", d);
+    
+    if (sgn_d != 0) {
+        G_debug(0, "    general position");
+
+        det22(dda, bx1, ax1, bx1, bx2, by1, ay1, by1, by2);
+        det22(ddb, ax2, ax1, bx1, ax1, ay2, ay1, by1, ay1);
+        sgn_da = mpf_sgn(dda);
+        sgn_db = mpf_sgn(ddb);
+        
+        ra = da/d;
+        rb = db/d;
+        mpf_div(rra, dda, dd);
+        mpf_div(rrb, ddb, dd);
+
+        s = mpf_get_str(NULL, &exp, 10, 40, rra);
+        G_debug(0, "        rra = %sE%d", (s[0]==0)?"0":s, exp);
+        G_debug(0, "        ra = %.18E", ra);
+        s = mpf_get_str(NULL, &exp, 10, 40, rrb);
+        G_debug(0, "        rrb = %sE%d", (s[0]==0)?"0":s, exp);
+        G_debug(0, "        rb = %.18E", rb);
+        
+        if (sgn_d > 0) {
+            if ((sgn_da < 0) || (mpf_cmp(dda, dd) > 0)) {
+                G_debug(DLEVEL, "        no intersection");
+                return 0;
+            }
+            
+            if ((sgn_db < 0) || (mpf_cmp(ddb, dd) > 0)) {
+                G_debug(DLEVEL, "        no intersection");
+                return 0;
+            }
+        }
+        else { /* if sgn_d < 0 */
+            if ((sgn_da > 0) || (mpf_cmp(dda, dd) < 0)) {
+                G_debug(DLEVEL, "        no intersection");
+                return 0;
+            }
+            
+            if ((sgn_db > 0) || (mpf_cmp(ddb, dd) < 0)) {
+                G_debug(DLEVEL, "        no intersection");
+                return 0;
+            }
+        }
+
+        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(0, "        intersection at:");
+        G_debug(0, "            xx = %.18e", *x1);
+        G_debug(0, "             x = %.18e", ax1 + ra*(ax2-ax1));
+        G_debug(0, "            yy = %.18e", *y1);
+        G_debug(0, "             y = %.18e", ay1 + ra*(ay2-ay1));
+        return 1;
+    }
+    
+    G_debug(0, "    parallel/collinear...");
+    return -1;
+}
+#endif


Property changes on: grass/branches/develbranch_6/vector/v.buffer2/e_intersect.c
___________________________________________________________________
Name: svn:eol-style
   + native

Added: grass/branches/develbranch_6/vector/v.buffer2/e_intersect.h
===================================================================
--- grass/branches/develbranch_6/vector/v.buffer2/e_intersect.h	                        (rev 0)
+++ grass/branches/develbranch_6/vector/v.buffer2/e_intersect.h	2008-11-26 21:23:44 UTC (rev 34511)
@@ -0,0 +1,21 @@
+#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_test(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_tol(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 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);
+
+    
+int almost_equal(double a, double b, int bits);
+
+#endif


Property changes on: grass/branches/develbranch_6/vector/v.buffer2/e_intersect.h
___________________________________________________________________
Name: svn:eol-style
   + native

Added: grass/branches/develbranch_6/vector/v.buffer2/local_proto.h
===================================================================
--- grass/branches/develbranch_6/vector/v.buffer2/local_proto.h	                        (rev 0)
+++ grass/branches/develbranch_6/vector/v.buffer2/local_proto.h	2008-11-26 21:23:44 UTC (rev 34511)
@@ -0,0 +1,9 @@
+#ifndef LOCAL_PROTO_H
+#define LOCAL_PROTO_H
+
+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);
+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);
+void Vect_point_buffer2(double px, double py, double da, double db, double dalpha, int round, double tol, struct line_pnts **oPoints);
+void Vect_line_parallel2(struct line_pnts *InPoints, double da, double db, double dalpha, int side, int round, double tol, struct line_pnts *OutPoints);
+
+#endif


Property changes on: grass/branches/develbranch_6/vector/v.buffer2/local_proto.h
___________________________________________________________________
Name: svn:eol-style
   + native

Added: grass/branches/develbranch_6/vector/v.buffer2/main.c
===================================================================
--- grass/branches/develbranch_6/vector/v.buffer2/main.c	                        (rev 0)
+++ grass/branches/develbranch_6/vector/v.buffer2/main.c	2008-11-26 21:23:44 UTC (rev 34511)
@@ -0,0 +1,552 @@
+/****************************************************************
+ *
+ * MODULE:       v.buffer
+ * 
+ * AUTHOR(S):    Radim Blazek, Rosen Matev
+ *               
+ * PURPOSE:      Vector buffer
+ *               
+ * COPYRIGHT:    (C) 2001 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.
+ *
+ **************************************************************/
+#include <stdlib.h>
+#include <string.h>
+#include <grass/gis.h>
+#include <grass/Vect.h>
+#include <grass/dbmi.h>
+#include <grass/glocale.h>
+#include "local_proto.h"
+
+#define PI M_PI
+#ifndef MIN
+    #define MIN(X,Y) ((X<Y)?X:Y)
+#endif
+#ifndef MAX
+    #define MAX(X,Y) ((X>Y)?X:Y)
+#endif    
+
+
+void stop(struct Map_info *In, struct Map_info *Out) {
+    Vect_close (In);
+
+    G_message ( _("Rebuilding topology...") );
+    Vect_build_partial ( Out, GV_BUILD_NONE);
+    Vect_build (Out);
+    Vect_close (Out);
+}
+
+/* returns 1 if unit_tolerance is adjusted, 0 otherwise */
+int adjust_tolerance(double *tolerance) {
+    double t = 0.999 * ( 1 - cos ( 2 * PI / 8 / 2 ) );
+    G_debug(2, "Maximum tolerance = %f", t);
+    if (*tolerance > t) {
+        *tolerance = t;
+        return 1;
+    }
+    return 0;
+}
+
+int db_CatValArray_get_value_di(dbCatValArray *cvarr, int cat, double *value) {
+    int t;
+    int ctype = cvarr->ctype;
+    int ret;
+    
+    if (ctype == DB_C_TYPE_INT) {
+        ret = db_CatValArray_get_value_int(cvarr, cat, &t);
+        if (ret != DB_OK)
+            return ret;
+        *value = (double)t;
+        return DB_OK;
+    }
+
+    if (ctype == DB_C_TYPE_DOUBLE) {
+        ret = db_CatValArray_get_value_double(cvarr, cat, value);
+        return ret;
+    }
+    
+    return DB_FAILED;
+}
+
+struct buf_contours {
+    int inner_count;
+    struct line_pnts *oPoints;
+    struct line_pnts **iPoints;
+};
+
+int point_in_buffer(struct buf_contours *arr_bc, int buffers_count, struct Map_info *Out, double x, double y) {
+    int i, j, ret, flag;
+    
+    for (i = 0; i < buffers_count; i++) {
+        ret = Vect_point_in_poly(x, y, arr_bc[i].oPoints);
+        if (ret == 0)
+            continue;
+            
+        flag = 1;
+        for (j = 0; j < arr_bc[i].inner_count; j++) {
+            ret = Vect_point_in_poly(x, y, arr_bc[i].iPoints[j]);
+            if (ret != 0) { /* inside inner contour */
+                flag = 0;
+                break;
+            }
+        }
+        
+        if (flag) {
+            /* (x,y) is inside outer contour and outside inner contours of arr_bc[i] */
+            return 1;
+        }
+    }
+    
+    return 0;
+}
+
+int main (int argc, char *argv[])
+{
+    struct Map_info In, Out;
+    struct line_pnts *Points;
+    struct line_cats *Cats, *BCats;
+    char *mapset;
+    struct GModule *module;
+    struct Option *in_opt, *out_opt, *type_opt, *dista_opt, *distb_opt, *angle_opt;
+    struct Flag *straight_flag, *nocaps_flag;
+    struct Option *tol_opt, *bufcol_opt, *scale_opt, *field_opt;
+    double da, db, dalpha, tolerance, unit_tolerance;
+    int type;
+    int i, j, ret, nareas, area, nlines, line;
+    char *Areas, *Lines;
+    int field;
+    struct buf_contours *arr_bc;
+    int buffers_count;
+
+ /* Attributes if sizecol is used */
+    int nrec, ctype;
+    struct field_info *Fi;
+    dbDriver *Driver;
+    dbCatValArray cvarr;
+    double size_val, scale;
+
+
+    module = G_define_module();
+    module->keywords = _("vector");
+    module->description = _("Creates a buffer around features of given type (areas must contain centroid).");
+    
+    in_opt = G_define_standard_option(G_OPT_V_INPUT);
+    out_opt = G_define_standard_option(G_OPT_V_OUTPUT);
+    
+    type_opt = G_define_standard_option(G_OPT_V_TYPE) ;
+    type_opt->options = "point,line,boundary,centroid,area";
+    type_opt->answer = "point,line,area";
+    
+    field_opt = G_define_standard_option(G_OPT_V_FIELD);
+    
+    dista_opt = G_define_option();
+    dista_opt->key = "distance";
+    dista_opt->type = TYPE_DOUBLE;
+    dista_opt->required = NO;
+    dista_opt->description = _("Buffer distance along major axis in map units");
+    
+    distb_opt = G_define_option();
+    distb_opt->key = "minordistance";
+    distb_opt->type = TYPE_DOUBLE;
+    distb_opt->required = NO;
+    distb_opt->description = _("Buffer distance along minor axis in map units");
+    distb_opt->guisection  = _("Advanced");
+    
+    angle_opt = G_define_option();
+    angle_opt->key = "angle";
+    angle_opt->type =  TYPE_DOUBLE;
+    angle_opt->required = NO;
+    angle_opt->answer = "0";
+    angle_opt->description = _("Angle of major axis in degrees");
+    angle_opt->guisection  = _("Advanced");
+    
+    bufcol_opt = G_define_option();
+    bufcol_opt->key = "bufcol";
+    bufcol_opt->type = TYPE_STRING;
+    bufcol_opt->required = NO;
+    bufcol_opt->description = _("Attribute column to use for buffer distances");
+    bufcol_opt->guisection  = _("Advanced");
+    
+    scale_opt = G_define_option();
+    scale_opt->key = "scale";
+    scale_opt->type = TYPE_DOUBLE;
+    scale_opt->required = NO;
+    scale_opt->answer = "1.0";
+    scale_opt->description = _("Scaling factor for attribute column values");
+    scale_opt->guisection  = _("Advanced");
+
+    tol_opt = G_define_option();
+    tol_opt->key = "tolerance";
+    tol_opt->type = TYPE_DOUBLE;
+    tol_opt->required = NO;
+    tol_opt->answer = "0.01";
+    tol_opt->guisection = _("Advanced");    
+    tol_opt->description = _("Maximum distance between theoretical arc and polygon segments as multiple of buffer");
+
+    straight_flag = G_define_flag();
+    straight_flag->key = 's';
+    straight_flag->description = _("Make outside corners straight");
+
+    nocaps_flag = G_define_flag();
+    nocaps_flag->key = 'c';
+    nocaps_flag->description = _("Don't make caps at the ends of polylines");
+
+    G_gisinit(argv[0]);
+    if (G_parser (argc, argv))
+        exit(EXIT_FAILURE);
+    
+    type = Vect_option_to_types ( type_opt );
+    field = atoi( field_opt->answer );
+    
+    if ((dista_opt->answer && bufcol_opt->answer) || (!(dista_opt->answer || bufcol_opt->answer)))
+        G_fatal_error("Select a buffer distance/minordistance/angle or column, but not both.");
+    
+    if (bufcol_opt->answer)
+        G_warning(_("The bufcol option may contain bugs during the cleaning "
+            "step. If you encounter problems, use the debug "
+            "option or clean manually with v.clean tool=break; "
+            "v.category step=0; v.extract -d type=area"));
+    
+    tolerance = atof(tol_opt->answer);
+    if (adjust_tolerance(&tolerance))
+        G_warning(_("The tolerance was reset to %g."), tolerance);
+    
+    scale = atof( scale_opt->answer );
+    if (scale <= 0.0)
+        G_fatal_error("Illegal scale value");
+
+
+    if (dista_opt->answer) {
+        da = atof(dista_opt->answer);
+        
+        if (distb_opt->answer)
+            db = atof(distb_opt->answer);
+        else
+            db = da;
+            
+        if (angle_opt->answer)
+            dalpha = atof(angle_opt->answer);
+        else
+            dalpha = 0;
+        
+        unit_tolerance = tolerance * MIN(da, db);
+        G_message(_("The tolerance in map units = %g"), unit_tolerance);
+    }
+
+    Vect_check_input_output_name(in_opt->answer, out_opt->answer, GV_FATAL_EXIT);
+    
+    Points = Vect_new_line_struct();
+    Cats = Vect_new_cats_struct();
+    BCats = Vect_new_cats_struct();
+
+    /* open input vector */
+    if ((mapset = G_find_vector2(in_opt->answer, "")) == NULL)
+        G_fatal_error(_("Vector map <%s> not found"), in_opt->answer);
+
+    Vect_set_open_level(2);
+
+    if (1 > Vect_open_old(&In, in_opt->answer, mapset))
+        G_fatal_error(_("Unable to open vector map <%s>"), in_opt->answer);
+
+    if (0 > Vect_open_new(&Out, out_opt->answer, WITHOUT_Z)) {
+        Vect_close(&In);
+        G_fatal_error(_("Unable to create vector map <%s>"), out_opt->answer);
+    }
+    
+    /* check and load attribute column data */
+    if (bufcol_opt->answer) {
+        db_CatValArray_init ( &cvarr );
+        
+        Fi = Vect_get_field(&In, field);
+        if ( Fi == NULL )
+            G_fatal_error(_("Unable to get layer info for vector map"));
+        
+        Driver = db_start_driver_open_database(Fi->driver, Fi->database);
+        if (Driver == NULL)
+            G_fatal_error(_("Unable to open database <%s> by driver <%s>"), Fi->database, Fi->driver);
+        
+        /* Note do not check if the column exists in the table because it may be expression */
+
+        /* TODO: only select values we need instead of all in column */
+        nrec = db_select_CatValArray(Driver, Fi->table, Fi->key, bufcol_opt->answer, NULL, &cvarr);
+        if (nrec < 0)
+            G_fatal_error(_("Unable to select data from table"));
+        G_debug(2, "%d records selected from table", nrec);
+
+        ctype = cvarr.ctype;
+        if (ctype != DB_C_TYPE_INT && ctype != DB_C_TYPE_DOUBLE)
+            G_fatal_error(_("Column type not supported"));
+        
+        db_close_database_shutdown_driver(Driver);
+        
+        /* Output cats/values list */
+        for (i = 0; i < cvarr.n_values; i++) {
+            if (ctype == DB_C_TYPE_INT) {
+                G_debug(4, "cat = %d val = %d", cvarr.value[i].cat, cvarr.value[i].val.i);
+            } else if (ctype == DB_C_TYPE_DOUBLE) {
+                G_debug(4, "cat = %d val = %f", cvarr.value[i].cat, cvarr.value[i].val.d);
+            }
+        }
+    }
+        
+    Vect_copy_head_data(&In, &Out);
+    Vect_hist_copy(&In, &Out);
+    Vect_hist_command(&Out);
+
+
+    /* Create buffers' boundaries */
+    nlines = Vect_get_num_lines(&In);
+    nareas = Vect_get_num_areas (&In);
+    /* TODO: don't allocate so much space */
+    buffers_count = 0;
+    arr_bc = G_malloc((nlines+nareas)*sizeof(struct buf_contours));
+    
+    /* Lines (and Points) */
+    if ((type & GV_POINTS) || (type & GV_LINES)) {
+        int ltype;
+        
+        G_message(_("Lines buffers... "));
+        for (line = 1; line <= nlines; line++) {
+            int cat;
+            
+            G_debug(2, "line = %d", line);
+            G_percent(line, nlines, 2);
+
+            ltype = Vect_read_line(&In, Points, Cats, line);
+            if (!(ltype & type))
+                continue;
+
+            if (!Vect_cat_get(Cats, field, &cat))
+                continue;
+
+            if (bufcol_opt->answer) {
+                ret = db_CatValArray_get_value_di(&cvarr, cat, &size_val);
+                if (ret != DB_OK) {
+                    G_warning(_("No record for category %d in table <%s>"), cat, Fi->table);
+                    continue;
+                }
+                
+                if (size_val < 0.0) {
+                    G_warning(_("Attribute is of invalid size (%.3f) for category %d"), size_val, cat);
+                    continue;
+                }
+                
+                if (size_val == 0.0)
+                    continue;
+                
+                da = size_val * scale;
+                db = da;
+                dalpha = 0;
+                unit_tolerance = tolerance * MIN(da, db);
+                
+                G_debug(2, "    dynamic buffer size = %.2f", da);                
+                G_debug(2, _("The tolerance in map units: %g"), unit_tolerance);
+            }
+            
+            if (ltype & GV_POINTS) {
+                Vect_point_buffer2(Points->x[0], Points->y[0], da, db, dalpha, !(straight_flag->answer), unit_tolerance, &(arr_bc[buffers_count].oPoints));
+                arr_bc[buffers_count].iPoints = NULL;
+                arr_bc[buffers_count].inner_count = 0;
+                buffers_count++;
+                    
+            }
+            else {
+                Vect_line_buffer2(Points, da, db, dalpha, !(straight_flag->answer), !(nocaps_flag->answer), unit_tolerance, &(arr_bc[buffers_count].oPoints), &(arr_bc[buffers_count].iPoints), &(arr_bc[buffers_count].inner_count));
+                buffers_count++;
+            }
+        }
+    }
+
+    /* Areas */
+    if (type & GV_AREA) {
+        int centroid;
+        
+        G_message ( _("Areas buffers... "));
+        for (area = 1; area <= nareas; area++) {
+            int cat;
+            G_percent(area, nareas, 2);
+            centroid = Vect_get_area_centroid(&In, area);
+            if (centroid == 0)
+                continue;
+            
+            Vect_read_line(&In, NULL, Cats, centroid);
+            if (!Vect_cat_get(Cats, field, &cat))
+                continue;
+            
+            if (bufcol_opt->answer) {
+                ret = db_CatValArray_get_value_di(&cvarr, cat, &size_val);
+                if (ret != DB_OK) {
+                    G_warning(_("No record for category %d in table <%s>"), cat, Fi->table);
+                    continue;
+                }
+                
+                if (size_val < 0.0) {
+                    G_warning(_("Attribute is of invalid size (%.3f) for category %d"), size_val, cat);
+                    continue;
+                }
+                
+                if (size_val == 0.0)
+                    continue;
+                
+                da = size_val * scale;
+                db = da;
+                dalpha = 0;
+                unit_tolerance = tolerance * MIN(da, db);
+                
+                G_debug(2, "    dynamic buffer size = %.2f", da);                
+                G_debug(2, _("The tolerance in map units: %g"), unit_tolerance);
+            }
+            
+            Vect_area_buffer2(&In, area, da, db, dalpha, !(straight_flag->answer), !(nocaps_flag->answer), unit_tolerance, &(arr_bc[buffers_count].oPoints), &(arr_bc[buffers_count].iPoints), &(arr_bc[buffers_count].inner_count));
+            buffers_count++;
+        }
+    }
+
+    /* write all buffer contours */
+    for (i = 0; i < buffers_count; i++) {
+        Vect_write_line(&Out, GV_BOUNDARY, arr_bc[i].oPoints, BCats);
+        for (j = 0; j < arr_bc[i].inner_count; j++)
+            Vect_write_line(&Out, GV_BOUNDARY, arr_bc[i].iPoints[j], BCats);
+    }
+
+    /* Create areas */
+    
+    /* Break lines */
+    G_message (_("Building parts of topology...") );
+    Vect_build_partial(&Out, GV_BUILD_BASE);
+
+    G_message(_("Snapping boundaries...") );
+    Vect_snap_lines ( &Out, GV_BOUNDARY, 1e-7, NULL);
+
+    G_message(_( "Breaking boundaries...") );
+    Vect_break_lines ( &Out, GV_BOUNDARY, NULL );
+
+    G_message(_( "Removing duplicates...") );
+    Vect_remove_duplicates ( &Out, GV_BOUNDARY, NULL );
+
+    /* Dangles and bridges don't seem to be necessary if snapping is small enough. */
+    /*
+    G_message (  "Removing dangles..." );
+    Vect_remove_dangles ( &Out, GV_BOUNDARY, -1, NULL, stderr );
+
+    G_message (  "Removing bridges..." );
+    Vect_remove_bridges ( &Out, NULL, stderr );
+    */
+
+    G_message(_( "Attaching islands..."));
+    Vect_build_partial(&Out, GV_BUILD_ATTACH_ISLES);
+    
+    /* Calculate new centroids for all areas */
+    nareas = Vect_get_num_areas ( &Out );
+    Areas = (char *) G_calloc ( nareas+1, sizeof(char) );
+    for (area = 1; area <= nareas; area++) {
+        double x, y;
+        
+        G_debug ( 3, "area = %d", area );
+
+        if (!Vect_area_alive(&Out, area))
+            continue;
+                    
+        ret = Vect_get_point_in_area(&Out, area, &x, &y);
+        if ( ret < 0 ) {
+            G_warning (_("Cannot calculate area centroid"));
+            continue;
+        }
+        
+        ret = point_in_buffer(arr_bc, buffers_count, &Out, x, y);
+        
+        if (ret) {
+            G_debug (3, "  -> in buffer");
+            Areas[area] = 1;
+        }
+    }
+    
+    /* Make a list of boundaries to be deleted (both sides inside) */
+    nlines = Vect_get_num_lines(&Out);
+    G_debug ( 3, "nlines = %d", nlines );
+    Lines = (char *) G_calloc ( nlines+1, sizeof(char) );
+
+    for ( line = 1; line <= nlines; line++ ) {
+        int j, side[2], areas[2];
+        
+        G_debug ( 3, "line = %d", line );
+        
+        if ( !Vect_line_alive ( &Out, line ) )
+            continue;
+        
+        Vect_get_line_areas ( &Out, line, &side[0], &side[1] );
+        
+        for ( j = 0; j < 2; j++ ) { 
+            if ( side[j] == 0 ) { /* area/isle not build */
+                areas[j] = 0;
+            } else if ( side[j] > 0 ) { /* area */
+                areas[j] = side[j];
+            } else { /* < 0 -> island */
+                areas[j] = Vect_get_isle_area ( &Out, abs ( side[j] ) );
+            }
+        }
+        
+        G_debug ( 3, " areas = %d , %d -> Areas = %d, %d", areas[0], areas[1], Areas[areas[0]], Areas[areas[1]] );
+        if ( Areas[areas[0]] && Areas[areas[1]] )
+            Lines[line] = 1;
+    }
+    G_free(Areas);
+    
+    /* Delete boundaries */
+    for ( line = 1; line <= nlines; line++ ) {
+        if ( Lines[line] ) {
+            G_debug ( 3, " delete line %d", line );
+            Vect_delete_line ( &Out, line );
+        }
+    }
+
+    G_free( Lines );
+
+    /* Create new centroids */
+    Vect_reset_cats ( Cats );
+    Vect_cat_set (Cats, 1, 1 );
+    nareas = Vect_get_num_areas(&Out);
+
+    for (area = 1; area <= nareas; area++) {
+        double x, y;
+        
+        G_debug ( 3, "area = %d", area );
+
+        if (!Vect_area_alive(&Out, area))
+            continue;
+                    
+        ret = Vect_get_point_in_area(&Out, area, &x, &y);
+        if ( ret < 0 ) {
+            G_warning (_("Cannot calculate area centroid"));
+            continue;
+        }
+        
+        ret = point_in_buffer(arr_bc, buffers_count, &Out, x, y);
+        
+        if (ret) {
+            Vect_reset_line(Points);
+            Vect_append_point(Points, x, y, 0.);
+            Vect_write_line(&Out, GV_CENTROID, Points, Cats);
+        }
+    }
+
+    /* free arr_bc[] */
+    /* will only slow down the module
+    for (i = 0; i < buffers_count; i++) {
+        Vect_destroy_line_struct(arr_bc[i].oPoints);
+        for (j = 0; j < arr_bc[i].inner_count; j++)
+            Vect_destroy_line_struct(arr_bc[i].iPoints[j]);
+        G_free(arr_bc[i].iPoints);
+    } */
+    
+    G_message(_("Attaching centroids...") );
+    Vect_build_partial ( &Out, GV_BUILD_CENTROIDS);
+    
+    stop ( &In, &Out );
+    exit(EXIT_SUCCESS);
+}
+


Property changes on: grass/branches/develbranch_6/vector/v.buffer2/main.c
___________________________________________________________________
Name: svn:eol-style
   + native

Added: grass/branches/develbranch_6/vector/v.buffer2/vlib_buffer.c
===================================================================
--- grass/branches/develbranch_6/vector/v.buffer2/vlib_buffer.c	                        (rev 0)
+++ grass/branches/develbranch_6/vector/v.buffer2/vlib_buffer.c	2008-11-26 21:23:44 UTC (rev 34511)
@@ -0,0 +1,1012 @@
+/* Functions: ...
+**
+** Author: Radim Blazek, Rosen Matev; July 2008
+**
+**
+*/
+#include <stdlib.h>
+#include <grass/Vect.h>
+#include <grass/gis.h>
+#include "dgraph.h"
+
+#define LENGTH(DX, DY) (sqrt((DX*DX)+(DY*DY)))
+#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
+#define RIGHT_SIDE 1
+#define LEFT_SIDE -1
+#define LOOPED_LINE 1
+#define NON_LOOPED_LINE 0
+
+/* norm_vector() calculates normalized vector form two points */
+static void norm_vector(double x1, double y1, double x2, double y2, double *x, double *y) {
+    double dx, dy, l;
+    dx  = x2 - x1;
+    dy  = y2 - y1;
+    if ((dx == 0) && (dy == 0)) {
+        /* assume that dx == dy == 0, which should give (NaN,NaN) */
+        /* without this, very small dx or dy could result in Infinity */
+        *x = 0;
+        *y = 0;
+        return;
+    }
+    l = LENGTH(dx, dy);
+    *x = dx/l;
+    *y = dy/l;
+    return;
+}
+
+static void rotate_vector(double x, double y, double cosa, double sina, double *nx, double *ny) {
+    *nx = x*cosa - y*sina;
+    *ny = x*sina + y*cosa;
+    return;   
+}
+
+/*
+* (x,y) shoud be normalized vector for common transforms; This func transforms (x,y) to a vector corresponding to da, db, dalpha params
+* dalpha is in radians
+*/
+static void elliptic_transform(double x, double y, double da, double db, double dalpha, double *nx, double *ny) {
+    double cosa = cos(dalpha);
+    double sina = sin(dalpha);
+/*    double cc = cosa*cosa;
+    double ss = sina*sina;
+    double t = (da-db)*sina*cosa;
+    
+    *nx = (da*cc + db*ss)*x + t*y;
+    *ny = (da*ss + db*cc)*y + t*x;
+    return;*/
+    
+    double va, vb;
+    va = (x*cosa + y*sina)*da;
+    vb = (x*(-sina) + y*cosa)*db;
+    *nx = va*cosa + vb*(-sina);
+    *ny = va*sina + vb*cosa;
+    return; 
+}
+
+/*
+* vect(x,y) must be normalized
+* gives the tangent point of the tangent to ellpise(da,db,dalpha) parallel to vect(x,y)
+* dalpha is in radians
+* ellipse center is in (0,0)
+*/
+static void elliptic_tangent(double x, double y, double da, double db, double dalpha, double *px, double *py) {
+    double cosa = cos(dalpha);
+    double sina = sin(dalpha);
+    double u, v, len;
+    /* rotate (x,y) -dalpha radians */
+    rotate_vector(x, y, cosa, -sina, &x, &y);
+    /*u = (x + da*y/db)/2;
+    v = (y - db*x/da)/2;*/
+    u = da*da*y;
+    v = -db*db*x;
+    len = da*db/sqrt(da*da*v*v + db*db*u*u);
+    u *= len;
+    v *= len;
+    rotate_vector(u, v, cosa, sina, px, py);
+    return; 
+}
+
+
+/*
+* !!! This is not line in GRASS' sense. See http://en.wikipedia.org/wiki/Line_%28mathematics%29
+*/
+static void line_coefficients(double x1, double y1, double x2, double y2, double *a, double *b, double *c) {
+    *a = y2 - y1;
+    *b = x1 - x2;
+    *c = x2*y1 - x1*y2;
+    return;
+}
+
+/*
+* Finds intersection of two straight lines. Returns 0 if the lines are parallel, 1 if they cross,
+* 2 if they are the same line.
+* !!!!!!!!!!!!!!!! FIX THIS TOLLERANCE CONSTANTS BAD (and UGLY) CODE !!!!!!!!!
+*/
+static int line_intersection(double a1, double b1, double c1, double a2, double b2, double c2, double *x, double *y) {
+    double d;
+    
+    if (fabs(a2*b1 - a1*b2) == 0) {
+        if (fabs(a2*c1 - a1*c2) == 0)
+            return 2;
+        else                   
+            return 0;
+    }
+    else {
+        d = a1*b2 - a2*b1;
+        *x = (b1*c2 - b2*c1)/d;
+        *y = (c1*a2 - c2*a1)/d;
+        return 1;
+    }
+}
+
+static double angular_tolerance(double tol, double da, double db) {
+    double a = MAX(da, db);
+    if (tol > a)
+        tol = a;
+    return 2*acos(1-tol/a);
+}
+
+/*
+* This function generates parallel line (with loops, but not like the old ones).
+* It is not to be used directly for creating buffers.
+* + added elliptical buffers/par.lines support
+*
+* dalpha - direction of elliptical buffer major axis in degrees
+* da - distance along major axis
+* db: distance along minor (perp.) axis
+* side: side >= 0 - right side, side < 0 - left side
+* when (da == db) we have plain distances (old case)
+* round - 1 for round corners, 0 for sharp corners. (tol is used only if round == 1)
+*/
+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, j, res, np;
+    double *x, *y;
+    double tx, ty, vx, vy, wx, wy, nx, ny, mx, my, rx, ry;
+    double vx1, vy1, wx1, wy1;
+    double a0, b0, c0, a1, b1, c1;
+    double phi1, phi2, delta_phi;
+    double nsegments, angular_tol, angular_step;
+    int inner_corner, turns360;
+    
+    G_debug(3, "parallel_line()");
+    
+    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;
+
+    if ((da == 0) || (db == 0)) {
+        Vect_copy_xyz_to_pnts(nPoints, x, y, NULL, np);
+        return;
+    }
+
+    side = (side >= 0)?(1):(-1); /* normalize variable */
+    dalpha *= PI/180; /* convert dalpha from degrees to radians */
+    angular_tol = angular_tolerance(tol, da, db);
+    
+    for (i = 0; i < np-1; i++)
+    {
+        /* 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);
+        if ((tx == 0) && (ty == 0))
+            continue;
+        
+        elliptic_tangent(side*tx, side*ty, da, db, dalpha, &vx, &vy);
+        
+        nx = x[i] + vx;
+        ny = y[i] + vy;
+        
+        mx = x[i+1] + vx;
+        my = y[i+1] + vy;
+
+        line_coefficients(nx, ny, mx, my, &a1, &b1, &c1);
+
+        if (i == 0) {
+            if (!looped)
+                Vect_append_point(nPoints, nx, ny, 0);
+            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] */
+        turns360 = (fabs(fabs(delta_phi) - PI) < 1e-15);
+        inner_corner = (side*delta_phi <= 0) && (!turns360);
+        
+        if ((turns360) && (!(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) {
+                if (!round)
+                    Vect_append_point(nPoints, rx, ry, 0);
+                else {
+/*                    d = dig_distance2_point_to_line(rx, ry, 0, x[i-1], y[i-1], 0, x[i], y[i], 0,
+                        0, NULL, NULL, NULL, NULL, NULL);
+                    if (*/
+                    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 ((!looped) && (i == np-2)) {
+            Vect_append_point(nPoints, mx, my, 0);
+        }
+    }
+    
+    if (looped) {
+        Vect_append_point(nPoints, nPoints->x[0], nPoints->y[0], nPoints->z[0]);
+    }
+    
+    Vect_line_prune ( nPoints );
+    
+    if (looped) {
+        Vect_line_delete_point(Points, Points->n_points-1);
+    }
+}
+
+/* input line must be looped */
+void convolution_line(struct line_pnts *Points, double da, double db, double dalpha, int side, int round, int caps, double tol, struct line_pnts *nPoints)
+{
+    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;
+    double a0, b0, c0, a1, b1, c1;
+    double phi1, phi2, delta_phi;
+    double nsegments, angular_tol, angular_step;
+    double angle0, angle1;
+    int inner_corner, turns360;
+    
+    G_debug(3, "convolution_line() side = %d", side);
+
+    np = Points->n_points;
+    x = Points->x;
+    y = Points->y;
+    if ((np == 0) || (np == 1))
+        return;
+    if ((x[0] != x[np-1]) || (y[0] != y[np-1])) {
+        G_fatal_error("line is not looped");
+        return;
+    }
+    
+    Vect_reset_line(nPoints);
+    
+    if ((da == 0) || (db == 0)) {
+        Vect_copy_xyz_to_pnts(nPoints, x, y, NULL, np);
+        return;
+    }
+
+    side = (side >= 0)?(1):(-1); /* normalize variable */
+    dalpha *= PI/180; /* convert dalpha from degrees to radians */
+    angular_tol = angular_tolerance(tol, da, db);
+    
+    i = np-2;
+    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);
+    angle1 = atan2(ty, tx);
+    nx = x[i] + vx;
+    ny = y[i] + vy;
+    mx = x[i+1] + vx;
+    my = y[i+1] + vy;
+    if (!round)
+        line_coefficients(nx, ny, mx, my, &a1, &b1, &c1);
+    
+    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;
+            b0 = b1;
+            c0 = c1;
+        }
+        wx = vx;
+        wy = vy;
+        angle0 = angle1;
+        
+        norm_vector(x[i], y[i], x[i+1], y[i+1], &tx, &ty);
+        if ((tx == 0) && (ty == 0))
+            continue;
+        elliptic_tangent(side*tx, side*ty, da, db, dalpha, &vx, &vy);
+        angle1 = atan2(ty, tx);
+        nx = x[i] + vx;
+        ny = y[i] + vy;
+        mx = x[i+1] + vx;
+        my = y[i+1] + vy;
+        if (!round)
+            line_coefficients(nx, ny, mx, my, &a1, &b1, &c1);
+
+        
+        delta_phi = angle1 - angle0;
+        if (delta_phi > PI)
+            delta_phi -= 2*PI;
+        else if (delta_phi <= -PI)
+            delta_phi += 2*PI;
+        /* now delta_phi is in [-pi;pi] */
+        turns360 = (fabs(fabs(delta_phi) - PI) < 1e-15);
+        inner_corner = (side*delta_phi <= 0) && (!turns360);
+
+        
+        /* if <line turns 360> and (<caps> and <not round>) */
+        if (turns360 && caps && (!round)) {
+            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 ((!turns360) && (!round) && (!inner_corner)) {
+            res = line_intersection(a0, b0, c0, a1, b1, c1, &rx, &ry);
+            if (res == 1) {
+                Vect_append_point(nPoints, rx, ry, 0);
+                G_debug(4, " append point (o) x=%.16f y=%.16f", rx, ry);
+            }
+            else if (res == 2) {
+                /* no need to append point in this case */
+            }
+            else
+                G_fatal_error("unexpected result of line_intersection() res = %d", res);
+        }
+        
+        if (round && (!inner_corner) && (!turns360 || caps)) {
+            /* 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);
+            
+            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 */
+    Vect_append_point(nPoints, nPoints->x[0], nPoints->y[0], nPoints->z[0]);
+/*    Vect_line_prune ( nPoints ); */
+}
+
+/*
+* side: side >= 0 - extracts contour on right side of edge, side < 0 - extracts contour on left side of edge
+* if the extracted contour is the outer contour, it is returned in ccw order
+* else if it is inner contour, it is returned in cw order
+*/
+static void extract_contour(struct planar_graph *pg, struct pg_edge *first, int side, int winding, int stop_at_line_end, struct line_pnts *nPoints) {
+    int j;
+    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 */
+    /* we will always turn right and dont need that one */
+    double opt_angle, tangle;
+    int opt_j, opt_side, opt_flag;
+    
+    G_debug(3, "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);
+    
+    edge = first;
+    if (side >= 0) {
+        eside = 1;
+        v0 = edge->v1;
+        v = edge->v2;
+    }
+    else {
+        eside = -1;
+        v0 = edge->v2;
+        v = edge->v1;
+    }
+    vert0 = &(pg->v[v0]);    
+    vert = &(pg->v[v]);
+    eangle = atan2(vert->y - vert0->y, vert->x - vert0->x);
+    
+    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: 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;
+            edge->winding_right = winding;
+        }
+        else {
+            edge->visited_left = 1;
+            edge->winding_left = winding;
+        }
+        
+        opt_flag = 1;
+        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 || (tangle < opt_angle)) {
+                    opt_j = j;
+                    opt_side = (vert->edges[j]->v1 == v)?(1):(-1);
+                    opt_angle = tangle;
+                    opt_flag = 0;
+                }
+            }
+        }
+        
+//        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) {
+                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;
+        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;
+}
+
+/*
+* This function extracts the outer contour of a (self crossing) line.
+* It can generate left/right contour if none of the line ends are in a loop.
+* If one or both of them is in a loop, then there's only one contour
+* 
+* 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
+* 
+* TODO: Implement side != 0 feature;
+*/
+void 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;
+
+    G_debug(3, "extract_outer_contour()");
+
+    if (side != 0) {
+        G_fatal_error("    side != 0 feature not implemented");
+        return;
+    }
+    
+    /* find a line segment which is on the outer contour */
+    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;
+        }
+    }
+    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;
+        }
+    }
+    
+    /* the winding on the outer contour is 0 */
+    extract_contour(pg, edge, (edge->v1 == v)?RIGHT_SIDE:LEFT_SIDE, 0, 0, nPoints);
+    
+    return;
+}
+
+/*
+* 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, 1
+*/
+int extract_inner_contour(struct planar_graph *pg, int *winding, struct line_pnts *nPoints) {
+    int i, w;
+    struct pg_edge *edge;
+    
+    G_debug(3, "extract_inner_contour()");
+
+    for (i = 0; i < pg->ecount; i++) {
+        edge = &(pg->e[i]);
+        if (edge->visited_left) {
+            if (!(pg->e[i].visited_right)) {
+                w = edge->winding_left - 1;
+                extract_contour(pg, &(pg->e[i]), RIGHT_SIDE, w, 0, nPoints);
+                *winding = w;
+                return 1;
+            }
+        }
+        else {
+            if (pg->e[i].visited_right) {
+                w = edge->winding_right + 1;
+                extract_contour(pg, &(pg->e[i]), LEFT_SIDE, w, 0, nPoints);
+                *winding = w;
+                return 1;
+            }
+        }
+    }
+    
+    return 0;
+}
+
+/* point_in_buf - test if point px,py is in d buffer of Points
+** dalpha is in degrees
+** returns:  1 in buffer
+**           0 not in buffer
+*/
+int point_in_buf(struct line_pnts *Points, double px, double py, double da, double db, double dalpha) {
+    int i, np;
+    double cx, cy;
+    double delta, delta_k, k;
+    double vx, vy, wx, wy, mx, my, nx, ny;
+    double len, tx, ty, d, da2;
+    
+    G_debug(3, "point_in_buf()");
+    
+    dalpha *= PI/180; /* convert dalpha from degrees to radians */
+    
+    np = Points->n_points;
+    da2 = da*da;
+    for (i = 0; i < np-1; i++) {
+        vx = Points->x[i];
+        vy = Points->y[i];
+        wx = Points->x[i+1];
+        wy = Points->y[i+1];
+        
+        if (da != db) {
+            mx = wx - vx;
+            my = wy - vy;
+            len = LENGTH(mx, my);
+            elliptic_tangent(mx/len, my/len, da, db, dalpha, &cx, &cy);
+            
+            delta = mx*cy - my*cx;
+            delta_k = (px-vx)*cy - (py-vy)*cx;
+            k = delta_k/delta;
+/*            G_debug(4, "k = %g, k1 = %g", k, (mx * (px - vx) + my * (py - vy)) / (mx * mx + my * my)); */
+            if (k <= 0) {
+                nx = vx;
+                ny = vy;
+            }
+            else if (k >= 1) {
+                nx = wx;
+                ny = wy;
+            }
+            else {
+                nx = vx + k*mx;
+                ny = vy + k*my;
+            }
+            
+            /* inverse transform */
+            elliptic_transform(px - nx, py - ny, 1/da, 1/db, dalpha, &tx, &ty);
+            
+            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, "sqrt(d)*da = %g, len' = %g, olen = %g", sqrt(d)*da, da*LENGTH(tx,ty), LENGTH((px-nx),(py-ny)));*/
+            if (d <= 1) {
+                //G_debug(1, "d=%g", d);
+                return 1;
+            }
+        }
+        else { 
+            d = dig_distance2_point_to_line(px, py, 0, vx, vy, 0, wx, wy, 0,
+                0, NULL, NULL, NULL, NULL, NULL);
+/*            G_debug(4, "sqrt(d)     = %g", sqrt(d));*/
+            if (d <= da2) {
+                return 1;
+            }
+        }
+    }
+    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 (area > 0);
+}
+
+/* internal */
+void add_line_to_array(struct line_pnts *Points, struct line_pnts ***arrPoints, int *count, int *allocated, int more) {
+    if (*allocated == *count) {
+        *allocated += more;
+        *arrPoints = G_realloc(*arrPoints, (*allocated)*sizeof(struct line_pnts *));
+    }
+    (*arrPoints)[*count] = Points;
+    (*count)++;
+    return;
+}
+
+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(3, "buffer_lines()");
+    
+    auto_side = (side == 0);
+    
+    /* initializations */
+    sPoints = Vect_new_line_struct();
+    cPoints = Vect_new_line_struct();
+    arrPoints = NULL;
+    
+    /* outer contour */
+    G_debug(3, "    processing outer contour");
+    *oPoints = Vect_new_line_struct();
+    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);
+    while (res != 0) {
+        if (winding == 0) {
+            if (!Vect_point_in_poly(cPoints->x[0], cPoints->y[0], area_outer)) {
+                if (Vect_get_point_in_poly(cPoints, &px, &py) != 0)
+                    G_fatal_error("Vect_get_point_in_poly() failed.");
+                if (!point_in_buf(area_outer, px, py, da, db, dalpha)) {
+                    add_line_to_array(cPoints, &arrPoints, &count, &allocated, more);
+                    cPoints = Vect_new_line_struct();
+                }
+            }
+        }
+        res = extract_inner_contour(pg2, &winding, cPoints);
+    }
+    pg_destroy_struct(pg2);
+    
+    /* inner contours */
+    G_debug(3, "    processing inner contours");
+    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'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], 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(area_isles[i], px, py, da, db, dalpha)) {
+                        add_line_to_array(cPoints, &arrPoints, &count, &allocated, more);
+                        cPoints = Vect_new_line_struct();
+                    }
+                }
+            }
+            res = extract_inner_contour(pg2, &winding, cPoints);
+        }
+        pg_destroy_struct(pg2);
+    }
+
+    arrPoints = G_realloc(arrPoints, count*sizeof(struct line_pnts *));
+    *inner_count = count;
+    *iPoints = arrPoints;
+    
+    Vect_destroy_line_struct(sPoints);
+    Vect_destroy_line_struct(cPoints);
+    
+    G_debug(3, "buffer_lines() ... done");
+    
+    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(2, "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, isle;
+    int more = 8;
+    int isles_allocated = 0;
+    
+    G_debug(2, "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 **oPoints) {
+  \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 **oPoints) {
+    double tx, ty;
+    double angular_tol, angular_step, phi1;
+    int j, nsegments;
+    
+    G_debug(2, "Vect_point_buffer()");
+    
+    *oPoints = Vect_new_line_struct();
+    
+    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(*oPoints, px + tx, py + ty, 0);
+            phi1 += angular_step;
+        }
+    }
+    else {
+        
+    }
+    
+    /* close the output line */
+    Vect_append_point(*oPoints, (*oPoints)->x[0], (*oPoints)->y[0], (*oPoints)->z[0]);    
+    
+    return;
+}
+
+
+/*
+  \fn void Vect_line_parallel2 ( struct line_pnts *InPoints, double distance, double tolerance, int rm_end,
+                       struct line_pnts *OutPoints )
+  \brief Create parrallel line
+  \param InPoints input line
+  \param distance
+  \param tolerance maximum distance between theoretical arc and polygon segments
+  \param rm_end remove end points falling into distance
+  \param OutPoints output line
+*/
+void Vect_line_parallel2(struct line_pnts *InPoints, double da, double db, double dalpha, int side, int round, double tol, struct line_pnts *OutPoints )
+{
+    G_debug(2, "Vect_line_parallel(): npoints = %d, da = %f, db = %f, dalpha = %f, side = %d, round_corners = %d, tol = %f",
+            InPoints->n_points, da, db, dalpha, side, round, tol);
+
+    parallel_line(InPoints, da, db, dalpha, side, round, 1, NON_LOOPED_LINE, tol, OutPoints);
+    
+/*    if (!loops)
+        clean_parallel(OutPoints, InPoints, distance, rm_end);
+*/
+    return;
+}


Property changes on: grass/branches/develbranch_6/vector/v.buffer2/vlib_buffer.c
___________________________________________________________________
Name: svn:eol-style
   + native



More information about the grass-commit mailing list