www.pudn.com > geosteiner-3.1.zip > efuncs.h
/***********************************************************************
File: efuncs.h
Rev: a-1
Date: 01/22/00
Copyright (c) 2001 by Pawel Winter, Martin Zachariasen
& David M. Warme
************************************************************************
Basic plane geometry functions used by Euclidean FST generator
and pruning algorithms.
This file is merely intended for inlining these speed-critical
functions in various c-files.
************************************************************************
Modification Log:
a-1: 01/22/00 martinz
: Created. Split off from efst.c
************************************************************************/
#ifndef EFUNCS_H
#define EFUNCS_H
#include "steiner.h"
/*
* Local Constants
*/
#define PI (3.1415926535897932384626433832795028841971693)
#define RAD120 (2.0*PI/3.0)
#define SQRT3 (1.73205080756887729352)
#define HALF_SQRT3 (0.86602540378443864676)
/*
* Local Macros
*/
#undef MIN
#define MIN(a,b) (((a) < (b)) ? (a) : (b))
#undef MAX
#define MAX(a,b) (((a) > (b)) ? (a) : (b))
/*
* Does the sequence p1, p2, p3 make a right turn?
*/
static
bool
right_turn (
struct point * p1, /* IN - first point */
struct point * p2, /* IN - second point */
struct point * p3 /* IN - third point */
)
{
return ( (p1 -> x - p2 -> x) * (p1 -> y - p3 -> y) <
(p1 -> y - p2 -> y) * (p1 -> x - p3 -> x) );
}
/*
* Does the sequence p1, p2, p3 make a left turn?
*/
static
bool
left_turn (
struct point * p1, /* IN - first point */
struct point * p2, /* IN - second point */
struct point * p3 /* IN - third point */
)
{
return ( (p1 -> x - p2 -> x) * (p1 -> y - p3 -> y) >
(p1 -> y - p2 -> y) * (p1 -> x - p3 -> x) );
}
/*
* Square of distance between points p1 and p2
*/
static
dist_t
sqr_dist (
struct point * p1, /* IN - first point */
struct point * p2 /* IN - second point */
)
{
double dx, dy;
dx = p1 -> x - p2 -> x;
dy = p1 -> y - p2 -> y;
return (dx*dx + dy*dy);
}
/*
* Rotate vector 60 degrees counter-clockwise
*/
static
void
rotate60 (
double dx, /* IN - x-coordinate of vector */
double dy, /* IN - y-coordinate of vector */
double * rdx, /* OUT - x-coordinate of rotated vector */
double * rdy /* OUT - y-coordinate of rotated vector */
)
{
*rdx = 0.5*dx - HALF_SQRT3*dy;
*rdy = 0.5*dy + HALF_SQRT3*dx;
}
/*
* Find equilateral point for points p1 and p2
*/
static
void
eq_point (
struct point * p1, /* IN - first point */
struct point * p2, /* IN - second point */
struct point * e /* OUT - equilateral point */
)
{
double dx, dy;
dx = p1 -> x - p2 -> x;
dy = p1 -> y - p2 -> y;
e -> x = p2 -> x + 0.5*dx - HALF_SQRT3*dy;
e -> y = p2 -> y + 0.5*dy + HALF_SQRT3*dx;
}
/*
* Find equilateral point for points p1 and p2 (alternative version)
*/
#if 0
static
void
eq_point (
struct point * p1, /* IN - first point */
struct point * p2, /* IN - second point */
struct point * e /* OUT - equilateral point */
)
{
e -> x = 0.5*(p1 -> x + p2 -> x) - HALF_SQRT3*(p1 -> y - p2 -> y);
e -> y = 0.5*(p1 -> y + p2 -> y) + HALF_SQRT3*(p1 -> x - p2 -> x);
}
#endif
/*
* Find center of equilateral circle
*/
static
void
eq_circle_center (
struct point * p1, /* IN - first point on circle */
struct point * p2, /* IN - second second on circle */
struct point * e, /* IN - third point on circle */
struct point * c /* OUT - center of circle */
)
{
c -> x = (p1 -> x + p2 -> x + e -> x) / 3.0;
c -> y = (p1 -> y + p2 -> y + e -> y) / 3.0;
}
/*
* Compute vector corresponding to angle between points a, b and c
*/
static
void
get_angle_vector (
struct point * a, /* IN - first point */
struct point * b, /* IN - second point */
struct point * c, /* IN - third point */
double * dx, /* OUT - x-coordinate of vector */
double * dy /* OUT - y-coordinate of vector */
)
{
double abx, aby, cbx, cby;
abx = a -> x - b -> x;
aby = a -> y - b -> y;
cbx = c -> x - b -> x;
cby = c -> y - b -> y;
*dx = abx*cbx + aby*cby;
*dy = abx*cby - cbx*aby;
}
/*
* Rotate point P around C at an angle v whose sin(v)
* and cos(v) are given.
*/
static
void
rotate (
struct point * P, /* IN - point to be rotated */
struct point * C, /* IN - center of rotation */
double sinv, /* IN - sin(v) of angle */
double cosv, /* IN - cos(v) of angle */
struct point * PN /* OUT - rotated point */
)
{
double ox, oy, dx, dy;
ox = C -> x;
oy = C -> y;
dx = P -> x - ox;
dy = P -> y - oy;
PN -> x = ox + dx*cosv - dy*sinv;
PN -> y = oy + dx*sinv + dy*cosv;
}
/*
* Rotate point P around C by an angle 2v, where
* sin(v) and cos(v) are given.
*/
static
void
rotate2 (
struct point * P, /* IN - point to be rotaed */
struct point * C, /* IN - center of rotation */
double sinv, /* IN - sin(v) of angle */
double cosv, /* IN - cos(v) of angle */
struct point * PN /* OUT - rotated point */
)
{
double cosv2, sin2v, cos2v, ox, oy, dx, dy;
if (cosv EQ 2.0) { /* not defined */
cosv2 = 1.0 - sinv*sinv;
cosv = sqrt(cosv2);
}
else {
cosv2 = cosv*cosv;
}
sin2v = 2.0*sinv*cosv;
cos2v = 2.0*cosv2 - 1.0;
ox = C -> x;
oy = C -> y;
dx = P -> x - ox;
dy = P -> y - oy;
PN -> x = ox + dx*cos2v - dy*sin2v;
PN -> y = oy + dx*sin2v + dy*cos2v;
}
/*
* Given a point E on a circle with center C the procedure returns
* the projection EN of E onto the same circle in direction ER.
*/
static
void
project_point (
struct point * E, /* IN - point to be projected */
struct point * C, /* IN - center of circle */
struct point * R, /* IN - direction of projection */
struct point * EN /* OUT - projected point */
)
{
double Ax, Ay, Bx, By, AdotB, BdotB, lambda;
Ax = C -> x - E -> x;
Ay = C -> y - E -> y;
Bx = R -> x - E -> x;
By = R -> y - E -> y;
AdotB = Ax * Bx + Ay * By;
BdotB = Bx * Bx + By * By;
if (BdotB EQ 0.0) {
*EN = *E;
EN -> pnum = -1;
}
else {
lambda = AdotB / BdotB;
lambda += lambda; /* multiply by 2 */
EN -> x = E -> x + lambda * Bx;
EN -> y = E -> y + lambda * By;
EN -> pnum = -1;
}
}
/*
* Given a point E on a circle with center C the procedure returns
* the projection of E onto the same circle in a direction perpendicular
* to CR.
*/
static
void
project_point_perp (
struct point * E, /* IN - point to be projected */
struct point * C, /* IN - center of circle */
struct point * R, /* IN - perpendicular direction of projection */
struct point * EN /* OUT - projected point */
)
{
double Ax, Ay, Bx, By, AdotB, BdotB, lambda;
Ax = C -> x - E -> x;
Ay = C -> y - E -> y;
Bx = C -> y - R -> y;
By = R -> x - C -> x;
AdotB = Ax * Bx + Ay * By;
BdotB = Bx * Bx + By * By;
if (BdotB EQ 0.0) {
*EN = *E;
EN -> pnum = -1;
}
else {
lambda = AdotB / BdotB;
lambda += lambda; /* multiply by 2 */
EN -> x = E -> x + lambda * Bx;
EN -> y = E -> y + lambda * By;
EN -> pnum = -1;
}
}
/*
* Test if vector corresponds to an angle greater than 120 degrees.
* Special fast version which avoids using trigonometric functions.
*/
static
bool
angle_geq_120 (
double dx, /* IN - x-coordinate of vector */
double dy /* IN - y-coordinate of vector */
)
{
if (dy < 0.0) return TRUE;
if ((dy EQ 0.0) AND (dx <= 0.0)) return TRUE;
if ((dy > 0.0) AND (dx < 0.0) AND ( dy <= -SQRT3*dx)) return TRUE;
return FALSE;
}
/*
* Test if vector corresponds to an angle greater than 60 degrees.
* Special fast version which avoids using trigonometric functions.
*/
static
bool
angle_geq_60 (
double dx, /* IN - x-coordinate of vector */
double dy /* IN - y-coordinate of vector */
)
{
if ((dx < 0.0) OR (dy < 0.0)) return TRUE;
if ((dy EQ 0.0) AND (dx <= 0.0)) return TRUE;
if ((dy > 0.0) AND (dy >= SQRT3*dx)) return TRUE;
return FALSE;
}
/*
* Test if the points p1,p2,p3 make up a wedge
* that is greater than 120 degrees
*/
static
bool
wedge120 (
struct point * p1, /* IN - first point */
struct point * p2, /* IN - second point */
struct point * p3 /* IN - third point */
)
{
struct point e;
struct point c;
if (left_turn (p1,p3,p2)) {
eq_point (p1, p3, &e);
eq_circle_center (p1, p3, &e, &c);
}
else {
eq_point (p3, p1, &e);
eq_circle_center (p3, p1, &e, &c);
}
if (sqr_dist (&c, p2) <= sqr_dist (&c, p1)) {
/* p2 is inside circle -> angle greater than 120 deg */
return (TRUE);
}
/* less than 120 deg */
return (FALSE);
}
/*
* Intersection between two line segments pq and rs.
* Adapted from Dave's more generic intersection code.
*/
static
bool
segment_intersection (
struct point * p, /* IN - first endpoint of first segment */
struct point * q, /* IN - second endpoint of first segment */
struct point * r, /* IN - first endpoint of second segment */
struct point * s, /* IN - second endpoint of second segment */
struct point * is /* OUT - intersection point */
)
{
/* This routine takes two line segments PQ and RS and determines
their intersection (if any). It makes use of a coordinate
transformation to reduce the number of special cases to a
minimum. Line segment PQ runs from point P to point Q, and
line segment RS runs from point R to point S. We create a
parametric form for line segment PQ as:
XPQ(alpha) = (1 - alpha) * P.x + alpha * Q.x
YPQ(alpha) = (1 - alpha) * P.y + alpha * Q.y
As alpha varies from 0.0 to 1.0, the parametric point follows
line segment PQ from P to Q. Similarly, we parameterize line
segment RS as a function of beta:
XRS(beta) = (1 - beta) * R.x + beta * S.x
YRS(beta) = (1 - beta) * R.y + beta * S.y
The (infinite) lines containing segments PQ and RS intersect when:
[XPQ(alpha) = XRS(beta)] AND [YPQ(alpha) = YRS(beta)]
Note that if (0 <= alpha <= 1) AND (0 <= beta <= 1) then the
line segments actually intersect, not just the lines. This is
the really nifty property of this coordinate transformation!
Expanding out the simultaneous equations gives:
(1 - alpha) * P.x + alpha * Q.x = (1 - beta) * R.x + beta * S.x
(1 - alpha) * P.y + alpha * Q.y = (1 - beta) * R.y + beta * S.y
or:
(Q.x - P.x) * alpha + (R.x - S.x) * beta = (R.x - P.x)
(Q.y - P.y) * alpha + (R.y - S.y) * beta = (R.y - P.y)
which we write as:
A * alpha + B * beta = C
D * alpha + E * beta = F
Unless this system is singular, its solutions are:
B*F - C*E C*D - A*F
alpha = ----------- beta = -----------
B*D - A*E B*D - A*E
The system is singular if either line segment is zero-length, or
if the lines have the same slope (parallel or colinear). These
cases are handled separately.
*/
double alpha, beta;
double A, B, C, D, E, F;
double denom;
struct point * rp;
struct point * sp;
struct point * tmp;
A = q -> x - p -> x;
B = r -> x - s -> x;
C = r -> x - p -> x;
D = q -> y - p -> y;
E = r -> y - s -> y;
F = r -> y - p -> y;
rp = r;
sp = s;
denom = B * D - A * E;
if (denom < 0) {
/* switch points R and S... */
tmp = rp;
rp = sp;
sp = tmp;
denom = - denom;
B = - B;
C = rp -> x - p -> x;
E = - E;
F = rp -> y - p -> y;
}
if (denom NE 0.0) {
/* Easy case -- lines are NOT parallel. */
alpha = (B * F - C * E);
if ((alpha < 0.0) OR (denom < alpha)) {
/* First segment does not reach to second */
/* line segment. */
return (FALSE);
}
beta = C * D - A * F;
if ((beta < 0.0) OR (denom < beta)) {
/* Second line segment does not reach to first */
/* line segment. */
return (FALSE);
}
/* We have a single point of intersection! */
alpha /= denom;
is -> x = (1.0 - alpha) * p -> x + alpha * q -> x;
is -> y = (1.0 - alpha) * p -> y + alpha * q -> y;
return (TRUE);
}
/* The two line segments are parallel.
Report no (unique) intersection point. */
return (FALSE);
}
#endif