Math-Geometry-Delaunay
view release on metacpan or search on metacpan
src/triangle.c view on Meta::CPAN
/* purpose of point location. Splay trees are described by Daniel */
/* Dominic Sleator and Robert Endre Tarjan, "Self-Adjusting Binary Search */
/* Trees," Journal of the ACM 32(3):652-686, July 1985, */
/* http://portal.acm.org/citation.cfm?id=3835 . */
/* */
/* The algorithms for exact computation of the signs of determinants are */
/* described in Jonathan Richard Shewchuk, "Adaptive Precision Floating- */
/* Point Arithmetic and Fast Robust Geometric Predicates," Discrete & */
/* Computational Geometry 18(3):305-363, October 1997. (Also available */
/* as Technical Report CMU-CS-96-140, School of Computer Science, */
/* Carnegie Mellon University, Pittsburgh, Pennsylvania, May 1996.) [*] */
/* An abbreviated version appears as Jonathan Richard Shewchuk, "Robust */
/* Adaptive Floating-Point Geometric Predicates," Proceedings of the */
/* Twelfth Annual Symposium on Computational Geometry, ACM, May 1996. [*] */
/* Many of the ideas for my exact arithmetic routines originate with */
/* Douglas M. Priest, "Algorithms for Arbitrary Precision Floating Point */
/* Arithmetic," Tenth Symposium on Computer Arithmetic, pp. 132-143, IEEE */
/* Computer Society Press, 1991. [*] Many of the ideas for the correct */
/* evaluation of the signs of determinants are taken from Steven Fortune */
/* and Christopher J. Van Wyk, "Efficient Exact Arithmetic for Computa- */
/* tional Geometry," Proceedings of the Ninth Annual Symposium on */
/* Computational Geometry, ACM, pp. 163-172, May 1993, and from Steven */
/* Fortune, "Numerical Stability of Algorithms for 2D Delaunay Triangu- */
/* lations," International Journal of Computational Geometry & Applica- */
/* tions 5(1-2):193-213, March-June 1995. */
/* */
/* The method of inserting new vertices off-center (not precisely at the */
/* circumcenter of every poor-quality triangle) is from Alper Ungor, */
/* "Off-centers: A New Type of Steiner Points for Computing Size-Optimal */
/* Quality-Guaranteed Delaunay Triangulations," Proceedings of LATIN */
/* 2004 (Buenos Aires, Argentina), April 2004. */
/* */
/* For definitions of and results involving Delaunay triangulations, */
/* constrained and conforming versions thereof, and other aspects of */
/* triangular mesh generation, see the excellent survey by Marshall Bern */
/* and David Eppstein, "Mesh Generation and Optimal Triangulation," in */
/* Computing and Euclidean Geometry, Ding-Zhu Du and Frank Hwang, */
/* editors, World Scientific, Singapore, pp. 23-90, 1992. [*] */
/* */
/* The time for incrementally adding PSLG (planar straight line graph) */
/* segments to create a constrained Delaunay triangulation is probably */
/* O(t^2) per segment in the worst case and O(t) per segment in the */
/* common case, where t is the number of triangles that intersect the */
/* segment before it is inserted. This doesn't count point location, */
/* which can be much more expensive. I could improve this to O(d log d) */
/* time, but d is usually quite small, so it's not worth the bother. */
/* (This note does not apply when the -s switch is used, invoking a */
/* different method is used to insert segments.) */
/* */
/* The time for deleting a vertex from a Delaunay triangulation is O(d^2) */
/* in the worst case and O(d) in the common case, where d is the degree */
/* of the vertex being deleted. I could improve this to O(d log d) time, */
/* but d is usually quite small, so it's not worth the bother. */
/* */
/* Ruppert's Delaunay refinement algorithm typically generates triangles */
/* at a linear rate (constant time per triangle) after the initial */
/* triangulation is formed. There may be pathological cases where */
/* quadratic time is required, but these never arise in practice. */
/* */
/* The geometric predicates (circumcenter calculations, segment */
/* intersection formulae, etc.) appear in my "Lecture Notes on Geometric */
/* Robustness" at http://www.cs.berkeley.edu/~jrs/mesh . */
/* */
/* If you make any improvements to this code, please please please let me */
/* know, so that I may obtain the improvements. Even if you don't change */
/* the code, I'd still love to hear what it's being used for. */
/* */
/*****************************************************************************/
/* For single precision (which will save some memory and reduce paging), */
/* define the symbol SINGLE by using the -DSINGLE compiler switch or by */
/* writing "#define SINGLE" below. */
/* */
/* For double precision (which will allow you to refine meshes to a smaller */
/* edge length), leave SINGLE undefined. */
/* */
/* Double precision uses more memory, but improves the resolution of the */
/* meshes you can generate with Triangle. It also reduces the likelihood */
/* of a floating exception due to overflow. Finally, it is much faster */
/* than single precision on 64-bit architectures like the DEC Alpha. I */
/* recommend double precision unless you want to generate a mesh for which */
/* you do not have enough memory. */
/* #define SINGLE */
#ifdef SINGLE
#define REAL float
#else /* not SINGLE */
#define REAL double
#endif /* not SINGLE */
/* If yours is not a Unix system, define the NO_TIMER compiler switch to */
/* remove the Unix-specific timing code. */
/* #define NO_TIMER */
/* To insert lots of self-checks for internal errors, define the SELF_CHECK */
/* symbol. This will slow down the program significantly. It is best to */
/* define the symbol using the -DSELF_CHECK compiler switch, but you could */
/* write "#define SELF_CHECK" below. If you are modifying this code, I */
/* recommend you turn self-checks on until your work is debugged. */
/* #define SELF_CHECK */
/* To compile Triangle as a callable object library (triangle.o), define the */
/* TRILIBRARY symbol. Read the file triangle.h for details on how to call */
/* the procedure triangulate() that results. */
/* #define TRILIBRARY */
/* It is possible to generate a smaller version of Triangle using one or */
/* both of the following symbols. Define the REDUCED symbol to eliminate */
/* all features that are primarily of research interest; specifically, the */
/* -i, -F, -s, and -C switches. Define the CDT_ONLY symbol to eliminate */
/* all meshing algorithms above and beyond constrained Delaunay */
/* triangulation; specifically, the -r, -q, -a, -u, -D, -S, and -s */
/* switches. These reductions are most likely to be useful when */
/* generating an object library (triangle.o) by defining the TRILIBRARY */
/* symbol. */
/* #define REDUCED */
src/triangle.c view on Meta::CPAN
}
if (leftccw == 0.0) {
return LEFTCOLLINEAR;
} else if (rightccw == 0.0) {
return RIGHTCOLLINEAR;
} else {
return WITHIN;
}
}
/*****************************************************************************/
/* */
/* segmentintersection() Find the intersection of an existing segment */
/* and a segment that is being inserted. Insert */
/* a vertex at the intersection, splitting an */
/* existing subsegment. */
/* */
/* The segment being inserted connects the apex of splittri to endpoint2. */
/* splitsubseg is the subsegment being split, and MUST adjoin splittri. */
/* Hence, endpoints of the subsegment being split are the origin and */
/* destination of splittri. */
/* */
/* On completion, splittri is a handle having the newly inserted */
/* intersection point as its origin, and endpoint1 as its destination. */
/* */
/*****************************************************************************/
#ifdef ANSI_DECLARATORS
void segmentintersection(struct mesh *m, struct behavior *b,
struct otri *splittri, struct osub *splitsubseg,
vertex endpoint2)
#else /* not ANSI_DECLARATORS */
void segmentintersection(m, b, splittri, splitsubseg, endpoint2)
struct mesh *m;
struct behavior *b;
struct otri *splittri;
struct osub *splitsubseg;
vertex endpoint2;
#endif /* not ANSI_DECLARATORS */
{
struct osub opposubseg;
vertex endpoint1;
vertex torg, tdest;
vertex leftvertex, rightvertex;
vertex newvertex;
enum insertvertexresult success;
enum finddirectionresult collinear;
REAL ex, ey;
REAL tx, ty;
REAL etx, ety;
REAL split, denom;
int i;
triangle ptr; /* Temporary variable used by onext(). */
subseg sptr; /* Temporary variable used by snext(). */
/* Find the other three segment endpoints. */
apex(*splittri, endpoint1);
org(*splittri, torg);
dest(*splittri, tdest);
/* Segment intersection formulae; see the Antonio reference. */
tx = tdest[0] - torg[0];
ty = tdest[1] - torg[1];
ex = endpoint2[0] - endpoint1[0];
ey = endpoint2[1] - endpoint1[1];
etx = torg[0] - endpoint2[0];
ety = torg[1] - endpoint2[1];
denom = ty * ex - tx * ey;
if (denom == 0.0) {
printf("Internal error in segmentintersection():");
printf(" Attempt to find intersection of parallel segments.\n");
internalerror();
}
split = (ey * etx - ex * ety) / denom;
/* Create the new vertex. */
newvertex = (vertex) poolalloc(&m->vertices);
/* Interpolate its coordinate and attributes. */
for (i = 0; i < 2 + m->nextras; i++) {
newvertex[i] = torg[i] + split * (tdest[i] - torg[i]);
}
setvertexmark(newvertex, mark(*splitsubseg));
setvertextype(newvertex, INPUTVERTEX);
if (b->verbose > 1) {
printf(
" Splitting subsegment (%.12g, %.12g) (%.12g, %.12g) at (%.12g, %.12g).\n",
torg[0], torg[1], tdest[0], tdest[1], newvertex[0], newvertex[1]);
}
/* Insert the intersection vertex. This should always succeed. */
success = insertvertex(m, b, newvertex, splittri, splitsubseg, 0, 0);
if (success != SUCCESSFULVERTEX) {
printf("Internal error in segmentintersection():\n");
printf(" Failure to split a segment.\n");
internalerror();
}
/* Record a triangle whose origin is the new vertex. */
setvertex2tri(newvertex, encode(*splittri));
if (m->steinerleft > 0) {
m->steinerleft--;
}
/* Divide the segment into two, and correct the segment endpoints. */
ssymself(*splitsubseg);
spivot(*splitsubseg, opposubseg);
sdissolve(*splitsubseg);
sdissolve(opposubseg);
do {
setsegorg(*splitsubseg, newvertex);
snextself(*splitsubseg);
} while (splitsubseg->ss != m->dummysub);
do {
setsegorg(opposubseg, newvertex);
snextself(opposubseg);
} while (opposubseg.ss != m->dummysub);
/* Inserting the vertex may have caused edge flips. We wish to rediscover */
/* the edge connecting endpoint1 to the new intersection vertex. */
collinear = finddirection(m, b, splittri, endpoint1);
dest(*splittri, rightvertex);
apex(*splittri, leftvertex);
if ((leftvertex[0] == endpoint1[0]) && (leftvertex[1] == endpoint1[1])) {
onextself(*splittri);
( run in 0.518 second using v1.01-cache-2.11-cpan-39bf76dae61 )