3D图形编程

开发平台：

C/C++

triangle.c：源码内容

REAL bheight;
REAL cheight;
REAL dheight;
#endif /* not ANSI_DECLARATORS */
{
REAL adx, bdx, cdx, ady, bdy, cdy, adheight, bdheight, cdheight;
REAL bdxcdy, cdxbdy, cdxady, adxcdy, adxbdy, bdxady;
REAL det;
REAL permanent, errbound;
m->orient3dcount++;
adx = pa[0] - pd[0];
bdx = pb[0] - pd[0];
cdx = pc[0] - pd[0];
ady = pa[1] - pd[1];
bdy = pb[1] - pd[1];
cdy = pc[1] - pd[1];
adheight = aheight - dheight;
bdheight = bheight - dheight;
cdheight = cheight - dheight;
bdxcdy = bdx * cdy;
cdxbdy = cdx * bdy;
cdxady = cdx * ady;
adxcdy = adx * cdy;
adxbdy = adx * bdy;
bdxady = bdx * ady;
det = adheight * (bdxcdy - cdxbdy)
+ bdheight * (cdxady - adxcdy)
+ cdheight * (adxbdy - bdxady);
if (b->noexact) {
return det;
}
permanent = (Absolute(bdxcdy) + Absolute(cdxbdy)) * Absolute(adheight)
+ (Absolute(cdxady) + Absolute(adxcdy)) * Absolute(bdheight)
+ (Absolute(adxbdy) + Absolute(bdxady)) * Absolute(cdheight);
errbound = o3derrboundA * permanent;
if ((det > errbound) || (-det > errbound)) {
return det;
}
return orient3dadapt(pa, pb, pc, pd, aheight, bheight, cheight, dheight,
permanent);
}
/*****************************************************************************/
/* */
/* nonregular() Return a positive value if the point pd is incompatible */
/* with the circle or plane passing through pa, pb, and pc */
/* (meaning that pd is inside the circle or below the */
/* plane); a negative value if it is compatible; and zero if */
/* the four points are cocircular/coplanar. The points pa, */
/* pb, and pc must be in counterclockwise order, or the sign */
/* of the result will be reversed. */
/* */
/* If the -w switch is used, the points are lifted onto the parabolic */
/* lifting map, then they are dropped according to their weights, then the */
/* 3D orientation test is applied. If the -W switch is used, the points' */
/* heights are already provided, so the 3D orientation test is applied */
/* directly. If neither switch is used, the incircle test is applied. */
/* */
/*****************************************************************************/
#ifdef ANSI_DECLARATORS
REAL nonregular(struct mesh *m, struct behavior *b,
vertex pa, vertex pb, vertex pc, vertex pd)
#else /* not ANSI_DECLARATORS */
REAL nonregular(m, b, pa, pb, pc, pd)
struct mesh *m;
struct behavior *b;
vertex pa;
vertex pb;
vertex pc;
vertex pd;
#endif /* not ANSI_DECLARATORS */
{
if (b->weighted == 0) {
return incircle(m, b, pa, pb, pc, pd);
} else if (b->weighted == 1) {
return orient3d(m, b, pa, pb, pc, pd,
pa[0] * pa[0] + pa[1] * pa[1] - pa[2],
pb[0] * pb[0] + pb[1] * pb[1] - pb[2],
pc[0] * pc[0] + pc[1] * pc[1] - pc[2],
pd[0] * pd[0] + pd[1] * pd[1] - pd[2]);
} else {
return orient3d(m, b, pa, pb, pc, pd, pa[2], pb[2], pc[2], pd[2]);
}
}
/*****************************************************************************/
/* */
/* findcircumcenter() Find the circumcenter of a triangle. */
/* */
/* The result is returned both in terms of x-y coordinates and xi-eta */
/* (barycentric) coordinates. The xi-eta coordinate system is defined in */
/* terms of the triangle: the origin of the triangle is the origin of the */
/* coordinate system; the destination of the triangle is one unit along the */
/* xi axis; and the apex of the triangle is one unit along the eta axis. */
/* This procedure also returns the square of the length of the triangle's */
/* shortest edge. */
/* */
/*****************************************************************************/
#ifdef ANSI_DECLARATORS
void findcircumcenter(struct mesh *m, struct behavior *b,
vertex torg, vertex tdest, vertex tapex,
vertex circumcenter, REAL *xi, REAL *eta, int offcenter)
#else /* not ANSI_DECLARATORS */
void findcircumcenter(m, b, torg, tdest, tapex, circumcenter, xi, eta,
offcenter)
struct mesh *m;
struct behavior *b;
vertex torg;
vertex tdest;
vertex tapex;
vertex circumcenter;
REAL *xi;
REAL *eta;
int offcenter;
#endif /* not ANSI_DECLARATORS */
{
REAL xdo, ydo, xao, yao;
REAL dodist, aodist, dadist;
REAL denominator;
REAL dx, dy, dxoff, dyoff;
m->circumcentercount++;
/* Compute the circumcenter of the triangle. */
xdo = tdest[0] - torg[0];
ydo = tdest[1] - torg[1];
xao = tapex[0] - torg[0];
yao = tapex[1] - torg[1];
dodist = xdo * xdo + ydo * ydo;
aodist = xao * xao + yao * yao;
dadist = (tdest[0] - tapex[0]) * (tdest[0] - tapex[0]) +
(tdest[1] - tapex[1]) * (tdest[1] - tapex[1]);
if (b->noexact) {
denominator = 0.5 / (xdo * yao - xao * ydo);
} else {
/* Use the counterclockwise() routine to ensure a positive (and */
/* reasonably accurate) result, avoiding any possibility of */
/* division by zero. */
denominator = 0.5 / counterclockwise(m, b, tdest, tapex, torg);
/* Don't count the above as an orientation test. */
m->counterclockcount--;
}
dx = (yao * dodist - ydo * aodist) * denominator;
dy = (xdo * aodist - xao * dodist) * denominator;
/* Find the (squared) length of the triangle's shortest edge. This */
/* serves as a conservative estimate of the insertion radius of the */
/* circumcenter's parent. The estimate is used to ensure that */
/* the algorithm terminates even if very small angles appear in */
/* the input PSLG. */
if ((dodist < aodist) && (dodist < dadist)) {
if (offcenter && (b->offconstant > 0.0)) {
/* Find the position of the off-center, as described by Alper Ungor. */
dxoff = 0.5 * xdo - b->offconstant * ydo;
dyoff = 0.5 * ydo + b->offconstant * xdo;
/* If the off-center is closer to the origin than the */
/* circumcenter, use the off-center instead. */
if (dxoff * dxoff + dyoff * dyoff < dx * dx + dy * dy) {
dx = dxoff;
dy = dyoff;
}
}
} else if (aodist < dadist) {
if (offcenter && (b->offconstant > 0.0)) {
dxoff = 0.5 * xao + b->offconstant * yao;
dyoff = 0.5 * yao - b->offconstant * xao;
/* If the off-center is closer to the origin than the */
/* circumcenter, use the off-center instead. */
if (dxoff * dxoff + dyoff * dyoff < dx * dx + dy * dy) {
dx = dxoff;
dy = dyoff;
}
}
} else {
if (offcenter && (b->offconstant > 0.0)) {
dxoff = 0.5 * (tapex[0] - tdest[0]) -
b->offconstant * (tapex[1] - tdest[1]);
dyoff = 0.5 * (tapex[1] - tdest[1]) +
b->offconstant * (tapex[0] - tdest[0]);
/* If the off-center is closer to the destination than the */
/* circumcenter, use the off-center instead. */
if (dxoff * dxoff + dyoff * dyoff <
(dx - xdo) * (dx - xdo) + (dy - ydo) * (dy - ydo)) {
dx = xdo + dxoff;
dy = ydo + dyoff;
}
}
}
circumcenter[0] = torg[0] + dx;
circumcenter[1] = torg[1] + dy;
/* To interpolate vertex attributes for the new vertex inserted at */
/* the circumcenter, define a coordinate system with a xi-axis, */
/* directed from the triangle's origin to its destination, and */
/* an eta-axis, directed from its origin to its apex. */
/* Calculate the xi and eta coordinates of the circumcenter. */
*xi = (yao * dx - xao * dy) * (2.0 * denominator);
*eta = (xdo * dy - ydo * dx) * (2.0 * denominator);
}
/** **/
/** **/
/********* Geometric primitives end here *********/
/*****************************************************************************/
/* */
/* triangleinit() Initialize some variables. */
/* */
/*****************************************************************************/
#ifdef ANSI_DECLARATORS
void triangleinit(struct mesh *m)
#else /* not ANSI_DECLARATORS */
void triangleinit(m)
struct mesh *m;
#endif /* not ANSI_DECLARATORS */
{
poolzero(&m->vertices);
poolzero(&m->triangles);
poolzero(&m->subsegs);
poolzero(&m->viri);
poolzero(&m->badsubsegs);
poolzero(&m->badtriangles);
poolzero(&m->flipstackers);
poolzero(&m->splaynodes);
m->recenttri.tri = (triangle *) NULL; /* No triangle has been visited yet. */
m->undeads = 0; /* No eliminated input vertices yet. */
m->samples = 1; /* Point location should take at least one sample. */
m->checksegments = 0; /* There are no segments in the triangulation yet. */
m->checkquality = 0; /* The quality triangulation stage has not begun. */
m->incirclecount = m->counterclockcount = m->orient3dcount = 0;
m->hyperbolacount = m->circletopcount = m->circumcentercount = 0;
randomseed = 1;
exactinit(); /* Initialize exact arithmetic constants. */
}
/*****************************************************************************/
/* */
/* randomnation() Generate a random number between 0 and `choices' - 1. */
/* */
/* This is a simple linear congruential random number generator. Hence, it */
/* is a bad random number generator, but good enough for most randomized */
/* geometric algorithms. */
/* */
/*****************************************************************************/
#ifdef ANSI_DECLARATORS
unsigned long randomnation(unsigned int choices)
#else /* not ANSI_DECLARATORS */
unsigned long randomnation(choices)
unsigned int choices;
#endif /* not ANSI_DECLARATORS */
{
randomseed = (randomseed * 1366l + 150889l) % 714025l;
return randomseed / (714025l / choices + 1);
}
/********* Mesh quality testing routines begin here *********/
/** **/
/** **/
/*****************************************************************************/
/* */
/* checkmesh() Test the mesh for topological consistency. */
/* */
/*****************************************************************************/
#ifndef REDUCED
#ifdef ANSI_DECLARATORS
void checkmesh(struct mesh *m, struct behavior *b)
#else /* not ANSI_DECLARATORS */
void checkmesh(m, b)
struct mesh *m;
struct behavior *b;
#endif /* not ANSI_DECLARATORS */
{
struct otri triangleloop;
struct otri oppotri, oppooppotri;
vertex triorg, tridest, triapex;
vertex oppoorg, oppodest;
int horrors;
int saveexact;
triangle ptr; /* Temporary variable used by sym(). */
/* Temporarily turn on exact arithmetic if it's off. */
saveexact = b->noexact;
b->noexact = 0;
if (!b->quiet) {
printf(" Checking consistency of mesh...n");
}
horrors = 0;
/* Run through the list of triangles, checking each one. */
traversalinit(&m->triangles);
triangleloop.tri = triangletraverse(m);
while (triangleloop.tri != (triangle *) NULL) {
/* Check all three edges of the triangle. */
for (triangleloop.orient = 0; triangleloop.orient < 3;
triangleloop.orient++) {
org(triangleloop, triorg);
dest(triangleloop, tridest);
if (triangleloop.orient == 0) { /* Only test for inversion once. */
/* Test if the triangle is flat or inverted. */
apex(triangleloop, triapex);
if (counterclockwise(m, b, triorg, tridest, triapex) <= 0.0) {
printf(" !! !! Inverted ");
printtriangle(m, b, &triangleloop);
horrors++;
}
}
/* Find the neighboring triangle on this edge. */
sym(triangleloop, oppotri);
if (oppotri.tri != m->dummytri) {
/* Check that the triangle's neighbor knows it's a neighbor. */
sym(oppotri, oppooppotri);
if ((triangleloop.tri != oppooppotri.tri)
|| (triangleloop.orient != oppooppotri.orient)) {
printf(" !! !! Asymmetric triangle-triangle bond:n");
if (triangleloop.tri == oppooppotri.tri) {
printf(" (Right triangle, wrong orientation)n");
}
printf(" First ");
printtriangle(m, b, &triangleloop);
printf(" Second (nonreciprocating) ");
printtriangle(m, b, &oppotri);
horrors++;
}
/* Check that both triangles agree on the identities */
/* of their shared vertices. */
org(oppotri, oppoorg);
dest(oppotri, oppodest);
if ((triorg != oppodest) || (tridest != oppoorg)) {
printf(" !! !! Mismatched edge coordinates between two triangles:n"
);
printf(" First mismatched ");
printtriangle(m, b, &triangleloop);
printf(" Second mismatched ");
printtriangle(m, b, &oppotri);
horrors++;
}
}
}
triangleloop.tri = triangletraverse(m);
}
if (horrors == 0) {
if (!b->quiet) {
printf(" In my studied opinion, the mesh appears to be consistent.n");
}
} else if (horrors == 1) {
printf(" !! !! !! !! Precisely one festering wound discovered.n");
} else {
printf(" !! !! !! !! %d abominations witnessed.n", horrors);
}
/* Restore the status of exact arithmetic. */
b->noexact = saveexact;
}
#endif /* not REDUCED */
/*****************************************************************************/
/* */
/* checkdelaunay() Ensure that the mesh is (constrained) Delaunay. */
/* */
/*****************************************************************************/
#ifndef REDUCED
#ifdef ANSI_DECLARATORS
void checkdelaunay(struct mesh *m, struct behavior *b)
#else /* not ANSI_DECLARATORS */
void checkdelaunay(m, b)
struct mesh *m;
struct behavior *b;
#endif /* not ANSI_DECLARATORS */
{
struct otri triangleloop;
struct otri oppotri;
struct osub opposubseg;
vertex triorg, tridest, triapex;
vertex oppoapex;
int shouldbedelaunay;
int horrors;
int saveexact;
triangle ptr; /* Temporary variable used by sym(). */
subseg sptr; /* Temporary variable used by tspivot(). */
/* Temporarily turn on exact arithmetic if it's off. */
saveexact = b->noexact;
b->noexact = 0;
if (!b->quiet) {
printf(" Checking Delaunay property of mesh...n");
}
horrors = 0;
/* Run through the list of triangles, checking each one. */
traversalinit(&m->triangles);
triangleloop.tri = triangletraverse(m);
while (triangleloop.tri != (triangle *) NULL) {
/* Check all three edges of the triangle. */
for (triangleloop.orient = 0; triangleloop.orient < 3;
triangleloop.orient++) {
org(triangleloop, triorg);
dest(triangleloop, tridest);
apex(triangleloop, triapex);
sym(triangleloop, oppotri);
apex(oppotri, oppoapex);
/* Only test that the edge is locally Delaunay if there is an */
/* adjoining triangle whose pointer is larger (to ensure that */
/* each pair isn't tested twice). */
shouldbedelaunay = (oppotri.tri != m->dummytri) &&
!deadtri(oppotri.tri) && (triangleloop.tri < oppotri.tri) &&
(triorg != m->infvertex1) && (triorg != m->infvertex2) &&
(triorg != m->infvertex3) &&
(tridest != m->infvertex1) && (tridest != m->infvertex2) &&
(tridest != m->infvertex3) &&
(triapex != m->infvertex1) && (triapex != m->infvertex2) &&
(triapex != m->infvertex3) &&
(oppoapex != m->infvertex1) && (oppoapex != m->infvertex2) &&
(oppoapex != m->infvertex3);
if (m->checksegments && shouldbedelaunay) {
/* If a subsegment separates the triangles, then the edge is */
/* constrained, so no local Delaunay test should be done. */
tspivot(triangleloop, opposubseg);
if (opposubseg.ss != m->dummysub){
shouldbedelaunay = 0;
}
}
if (shouldbedelaunay) {
if (nonregular(m, b, triorg, tridest, triapex, oppoapex) > 0.0) {
if (!b->weighted) {
printf(" !! !! Non-Delaunay pair of triangles:n");
printf(" First non-Delaunay ");
printtriangle(m, b, &triangleloop);
printf(" Second non-Delaunay ");
} else {
printf(" !! !! Non-regular pair of triangles:n");
printf(" First non-regular ");
printtriangle(m, b, &triangleloop);
printf(" Second non-regular ");
}
printtriangle(m, b, &oppotri);
horrors++;
}
}
}
triangleloop.tri = triangletraverse(m);
}
if (horrors == 0) {
if (!b->quiet) {
printf(
" By virtue of my perceptive intelligence, I declare the mesh Delaunay.n");
}
} else if (horrors == 1) {
printf(
" !! !! !! !! Precisely one terrifying transgression identified.n");
} else {
printf(" !! !! !! !! %d obscenities viewed with horror.n", horrors);
}
/* Restore the status of exact arithmetic. */
b->noexact = saveexact;
}
#endif /* not REDUCED */
/*****************************************************************************/
/* */
/* enqueuebadtriang() Add a bad triangle data structure to the end of a */
/* queue. */
/* */
/* The queue is actually a set of 4096 queues. I use multiple queues to */
/* give priority to smaller angles. I originally implemented a heap, but */
/* the queues are faster by a larger margin than I'd suspected. */
/* */
/*****************************************************************************/
#ifndef CDT_ONLY
#ifdef ANSI_DECLARATORS
void enqueuebadtriang(struct mesh *m, struct behavior *b,
struct badtriang *badtri)
#else /* not ANSI_DECLARATORS */
void enqueuebadtriang(m, b, badtri)
struct mesh *m;
struct behavior *b;
struct badtriang *badtri;
#endif /* not ANSI_DECLARATORS */
{
REAL length, multiplier;
int exponent, expincrement;
int queuenumber;
int posexponent;
int i;
if (b->verbose > 2) {
printf(" Queueing bad triangle:n");
printf(" (%.12g, %.12g) (%.12g, %.12g) (%.12g, %.12g)n",
badtri->triangorg[0], badtri->triangorg[1],
badtri->triangdest[0], badtri->triangdest[1],
badtri->triangapex[0], badtri->triangapex[1]);
}
/* Determine the appropriate queue to put the bad triangle into. */
/* Recall that the key is the square of its shortest edge length. */
if (badtri->key >= 1.0) {
length = badtri->key;
posexponent = 1;
} else {
/* `badtri->key' is 2.0 to a negative exponent, so we'll record that */
/* fact and use the reciprocal of `badtri->key', which is > 1.0. */
length = 1.0 / badtri->key;
posexponent = 0;
}
/* `length' is approximately 2.0 to what exponent? The following code */
/* determines the answer in time logarithmic in the exponent. */
exponent = 0;
while (length > 2.0) {
/* Find an approximation by repeated squaring of two. */
expincrement = 1;
multiplier = 0.5;
while (length * multiplier * multiplier > 1.0) {
expincrement *= 2;
multiplier *= multiplier;
}
/* Reduce the value of `length', then iterate if necessary. */
exponent += expincrement;
length *= multiplier;
}
/* `length' is approximately squareroot(2.0) to what exponent? */
exponent = 2.0 * exponent + (length > SQUAREROOTTWO);
/* `exponent' is now in the range 0...2047 for IEEE double precision. */
/* Choose a queue in the range 0...4095. The shortest edges have the */
/* highest priority (queue 4095). */
if (posexponent) {
queuenumber = 2047 - exponent;
} else {
queuenumber = 2048 + exponent;
}
/* Are we inserting into an empty queue? */
if (m->queuefront[queuenumber] == (struct badtriang *) NULL) {
/* Yes, we are inserting into an empty queue. */
/* Will this become the highest-priority queue? */
if (queuenumber > m->firstnonemptyq) {
/* Yes, this is the highest-priority queue. */
m->nextnonemptyq[queuenumber] = m->firstnonemptyq;
m->firstnonemptyq = queuenumber;
} else {
/* No, this is not the highest-priority queue. */
/* Find the queue with next higher priority. */
i = queuenumber + 1;
while (m->queuefront[i] == (struct badtriang *) NULL) {
i++;
}
/* Mark the newly nonempty queue as following a higher-priority queue. */
m->nextnonemptyq[queuenumber] = m->nextnonemptyq[i];
m->nextnonemptyq[i] = queuenumber;
}
/* Put the bad triangle at the beginning of the (empty) queue. */
m->queuefront[queuenumber] = badtri;
} else {
/* Add the bad triangle to the end of an already nonempty queue. */
m->queuetail[queuenumber]->nexttriang = badtri;
}
/* Maintain a pointer to the last triangle of the queue. */
m->queuetail[queuenumber] = badtri;
/* Newly enqueued bad triangle has no successor in the queue. */
badtri->nexttriang = (struct badtriang *) NULL;
}
#endif /* not CDT_ONLY */
/*****************************************************************************/
/* */
/* enqueuebadtri() Add a bad triangle to the end of a queue. */
/* */
/* Allocates a badtriang data structure for the triangle, then passes it to */
/* enqueuebadtriang(). */
/* */
/*****************************************************************************/
#ifndef CDT_ONLY
#ifdef ANSI_DECLARATORS
void enqueuebadtri(struct mesh *m, struct behavior *b, struct otri *enqtri,
REAL minedge, vertex enqapex, vertex enqorg, vertex enqdest)
#else /* not ANSI_DECLARATORS */
void enqueuebadtri(m, b, enqtri, minedge, enqapex, enqorg, enqdest)
struct mesh *m;
struct behavior *b;
struct otri *enqtri;
REAL minedge;
vertex enqapex;
vertex enqorg;
vertex enqdest;
#endif /* not ANSI_DECLARATORS */
{
struct badtriang *newbad;
/* Allocate space for the bad triangle. */
newbad = (struct badtriang *) poolalloc(&m->badtriangles);
newbad->poortri = encode(*enqtri);
newbad->key = minedge;
newbad->triangapex = enqapex;
newbad->triangorg = enqorg;
newbad->triangdest = enqdest;
enqueuebadtriang(m, b, newbad);
}
#endif /* not CDT_ONLY */
/*****************************************************************************/
/* */
/* dequeuebadtriang() Remove a triangle from the front of the queue. */
/* */
/*****************************************************************************/
#ifndef CDT_ONLY
#ifdef ANSI_DECLARATORS
struct badtriang *dequeuebadtriang(struct mesh *m)
#else /* not ANSI_DECLARATORS */
struct badtriang *dequeuebadtriang(m)
struct mesh *m;
#endif /* not ANSI_DECLARATORS */
{
struct badtriang *result;
/* If no queues are nonempty, return NULL. */
if (m->firstnonemptyq < 0) {
return (struct badtriang *) NULL;
}
/* Find the first triangle of the highest-priority queue. */
result = m->queuefront[m->firstnonemptyq];
/* Remove the triangle from the queue. */
m->queuefront[m->firstnonemptyq] = result->nexttriang;
/* If this queue is now empty, note the new highest-priority */
/* nonempty queue. */
if (result == m->queuetail[m->firstnonemptyq]) {
m->firstnonemptyq = m->nextnonemptyq[m->firstnonemptyq];
}
return result;
}
#endif /* not CDT_ONLY */
/*****************************************************************************/
/* */
/* checkseg4encroach() Check a subsegment to see if it is encroached; add */
/* it to the list if it is. */
/* */
/* A subsegment is encroached if there is a vertex in its diametral lens. */
/* For Ruppert's algorithm (-D switch), the "diametral lens" is the */
/* diametral circle. For Chew's algorithm (default), the diametral lens is */
/* just big enough to enclose two isosceles triangles whose bases are the */
/* subsegment. Each of the two isosceles triangles has two angles equal */
/* to `b->minangle'. */
/* */
/* Chew's algorithm does not require diametral lenses at all--but they save */
/* time. Any vertex inside a subsegment's diametral lens implies that the */
/* triangle adjoining the subsegment will be too skinny, so it's only a */
/* matter of time before the encroaching vertex is deleted by Chew's */
/* algorithm. It's faster to simply not insert the doomed vertex in the */
/* first place, which is why I use diametral lenses with Chew's algorithm. */
/* */
/* Returns a nonzero value if the subsegment is encroached. */
/* */
/*****************************************************************************/
#ifndef CDT_ONLY
#ifdef ANSI_DECLARATORS
int checkseg4encroach(struct mesh *m, struct behavior *b,
struct osub *testsubseg)
#else /* not ANSI_DECLARATORS */
int checkseg4encroach(m, b, testsubseg)
struct mesh *m;
struct behavior *b;
struct osub *testsubseg;
#endif /* not ANSI_DECLARATORS */
{
struct otri neighbortri;
struct osub testsym;
struct badsubseg *encroachedseg;
REAL dotproduct;
int encroached;
int sides;
vertex eorg, edest, eapex;
triangle ptr; /* Temporary variable used by stpivot(). */
encroached = 0;
sides = 0;
sorg(*testsubseg, eorg);
sdest(*testsubseg, edest);
/* Check one neighbor of the subsegment. */
stpivot(*testsubseg, neighbortri);
/* Does the neighbor exist, or is this a boundary edge? */
if (neighbortri.tri != m->dummytri) {
sides++;
/* Find a vertex opposite this subsegment. */
apex(neighbortri, eapex);
/* Check whether the apex is in the diametral lens of the subsegment */
/* (the diametral circle if `conformdel' is set). A dot product */
/* of two sides of the triangle is used to check whether the angle */
/* at the apex is greater than (180 - 2 `minangle') degrees (for */
/* lenses; 90 degrees for diametral circles). */
dotproduct = (eorg[0] - eapex[0]) * (edest[0] - eapex[0]) +
(eorg[1] - eapex[1]) * (edest[1] - eapex[1]);
if (dotproduct < 0.0) {
if (b->conformdel ||
(dotproduct * dotproduct >=
(2.0 * b->goodangle - 1.0) * (2.0 * b->goodangle - 1.0) *
((eorg[0] - eapex[0]) * (eorg[0] - eapex[0]) +
(eorg[1] - eapex[1]) * (eorg[1] - eapex[1])) *
((edest[0] - eapex[0]) * (edest[0] - eapex[0]) +
(edest[1] - eapex[1]) * (edest[1] - eapex[1])))) {
encroached = 1;
}
}
}
/* Check the other neighbor of the subsegment. */
ssym(*testsubseg, testsym);
stpivot(testsym, neighbortri);
/* Does the neighbor exist, or is this a boundary edge? */
if (neighbortri.tri != m->dummytri) {
sides++;
/* Find the other vertex opposite this subsegment. */
apex(neighbortri, eapex);
/* Check whether the apex is in the diametral lens of the subsegment */
/* (or the diametral circle, if `conformdel' is set). */
dotproduct = (eorg[0] - eapex[0]) * (edest[0] - eapex[0]) +
(eorg[1] - eapex[1]) * (edest[1] - eapex[1]);
if (dotproduct < 0.0) {
if (b->conformdel ||
(dotproduct * dotproduct >=
(2.0 * b->goodangle - 1.0) * (2.0 * b->goodangle - 1.0) *
((eorg[0] - eapex[0]) * (eorg[0] - eapex[0]) +
(eorg[1] - eapex[1]) * (eorg[1] - eapex[1])) *
((edest[0] - eapex[0]) * (edest[0] - eapex[0]) +
(edest[1] - eapex[1]) * (edest[1] - eapex[1])))) {
encroached += 2;
}
}
}
if (encroached && (!b->nobisect || ((b->nobisect == 1) && (sides == 2)))) {
if (b->verbose > 2) {
printf(
" Queueing encroached subsegment (%.12g, %.12g) (%.12g, %.12g).n",
eorg[0], eorg[1], edest[0], edest[1]);
}
/* Add the subsegment to the list of encroached subsegments. */
/* Be sure to get the orientation right. */
encroachedseg = (struct badsubseg *) poolalloc(&m->badsubsegs);
if (encroached == 1) {
encroachedseg->encsubseg = sencode(*testsubseg);
encroachedseg->subsegorg = eorg;
encroachedseg->subsegdest = edest;
} else {
encroachedseg->encsubseg = sencode(testsym);
encroachedseg->subsegorg = edest;
encroachedseg->subsegdest = eorg;
}
}
return encroached;
}
#endif /* not CDT_ONLY */
/*****************************************************************************/
/* */
/* testtriangle() Test a triangle for quality and size. */
/* */
/* Tests a triangle to see if it satisfies the minimum angle condition and */
/* the maximum area condition. Triangles that aren't up to spec are added */
/* to the bad triangle queue. */
/* */
/*****************************************************************************/
#ifndef CDT_ONLY
#ifdef ANSI_DECLARATORS
void testtriangle(struct mesh *m, struct behavior *b, struct otri *testtri)
#else /* not ANSI_DECLARATORS */
void testtriangle(m, b, testtri)
struct mesh *m;
struct behavior *b;
struct otri *testtri;
#endif /* not ANSI_DECLARATORS */
{
struct otri tri1, tri2;
struct osub testsub;
vertex torg, tdest, tapex;
vertex base1, base2;
vertex org1, dest1, org2, dest2;
vertex joinvertex;
REAL dxod, dyod, dxda, dyda, dxao, dyao;
REAL dxod2, dyod2, dxda2, dyda2, dxao2, dyao2;
REAL apexlen, orglen, destlen, minedge;
REAL angle;
REAL area;
REAL dist1, dist2;
subseg sptr; /* Temporary variable used by tspivot(). */
triangle ptr; /* Temporary variable used by oprev() and dnext(). */
org(*testtri, torg);
dest(*testtri, tdest);
apex(*testtri, tapex);
dxod = torg[0] - tdest[0];
dyod = torg[1] - tdest[1];
dxda = tdest[0] - tapex[0];
dyda = tdest[1] - tapex[1];
dxao = tapex[0] - torg[0];
dyao = tapex[1] - torg[1];
dxod2 = dxod * dxod;
dyod2 = dyod * dyod;
dxda2 = dxda * dxda;
dyda2 = dyda * dyda;
dxao2 = dxao * dxao;
dyao2 = dyao * dyao;
/* Find the lengths of the triangle's three edges. */
apexlen = dxod2 + dyod2;
orglen = dxda2 + dyda2;
destlen = dxao2 + dyao2;
if ((apexlen < orglen) && (apexlen < destlen)) {
/* The edge opposite the apex is shortest. */
minedge = apexlen;
/* Find the square of the cosine of the angle at the apex. */
angle = dxda * dxao + dyda * dyao;
angle = angle * angle / (orglen * destlen);
base1 = torg;
base2 = tdest;
otricopy(*testtri, tri1);
} else if (orglen < destlen) {
/* The edge opposite the origin is shortest. */
minedge = orglen;
/* Find the square of the cosine of the angle at the origin. */
angle = dxod * dxao + dyod * dyao;
angle = angle * angle / (apexlen * destlen);
base1 = tdest;
base2 = tapex;
lnext(*testtri, tri1);
} else {
/* The edge opposite the destination is shortest. */
minedge = destlen;
/* Find the square of the cosine of the angle at the destination. */
angle = dxod * dxda + dyod * dyda;
angle = angle * angle / (apexlen * orglen);
base1 = tapex;
base2 = torg;
lprev(*testtri, tri1);
}
if (b->vararea || b->fixedarea || b->usertest) {
/* Check whether the area is larger than permitted. */
area = 0.5 * (dxod * dyda - dyod * dxda);
if (b->fixedarea && (area > b->maxarea)) {
/* Add this triangle to the list of bad triangles. */
enqueuebadtri(m, b, testtri, minedge, tapex, torg, tdest);
return;
}
/* Nonpositive area constraints are treated as unconstrained. */
if ((b->vararea) && (area > areabound(*testtri)) &&
(areabound(*testtri) > 0.0)) {
/* Add this triangle to the list of bad triangles. */
enqueuebadtri(m, b, testtri, minedge, tapex, torg, tdest);
return;
}
if (b->usertest) {
/* Check whether the user thinks this triangle is too large. */
if (triunsuitable(torg, tdest, tapex, area)) {
enqueuebadtri(m, b, testtri, minedge, tapex, torg, tdest);
return;
}
}
}
/* Check whether the angle is smaller than permitted. */
if (angle > b->goodangle) {
/* Use the rules of Miller, Pav, and Walkington to decide that certain */
/* triangles should not be split, even if they have bad angles. */
/* A skinny triangle is not split if its shortest edge subtends a */
/* small input angle, and both endpoints of the edge lie on a */
/* concentric circular shell. For convenience, I make a small */
/* adjustment to that rule: I check if the endpoints of the edge */
/* both lie in segment interiors, equidistant from the apex where */
/* the two segments meet. */
/* First, check if both points lie in segment interiors. */
if ((vertextype(base1) == SEGMENTVERTEX) &&
(vertextype(base2) == SEGMENTVERTEX)) {
/* Check if both points lie in a common segment. If they do, the */
/* skinny triangle is enqueued to be split as usual. */
tspivot(tri1, testsub);
if (testsub.ss == m->dummysub) {
/* No common segment. Find a subsegment that contains `torg'. */
otricopy(tri1, tri2);
do {
oprevself(tri1);
tspivot(tri1, testsub);
} while (testsub.ss == m->dummysub);
/* Find the endpoints of the containing segment. */
segorg(testsub, org1);
segdest(testsub, dest1);
/* Find a subsegment that contains `tdest'. */
do {
dnextself(tri2);
tspivot(tri2, testsub);
} while (testsub.ss == m->dummysub);
/* Find the endpoints of the containing segment. */
segorg(testsub, org2);
segdest(testsub, dest2);
/* Check if the two containing segments have an endpoint in common. */
joinvertex = (vertex) NULL;
if ((dest1[0] == org2[0]) && (dest1[1] == org2[1])) {
joinvertex = dest1;
} else if ((org1[0] == dest2[0]) && (org1[1] == dest2[1])) {
joinvertex = org1;
}
if (joinvertex != (vertex) NULL) {
/* Compute the distance from the common endpoint (of the two */
/* segments) to each of the endpoints of the shortest edge. */
dist1 = ((base1[0] - joinvertex[0]) * (base1[0] - joinvertex[0]) +
(base1[1] - joinvertex[1]) * (base1[1] - joinvertex[1]));
dist2 = ((base2[0] - joinvertex[0]) * (base2[0] - joinvertex[0]) +
(base2[1] - joinvertex[1]) * (base2[1] - joinvertex[1]));
/* If the two distances are equal, don't split the triangle. */
if ((dist1 < 1.001 * dist2) && (dist1 > 0.999 * dist2)) {
/* Return now to avoid enqueueing the bad triangle. */
return;
}
}
}
}
/* Add this triangle to the list of bad triangles. */
enqueuebadtri(m, b, testtri, minedge, tapex, torg, tdest);
}
}
#endif /* not CDT_ONLY */
/** **/
/** **/
/********* Mesh quality testing routines end here *********/
/********* Point location routines begin here *********/
/** **/
/** **/
/*****************************************************************************/
/* */
/* makevertexmap() Construct a mapping from vertices to triangles to */
/* improve the speed of point location for segment */
/* insertion. */
/* */
/* Traverses all the triangles, and provides each corner of each triangle */
/* with a pointer to that triangle. Of course, pointers will be */
/* overwritten by other pointers because (almost) each vertex is a corner */
/* of several triangles, but in the end every vertex will point to some */
/* triangle that contains it. */
/* */
/*****************************************************************************/
#ifdef ANSI_DECLARATORS
void makevertexmap(struct mesh *m, struct behavior *b)
#else /* not ANSI_DECLARATORS */
void makevertexmap(m, b)
struct mesh *m;
struct behavior *b;
#endif /* not ANSI_DECLARATORS */
{
struct otri triangleloop;
vertex triorg;
if (b->verbose) {
printf(" Constructing mapping from vertices to triangles.n");
}
traversalinit(&m->triangles);
triangleloop.tri = triangletraverse(m);
while (triangleloop.tri != (triangle *) NULL) {
/* Check all three vertices of the triangle. */
for (triangleloop.orient = 0; triangleloop.orient < 3;
triangleloop.orient++) {
org(triangleloop, triorg);
setvertex2tri(triorg, encode(triangleloop));
}
triangleloop.tri = triangletraverse(m);
}
}
/*****************************************************************************/
/* */
/* preciselocate() Find a triangle or edge containing a given point. */
/* */
/* Begins its search from `searchtri'. It is important that `searchtri' */
/* be a handle with the property that `searchpoint' is strictly to the left */
/* of the edge denoted by `searchtri', or is collinear with that edge and */
/* does not intersect that edge. (In particular, `searchpoint' should not */
/* be the origin or destination of that edge.) */
/* */
/* These conditions are imposed because preciselocate() is normally used in */
/* one of two situations: */
/* */
/* (1) To try to find the location to insert a new point. Normally, we */
/* know an edge that the point is strictly to the left of. In the */
/* incremental Delaunay algorithm, that edge is a bounding box edge. */
/* In Ruppert's Delaunay refinement algorithm for quality meshing, */
/* that edge is the shortest edge of the triangle whose circumcenter */
/* is being inserted. */
/* */
/* (2) To try to find an existing point. In this case, any edge on the */
/* convex hull is a good starting edge. You must screen out the */
/* possibility that the vertex sought is an endpoint of the starting */
/* edge before you call preciselocate(). */
/* */
/* On completion, `searchtri' is a triangle that contains `searchpoint'. */
/* */
/* This implementation differs from that given by Guibas and Stolfi. It */
/* walks from triangle to triangle, crossing an edge only if `searchpoint' */
/* is on the other side of the line containing that edge. After entering */
/* a triangle, there are two edges by which one can leave that triangle. */
/* If both edges are valid (`searchpoint' is on the other side of both */
/* edges), one of the two is chosen by drawing a line perpendicular to */
/* the entry edge (whose endpoints are `forg' and `fdest') passing through */
/* `fapex'. Depending on which side of this perpendicular `searchpoint' */
/* falls on, an exit edge is chosen. */
/* */
/* This implementation is empirically faster than the Guibas and Stolfi */
/* point location routine (which I originally used), which tends to spiral */
/* in toward its target. */
/* */
/* Returns ONVERTEX if the point lies on an existing vertex. `searchtri' */
/* is a handle whose origin is the existing vertex. */
/* */
/* Returns ONEDGE if the point lies on a mesh edge. `searchtri' is a */
/* handle whose primary edge is the edge on which the point lies. */
/* */
/* Returns INTRIANGLE if the point lies strictly within a triangle. */
/* `searchtri' is a handle on the triangle that contains the point. */
/* */
/* Returns OUTSIDE if the point lies outside the mesh. `searchtri' is a */
/* handle whose primary edge the point is to the right of. This might */
/* occur when the circumcenter of a triangle falls just slightly outside */
/* the mesh due to floating-point roundoff error. It also occurs when */
/* seeking a hole or region point that a foolish user has placed outside */
/* the mesh. */
/* */
/* If `stopatsubsegment' is nonzero, the search will stop if it tries to */
/* walk through a subsegment, and will return OUTSIDE. */
/* */
/* WARNING: This routine is designed for convex triangulations, and will */
/* not generally work after the holes and concavities have been carved. */
/* However, it can still be used to find the circumcenter of a triangle, as */
/* long as the search is begun from the triangle in question. */
/* */
/*****************************************************************************/
#ifdef ANSI_DECLARATORS
enum locateresult preciselocate(struct mesh *m, struct behavior *b,
vertex searchpoint, struct otri *searchtri,
int stopatsubsegment)
#else /* not ANSI_DECLARATORS */
enum locateresult preciselocate(m, b, searchpoint, searchtri, stopatsubsegment)
struct mesh *m;
struct behavior *b;
vertex searchpoint;
struct otri *searchtri;
int stopatsubsegment;
#endif /* not ANSI_DECLARATORS */
{
struct otri backtracktri;
struct osub checkedge;
vertex forg, fdest, fapex;
REAL orgorient, destorient;
int moveleft;
triangle ptr; /* Temporary variable used by sym(). */
subseg sptr; /* Temporary variable used by tspivot(). */
if (b->verbose > 2) {
printf(" Searching for point (%.12g, %.12g).n",
searchpoint[0], searchpoint[1]);
}
/* Where are we? */
org(*searchtri, forg);
dest(*searchtri, fdest);
apex(*searchtri, fapex);
while (1) {
if (b->verbose > 2) {
printf(" At (%.12g, %.12g) (%.12g, %.12g) (%.12g, %.12g)n",
forg[0], forg[1], fdest[0], fdest[1], fapex[0], fapex[1]);
}
/* Check whether the apex is the point we seek. */
if ((fapex[0] == searchpoint[0]) && (fapex[1] == searchpoint[1])) {
lprevself(*searchtri);
return ONVERTEX;
}
/* Does the point lie on the other side of the line defined by the */
/* triangle edge opposite the triangle's destination? */
destorient = counterclockwise(m, b, forg, fapex, searchpoint);
/* Does the point lie on the other side of the line defined by the */
/* triangle edge opposite the triangle's origin? */
orgorient = counterclockwise(m, b, fapex, fdest, searchpoint);
if (destorient > 0.0) {
if (orgorient > 0.0) {
/* Move left if the inner product of (fapex - searchpoint) and */
/* (fdest - forg) is positive. This is equivalent to drawing */
/* a line perpendicular to the line (forg, fdest) and passing */
/* through `fapex', and determining which side of this line */
/* `searchpoint' falls on. */
moveleft = (fapex[0] - searchpoint[0]) * (fdest[0] - forg[0]) +
(fapex[1] - searchpoint[1]) * (fdest[1] - forg[1]) > 0.0;
} else {
moveleft = 1;
}
} else {
if (orgorient > 0.0) {
moveleft = 0;
} else {
/* The point we seek must be on the boundary of or inside this */
/* triangle. */
if (destorient == 0.0) {
lprevself(*searchtri);
return ONEDGE;
}
if (orgorient == 0.0) {
lnextself(*searchtri);
return ONEDGE;
}
return INTRIANGLE;
}
}
/* Move to another triangle. Leave a trace `backtracktri' in case */
/* floating-point roundoff or some such bogey causes us to walk */
/* off a boundary of the triangulation. */
if (moveleft) {
lprev(*searchtri, backtracktri);
fdest = fapex;
} else {
lnext(*searchtri, backtracktri);
forg = fapex;
}
sym(backtracktri, *searchtri);
if (m->checksegments && stopatsubsegment) {
/* Check for walking through a subsegment. */
tspivot(backtracktri, checkedge);
if (checkedge.ss != m->dummysub) {
/* Go back to the last triangle. */
otricopy(backtracktri, *searchtri);
return OUTSIDE;
}
}
/* Check for walking right out of the triangulation. */
if (searchtri->tri == m->dummytri) {
/* Go back to the last triangle. */
otricopy(backtracktri, *searchtri);
return OUTSIDE;
}
apex(*searchtri, fapex);
}
}
/*****************************************************************************/
/* */
/* locate() Find a triangle or edge containing a given point. */
/* */
/* Searching begins from one of: the input `searchtri', a recently */
/* encountered triangle `recenttri', or from a triangle chosen from a */
/* random sample. The choice is made by determining which triangle's */
/* origin is closest to the point we are searching for. Normally, */
/* `searchtri' should be a handle on the convex hull of the triangulation. */
/* */
/* Details on the random sampling method can be found in the Mucke, Saias, */
/* and Zhu paper cited in the header of this code. */
/* */
/* On completion, `searchtri' is a triangle that contains `searchpoint'. */
/* */
/* Returns ONVERTEX if the point lies on an existing vertex. `searchtri' */
/* is a handle whose origin is the existing vertex. */
/* */
/* Returns ONEDGE if the point lies on a mesh edge. `searchtri' is a */
/* handle whose primary edge is the edge on which the point lies. */
/* */
/* Returns INTRIANGLE if the point lies strictly within a triangle. */
/* `searchtri' is a handle on the triangle that contains the point. */
/* */
/* Returns OUTSIDE if the point lies outside the mesh. `searchtri' is a */
/* handle whose primary edge the point is to the right of. This might */
/* occur when the circumcenter of a triangle falls just slightly outside */
/* the mesh due to floating-point roundoff error. It also occurs when */
/* seeking a hole or region point that a foolish user has placed outside */
/* the mesh. */
/* */
/* WARNING: This routine is designed for convex triangulations, and will */
/* not generally work after the holes and concavities have been carved. */
/* */
/*****************************************************************************/
#ifdef ANSI_DECLARATORS
enum locateresult locate(struct mesh *m, struct behavior *b,
vertex searchpoint, struct otri *searchtri)
#else /* not ANSI_DECLARATORS */
enum locateresult locate(m, b, searchpoint, searchtri)
struct mesh *m;
struct behavior *b;
vertex searchpoint;
struct otri *searchtri;
#endif /* not ANSI_DECLARATORS */
{
VOID **sampleblock;
char *firsttri;
struct otri sampletri;
vertex torg, tdest;
unsigned long alignptr;
REAL searchdist, dist;
REAL ahead;
long samplesperblock, totalsamplesleft, samplesleft;
long population, totalpopulation;
triangle ptr; /* Temporary variable used by sym(). */
if (b->verbose > 2) {
printf(" Randomly sampling for a triangle near point (%.12g, %.12g).n",
searchpoint[0], searchpoint[1]);
}
/* Record the distance from the suggested starting triangle to the */
/* point we seek. */
org(*searchtri, torg);
searchdist = (searchpoint[0] - torg[0]) * (searchpoint[0] - torg[0]) +
(searchpoint[1] - torg[1]) * (searchpoint[1] - torg[1]);
if (b->verbose > 2) {
printf(" Boundary triangle has origin (%.12g, %.12g).n",
torg[0], torg[1]);
}
/* If a recently encountered triangle has been recorded and has not been */
/* deallocated, test it as a good starting point. */
if (m->recenttri.tri != (triangle *) NULL) {
if (!deadtri(m->recenttri.tri)) {
org(m->recenttri, torg);
if ((torg[0] == searchpoint[0]) && (torg[1] == searchpoint[1])) {
otricopy(m->recenttri, *searchtri);
return ONVERTEX;
}
dist = (searchpoint[0] - torg[0]) * (searchpoint[0] - torg[0]) +
(searchpoint[1] - torg[1]) * (searchpoint[1] - torg[1]);
if (dist < searchdist) {
otricopy(m->recenttri, *searchtri);
searchdist = dist;
if (b->verbose > 2) {
printf(" Choosing recent triangle with origin (%.12g, %.12g).n",
torg[0], torg[1]);
}
}
}
}
/* The number of random samples taken is proportional to the cube root of */
/* the number of triangles in the mesh. The next bit of code assumes */
/* that the number of triangles increases monotonically (or at least */
/* doesn't decrease enough to matter). */
while (SAMPLEFACTOR * m->samples * m->samples * m->samples <
m->triangles.items) {
m->samples++;
}
/* We'll draw ceiling(samples * TRIPERBLOCK / maxitems) random samples */
/* from each block of triangles (except the first)--until we meet the */
/* sample quota. The ceiling means that blocks at the end might be */
/* neglected, but I don't care. */
samplesperblock = (m->samples * TRIPERBLOCK - 1) / m->triangles.maxitems + 1;
/* We'll draw ceiling(samples * itemsfirstblock / maxitems) random samples */
/* from the first block of triangles. */
samplesleft = (m->samples * m->triangles.itemsfirstblock - 1) /
m->triangles.maxitems + 1;
totalsamplesleft = m->samples;
population = m->triangles.itemsfirstblock;
totalpopulation = m->triangles.maxitems;
sampleblock = m->triangles.firstblock;
sampletri.orient = 0;
while (totalsamplesleft > 0) {
/* If we're in the last block, `population' needs to be corrected. */
if (population > totalpopulation) {
population = totalpopulation;
}
/* Find a pointer to the first triangle in the block. */
alignptr = (unsigned long) (sampleblock + 1);
firsttri = (char *) (alignptr +
(unsigned long) m->triangles.alignbytes -
(alignptr %
(unsigned long) m->triangles.alignbytes));
/* Choose `samplesleft' randomly sampled triangles in this block. */
do {
sampletri.tri = (triangle *) (firsttri +
(randomnation((unsigned int) population) *
m->triangles.itembytes));
if (!deadtri(sampletri.tri)) {
org(sampletri, torg);
dist = (searchpoint[0] - torg[0]) * (searchpoint[0] - torg[0]) +
(searchpoint[1] - torg[1]) * (searchpoint[1] - torg[1]);
if (dist < searchdist) {
otricopy(sampletri, *searchtri);
searchdist = dist;
if (b->verbose > 2) {
printf(" Choosing triangle with origin (%.12g, %.12g).n",
torg[0], torg[1]);
}
}
}
samplesleft--;
totalsamplesleft--;
} while ((samplesleft > 0) && (totalsamplesleft > 0));
if (totalsamplesleft > 0) {
sampleblock = (VOID **) *sampleblock;
samplesleft = samplesperblock;
totalpopulation -= population;
population = TRIPERBLOCK;
}
}
/* Where are we? */
org(*searchtri, torg);
dest(*searchtri, tdest);
/* Check the starting triangle's vertices. */
if ((torg[0] == searchpoint[0]) && (torg[1] == searchpoint[1])) {
return ONVERTEX;
}
if ((tdest[0] == searchpoint[0]) && (tdest[1] == searchpoint[1])) {
lnextself(*searchtri);
return ONVERTEX;
}
/* Orient `searchtri' to fit the preconditions of calling preciselocate(). */
ahead = counterclockwise(m, b, torg, tdest, searchpoint);
if (ahead < 0.0) {
/* Turn around so that `searchpoint' is to the left of the */
/* edge specified by `searchtri'. */
symself(*searchtri);
} else if (ahead == 0.0) {
/* Check if `searchpoint' is between `torg' and `tdest'. */
if (((torg[0] < searchpoint[0]) == (searchpoint[0] < tdest[0])) &&
((torg[1] < searchpoint[1]) == (searchpoint[1] < tdest[1]))) {
return ONEDGE;
}
}
return preciselocate(m, b, searchpoint, searchtri, 0);
}
/** **/
/** **/
/********* Point location routines end here *********/
/********* Mesh transformation routines begin here *********/
/** **/
/** **/
/*****************************************************************************/
/* */
/* insertsubseg() Create a new subsegment and insert it between two */
/* triangles. */
/* */
/* The new subsegment is inserted at the edge described by the handle */
/* `tri'. Its vertices are properly initialized. The marker `subsegmark' */
/* is applied to the subsegment and, if appropriate, its vertices. */
/* */
/*****************************************************************************/
#ifdef ANSI_DECLARATORS
void insertsubseg(struct mesh *m, struct behavior *b, struct otri *tri,
int subsegmark)
#else /* not ANSI_DECLARATORS */
void insertsubseg(m, b, tri, subsegmark)
struct mesh *m;
struct behavior *b;
struct otri *tri; /* Edge at which to insert the new subsegment. */
int subsegmark; /* Marker for the new subsegment. */
#endif /* not ANSI_DECLARATORS */
{
struct otri oppotri;
struct osub newsubseg;
vertex triorg, tridest;
triangle ptr; /* Temporary variable used by sym(). */
subseg sptr; /* Temporary variable used by tspivot(). */
org(*tri, triorg);
dest(*tri, tridest);
/* Mark vertices if possible. */
if (vertexmark(triorg) == 0) {
setvertexmark(triorg, subsegmark);
}
if (vertexmark(tridest) == 0) {
setvertexmark(tridest, subsegmark);
}
/* Check if there's already a subsegment here. */
tspivot(*tri, newsubseg);
if (newsubseg.ss == m->dummysub) {
/* Make new subsegment and initialize its vertices. */
makesubseg(m, &newsubseg);
setsorg(newsubseg, tridest);
setsdest(newsubseg, triorg);
setsegorg(newsubseg, tridest);
setsegdest(newsubseg, triorg);
/* Bond new subsegment to the two triangles it is sandwiched between. */
/* Note that the facing triangle `oppotri' might be equal to */
/* `dummytri' (outer space), but the new subsegment is bonded to it */
/* all the same. */
tsbond(*tri, newsubseg);
sym(*tri, oppotri);
ssymself(newsubseg);
tsbond(oppotri, newsubseg);
setmark(newsubseg, subsegmark);
if (b->verbose > 2) {
printf(" Inserting new ");
printsubseg(m, b, &newsubseg);
}
} else {
if (mark(newsubseg) == 0) {
setmark(newsubseg, subsegmark);
}
}
}
/*****************************************************************************/
/* */
/* Terminology */
/* */
/* A "local transformation" replaces a small set of triangles with another */
/* set of triangles. This may or may not involve inserting or deleting a */
/* vertex. */
/* */
/* The term "casing" is used to describe the set of triangles that are */
/* attached to the triangles being transformed, but are not transformed */
/* themselves. Think of the casing as a fixed hollow structure inside */
/* which all the action happens. A "casing" is only defined relative to */
/* a single transformation; each occurrence of a transformation will */
/* involve a different casing. */
/* */
/*****************************************************************************/
/*****************************************************************************/
/* */
/* flip() Transform two triangles to two different triangles by flipping */
/* an edge counterclockwise within a quadrilateral. */
/* */
/* Imagine the original triangles, abc and bad, oriented so that the */
/* shared edge ab lies in a horizontal plane, with the vertex b on the left */
/* and the vertex a on the right. The vertex c lies below the edge, and */
/* the vertex d lies above the edge. The `flipedge' handle holds the edge */
/* ab of triangle abc, and is directed left, from vertex a to vertex b. */
/* */
/* The triangles abc and bad are deleted and replaced by the triangles cdb */
/* and dca. The triangles that represent abc and bad are NOT deallocated; */
/* they are reused for dca and cdb, respectively. Hence, any handles that */
/* may have held the original triangles are still valid, although not */
/* directed as they were before. */
/* */
/* Upon completion of this routine, the `flipedge' handle holds the edge */
/* dc of triangle dca, and is directed down, from vertex d to vertex c. */
/* (Hence, the two triangles have rotated counterclockwise.) */
/* */
/* WARNING: This transformation is geometrically valid only if the */
/* quadrilateral adbc is convex. Furthermore, this transformation is */
/* valid only if there is not a subsegment between the triangles abc and */
/* bad. This routine does not check either of these preconditions, and */
/* it is the responsibility of the calling routine to ensure that they are */
/* met. If they are not, the streets shall be filled with wailing and */
/* gnashing of teeth. */
/* */
/*****************************************************************************/
#ifdef ANSI_DECLARATORS
void flip(struct mesh *m, struct behavior *b, struct otri *flipedge)
#else /* not ANSI_DECLARATORS */
void flip(m, b, flipedge)
struct mesh *m;
struct behavior *b;
struct otri *flipedge; /* Handle for the triangle abc. */
#endif /* not ANSI_DECLARATORS */
{
struct otri botleft, botright;
struct otri topleft, topright;
struct otri top;
struct otri botlcasing, botrcasing;
struct otri toplcasing, toprcasing;
struct osub botlsubseg, botrsubseg;
struct osub toplsubseg, toprsubseg;
vertex leftvertex, rightvertex, botvertex;
vertex farvertex;
triangle ptr; /* Temporary variable used by sym(). */
subseg sptr; /* Temporary variable used by tspivot(). */
/* Identify the vertices of the quadrilateral. */
org(*flipedge, rightvertex);
dest(*flipedge, leftvertex);
apex(*flipedge, botvertex);
sym(*flipedge, top);
#ifdef SELF_CHECK
if (top.tri == m->dummytri) {
printf("Internal error in flip(): Attempt to flip on boundary.n");
lnextself(*flipedge);
return;
}
if (m->checksegments) {
tspivot(*flipedge, toplsubseg);
if (toplsubseg.ss != m->dummysub) {
printf("Internal error in flip(): Attempt to flip a segment.n");
lnextself(*flipedge);
return;
}
}
#endif /* SELF_CHECK */
apex(top, farvertex);
/* Identify the casing of the quadrilateral. */
lprev(top, topleft);
sym(topleft, toplcasing);
lnext(top, topright);
sym(topright, toprcasing);
lnext(*flipedge, botleft);
sym(botleft, botlcasing);
lprev(*flipedge, botright);
sym(botright, botrcasing);
/* Rotate the quadrilateral one-quarter turn counterclockwise. */
bond(topleft, botlcasing);
bond(botleft, botrcasing);
bond(botright, toprcasing);
bond(topright, toplcasing);
if (m->checksegments) {
/* Check for subsegments and rebond them to the quadrilateral. */
tspivot(topleft, toplsubseg);
tspivot(botleft, botlsubseg);
tspivot(botright, botrsubseg);
tspivot(topright, toprsubseg);
if (toplsubseg.ss == m->dummysub) {
tsdissolve(topright);
} else {
tsbond(topright, toplsubseg);
}
if (botlsubseg.ss == m->dummysub) {
tsdissolve(topleft);
} else {
tsbond(topleft, botlsubseg);
}
if (botrsubseg.ss == m->dummysub) {
tsdissolve(botleft);
} else {
tsbond(botleft, botrsubseg);
}
if (toprsubseg.ss == m->dummysub) {
tsdissolve(botright);
} else {
tsbond(botright, toprsubseg);
}
}
/* New vertex assignments for the rotated quadrilateral. */
setorg(*flipedge, farvertex);
setdest(*flipedge, botvertex);
setapex(*flipedge, rightvertex);
setorg(top, botvertex);
setdest(top, farvertex);
setapex(top, leftvertex);
if (b->verbose > 2) {
printf(" Edge flip results in left ");
printtriangle(m, b, &top);
printf(" and right ");
printtriangle(m, b, flipedge);
}
}
/*****************************************************************************/
/* */
/* unflip() Transform two triangles to two different triangles by */
/* flipping an edge clockwise within a quadrilateral. Reverses */
/* the flip() operation so that the data structures representing */
/* the triangles are back where they were before the flip(). */
/* */
/* Imagine the original triangles, abc and bad, oriented so that the */
/* shared edge ab lies in a horizontal plane, with the vertex b on the left */
/* and the vertex a on the right. The vertex c lies below the edge, and */
/* the vertex d lies above the edge. The `flipedge' handle holds the edge */
/* ab of triangle abc, and is directed left, from vertex a to vertex b. */
/* */
/* The triangles abc and bad are deleted and replaced by the triangles cdb */
/* and dca. The triangles that represent abc and bad are NOT deallocated; */
/* they are reused for cdb and dca, respectively. Hence, any handles that */
/* may have held the original triangles are still valid, although not */
/* directed as they were before. */
/* */
/* Upon completion of this routine, the `flipedge' handle holds the edge */
/* cd of triangle cdb, and is directed up, from vertex c to vertex d. */
/* (Hence, the two triangles have rotated clockwise.) */
/* */
/* WARNING: This transformation is geometrically valid only if the */
/* quadrilateral adbc is convex. Furthermore, this transformation is */
/* valid only if there is not a subsegment between the triangles abc and */
/* bad. This routine does not check either of these preconditions, and */
/* it is the responsibility of the calling routine to ensure that they are */
/* met. If they are not, the streets shall be filled with wailing and */
/* gnashing of teeth. */
/* */
/*****************************************************************************/
#ifdef ANSI_DECLARATORS
void unflip(struct mesh *m, struct behavior *b, struct otri *flipedge)
#else /* not ANSI_DECLARATORS */
void unflip(m, b, flipedge)
struct mesh *m;
struct behavior *b;
struct otri *flipedge; /* Handle for the triangle abc. */
#endif /* not ANSI_DECLARATORS */
{
struct otri botleft, botright;
struct otri topleft, topright;
struct otri top;
struct otri botlcasing, botrcasing;
struct otri toplcasing, toprcasing;
struct osub botlsubseg, botrsubseg;
struct osub toplsubseg, toprsubseg;
vertex leftvertex, rightvertex, botvertex;
vertex farvertex;
triangle ptr; /* Temporary variable used by sym(). */
subseg sptr; /* Temporary variable used by tspivot(). */
/* Identify the vertices of the quadrilateral. */
org(*flipedge, rightvertex);
dest(*flipedge, leftvertex);
apex(*flipedge, botvertex);
sym(*flipedge, top);
#ifdef SELF_CHECK
if (top.tri == m->dummytri) {
printf("Internal error in unflip(): Attempt to flip on boundary.n");
lnextself(*flipedge);
return;
}
if (m->checksegments) {
tspivot(*flipedge, toplsubseg);
if (toplsubseg.ss != m->dummysub) {
printf("Internal error in unflip(): Attempt to flip a subsegment.n");
lnextself(*flipedge);
return;
}
}
#endif /* SELF_CHECK */
apex(top, farvertex);
/* Identify the casing of the quadrilateral. */
lprev(top, topleft);
sym(topleft, toplcasing);
lnext(top, topright);
sym(topright, toprcasing);
lnext(*flipedge, botleft);
sym(botleft, botlcasing);
lprev(*flipedge, botright);
sym(botright, botrcasing);
/* Rotate the quadrilateral one-quarter turn clockwise. */
bond(topleft, toprcasing);
bond(botleft, toplcasing);
bond(botright, botlcasing);
bond(topright, botrcasing);
if (m->checksegments) {
/* Check for subsegments and rebond them to the quadrilateral. */
tspivot(topleft, toplsubseg);
tspivot(botleft, botlsubseg);
tspivot(botright, botrsubseg);
tspivot(topright, toprsubseg);
if (toplsubseg.ss == m->dummysub) {
tsdissolve(botleft);
} else {
tsbond(botleft, toplsubseg);
}
if (botlsubseg.ss == m->dummysub) {
tsdissolve(botright);
} else {
tsbond(botright, botlsubseg);
}
if (botrsubseg.ss == m->dummysub) {
tsdissolve(topright);
} else {
tsbond(topright, botrsubseg);
}
if (toprsubseg.ss == m->dummysub) {
tsdissolve(topleft);
} else {
tsbond(topleft, toprsubseg);
}
}
/* New vertex assignments for the rotated quadrilateral. */
setorg(*flipedge, botvertex);
setdest(*flipedge, farvertex);
setapex(*flipedge, leftvertex);
setorg(top, farvertex);
setdest(top, botvertex);
setapex(top, rightvertex);
if (b->verbose > 2) {
printf(" Edge unflip results in left ");
printtriangle(m, b, flipedge);
printf(" and right ");
printtriangle(m, b, &top);
}
}
/*****************************************************************************/
/* */
/* insertvertex() Insert a vertex into a Delaunay triangulation, */
/* performing flips as necessary to maintain the Delaunay */
/* property. */
/* */
/* The point `insertvertex' is located. If `searchtri.tri' is not NULL, */
/* the search for the containing triangle begins from `searchtri'. If */
/* `searchtri.tri' is NULL, a full point location procedure is called. */
/* If `insertvertex' is found inside a triangle, the triangle is split into */
/* three; if `insertvertex' lies on an edge, the edge is split in two, */
/* thereby splitting the two adjacent triangles into four. Edge flips are */
/* used to restore the Delaunay property. If `insertvertex' lies on an */
/* existing vertex, no action is taken, and the value DUPLICATEVERTEX is */
/* returned. On return, `searchtri' is set to a handle whose origin is the */
/* existing vertex. */
/* */
/* Normally, the parameter `splitseg' is set to NULL, implying that no */
/* subsegment should be split. In this case, if `insertvertex' is found to */
/* lie on a segment, no action is taken, and the value VIOLATINGVERTEX is */
/* returned. On return, `searchtri' is set to a handle whose primary edge */
/* is the violated subsegment. */
/* */
/* If the calling routine wishes to split a subsegment by inserting a */
/* vertex in it, the parameter `splitseg' should be that subsegment. In */
/* this case, `searchtri' MUST be the triangle handle reached by pivoting */
/* from that subsegment; no point location is done. */
/* */
/* `segmentflaws' and `triflaws' are flags that indicate whether or not */
/* there should be checks for the creation of encroached subsegments or bad */
/* quality triangles. If a newly inserted vertex encroaches upon */
/* subsegments, these subsegments are added to the list of subsegments to */
/* be split if `segmentflaws' is set. If bad triangles are created, these */
/* are added to the queue if `triflaws' is set. */
/* */
/* If a duplicate vertex or violated segment does not prevent the vertex */
/* from being inserted, the return value will be ENCROACHINGVERTEX if the */
/* vertex encroaches upon a subsegment (and checking is enabled), or */
/* SUCCESSFULVERTEX otherwise. In either case, `searchtri' is set to a */
/* handle whose origin is the newly inserted vertex. */
/* */
/* insertvertex() does not use flip() for reasons of speed; some */
/* information can be reused from edge flip to edge flip, like the */
/* locations of subsegments. */
/* */
/*****************************************************************************/
#ifdef ANSI_DECLARATORS
enum insertvertexresult insertvertex(struct mesh *m, struct behavior *b,
vertex newvertex, struct otri *searchtri,
struct osub *splitseg,
int segmentflaws, int triflaws)
#else /* not ANSI_DECLARATORS */
enum insertvertexresult insertvertex(m, b, newvertex, searchtri, splitseg,
segmentflaws, triflaws)
struct mesh *m;
struct behavior *b;
vertex newvertex;
struct otri *searchtri;
struct osub *splitseg;
int segmentflaws;
int triflaws;
#endif /* not ANSI_DECLARATORS */
{
struct otri horiz;
struct otri top;
struct otri botleft, botright;
struct otri topleft, topright;
struct otri newbotleft, newbotright;
struct otri newtopright;
struct otri botlcasing, botrcasing;
struct otri toplcasing, toprcasing;
struct otri testtri;
struct osub botlsubseg, botrsubseg;
struct osub toplsubseg, toprsubseg;
struct osub brokensubseg;
struct osub checksubseg;
struct osub rightsubseg;
struct osub newsubseg;
struct badsubseg *encroached;
struct flipstacker *newflip;
vertex first;
vertex leftvertex, rightvertex, botvertex, topvertex, farvertex;
vertex segmentorg, segmentdest;
REAL attrib;
REAL area;
enum insertvertexresult success;
enum locateresult intersect;
int doflip;
int mirrorflag;
int enq;
int i;
triangle ptr; /* Temporary variable used by sym(). */
subseg sptr; /* Temporary variable used by spivot() and tspivot(). */
if (b->verbose > 1) {
printf(" Inserting (%.12g, %.12g).n", newvertex[0], newvertex[1]);
}
if (splitseg == (struct osub *) NULL) {
/* Find the location of the vertex to be inserted. Check if a good */
/* starting triangle has already been provided by the caller. */
if (searchtri->tri == m->dummytri) {
/* Find a boundary triangle. */
horiz.tri = m->dummytri;
horiz.orient = 0;
symself(horiz);
/* Search for a triangle containing `newvertex'. */
intersect = locate(m, b, newvertex, &horiz);
} else {
/* Start searching from the triangle provided by the caller. */
otricopy(*searchtri, horiz);
intersect = preciselocate(m, b, newvertex, &horiz, 1);
}
} else {
/* The calling routine provides the subsegment in which */
/* the vertex is inserted. */
otricopy(*searchtri, horiz);
intersect = ONEDGE;
}
if (intersect == ONVERTEX) {
/* There's already a vertex there. Return in `searchtri' a triangle */
/* whose origin is the existing vertex. */
otricopy(horiz, *searchtri);
otricopy(horiz, m->recenttri);
return DUPLICATEVERTEX;
}
if ((intersect == ONEDGE) || (intersect == OUTSIDE)) {
/* The vertex falls on an edge or boundary. */
if (m->checksegments && (splitseg == (struct osub *) NULL)) {
/* Check whether the vertex falls on a subsegment. */
tspivot(horiz, brokensubseg);
if (brokensubseg.ss != m->dummysub) {
/* The vertex falls on a subsegment, and hence will not be inserted. */
if (segmentflaws) {
enq = b->nobisect != 2;
if (enq && (b->nobisect == 1)) {
/* This subsegment may be split only if it is an */
/* internal boundary. */
sym(horiz, testtri);
enq = testtri.tri != m->dummytri;
}
if (enq) {
/* Add the subsegment to the list of encroached subsegments. */
encroached = (struct badsubseg *) poolalloc(&m->badsubsegs);
encroached->encsubseg = sencode(brokensubseg);
sorg(brokensubseg, encroached->subsegorg);
sdest(brokensubseg, encroached->subsegdest);
if (b->verbose > 2) {
printf(
" Queueing encroached subsegment (%.12g, %.12g) (%.12g, %.12g).n",
encroached->subsegorg[0], encroached->subsegorg[1],
encroached->subsegdest[0], encroached->subsegdest[1]);
}
}
}
/* Return a handle whose primary edge contains the vertex, */
/* which has not been inserted. */
otricopy(horiz, *searchtri);
otricopy(horiz, m->recenttri);
return VIOLATINGVERTEX;
}
}
/* Insert the vertex on an edge, dividing one triangle into two (if */
/* the edge lies on a boundary) or two triangles into four. */
lprev(horiz, botright);
sym(botright, botrcasing);
sym(horiz, topright);
/* Is there a second triangle? (Or does this edge lie on a boundary?) */
mirrorflag = topright.tri != m->dummytri;
if (mirrorflag) {
lnextself(topright);
sym(topright, toprcasing);
maketriangle(m, b, &newtopright);
} else {
/* Splitting a boundary edge increases the number of boundary edges. */
m->hullsize++;
}
maketriangle(m, b, &newbotright);
/* Set the vertices of changed and new triangles. */
org(horiz, rightvertex);
dest(horiz, leftvertex);
apex(horiz, botvertex);
setorg(newbotright, botvertex);
setdest(newbotright, rightvertex);
setapex(newbotright, newvertex);
setorg(horiz, newvertex);
for (i = 0; i < m->eextras; i++) {
/* Set the element attributes of a new triangle. */
setelemattribute(newbotright, i, elemattribute(botright, i));
}
if (b->vararea) {
/* Set the area constraint of a new triangle. */
setareabound(newbotright, areabound(botright));
}
if (mirrorflag) {
dest(topright, topvertex);
setorg(newtopright, rightvertex);
setdest(newtopright, topvertex);
setapex(newtopright, newvertex);
setorg(topright, newvertex);
for (i = 0; i < m->eextras; i++) {
/* Set the element attributes of another new triangle. */
setelemattribute(newtopright, i, elemattribute(topright, i));
}
if (b->vararea) {
/* Set the area constraint of another new triangle. */
setareabound(newtopright, areabound(topright));
}
}
/* There may be subsegments that need to be bonded */
/* to the new triangle(s). */
if (m->checksegments) {
tspivot(botright, botrsubseg);
if (botrsubseg.ss != m->dummysub) {
tsdissolve(botright);
tsbond(newbotright, botrsubseg);
}
if (mirrorflag) {
tspivot(topright, toprsubseg);
if (toprsubseg.ss != m->dummysub) {
tsdissolve(topright);
tsbond(newtopright, toprsubseg);
}
}
}
/* Bond the new triangle(s) to the surrounding triangles. */
bond(newbotright, botrcasing);
lprevself(newbotright);
bond(newbotright, botright);
lprevself(newbotright);
if (mirrorflag) {
bond(newtopright, toprcasing);
lnextself(newtopright);
bond(newtopright, topright);
lnextself(newtopright);
bond(newtopright, newbotright);
}
if (splitseg != (struct osub *) NULL) {
/* Split the subsegment into two. */
setsdest(*splitseg, newvertex);
segorg(*splitseg, segmentorg);
segdest(*splitseg, segmentdest);
ssymself(*splitseg);
spivot(*splitseg, rightsubseg);
insertsubseg(m, b, &newbotright, mark(*splitseg));
tspivot(newbotright, newsubseg);
setsegorg(newsubseg, segmentorg);
setsegdest(newsubseg, segmentdest);
sbond(*splitseg, newsubseg);
ssymself(newsubseg);
sbond(newsubseg, rightsubseg);
ssymself(*splitseg);
/* Transfer the subsegment's boundary marker to the vertex */
/* if required. */
if (vertexmark(newvertex) == 0) {
setvertexmark(newvertex, mark(*splitseg));
}
}
if (m->checkquality) {
poolrestart(&m->flipstackers);
m->lastflip = (struct flipstacker *) poolalloc(&m->flipstackers);
m->lastflip->flippedtri = encode(horiz);
m->lastflip->prevflip = (struct flipstacker *) &insertvertex;
}
#ifdef SELF_CHECK
if (counterclockwise(m, b, rightvertex, leftvertex, botvertex) < 0.0) {
printf("Internal error in insertvertex():n");
printf(
" Clockwise triangle prior to edge vertex insertion (bottom).n");
}
if (mirrorflag) {
if (counterclockwise(m, b, leftvertex, rightvertex, topvertex) < 0.0) {
printf("Internal error in insertvertex():n");
printf(" Clockwise triangle prior to edge vertex insertion (top).n");
}
if (counterclockwise(m, b, rightvertex, topvertex, newvertex) < 0.0) {
printf("Internal error in insertvertex():n");
printf(
" Clockwise triangle after edge vertex insertion (top right).n");
}
if (counterclockwise(m, b, topvertex, leftvertex, newvertex) < 0.0) {
printf("Internal error in insertvertex():n");
printf(
" Clockwise triangle after edge vertex insertion (top left).n");
}
}
if (counterclockwise(m, b, leftvertex, botvertex, newvertex) < 0.0) {
printf("Internal error in insertvertex():n");
printf(
" Clockwise triangle after edge vertex insertion (bottom left).n");
}
if (counterclockwise(m, b, botvertex, rightvertex, newvertex) < 0.0) {
printf("Internal error in insertvertex():n");
printf(
" Clockwise triangle after edge vertex insertion (bottom right).n");
}
#endif /* SELF_CHECK */
if (b->verbose > 2) {
printf(" Updating bottom left ");
printtriangle(m, b, &botright);
if (mirrorflag) {
printf(" Updating top left ");
printtriangle(m, b, &topright);
printf(" Creating top right ");
printtriangle(m, b, &newtopright);
}
printf(" Creating bottom right ");
printtriangle(m, b, &newbotright);
}
/* Position `horiz' on the first edge to check for */
/* the Delaunay property. */
lnextself(horiz);
} else {
/* Insert the vertex in a triangle, splitting it into three. */
lnext(horiz, botleft);
lprev(horiz, botright);
sym(botleft, botlcasing);
sym(botright, botrcasing);
maketriangle(m, b, &newbotleft);
maketriangle(m, b, &newbotright);
/* Set the vertices of changed and new triangles. */
org(horiz, rightvertex);
dest(horiz, leftvertex);
apex(horiz, botvertex);
setorg(newbotleft, leftvertex);
setdest(newbotleft, botvertex);
setapex(newbotleft, newvertex);
setorg(newbotright, botvertex);
setdest(newbotright, rightvertex);
setapex(newbotright, newvertex);
setapex(horiz, newvertex);
for (i = 0; i < m->eextras; i++) {
/* Set the element attributes of the new triangles. */
attrib = elemattribute(horiz, i);
setelemattribute(newbotleft, i, attrib);
setelemattribute(newbotright, i, attrib);
}
if (b->vararea) {
/* Set the area constraint of the new triangles. */
area = areabound(horiz);
setareabound(newbotleft, area);
setareabound(newbotright, area);
}
/* There may be subsegments that need to be bonded */
/* to the new triangles. */
if (m->checksegments) {
tspivot(botleft, botlsubseg);
if (botlsubseg.ss != m->dummysub) {
tsdissolve(botleft);
tsbond(newbotleft, botlsubseg);
}
tspivot(botright, botrsubseg);
if (botrsubseg.ss != m->dummysub) {
tsdissolve(botright);
tsbond(newbotright, botrsubseg);
}
}
/* Bond the new triangles to the surrounding triangles. */
bond(newbotleft, botlcasing);
bond(newbotright, botrcasing);
lnextself(newbotleft);
lprevself(newbotright);
bond(newbotleft, newbotright);
lnextself(newbotleft);
bond(botleft, newbotleft);
lprevself(newbotright);
bond(botright, newbotright);
if (m->checkquality) {
poolrestart(&m->flipstackers);
m->lastflip = (struct flipstacker *) poolalloc(&m->flipstackers);
m->lastflip->flippedtri = encode(horiz);
m->lastflip->prevflip = (struct flipstacker *) NULL;
}
#ifdef SELF_CHECK
if (counterclockwise(m, b, rightvertex, leftvertex, botvertex) < 0.0) {
printf("Internal error in insertvertex():n");
printf(" Clockwise triangle prior to vertex insertion.n");
}
if (counterclockwise(m, b, rightvertex, leftvertex, newvertex) < 0.0) {
printf("Internal error in insertvertex():n");
printf(" Clockwise triangle after vertex insertion (top).n");
}
if (counterclockwise(m, b, leftvertex, botvertex, newvertex) < 0.0) {
printf("Internal error in insertvertex():n");
printf(" Clockwise triangle after vertex insertion (left).n");
}
if (counterclockwise(m, b, botvertex, rightvertex, newvertex) < 0.0) {
printf("Internal error in insertvertex():n");
printf(" Clockwise triangle after vertex insertion (right).n");
}
#endif /* SELF_CHECK */
if (b->verbose > 2) {
printf(" Updating top ");
printtriangle(m, b, &horiz);
printf(" Creating left ");
printtriangle(m, b, &newbotleft);
printf(" Creating right ");
printtriangle(m, b, &newbotright);
}
}
/* The insertion is successful by default, unless an encroached */
/* subsegment is found. */
success = SUCCESSFULVERTEX;
/* Circle around the newly inserted vertex, checking each edge opposite */
/* it for the Delaunay property. Non-Delaunay edges are flipped. */
/* `horiz' is always the edge being checked. `first' marks where to */
/* stop circling. */
org(horiz, first);
rightvertex = first;
dest(horiz, leftvertex);
/* Circle until finished. */
while (1) {
/* By default, the edge will be flipped. */
doflip = 1;
if (m->checksegments) {
/* Check for a subsegment, which cannot be flipped. */
tspivot(horiz, checksubseg);
if (checksubseg.ss != m->dummysub) {
/* The edge is a subsegment and cannot be flipped. */
doflip = 0;
#ifndef CDT_ONLY
if (segmentflaws) {
/* Does the new vertex encroach upon this subsegment? */
if (checkseg4encroach(m, b, &checksubseg)) {
success = ENCROACHINGVERTEX;
}
}
#endif /* not CDT_ONLY */
}
}
if (doflip) {
/* Check if the edge is a boundary edge. */
sym(horiz, top);
if (top.tri == m->dummytri) {
/* The edge is a boundary edge and cannot be flipped. */
doflip = 0;
} else {
/* Find the vertex on the other side of the edge. */
apex(top, farvertex);
/* In the incremental Delaunay triangulation algorithm, any of */
/* `leftvertex', `rightvertex', and `farvertex' could be vertices */
/* of the triangular bounding box. These vertices must be */
/* treated as if they are infinitely distant, even though their */
/* "coordinates" are not. */
if ((leftvertex == m->infvertex1) || (leftvertex == m->infvertex2) ||
(leftvertex == m->infvertex3)) {
/* `leftvertex' is infinitely distant. Check the convexity of */
/* the boundary of the triangulation. 'farvertex' might be */
/* infinite as well, but trust me, this same condition should */
/* be applied. */
doflip = counterclockwise(m, b, newvertex, rightvertex, farvertex)
> 0.0;
} else if ((rightvertex == m->infvertex1) ||
(rightvertex == m->infvertex2) ||
(rightvertex == m->infvertex3)) {
/* `rightvertex' is infinitely distant. Check the convexity of */
/* the boundary of the triangulation. 'farvertex' might be */
/* infinite as well, but trust me, this same condition should */
/* be applied. */
doflip = counterclockwise(m, b, farvertex, leftvertex, newvertex)
> 0.0;
} else if ((farvertex == m->infvertex1) ||
(farvertex == m->infvertex2) ||
(farvertex == m->infvertex3)) {
/* `farvertex' is infinitely distant and cannot be inside */
/* the circumcircle of the triangle `horiz'. */
doflip = 0;
} else {
/* Test whether the edge is locally Delaunay. */
doflip = incircle(m, b, leftvertex, newvertex, rightvertex,
farvertex) > 0.0;
}
if (doflip) {
/* We made it! Flip the edge `horiz' by rotating its containing */
/* quadrilateral (the two triangles adjacent to `horiz'). */
/* Identify the casing of the quadrilateral. */
lprev(top, topleft);
sym(topleft, toplcasing);
lnext(top, topright);
sym(topright, toprcasing);
lnext(horiz, botleft);
sym(botleft, botlcasing);
lprev(horiz, botright);
sym(botright, botrcasing);
/* Rotate the quadrilateral one-quarter turn counterclockwise. */
bond(topleft, botlcasing);
bond(botleft, botrcasing);
bond(botright, toprcasing);
bond(topright, toplcasing);
if (m->checksegments) {
/* Check for subsegments and rebond them to the quadrilateral. */
tspivot(topleft, toplsubseg);
tspivot(botleft, botlsubseg);
tspivot(botright, botrsubseg);
tspivot(topright, toprsubseg);
if (toplsubseg.ss == m->dummysub) {
tsdissolve(topright);
} else {
tsbond(topright, toplsubseg);
}
if (botlsubseg.ss == m->dummysub) {
tsdissolve(topleft);
} else {
tsbond(topleft, botlsubseg);
}
if (botrsubseg.ss == m->dummysub) {
tsdissolve(botleft);
} else {
tsbond(botleft, botrsubseg);
}
if (toprsubseg.ss == m->dummysub) {
tsdissolve(botright);
} else {
tsbond(botright, toprsubseg);
}
}
/* New vertex assignments for the rotated quadrilateral. */
setorg(horiz, farvertex);
setdest(horiz, newvertex);
setapex(horiz, rightvertex);
setorg(top, newvertex);
setdest(top, farvertex);
setapex(top, leftvertex);
for (i = 0; i < m->eextras; i++) {
/* Take the average of the two triangles' attributes. */
attrib = 0.5 * (elemattribute(top, i) + elemattribute(horiz, i));
setelemattribute(top, i, attrib);
setelemattribute(horiz, i, attrib);
}
if (b->vararea) {
if ((areabound(top) <= 0.0) || (areabound(horiz) <= 0.0)) {
area = -1.0;
} else {
/* Take the average of the two triangles' area constraints. */
/* This prevents small area constraints from migrating a */
/* long, long way from their original location due to flips. */
area = 0.5 * (areabound(top) + areabound(horiz));
}
setareabound(top, area);
setareabound(horiz, area);
}
if (m->checkquality) {
newflip = (struct flipstacker *) poolalloc(&m->flipstackers);
newflip->flippedtri = encode(horiz);
newflip->prevflip = m->lastflip;
m->lastflip = newflip;
}
#ifdef SELF_CHECK
if (newvertex != (vertex) NULL) {
if (counterclockwise(m, b, leftvertex, newvertex, rightvertex) <
0.0) {
printf("Internal error in insertvertex():n");
printf(" Clockwise triangle prior to edge flip (bottom).n");
}
/* The following test has been removed because constrainededge() */
/* sometimes generates inverted triangles that insertvertex() */
/* removes. */
/*
if (counterclockwise(m, b, rightvertex, farvertex, leftvertex) <
0.0) {
printf("Internal error in insertvertex():n");
printf(" Clockwise triangle prior to edge flip (top).n");
}
*/
if (counterclockwise(m, b, farvertex, leftvertex, newvertex) <
0.0) {
printf("Internal error in insertvertex():n");
printf(" Clockwise triangle after edge flip (left).n");
}
if (counterclockwise(m, b, newvertex, rightvertex, farvertex) <
0.0) {
printf("Internal error in insertvertex():n");
printf(" Clockwise triangle after edge flip (right).n");
}
}
#endif /* SELF_CHECK */
if (b->verbose > 2) {
printf(" Edge flip results in left ");
lnextself(topleft);
printtriangle(m, b, &topleft);
printf(" and right ");
printtriangle(m, b, &horiz);
}
/* On the next iterations, consider the two edges that were */
/* exposed (this is, are now visible to the newly inserted */
/* vertex) by the edge flip. */
lprevself(horiz);
leftvertex = farvertex;
}
}
}
if (!doflip) {
/* The handle `horiz' is accepted as locally Delaunay. */
#ifndef CDT_ONLY
if (triflaws) {
/* Check the triangle `horiz' for quality. */
testtriangle(m, b, &horiz);
}
#endif /* not CDT_ONLY */
/* Look for the next edge around the newly inserted vertex. */
lnextself(horiz);
sym(horiz, testtri);
/* Check for finishing a complete revolution about the new vertex, or */
/* falling outside of the triangulation. The latter will happen */
/* when a vertex is inserted at a boundary. */
if ((leftvertex == first) || (testtri.tri == m->dummytri)) {
/* We're done. Return a triangle whose origin is the new vertex. */
lnext(horiz, *searchtri);
lnext(horiz, m->recenttri);
return success;
}
/* Finish finding the next edge around the newly inserted vertex. */
lnext(testtri, horiz);
rightvertex = leftvertex;
dest(horiz, leftvertex);
}
}
}
/*****************************************************************************/
/* */
/* triangulatepolygon() Find the Delaunay triangulation of a polygon that */
/* has a certain "nice" shape. This includes the */
/* polygons that result from deletion of a vertex or */
/* insertion of a segment. */
/* */
/* This is a conceptually difficult routine. The starting assumption is */
/* that we have a polygon with n sides. n - 1 of these sides are currently */
/* represented as edges in the mesh. One side, called the "base", need not */
/* be. */
/* */
/* Inside the polygon is a structure I call a "fan", consisting of n - 1 */
/* triangles that share a common origin. For each of these triangles, the */
/* edge opposite the origin is one of the sides of the polygon. The */
/* primary edge of each triangle is the edge directed from the origin to */
/* the destination; note that this is not the same edge that is a side of */
/* the polygon. `firstedge' is the primary edge of the first triangle. */
/* From there, the triangles follow in counterclockwise order about the */
/* polygon, until `lastedge', the primary edge of the last triangle. */
/* `firstedge' and `lastedge' are probably connected to other triangles */
/* beyond the extremes of the fan, but their identity is not important, as */
/* long as the fan remains connected to them. */
/* */
/* Imagine the polygon oriented so that its base is at the bottom. This */
/* puts `firstedge' on the far right, and `lastedge' on the far left. */
/* The right vertex of the base is the destination of `firstedge', and the */
/* left vertex of the base is the apex of `lastedge'. */
/* */
/* The challenge now is to find the right sequence of edge flips to */
/* transform the fan into a Delaunay triangulation of the polygon. Each */
/* edge flip effectively removes one triangle from the fan, committing it */
/* to the polygon. The resulting polygon has one fewer edge. If `doflip' */
/* is set, the final flip will be performed, resulting in a fan of one */
/* (useless?) triangle. If `doflip' is not set, the final flip is not */
/* performed, resulting in a fan of two triangles, and an unfinished */
/* triangular polygon that is not yet filled out with a single triangle. */
/* On completion of the routine, `lastedge' is the last remaining triangle, */
/* or the leftmost of the last two. */
/* */
/* Although the flips are performed in the order described above, the */
/* decisions about what flips to perform are made in precisely the reverse */
/* order. The recursive triangulatepolygon() procedure makes a decision, */
/* uses up to two recursive calls to triangulate the "subproblems" */
/* (polygons with fewer edges), and then performs an edge flip. */
/* */
/* The "decision" it makes is which vertex of the polygon should be */
/* connected to the base. This decision is made by testing every possible */
/* vertex. Once the best vertex is found, the two edges that connect this */
/* vertex to the base become the bases for two smaller polygons. These */
/* are triangulated recursively. Unfortunately, this approach can take */
/* O(n^2) time not only in the worst case, but in many common cases. It's */
/* rarely a big deal for vertex deletion, where n is rarely larger than */
/* ten, but it could be a big deal for segment insertion, especially if */
/* there's a lot of long segments that each cut many triangles. I ought to */
/* code a faster algorithm some day. */
/* */
/* The `edgecount' parameter is the number of sides of the polygon, */
/* including its base. `triflaws' is a flag that determines whether the */
/* new triangles should be tested for quality, and enqueued if they are */
/* bad. */
/* */
/*****************************************************************************/
#ifdef ANSI_DECLARATORS
void triangulatepolygon(struct mesh *m, struct behavior *b,
struct otri *firstedge, struct otri *lastedge,
int edgecount, int doflip, int triflaws)
#else /* not ANSI_DECLARATORS */
void triangulatepolygon(m, b, firstedge, lastedge, edgecount, doflip, triflaws)
struct mesh *m;
struct behavior *b;
struct otri *firstedge;
struct otri *lastedge;
int edgecount;
int doflip;
int triflaws;
#endif /* not ANSI_DECLARATORS */
{
struct otri testtri;
struct otri besttri;
struct otri tempedge;
vertex leftbasevertex, rightbasevertex;
vertex testvertex;
vertex bestvertex;
int bestnumber;
int i;
triangle ptr; /* Temporary variable used by sym(), onext(), and oprev(). */
/* Identify the base vertices. */
apex(*lastedge, leftbasevertex);
dest(*firstedge, rightbasevertex);
if (b->verbose > 2) {
printf(" Triangulating interior polygon at edgen");
printf(" (%.12g, %.12g) (%.12g, %.12g)n", leftbasevertex[0],
leftbasevertex[1], rightbasevertex[0], rightbasevertex[1]);
}
/* Find the best vertex to connect the base to. */
onext(*firstedge, besttri);
dest(besttri, bestvertex);
otricopy(besttri, testtri);
bestnumber = 1;
for (i = 2; i <= edgecount - 2; i++) {
onextself(testtri);
dest(testtri, testvertex);
/* Is this a better vertex? */
if (incircle(m, b, leftbasevertex, rightbasevertex, bestvertex,
testvertex) > 0.0) {
otricopy(testtri, besttri);
bestvertex = testvertex;
bestnumber = i;
}
}
if (b->verbose > 2) {
printf(" Connecting edge to (%.12g, %.12g)n", bestvertex[0],
bestvertex[1]);
}
if (bestnumber > 1) {
/* Recursively triangulate the smaller polygon on the right. */
oprev(besttri, tempedge);
triangulatepolygon(m, b, firstedge, &tempedge, bestnumber + 1, 1,
triflaws);
}
if (bestnumber < edgecount - 2) {
/* Recursively triangulate the smaller polygon on the left. */
sym(besttri, tempedge);
triangulatepolygon(m, b, &besttri, lastedge, edgecount - bestnumber, 1,
triflaws);
/* Find `besttri' again; it may have been lost to edge flips. */
sym(tempedge, besttri);
}
if (doflip) {
/* Do one final edge flip. */
flip(m, b, &besttri);
#ifndef CDT_ONLY
if (triflaws) {
/* Check the quality of the newly committed triangle. */
sym(besttri, testtri);
testtriangle(m, b, &testtri);
}
#endif /* not CDT_ONLY */
}
/* Return the base triangle. */
otricopy(besttri, *lastedge);
}
/*****************************************************************************/
/* */
/* deletevertex() Delete a vertex from a Delaunay triangulation, ensuring */
/* that the triangulation remains Delaunay. */
/* */
/* The origin of `deltri' is deleted. The union of the triangles adjacent */
/* to this vertex is a polygon, for which the Delaunay triangulation is */
/* found. Two triangles are removed from the mesh. */
/* */
/* Only interior vertices that do not lie on segments or boundaries may be */
/* deleted. */
/* */
/*****************************************************************************/
#ifndef CDT_ONLY
#ifdef ANSI_DECLARATORS
void deletevertex(struct mesh *m, struct behavior *b, struct otri *deltri)
#else /* not ANSI_DECLARATORS */
void deletevertex(m, b, deltri)
struct mesh *m;
struct behavior *b;
struct otri *deltri;
#endif /* not ANSI_DECLARATORS */
{
struct otri countingtri;
struct otri firstedge, lastedge;
struct otri deltriright;
struct otri lefttri, righttri;
struct otri leftcasing, rightcasing;
struct osub leftsubseg, rightsubseg;
vertex delvertex;
vertex neworg;
int edgecount;
triangle ptr; /* Temporary variable used by sym(), onext(), and oprev(). */
subseg sptr; /* Temporary variable used by tspivot(). */
org(*deltri, delvertex);
if (b->verbose > 1) {
printf(" Deleting (%.12g, %.12g).n", delvertex[0], delvertex[1]);
}
vertexdealloc(m, delvertex);
/* Count the degree of the vertex being deleted. */
onext(*deltri, countingtri);
edgecount = 1;
while (!otriequal(*deltri, countingtri)) {
#ifdef SELF_CHECK
if (countingtri.tri == m->dummytri) {
printf("Internal error in deletevertex():n");
printf(" Attempt to delete boundary vertex.n");
internalerror();
}
#endif /* SELF_CHECK */
edgecount++;
onextself(countingtri);
}
#ifdef SELF_CHECK
if (edgecount < 3) {
printf("Internal error in deletevertex():n Vertex has degree %d.n",
edgecount);
internalerror();
}
#endif /* SELF_CHECK */
if (edgecount > 3) {
/* Triangulate the polygon defined by the union of all triangles */
/* adjacent to the vertex being deleted. Check the quality of */
/* the resulting triangles. */
onext(*deltri, firstedge);
oprev(*deltri, lastedge);
triangulatepolygon(m, b, &firstedge, &lastedge, edgecount, 0,
!b->nobisect);
}
/* Splice out two triangles. */
lprev(*deltri, deltriright);
dnext(*deltri, lefttri);
sym(lefttri, leftcasing);
oprev(deltriright, righttri);
sym(righttri, rightcasing);
bond(*deltri, leftcasing);
bond(deltriright, rightcasing);
tspivot(lefttri, leftsubseg);
if (leftsubseg.ss != m->dummysub) {
tsbond(*deltri, leftsubseg);
}
tspivot(righttri, rightsubseg);
if (rightsubseg.ss != m->dummysub) {
tsbond(deltriright, rightsubseg);
}
/* Set the new origin of `deltri' and check its quality. */
org(lefttri, neworg);
setorg(*deltri, neworg);
if (!b->nobisect) {
testtriangle(m, b, deltri);
}
/* Delete the two spliced-out triangles. */
triangledealloc(m, lefttri.tri);
triangledealloc(m, righttri.tri);
}
#endif /* not CDT_ONLY */
/*****************************************************************************/
/* */
/* undovertex() Undo the most recent vertex insertion. */
/* */
/* Walks through the list of transformations (flips and a vertex insertion) */
/* in the reverse of the order in which they were done, and undoes them. */
/* The inserted vertex is removed from the triangulation and deallocated. */
/* Two triangles (possibly just one) are also deallocated. */
/* */
/*****************************************************************************/
#ifndef CDT_ONLY
#ifdef ANSI_DECLARATORS
void undovertex(struct mesh *m, struct behavior *b)
#else /* not ANSI_DECLARATORS */
void undovertex(m, b)
struct mesh *m;
struct behavior *b;
#endif /* not ANSI_DECLARATORS */
{
struct otri fliptri;
struct otri botleft, botright, topright;
struct otri botlcasing, botrcasing, toprcasing;
struct otri gluetri;
struct osub botlsubseg, botrsubseg, toprsubseg;
vertex botvertex, rightvertex;
triangle ptr; /* Temporary variable used by sym(). */
subseg sptr; /* Temporary variable used by tspivot(). */
/* Walk through the list of transformations (flips and a vertex insertion) */
/* in the reverse of the order in which they were done, and undo them. */
while (m->lastflip != (struct flipstacker *) NULL) {
/* Find a triangle involved in the last unreversed transformation. */
decode(m->lastflip->flippedtri, fliptri);
/* We are reversing one of three transformations: a trisection of one */
/* triangle into three (by inserting a vertex in the triangle), a */
/* bisection of two triangles into four (by inserting a vertex in an */
/* edge), or an edge flip. */
if (m->lastflip->prevflip == (struct flipstacker *) NULL) {
/* Restore a triangle that was split into three triangles, */
/* so it is again one triangle. */
dprev(fliptri, botleft);
lnextself(botleft);
onext(fliptri, botright);
lprevself(botright);
sym(botleft, botlcasing);
sym(botright, botrcasing);
dest(botleft, botvertex);
setapex(fliptri, botvertex);
lnextself(fliptri);
bond(fliptri, botlcasing);
tspivot(botleft, botlsubseg);
tsbond(fliptri, botlsubseg);
lnextself(fliptri);
bond(fliptri, botrcasing);
tspivot(botright, botrsubseg);
tsbond(fliptri, botrsubseg);
/* Delete the two spliced-out triangles. */
triangledealloc(m, botleft.tri);
triangledealloc(m, botright.tri);
} else if (m->lastflip->prevflip == (struct flipstacker *) &insertvertex) {
/* Restore two triangles that were split into four triangles, */
/* so they are again two triangles. */
lprev(fliptri, gluetri);
sym(gluetri, botright);
lnextself(botright);
sym(botright, botrcasing);
dest(botright, rightvertex);
setorg(fliptri, rightvertex);
bond(gluetri, botrcasing);
tspivot(botright, botrsubseg);
tsbond(gluetri, botrsubseg);
/* Delete the spliced-out triangle. */
triangledealloc(m, botright.tri);
sym(fliptri, gluetri);
if (gluetri.tri != m->dummytri) {
lnextself(gluetri);
dnext(gluetri, topright);
sym(topright, toprcasing);
setorg(gluetri, rightvertex);
bond(gluetri, toprcasing);
tspivot(topright, toprsubseg);
tsbond(gluetri, toprsubseg);
/* Delete the spliced-out triangle. */
triangledealloc(m, topright.tri);
}
/* This is the end of the list, sneakily encoded. */
m->lastflip->prevflip = (struct flipstacker *) NULL;
} else {
/* Undo an edge flip. */
unflip(m, b, &fliptri);
}
/* Go on and process the next transformation. */
m->lastflip = m->lastflip->prevflip;
}
}
#endif /* not CDT_ONLY */
/** **/
/** **/
/********* Mesh transformation routines end here *********/
/********* Divide-and-conquer Delaunay triangulation begins here *********/
/** **/
/** **/
/*****************************************************************************/
/* */
/* The divide-and-conquer bounding box */
/* */
/* I originally implemented the divide-and-conquer and incremental Delaunay */
/* triangulations using the edge-based data structure presented by Guibas */
/* and Stolfi. Switching to a triangle-based data structure doubled the */
/* speed. However, I had to think of a few extra tricks to maintain the */
/* elegance of the original algorithms. */
/* */
/* The "bounding box" used by my variant of the divide-and-conquer */
/* algorithm uses one triangle for each edge of the convex hull of the */
/* triangulation. These bounding triangles all share a common apical */
/* vertex, which is represented by NULL and which represents nothing. */
/* The bounding triangles are linked in a circular fan about this NULL */
/* vertex, and the edges on the convex hull of the triangulation appear */
/* opposite the NULL vertex. You might find it easiest to imagine that */
/* the NULL vertex is a point in 3D space behind the center of the */
/* triangulation, and that the bounding triangles form a sort of cone. */
/* */
/* This bounding box makes it easy to represent degenerate cases. For */
/* instance, the triangulation of two vertices is a single edge. This edge */
/* is represented by two bounding box triangles, one on each "side" of the */
/* edge. These triangles are also linked together in a fan about the NULL */
/* vertex. */
/* */
/* The bounding box also makes it easy to traverse the convex hull, as the */
/* divide-and-conquer algorithm needs to do. */
/* */
/*****************************************************************************/
/*****************************************************************************/
/* */
/* vertexsort() Sort an array of vertices by x-coordinate, using the */
/* y-coordinate as a secondary key. */
/* */
/* Uses quicksort. Randomized O(n log n) time. No, I did not make any of */
/* the usual quicksort mistakes. */
/* */
/*****************************************************************************/
#ifdef ANSI_DECLARATORS
void vertexsort(vertex *sortarray, int arraysize)
#else /* not ANSI_DECLARATORS */
void vertexsort(sortarray, arraysize)
vertex *sortarray;
int arraysize;
#endif /* not ANSI_DECLARATORS */
{
int left, right;
int pivot;
REAL pivotx, pivoty;
vertex temp;
if (arraysize == 2) {
/* Recursive base case. */
if ((sortarray[0][0] > sortarray[1][0]) ||
((sortarray[0][0] == sortarray[1][0]) &&
(sortarray[0][1] > sortarray[1][1]))) {
temp = sortarray[1];
sortarray[1] = sortarray[0];
sortarray[0] = temp;
}
return;
}
/* Choose a random pivot to split the array. */
pivot = (int) randomnation((unsigned int) arraysize);
pivotx = sortarray[pivot][0];
pivoty = sortarray[pivot][1];
/* Split the array. */
left = -1;
right = arraysize;
while (left < right) {
/* Search for a vertex whose x-coordinate is too large for the left. */
do {
left++;
} while ((left <= right) && ((sortarray[left][0] < pivotx) ||
((sortarray[left][0] == pivotx) &&
(sortarray[left][1] < pivoty))));
/* Search for a vertex whose x-coordinate is too small for the right. */
do {
right--;
} while ((left <= right) && ((sortarray[right][0] > pivotx) ||
((sortarray[right][0] == pivotx) &&
(sortarray[right][1] > pivoty))));
if (left < right) {
/* Swap the left and right vertices. */
temp = sortarray[left];
sortarray[left] = sortarray[right];
sortarray[right] = temp;
}
}
if (left > 1) {
/* Recursively sort the left subset. */
vertexsort(sortarray, left);
}
if (right < arraysize - 2) {
/* Recursively sort the right subset. */
vertexsort(&sortarray[right + 1], arraysize - right - 1);
}
}
/*****************************************************************************/
/* */
/* vertexmedian() An order statistic algorithm, almost. Shuffles an */
/* array of vertices so that the first `median' vertices */
/* occur lexicographically before the remaining vertices. */
/* */
/* Uses the x-coordinate as the primary key if axis == 0; the y-coordinate */
/* if axis == 1. Very similar to the vertexsort() procedure, but runs in */
/* randomized linear time. */
/* */
/*****************************************************************************/
#ifdef ANSI_DECLARATORS
void vertexmedian(vertex *sortarray, int arraysize, int median, int axis)
#else /* not ANSI_DECLARATORS */
void vertexmedian(sortarray, arraysize, median, axis)
vertex *sortarray;
int arraysize;
int median;
int axis;
#endif /* not ANSI_DECLARATORS */
{
int left, right;
int pivot;
REAL pivot1, pivot2;
vertex temp;
if (arraysize == 2) {
/* Recursive base case. */
if ((sortarray[0][axis] > sortarray[1][axis]) ||
((sortarray[0][axis] == sortarray[1][axis]) &&
(sortarray[0][1 - axis] > sortarray[1][1 - axis]))) {
temp = sortarray[1];
sortarray[1] = sortarray[0];
sortarray[0] = temp;
}
return;
}
/* Choose a random pivot to split the array. */
pivot = (int) randomnation((unsigned int) arraysize);
pivot1 = sortarray[pivot][axis];
pivot2 = sortarray[pivot][1 - axis];
/* Split the array. */
left = -1;
right = arraysize;
while (left < right) {
/* Search for a vertex whose x-coordinate is too large for the left. */
do {
left++;
} while ((left <= right) && ((sortarray[left][axis] < pivot1) ||
((sortarray[left][axis] == pivot1) &&
(sortarray[left][1 - axis] < pivot2))));
/* Search for a vertex whose x-coordinate is too small for the right. */
do {
right--;
} while ((left <= right) && ((sortarray[right][axis] > pivot1) ||
((sortarray[right][axis] == pivot1) &&
(sortarray[right][1 - axis] > pivot2))));
if (left < right) {
/* Swap the left and right vertices. */
temp = sortarray[left];
sortarray[left] = sortarray[right];
sortarray[right] = temp;
}
}
/* Unlike in vertexsort(), at most one of the following */
/* conditionals is true. */
if (left > median) {
/* Recursively shuffle the left subset. */
vertexmedian(sortarray, left, median, axis);
}
if (right < median - 1) {
/* Recursively shuffle the right subset. */
vertexmedian(&sortarray[right + 1], arraysize - right - 1,
median - right - 1, axis);
}
}
/*****************************************************************************/
/* */
/* alternateaxes() Sorts the vertices as appropriate for the divide-and- */
/* conquer algorithm with alternating cuts. */
/* */
/* Partitions by x-coordinate if axis == 0; by y-coordinate if axis == 1. */
/* For the base case, subsets containing only two or three vertices are */
/* always sorted by x-coordinate. */
/* */
/*****************************************************************************/
#ifdef ANSI_DECLARATORS
void alternateaxes(vertex *sortarray, int arraysize, int axis)
#else /* not ANSI_DECLARATORS */
void alternateaxes(sortarray, arraysize, axis)
vertex *sortarray;
int arraysize;
int axis;
#endif /* not ANSI_DECLARATORS */
{
int divider;
divider = arraysize >> 1;
if (arraysize <= 3) {
/* Recursive base case: subsets of two or three vertices will be */
/* handled specially, and should always be sorted by x-coordinate. */
axis = 0;
}
/* Partition with a horizontal or vertical cut. */
vertexmedian(sortarray, arraysize, divider, axis);
/* Recursively partition the subsets with a cross cut. */
if (arraysize - divider >= 2) {
if (divider >= 2) {
alternateaxes(sortarray, divider, 1 - axis);
}
alternateaxes(&sortarray[divider], arraysize - divider, 1 - axis);
}
}
/*****************************************************************************/
/* */
/* mergehulls() Merge two adjacent Delaunay triangulations into a */
/* single Delaunay triangulation. */
/* */
/* This is similar to the algorithm given by Guibas and Stolfi, but uses */
/* a triangle-based, rather than edge-based, data structure. */
/* */
/* The algorithm walks up the gap between the two triangulations, knitting */
/* them together. As they are merged, some of their bounding triangles */
/* are converted into real triangles of the triangulation. The procedure */
/* pulls each hull's bounding triangles apart, then knits them together */
/* like the teeth of two gears. The Delaunay property determines, at each */
/* step, whether the next "tooth" is a bounding triangle of the left hull */
/* or the right. When a bounding triangle becomes real, its apex is */
/* changed from NULL to a real vertex. */
/* */
/* Only two new triangles need to be allocated. These become new bounding */
/* triangles at the top and bottom of the seam. They are used to connect */
/* the remaining bounding triangles (those that have not been converted */
/* into real triangles) into a single fan. */
/* */
/* On entry, `farleft' and `innerleft' are bounding triangles of the left */
/* triangulation. The origin of `farleft' is the leftmost vertex, and */
/* the destination of `innerleft' is the rightmost vertex of the */
/* triangulation. Similarly, `innerright' and `farright' are bounding */
/* triangles of the right triangulation. The origin of `innerright' and */
/* destination of `farright' are the leftmost and rightmost vertices. */
/* */
/* On completion, the origin of `farleft' is the leftmost vertex of the */
/* merged triangulation, and the destination of `farright' is the rightmost */
/* vertex. */
/* */
/*****************************************************************************/
#ifdef ANSI_DECLARATORS
void mergehulls(struct mesh *m, struct behavior *b, struct otri *farleft,
struct otri *innerleft, struct otri *innerright,
struct otri *farright, int axis)
#else /* not ANSI_DECLARATORS */
void mergehulls(m, b, farleft, innerleft, innerright, farright, axis)
struct mesh *m;
struct behavior *b;
struct otri *farleft;
struct otri *innerleft;
struct otri *innerright;
struct otri *farright;
int axis;
#endif /* not ANSI_DECLARATORS */
{
struct otri leftcand, rightcand;
struct otri baseedge;
struct otri nextedge;
struct otri sidecasing, topcasing, outercasing;
struct otri checkedge;
vertex innerleftdest;
vertex innerrightorg;
vertex innerleftapex, innerrightapex;
vertex farleftpt, farrightpt;
vertex farleftapex, farrightapex;
vertex lowerleft, lowerright;
vertex upperleft, upperright;
vertex nextapex;
vertex checkvertex;
int changemade;
int badedge;
int leftfinished, rightfinished;
triangle ptr; /* Temporary variable used by sym(). */
dest(*innerleft, innerleftdest);
apex(*innerleft, innerleftapex);
org(*innerright, innerrightorg);
apex(*innerright, innerrightapex);
/* Special treatment for horizontal cuts. */
if (b->dwyer && (axis == 1)) {
org(*farleft, farleftpt);
apex(*farleft, farleftapex);
dest(*farright, farrightpt);
apex(*farright, farrightapex);
/* The pointers to the extremal vertices are shifted to point to the */
/* topmost and bottommost vertex of each hull, rather than the */
/* leftmost and rightmost vertices. */
while (farleftapex[1] < farleftpt[1]) {
lnextself(*farleft);
symself(*farleft);
farleftpt = farleftapex;
apex(*farleft, farleftapex);
}
sym(*innerleft, checkedge);
apex(checkedge, checkvertex);
while (checkvertex[1] > innerleftdest[1]) {
lnext(checkedge, *innerleft);
innerleftapex = innerleftdest;
innerleftdest = checkvertex;
sym(*innerleft, checkedge);
apex(checkedge, checkvertex);
}
while (innerrightapex[1] < innerrightorg[1]) {
lnextself(*innerright);
symself(*innerright);
innerrightorg = innerrightapex;
apex(*innerright, innerrightapex);
}
sym(*farright, checkedge);
apex(checkedge, checkvertex);
while (checkvertex[1] > farrightpt[1]) {
lnext(checkedge, *farright);
farrightapex = farrightpt;
farrightpt = checkvertex;
sym(*farright, checkedge);
apex(checkedge, checkvertex);
}
}
/* Find a line tangent to and below both hulls. */
do {
changemade = 0;
/* Make innerleftdest the "bottommost" vertex of the left hull. */
if (counterclockwise(m, b, innerleftdest, innerleftapex, innerrightorg) >
0.0) {
lprevself(*innerleft);
symself(*innerleft);
innerleftdest = innerleftapex;
apex(*innerleft, innerleftapex);
changemade = 1;
}
/* Make innerrightorg the "bottommost" vertex of the right hull. */
if (counterclockwise(m, b, innerrightapex, innerrightorg, innerleftdest) >
0.0) {
lnextself(*innerright);
symself(*innerright);
innerrightorg = innerrightapex;
apex(*innerright, innerrightapex);
changemade = 1;
}
} while (changemade);
/* Find the two candidates to be the next "gear tooth." */
sym(*innerleft, leftcand);
sym(*innerright, rightcand);
/* Create the bottom new bounding triangle. */
maketriangle(m, b, &baseedge);
/* Connect it to the bounding boxes of the left and right triangulations. */
bond(baseedge, *innerleft);
lnextself(baseedge);
bond(baseedge, *innerright);
lnextself(baseedge);
setorg(baseedge, innerrightorg);
setdest(baseedge, innerleftdest);
/* Apex is intentionally left NULL. */
if (b->verbose > 2) {
printf(" Creating base bounding ");
printtriangle(m, b, &baseedge);
}
/* Fix the extreme triangles if necessary. */
org(*farleft, farleftpt);
if (innerleftdest == farleftpt) {
lnext(baseedge, *farleft);
}
dest(*farright, farrightpt);
if (innerrightorg == farrightpt) {
lprev(baseedge, *farright);
}
/* The vertices of the current knitting edge. */
lowerleft = innerleftdest;
lowerright = innerrightorg;
/* The candidate vertices for knitting. */
apex(leftcand, upperleft);
apex(rightcand, upperright);
/* Walk up the gap between the two triangulations, knitting them together. */
while (1) {
/* Have we reached the top? (This isn't quite the right question, */
/* because even though the left triangulation might seem finished now, */
/* moving up on the right triangulation might reveal a new vertex of */
/* the left triangulation. And vice-versa.) */
leftfinished = counterclockwise(m, b, upperleft, lowerleft, lowerright) <=
0.0;
rightfinished = counterclockwise(m, b, upperright, lowerleft, lowerright)
<= 0.0;
if (leftfinished && rightfinished) {
/* Create the top new bounding triangle. */
maketriangle(m, b, &nextedge);
setorg(nextedge, lowerleft);
setdest(nextedge, lowerright);
/* Apex is intentionally left NULL. */
/* Connect it to the bounding boxes of the two triangulations. */
bond(nextedge, baseedge);
lnextself(nextedge);
bond(nextedge, rightcand);
lnextself(nextedge);
bond(nextedge, leftcand);
if (b->verbose > 2) {
printf(" Creating top bounding ");
printtriangle(m, b, &nextedge);
}
/* Special treatment for horizontal cuts. */
if (b->dwyer && (axis == 1)) {
org(*farleft, farleftpt);
apex(*farleft, farleftapex);
dest(*farright, farrightpt);
apex(*farright, farrightapex);
sym(*farleft, checkedge);
apex(checkedge, checkvertex);
/* The pointers to the extremal vertices are restored to the */
/* leftmost and rightmost vertices (rather than topmost and */
/* bottommost). */
while (checkvertex[0] < farleftpt[0]) {
lprev(checkedge, *farleft);
farleftapex = farleftpt;
farleftpt = checkvertex;
sym(*farleft, checkedge);
apex(checkedge, checkvertex);
}
while (farrightapex[0] > farrightpt[0]) {
lprevself(*farright);
symself(*farright);
farrightpt = farrightapex;
apex(*farright, farrightapex);
}
}
return;
}
/* Consider eliminating edges from the left triangulation. */
if (!leftfinished) {
/* What vertex would be exposed if an edge were deleted? */
lprev(leftcand, nextedge);