[GRASS-SVN] r34601 - in grass/branches/develbranch_6: include
lib/vector/Vlib vector/v.buffer2 vector/v.parallel2
svn_grass at osgeo.org
svn_grass at osgeo.org
Fri Nov 28 17:09:54 EST 2008
Author: martinl
Date: 2008-11-28 17:09:53 -0500 (Fri, 28 Nov 2008)
New Revision: 34601
Added:
grass/branches/develbranch_6/lib/vector/Vlib/buffer2.c
grass/branches/develbranch_6/lib/vector/Vlib/dgraph.c
grass/branches/develbranch_6/lib/vector/Vlib/dgraph.h
grass/branches/develbranch_6/lib/vector/Vlib/e_intersect.c
grass/branches/develbranch_6/lib/vector/Vlib/e_intersect.h
Removed:
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/vlib_buffer.c
grass/branches/develbranch_6/vector/v.parallel2/dgraph.c
grass/branches/develbranch_6/vector/v.parallel2/dgraph.h
grass/branches/develbranch_6/vector/v.parallel2/e_intersect.c
grass/branches/develbranch_6/vector/v.parallel2/e_intersect.h
grass/branches/develbranch_6/vector/v.parallel2/vlib_buffer.c
Modified:
grass/branches/develbranch_6/include/Vect.h
grass/branches/develbranch_6/lib/vector/Vlib/buffer.c
Log:
Vlib: Vect_line_buffe2(), Vect_area_buffer2(), Vect_point_buffer2() and
Vect_line_parallel() moved from v.buffer2 and v.parallel2 to
vector library (buffer2.c)
Modified: grass/branches/develbranch_6/include/Vect.h
===================================================================
--- grass/branches/develbranch_6/include/Vect.h 2008-11-28 17:55:19 UTC (rev 34600)
+++ grass/branches/develbranch_6/include/Vect.h 2008-11-28 22:09:53 UTC (rev 34601)
@@ -35,8 +35,23 @@
int Vect_line_box(struct line_pnts *, BOUND_BOX *);
void Vect_line_parallel(struct line_pnts *, double, double, int,
struct line_pnts *);
+void Vect_line_parallel2(struct line_pnts *, double, double,
+ double, int, int, double,
+ struct line_pnts *);
void Vect_line_buffer(struct line_pnts *, double, double, struct line_pnts *);
+void Vect_line_buffer2(struct line_pnts *, double, double,
+ double, int, int, double,
+ struct line_pnts **,
+ struct line_pnts ***, int *);
+void Vect_area_buffer2(struct Map_info *, int, double, double,
+ double, int, int, double,
+ struct line_pnts **,
+ struct line_pnts ***, int *);
+void Vect_point_buffer2(double, double, double, double,
+ double, int, double,
+ struct line_pnts **);
+
/* Categories */
struct line_cats *Vect_new_cats_struct(void);
int Vect_cat_set(struct line_cats *, int, int);
Modified: grass/branches/develbranch_6/lib/vector/Vlib/buffer.c
===================================================================
--- grass/branches/develbranch_6/lib/vector/Vlib/buffer.c 2008-11-28 17:55:19 UTC (rev 34600)
+++ grass/branches/develbranch_6/lib/vector/Vlib/buffer.c 2008-11-28 22:09:53 UTC (rev 34601)
@@ -1,7 +1,7 @@
/*!
\file buffer.c
- \brief Vector library - nearest, adjust_line, parallel_line
+ \brief Vector library - nearest, adjust, parallel lines
Higher level functions for reading/writing/manipulating vectors.
Added: grass/branches/develbranch_6/lib/vector/Vlib/buffer2.c
===================================================================
--- grass/branches/develbranch_6/lib/vector/Vlib/buffer2.c (rev 0)
+++ grass/branches/develbranch_6/lib/vector/Vlib/buffer2.c 2008-11-28 22:09:53 UTC (rev 34601)
@@ -0,0 +1,1203 @@
+/*!
+ \file buffer2.c
+
+ \brief Vector library - nearest, adjust, parallel lines
+
+ Higher level functions for reading/writing/manipulating vectors.
+
+ (C) 2001-2008 by the GRASS Development Team
+
+ This program is free software under the
+ GNU General Public License (>=v2).
+ Read the file COPYING that comes with GRASS
+ for details.
+
+ \author Original author Radim Blazek (see buffer.c)
+ \author Rewritten by Rosen Matev (Google Summer of Code 2008)
+
+ \date 2008
+ */
+
+#include <stdlib.h>
+#include <grass/Vect.h>
+#include <grass/gis.h>
+#include <grass/glocale.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)
+ */
+static 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 */
+static 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;
+ */
+static 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
+ */
+static 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
+ */
+static 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
+ */
+static 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 */
+static 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;
+}
+
+static 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
+ */
+static 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;
+}
+
+
+/*!
+ \brief Creates buffer around line.
+
+ See also Vect_line_buffer().
+
+ \param InPoints input line geometry
+ \param da distance along major axis
+ \param db 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[out] oPoints output polygon outer border (ccw order)
+ \param[out] inner_count number of holes
+ \param[out] 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;
+}
+
+/*!
+ \brief Creates buffer around area.
+
+ \param Map vector map
+ \param area area id
+ \param da distance along major axis
+ \param db 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[out] oPoints output polygon outer border (ccw order)
+ \param[out] inner_count number of holes
+ \param[out] iPoints array of output polygon's holes (cw order)
+ */
+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;
+}
+
+/*!
+ \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[out] 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;
+}
+
+
+/*
+ \brief Create parrallel line
+
+ See also Vect_line_parallel().
+
+ \param InPoints input line geometry
+ \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[out] 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;
+}
Copied: grass/branches/develbranch_6/lib/vector/Vlib/dgraph.c (from rev 34598, grass/branches/develbranch_6/vector/v.buffer2/dgraph.c)
===================================================================
--- grass/branches/develbranch_6/lib/vector/Vlib/dgraph.c (rev 0)
+++ grass/branches/develbranch_6/lib/vector/Vlib/dgraph.c 2008-11-28 22:09:53 UTC (rev 34601)
@@ -0,0 +1,521 @@
+#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/lib/vector/Vlib/dgraph.c
___________________________________________________________________
Name: svn:mergeinfo
+
Name: svn:eol-style
+ native
Copied: grass/branches/develbranch_6/lib/vector/Vlib/dgraph.h (from rev 34598, grass/branches/develbranch_6/vector/v.buffer2/dgraph.h)
===================================================================
--- grass/branches/develbranch_6/lib/vector/Vlib/dgraph.h (rev 0)
+++ grass/branches/develbranch_6/lib/vector/Vlib/dgraph.h 2008-11-28 22:09:53 UTC (rev 34601)
@@ -0,0 +1,38 @@
+#ifndef GRASS_DGRAPH_H
+#define GRASS_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/lib/vector/Vlib/dgraph.h
___________________________________________________________________
Name: svn:mergeinfo
+
Name: svn:eol-style
+ native
Copied: grass/branches/develbranch_6/lib/vector/Vlib/e_intersect.c (from rev 34598, grass/branches/develbranch_6/vector/v.buffer2/e_intersect.c)
===================================================================
--- grass/branches/develbranch_6/lib/vector/Vlib/e_intersect.c (rev 0)
+++ grass/branches/develbranch_6/lib/vector/Vlib/e_intersect.c 2008-11-28 22:09:53 UTC (rev 34601)
@@ -0,0 +1,1077 @@
+#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(3, "segment_intersection_2d_test()");
+ G_debug(3, " ax1 = %.18e, ay1 = %.18e", ax1, ay1);
+ G_debug(3, " ax2 = %.18e, ay2 = %.18e", ax2, ay2);
+ G_debug(3, " bx1 = %.18e, by1 = %.18e", bx1, by1);
+ G_debug(3, " 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(4, " identical segments");
+ *x1 = ax1;
+ *y1 = ay1;
+ *x2 = ax2;
+ *y2 = ay2;
+ return 5;
+ }
+ /* Check for identical endpoints */
+ if (f11 || f12) {
+ G_debug(4, " connected by endpoints");
+ *x1 = ax1;
+ *y1 = ay1;
+ return 1;
+ }
+ if (f21 || f22) {
+ G_debug(4, " connected by endpoints");
+ *x1 = ax2;
+ *y1 = ay2;
+ return 1;
+ }
+
+ if ((MAX(ax1, ax2) < MIN(bx1, bx2)) || (MAX(bx1, bx2) < MIN(ax1, ax2))) {
+ G_debug(4, " no intersection (disjoint bounding boxes)");
+ return 0;
+ }
+ if ((MAX(ay1, ay2) < MIN(by1, by2)) || (MAX(by1, by2) < MIN(ay1, ay2))) {
+ G_debug(4, " 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(3, " dd = %sE%d", (s[0] == 0) ? "0" : s, exp);
+ G_debug(3, " d = %.18E", d);
+
+ if (sgn_d != 0) {
+ G_debug(3, " 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(4, " rra = %sE%d", (s[0] == 0) ? "0" : s, exp);
+ G_debug(4, " ra = %.18E", ra);
+ s = mpf_get_str(NULL, &exp, 10, 40, rrb);
+ G_debug(4, " rrb = %sE%d", (s[0] == 0) ? "0" : s, exp);
+ G_debug(4, " 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(2, " intersection at:");
+ G_debug(2, " xx = %.18e", *x1);
+ G_debug(2, " x = %.18e", ax1 + ra * (ax2 - ax1));
+ G_debug(2, " yy = %.18e", *y1);
+ G_debug(2, " y = %.18e", ay1 + ra * (ay2 - ay1));
+ return 1;
+ }
+
+ G_debug(3, " parallel/collinear...");
+ return -1;
+}
+#endif
Property changes on: grass/branches/develbranch_6/lib/vector/Vlib/e_intersect.c
___________________________________________________________________
Name: svn:mergeinfo
+
Name: svn:eol-style
+ native
Copied: grass/branches/develbranch_6/lib/vector/Vlib/e_intersect.h (from rev 34598, grass/branches/develbranch_6/vector/v.buffer2/e_intersect.h)
===================================================================
--- grass/branches/develbranch_6/lib/vector/Vlib/e_intersect.h (rev 0)
+++ grass/branches/develbranch_6/lib/vector/Vlib/e_intersect.h 2008-11-28 22:09:53 UTC (rev 34601)
@@ -0,0 +1,21 @@
+#ifndef GRASS_E_INTERSECT_H
+#define GRASS_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/lib/vector/Vlib/e_intersect.h
___________________________________________________________________
Name: svn:mergeinfo
+
Name: svn:eol-style
+ native
Deleted: grass/branches/develbranch_6/vector/v.buffer2/dgraph.c
===================================================================
--- grass/branches/develbranch_6/vector/v.buffer2/dgraph.c 2008-11-28 17:55:19 UTC (rev 34600)
+++ grass/branches/develbranch_6/vector/v.buffer2/dgraph.c 2008-11-28 22:09:53 UTC (rev 34601)
@@ -1,521 +0,0 @@
-#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;
-}
Deleted: grass/branches/develbranch_6/vector/v.buffer2/dgraph.h
===================================================================
--- grass/branches/develbranch_6/vector/v.buffer2/dgraph.h 2008-11-28 17:55:19 UTC (rev 34600)
+++ grass/branches/develbranch_6/vector/v.buffer2/dgraph.h 2008-11-28 22:09:53 UTC (rev 34601)
@@ -1,38 +0,0 @@
-#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
Deleted: grass/branches/develbranch_6/vector/v.buffer2/e_intersect.c
===================================================================
--- grass/branches/develbranch_6/vector/v.buffer2/e_intersect.c 2008-11-28 17:55:19 UTC (rev 34600)
+++ grass/branches/develbranch_6/vector/v.buffer2/e_intersect.c 2008-11-28 22:09:53 UTC (rev 34601)
@@ -1,1077 +0,0 @@
-#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(3, "segment_intersection_2d_test()");
- G_debug(3, " ax1 = %.18e, ay1 = %.18e", ax1, ay1);
- G_debug(3, " ax2 = %.18e, ay2 = %.18e", ax2, ay2);
- G_debug(3, " bx1 = %.18e, by1 = %.18e", bx1, by1);
- G_debug(3, " 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(4, " identical segments");
- *x1 = ax1;
- *y1 = ay1;
- *x2 = ax2;
- *y2 = ay2;
- return 5;
- }
- /* Check for identical endpoints */
- if (f11 || f12) {
- G_debug(4, " connected by endpoints");
- *x1 = ax1;
- *y1 = ay1;
- return 1;
- }
- if (f21 || f22) {
- G_debug(4, " connected by endpoints");
- *x1 = ax2;
- *y1 = ay2;
- return 1;
- }
-
- if ((MAX(ax1, ax2) < MIN(bx1, bx2)) || (MAX(bx1, bx2) < MIN(ax1, ax2))) {
- G_debug(4, " no intersection (disjoint bounding boxes)");
- return 0;
- }
- if ((MAX(ay1, ay2) < MIN(by1, by2)) || (MAX(by1, by2) < MIN(ay1, ay2))) {
- G_debug(4, " 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(3, " dd = %sE%d", (s[0] == 0) ? "0" : s, exp);
- G_debug(3, " d = %.18E", d);
-
- if (sgn_d != 0) {
- G_debug(3, " 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(4, " rra = %sE%d", (s[0] == 0) ? "0" : s, exp);
- G_debug(4, " ra = %.18E", ra);
- s = mpf_get_str(NULL, &exp, 10, 40, rrb);
- G_debug(4, " rrb = %sE%d", (s[0] == 0) ? "0" : s, exp);
- G_debug(4, " 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(2, " intersection at:");
- G_debug(2, " xx = %.18e", *x1);
- G_debug(2, " x = %.18e", ax1 + ra * (ax2 - ax1));
- G_debug(2, " yy = %.18e", *y1);
- G_debug(2, " y = %.18e", ay1 + ra * (ay2 - ay1));
- return 1;
- }
-
- G_debug(3, " parallel/collinear...");
- return -1;
-}
-#endif
Deleted: grass/branches/develbranch_6/vector/v.buffer2/e_intersect.h
===================================================================
--- grass/branches/develbranch_6/vector/v.buffer2/e_intersect.h 2008-11-28 17:55:19 UTC (rev 34600)
+++ grass/branches/develbranch_6/vector/v.buffer2/e_intersect.h 2008-11-28 22:09:53 UTC (rev 34601)
@@ -1,21 +0,0 @@
-#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
Deleted: grass/branches/develbranch_6/vector/v.buffer2/vlib_buffer.c
===================================================================
--- grass/branches/develbranch_6/vector/v.buffer2/vlib_buffer.c 2008-11-28 17:55:19 UTC (rev 34600)
+++ grass/branches/develbranch_6/vector/v.buffer2/vlib_buffer.c 2008-11-28 22:09:53 UTC (rev 34601)
@@ -1,1168 +0,0 @@
-/* Functions: ...
- **
- ** Author: Radim Blazek, Rosen Matev; July 2008
- **
- **
- */
-#include <stdlib.h>
-#include <grass/Vect.h>
-#include <grass/gis.h>
-#include <grass/glocale.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;
-}
Deleted: grass/branches/develbranch_6/vector/v.parallel2/dgraph.c
===================================================================
--- grass/branches/develbranch_6/vector/v.parallel2/dgraph.c 2008-11-28 17:55:19 UTC (rev 34600)
+++ grass/branches/develbranch_6/vector/v.parallel2/dgraph.c 2008-11-28 22:09:53 UTC (rev 34601)
@@ -1,521 +0,0 @@
-#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;
-}
Deleted: grass/branches/develbranch_6/vector/v.parallel2/dgraph.h
===================================================================
--- grass/branches/develbranch_6/vector/v.parallel2/dgraph.h 2008-11-28 17:55:19 UTC (rev 34600)
+++ grass/branches/develbranch_6/vector/v.parallel2/dgraph.h 2008-11-28 22:09:53 UTC (rev 34601)
@@ -1,38 +0,0 @@
-#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
Deleted: grass/branches/develbranch_6/vector/v.parallel2/e_intersect.c
===================================================================
--- grass/branches/develbranch_6/vector/v.parallel2/e_intersect.c 2008-11-28 17:55:19 UTC (rev 34600)
+++ grass/branches/develbranch_6/vector/v.parallel2/e_intersect.c 2008-11-28 22:09:53 UTC (rev 34601)
@@ -1,1077 +0,0 @@
-#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(3, "segment_intersection_2d_test()");
- G_debug(3, " ax1 = %.18e, ay1 = %.18e", ax1, ay1);
- G_debug(3, " ax2 = %.18e, ay2 = %.18e", ax2, ay2);
- G_debug(3, " bx1 = %.18e, by1 = %.18e", bx1, by1);
- G_debug(3, " 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(4, " identical segments");
- *x1 = ax1;
- *y1 = ay1;
- *x2 = ax2;
- *y2 = ay2;
- return 5;
- }
- /* Check for identical endpoints */
- if (f11 || f12) {
- G_debug(4, " connected by endpoints");
- *x1 = ax1;
- *y1 = ay1;
- return 1;
- }
- if (f21 || f22) {
- G_debug(4, " connected by endpoints");
- *x1 = ax2;
- *y1 = ay2;
- return 1;
- }
-
- if ((MAX(ax1, ax2) < MIN(bx1, bx2)) || (MAX(bx1, bx2) < MIN(ax1, ax2))) {
- G_debug(4, " no intersection (disjoint bounding boxes)");
- return 0;
- }
- if ((MAX(ay1, ay2) < MIN(by1, by2)) || (MAX(by1, by2) < MIN(ay1, ay2))) {
- G_debug(4, " 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(3, " dd = %sE%d", (s[0] == 0) ? "0" : s, exp);
- G_debug(3, " d = %.18E", d);
-
- if (sgn_d != 0) {
- G_debug(3, " 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(3, " rra = %sE%d", (s[0] == 0) ? "0" : s, exp);
- G_debug(3, " ra = %.18E", ra);
- s = mpf_get_str(NULL, &exp, 10, 40, rrb);
- G_debug(3, " rrb = %sE%d", (s[0] == 0) ? "0" : s, exp);
- G_debug(3, " 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(3, " intersection at:");
- G_debug(3, " xx = %.18e", *x1);
- G_debug(3, " x = %.18e", ax1 + ra * (ax2 - ax1));
- G_debug(3, " yy = %.18e", *y1);
- G_debug(3, " y = %.18e", ay1 + ra * (ay2 - ay1));
- return 1;
- }
-
- G_debug(3, " parallel/collinear...");
- return -1;
-}
-#endif
Deleted: grass/branches/develbranch_6/vector/v.parallel2/e_intersect.h
===================================================================
--- grass/branches/develbranch_6/vector/v.parallel2/e_intersect.h 2008-11-28 17:55:19 UTC (rev 34600)
+++ grass/branches/develbranch_6/vector/v.parallel2/e_intersect.h 2008-11-28 22:09:53 UTC (rev 34601)
@@ -1,21 +0,0 @@
-#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
Deleted: grass/branches/develbranch_6/vector/v.parallel2/vlib_buffer.c
===================================================================
--- grass/branches/develbranch_6/vector/v.parallel2/vlib_buffer.c 2008-11-28 17:55:19 UTC (rev 34600)
+++ grass/branches/develbranch_6/vector/v.parallel2/vlib_buffer.c 2008-11-28 22:09:53 UTC (rev 34601)
@@ -1,1167 +0,0 @@
-/* Functions: ...
- **
- ** Author: Radim Blazek, Rosen Matev; July 2008
- **
- **
- */
-#include <stdlib.h>
-#include <grass/Vect.h>
-#include <grass/gis.h>
-#include <grass/glocale.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;
-}
More information about the grass-commit
mailing list