SaraP
491 lines
20 KiB
C
491 lines
20 KiB
C
/*****************************************************************************/
|
|
/* */
|
|
/* F I S T : Fast, Industrial-Strength Triangulation */
|
|
/* */
|
|
/*****************************************************************************/
|
|
/* */
|
|
/* (C) Martin Held */
|
|
/* (C) Universitaet Salzburg, Salzburg, Austria */
|
|
/* */
|
|
/* This code is not in the public domain. All rights reserved! Please make */
|
|
/* sure to read the full copyright statement contained in api_functions.cpp. */
|
|
/* */
|
|
/*****************************************************************************/
|
|
|
|
|
|
#define BaseLength(u, v, base, delta) { \
|
|
(base).x = (v).x - (u).x; \
|
|
(base).y = (v).y - (u).y; \
|
|
delta = Abs(TO_MDOUBLE((base).x)) + Abs(TO_MDOUBLE((base).y)); }
|
|
|
|
|
|
#define SideLength(u, v, base, delta) { \
|
|
(base).x = (v).x - (u).x; \
|
|
(base).y = (v).y - (u).y; \
|
|
delta = (base).x * (base).x + (base).y * (base).y; }
|
|
|
|
#define PntPntDist(u, v, base, delta) { \
|
|
(base).x = (v).x - (u).x; \
|
|
(base).y = (v).y - (u).y; \
|
|
delta = sqrt((base).x * (base).x + (base).y * (base).y); }
|
|
|
|
#define PntPntDistSqd(u, v, base, delta) { \
|
|
(base).x = TO_MDOUBLE((v).x) - TO_MDOUBLE((u).x); \
|
|
(base).y = TO_MDOUBLE((v).y) - TO_MDOUBLE((u).y); \
|
|
delta = (base).x * (base).x + (base).y * (base).y; }
|
|
|
|
|
|
/* */
|
|
/* this macro checks whether i3, which is collinear with i1, i2, is */
|
|
/* between i1, i2. note that we rely on the lexicographic sorting of the */
|
|
/* points! */
|
|
/* */
|
|
#define InBetween(i1, i2, i3) \
|
|
((i1 <= i3) && (i3 <= i2))
|
|
|
|
|
|
#define StrictlyInBetween(i1, i2, i3) \
|
|
((i1 < i3) && (i3 < i2))
|
|
|
|
|
|
/* */
|
|
/* this macro checks whether i3, which is collinear with i1, i2, is */
|
|
/* between i1, i2. note that we do _not_ rely on the lexicographic sorting */
|
|
/* of the points! */
|
|
/* */
|
|
#define InBetweenEps(data, i1, i2, i3) \
|
|
(((Min(data->points[i1].x, data->points[i2].x) - EPS) <= data->points[i3].x) && \
|
|
((Min(data->points[i1].y, data->points[i2].y) - EPS) <= data->points[i3].y) && \
|
|
((Max(data->points[i1].x, data->points[i2].x) + EPS) >= data->points[i3].x) && \
|
|
((Max(data->points[i1].y, data->points[i2].y) + EPS) >= data->points[i3].y))
|
|
|
|
|
|
/* */
|
|
/* this macro computes the determinant det(points[i],points[j],points[k]) */
|
|
/* in a consistent way. */
|
|
/* */
|
|
#ifdef JRC_PREDICATE
|
|
double orient2dfast(point pa, point pb, point pc);
|
|
double orient2d(point pa, point pb, point pc);
|
|
#define StableDet2D(i, j, k, det) \
|
|
{ assert(InPointsList(i)); \
|
|
assert(InPointsList(j)); \
|
|
assert(InPointsList(k)); \
|
|
\
|
|
if ((i == j) || (i == k) || (j == k)) { \
|
|
det = C_0_0; \
|
|
} \
|
|
else { \
|
|
det = orient2d(points[i], points[j], points[k]); \
|
|
} \
|
|
}
|
|
#else
|
|
#define StableDet2D(data, tmp, i, j, k, det) \
|
|
{ assert(InPointsList(data,i)); \
|
|
assert(InPointsList(data,j)); \
|
|
assert(InPointsList(data,k)); \
|
|
\
|
|
if ((i == j) || (i == k) || (j == k)) { \
|
|
det = C_0_0; \
|
|
} \
|
|
else { \
|
|
tmp.numerics_h_p = data->points[i]; \
|
|
tmp.numerics_h_q = data->points[j]; \
|
|
tmp.numerics_h_r = data->points[k]; \
|
|
\
|
|
if (i < j) { \
|
|
if (j < k) /* i < j < k */ \
|
|
det = Det2D(tmp.numerics_h_p, tmp.numerics_h_q, tmp.numerics_h_r); \
|
|
else if (i < k) /* i < k < j */ \
|
|
det = -Det2D(tmp.numerics_h_p, tmp.numerics_h_r, tmp.numerics_h_q); \
|
|
else /* k < i < j */ \
|
|
det = Det2D(tmp.numerics_h_r, tmp.numerics_h_p, tmp.numerics_h_q); \
|
|
} \
|
|
else { \
|
|
if (i < k) /* j < i < k */ \
|
|
det = -Det2D(tmp.numerics_h_q, tmp.numerics_h_p, tmp.numerics_h_r); \
|
|
else if (j < k) /* j < k < i */ \
|
|
det = Det2D(tmp.numerics_h_q, tmp.numerics_h_r, tmp.numerics_h_p); \
|
|
else /* k < j < i */ \
|
|
det = -Det2D(tmp.numerics_h_r, tmp.numerics_h_q, tmp.numerics_h_p); \
|
|
} \
|
|
} \
|
|
}
|
|
#endif
|
|
|
|
|
|
/* */
|
|
/* this macro sets ori to +1 if the points i, j, k are given in CCW order, */
|
|
/* -1 if the points i, j, k are given in CW order, */
|
|
/* 0 if the points i, j, k are collinear. */
|
|
/* */
|
|
#define Orientation(data, tmp, i, j, k, ori) \
|
|
{ StableDet2D(data, tmp, i, j, k, tmp.numerics_h_det); \
|
|
\
|
|
if (lt(tmp.numerics_h_det, EPS)) ori = -1; \
|
|
else if (gt(tmp.numerics_h_det, EPS)) ori = 1; \
|
|
else ori = 0; \
|
|
}
|
|
|
|
|
|
/* */
|
|
/* this macro sets ori to +1 if the points i, j, k are given in CCW order, */
|
|
/* -1 if the points i, j, k are given in CW order, */
|
|
/* 0 if the points i, j, k are collinear. */
|
|
/* */
|
|
#define OrientationThres(data, tmp, i, j, k, ori, thres) \
|
|
{ StableDet2D(data, tmp, i, j, k, tmp.numerics_h_det); \
|
|
\
|
|
if (lt(tmp.numerics_h_det, thres)) ori = -1; \
|
|
else if (gt(tmp.numerics_h_det, thres)) ori = 1; \
|
|
else ori = 0; \
|
|
}
|
|
|
|
|
|
/* */
|
|
/* this macro checks whether l is in the cone defined by i, j and j, k */
|
|
/* */
|
|
#define IsInConeStrict(data, tmp, i, j, k, l, convex, flag, thres) \
|
|
{ assert(InPointsList(data, i)); \
|
|
assert(InPointsList(data, j)); \
|
|
assert(InPointsList(data, k)); \
|
|
assert(InPointsList(data, l)); \
|
|
\
|
|
flag = true; \
|
|
if (convex) { \
|
|
if ((l != i) && (i != j)) { \
|
|
OrientationThres(data, tmp, i, j, l, tmp.numerics_h_ori1, thres); \
|
|
if (tmp.numerics_h_ori1 <= 0) flag = false; \
|
|
} \
|
|
if ((l != k) && (j != k) && (flag == true)) { \
|
|
OrientationThres(data, tmp, j, k, l, tmp.numerics_h_ori2, thres); \
|
|
if (tmp.numerics_h_ori2 <= 0) flag = false; \
|
|
} \
|
|
} \
|
|
else { \
|
|
OrientationThres(data, tmp, i, j, l, tmp.numerics_h_ori1, thres); \
|
|
if (tmp.numerics_h_ori1 <= 0) { \
|
|
OrientationThres(data, tmp, j, k, l, tmp.numerics_h_ori2, thres); \
|
|
if (tmp.numerics_h_ori2 <= 0) flag = false; \
|
|
} \
|
|
} \
|
|
}
|
|
|
|
|
|
/* */
|
|
/* this macro checks whether l is in the cone defined by i, j and j, k */
|
|
/* */
|
|
#define IsInCone(data, tmp, i, j, k, l, convex, flag) \
|
|
{ assert(InPointsList(data, i)); \
|
|
assert(InPointsList(data, j)); \
|
|
assert(InPointsList(data, k)); \
|
|
assert(InPointsList(data, l)); \
|
|
\
|
|
flag = true; \
|
|
if (convex) { \
|
|
if ((l != i) && (i != j)) { \
|
|
Orientation(data, tmp, i, j, l, tmp.numerics_h_ori1); \
|
|
if (tmp.numerics_h_ori1 < 0) flag = false; \
|
|
else if (tmp.numerics_h_ori1 == 0) { \
|
|
if (i < j) { \
|
|
if (!InBetween(i, j, l)) flag = false; \
|
|
} \
|
|
else { \
|
|
if (!InBetween(j, i, l)) flag = false; \
|
|
} \
|
|
} \
|
|
} \
|
|
if ((l != k) && (j != k) && (flag == true)) { \
|
|
Orientation(data, tmp, j, k, l, tmp.numerics_h_ori2); \
|
|
if (tmp.numerics_h_ori2 < 0) flag = false; \
|
|
else if (tmp.numerics_h_ori2 == 0) { \
|
|
if (j < k) { \
|
|
if (!InBetween(j, k, l)) flag = false; \
|
|
} \
|
|
else { \
|
|
if (!InBetween(k, j, l)) flag = false; \
|
|
} \
|
|
} \
|
|
} \
|
|
} \
|
|
else { \
|
|
Orientation(data, tmp, i, j, l, tmp.numerics_h_ori1); \
|
|
if (tmp.numerics_h_ori1 <= 0) { \
|
|
Orientation(data, tmp, j, k, l, tmp.numerics_h_ori2); \
|
|
if (tmp.numerics_h_ori2 <= 0) { \
|
|
if (tmp.numerics_h_ori1 == 0) { \
|
|
if (!InBetween(i, j, l)) flag = false; \
|
|
} \
|
|
else if (tmp.numerics_h_ori2 == 0) { \
|
|
if (!InBetween(j, k, l)) flag = false; \
|
|
} \
|
|
else flag = false; \
|
|
} \
|
|
} \
|
|
} \
|
|
}
|
|
|
|
|
|
|
|
/* */
|
|
/* this macro checks whether l is in the cone defined by i, j and j, k */
|
|
/* */
|
|
#define IsInConeEps(data, tmp, i, j, k, l, convex, flag) \
|
|
{ assert(InPointsList(data, i)); \
|
|
assert(InPointsList(data, j)); \
|
|
assert(InPointsList(data, k)); \
|
|
assert(InPointsList(data, l)); \
|
|
\
|
|
flag = true; \
|
|
if (convex) { \
|
|
if ((l != i) && (i != j)) { \
|
|
Orientation(data, tmp, i, j, l, tmp.numerics_h_ori1); \
|
|
if (tmp.numerics_h_ori1 < 0) flag = false; \
|
|
else if (tmp.numerics_h_ori1 == 0) { \
|
|
if (!InBetweenEps(data, i, j, l)) flag = false; \
|
|
} \
|
|
} \
|
|
if ((l != k) && (j != k) && (flag == true)) { \
|
|
Orientation(data, tmp, j, k, l, tmp.numerics_h_ori2); \
|
|
if (tmp.numerics_h_ori2 < 0) flag = false; \
|
|
else if (tmp.numerics_h_ori2 == 0) { \
|
|
if (!InBetweenEps(data, j, k, l)) flag = false; \
|
|
} \
|
|
} \
|
|
} \
|
|
else { \
|
|
Orientation(data, tmp, i, j, l, tmp.numerics_h_ori1); \
|
|
if (tmp.numerics_h_ori1 <= 0) { \
|
|
Orientation(data, tmp, j, k, l, tmp.numerics_h_ori2); \
|
|
if (tmp.numerics_h_ori2 <= 0) { \
|
|
if (tmp.numerics_h_ori1 == 0) { \
|
|
if (!InBetweenEps(data, i, j, l)) flag = false; \
|
|
} \
|
|
else if (tmp.numerics_h_ori2 == 0) { \
|
|
if (!InBetweenEps(data, j, k, l)) flag = false; \
|
|
} \
|
|
else flag = false; \
|
|
} \
|
|
} \
|
|
} \
|
|
}
|
|
|
|
|
|
|
|
|
|
/* */
|
|
/* this macro checks whether the diagonal j, l is in the cone defined by */
|
|
/* i, j and j, k. if left then j, k, l forms the ear to be tested. */
|
|
/* otherwise, i, j, k forms the ear to be tested. */
|
|
/* */
|
|
#define DiagIsInConeEps(data, tmp, i, j, k, l, convex, left, flag) \
|
|
{ assert(InPointsList(data, i)); \
|
|
assert(InPointsList(data, j)); \
|
|
assert(InPointsList(data, k)); \
|
|
assert(InPointsList(data, l)); \
|
|
\
|
|
flag = true; \
|
|
if (convex) { \
|
|
if (left) { \
|
|
if ((l != i) && (i != j)) { \
|
|
Orientation(data, tmp, i, j, l, tmp.numerics_h_ori1); \
|
|
if (tmp.numerics_h_ori1 < 0) flag = false; \
|
|
else if (tmp.numerics_h_ori1 == 0) { \
|
|
if (InBetweenEps(data, j, l, i)) flag = false; \
|
|
} \
|
|
} \
|
|
} \
|
|
else { \
|
|
if ((l != k) && (j != k)) { \
|
|
Orientation(data, tmp, j, k, l, tmp.numerics_h_ori2); \
|
|
if (tmp.numerics_h_ori2 < 0) flag = false; \
|
|
else if (tmp.numerics_h_ori2 == 0) { \
|
|
if (InBetweenEps(data, j, l, k)) flag = false; \
|
|
} \
|
|
} \
|
|
} \
|
|
} \
|
|
}
|
|
|
|
|
|
|
|
#define DiagIsInCone(data, tmp, i, j, k, l, convex, left, flag) \
|
|
{ assert(InPointsList(data, i)); \
|
|
assert(InPointsList(data, j)); \
|
|
assert(InPointsList(data, k)); \
|
|
assert(InPointsList(data, l)); \
|
|
\
|
|
flag = true; \
|
|
if (convex) { \
|
|
if (left) { \
|
|
if ((l != i) && (i != j)) { \
|
|
Orientation(data, tmp, i, j, l, tmp.numerics_h_ori1); \
|
|
if (tmp.numerics_h_ori1 < 0) flag = false; \
|
|
else if (tmp.numerics_h_ori1 == 0) { \
|
|
if (InBetween(j, l, i)) flag = false; \
|
|
} \
|
|
} \
|
|
} \
|
|
else { \
|
|
if ((l != k) && (j != k)) { \
|
|
Orientation(data, tmp, j, k, l, tmp.numerics_h_ori2); \
|
|
if (tmp.numerics_h_ori2 < 0) flag = false; \
|
|
else if (tmp.numerics_h_ori2 == 0) { \
|
|
if (InBetween(j, l, k)) flag = false; \
|
|
} \
|
|
} \
|
|
} \
|
|
} \
|
|
}
|
|
|
|
|
|
/* */
|
|
/* returns */
|
|
/* -2 ... if angle is 360 degrees */
|
|
/* -1 ... if angle between 180 and 360 degrees */
|
|
/* 0 ... if angle is 180 degrees */
|
|
/* 1 ... if angle between 0 and 180 degrees */
|
|
/* 2 ... if angle is 0 degrees */
|
|
/* 3 ... if angle is about 180 degrees and triangle is very small */
|
|
/* */
|
|
inline int IsConvexAngle(listdef *list, datadef *data, int i, int j, int k,
|
|
list_ind ind) {
|
|
assert(InPointsList(data, i));
|
|
assert(InPointsList(data, j));
|
|
assert(InPointsList(data, k));
|
|
tmp_data_def tmp;
|
|
|
|
if (i == j) {
|
|
if (j == k) {
|
|
/* */
|
|
/* all three vertices are identical. we set the angle to -2. */
|
|
/* using -2 means to err on the safe side, as all the */
|
|
/* incarnations of this vertex will be clipped right at the */
|
|
/* start of the ear-clipping algorithm. thus, eventually there */
|
|
/* will be no other duplicates at this vertex position, and the */
|
|
/* regular classification of angles will yield the correct */
|
|
/* answer for j. */
|
|
/* */
|
|
return -2;
|
|
}
|
|
else {
|
|
/* */
|
|
/* two of the three vertices are identical; we set the angle to */
|
|
/* to 2 order to enable clipping of j. */
|
|
/* */
|
|
return 2;
|
|
}
|
|
}
|
|
else if (j == k) {
|
|
/* */
|
|
/* two vertices are identical. we could either determine the angle */
|
|
/* by means of yet another lengthy analysis, or simply set the */
|
|
/* angle to -1. using -1 means to err on the safe side, as all the */
|
|
/* incarnations of this vertex will be clipped right at the start */
|
|
/* of the ear-clipping algorithm. thus, eventually there will be no */
|
|
/* other duplicates at this vertex position, and the regular */
|
|
/* classification of angles will yield the correct answer for j. */
|
|
/* */
|
|
return -1;
|
|
} else {
|
|
StableDet2D(data, tmp, i, j, k, tmp.numerics_h_det);
|
|
if (gt(tmp.numerics_h_det, EPS)) {
|
|
return 1;
|
|
}
|
|
else if (lt(tmp.numerics_h_det, EPS)) {
|
|
return -1;
|
|
}
|
|
else {
|
|
/* */
|
|
/* 0, 180, or 360 degrees. */
|
|
/* */
|
|
VectorSub2D(data->points[i], data->points[j], tmp.numerics_h_p);
|
|
VectorSub2D(data->points[k], data->points[j], tmp.numerics_h_q);
|
|
tmp.numerics_h_dot = DotProduct2D(tmp.numerics_h_p, tmp.numerics_h_q);
|
|
if (tmp.numerics_h_dot < C_0_0) {
|
|
/* */
|
|
/* 180 degrees. */
|
|
/* */
|
|
/*data->numerics_h_dot = sqrt(Length2D(data->numerics_h_p)); */
|
|
/*data->numerics_h_dot += sqrt(Length2D(data->numerics_h_q)); */
|
|
if (tmp.numerics_h_det > C_0_0) return 1;
|
|
else return 0;
|
|
}
|
|
else {
|
|
/* */
|
|
/* 0 or 360 degrees? this cannot be judged locally, and more */
|
|
/* work is needed. */
|
|
/* */
|
|
return SpikeAngle(list, data, i, j, k, ind);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
/* */
|
|
/* this macro checks whether point i4 is inside of or on the boundary */
|
|
/* of the triangle i1, i2, i3 */
|
|
/* */
|
|
#define PntInTriangle(data, tmp, i1, i2, i3, i4, inside) \
|
|
{ \
|
|
inside = false; \
|
|
Orientation(data, tmp, i2, i3, i4, tmp.numerics_h_ori1); \
|
|
if (tmp.numerics_h_ori1 >= 0) { \
|
|
Orientation(data, tmp, i1, i2, i4, tmp.numerics_h_ori1); \
|
|
if (tmp.numerics_h_ori1 >= 0) { \
|
|
Orientation(data, tmp, i3, i1, i4, tmp.numerics_h_ori1); \
|
|
if (tmp.numerics_h_ori1 >= 0) inside = true; \
|
|
} \
|
|
} \
|
|
}
|
|
|
|
|
|
/* */
|
|
/* this macro checks whether point i4 is inside of or on the boundary */
|
|
/* of the triangle i1, i2, i3. it also returns a classification if i4 is */
|
|
/* on the boundary of the triangle (except for the edge i2, i3). */
|
|
/* */
|
|
#define VtxInTriangle(data, tmp, i1, i2, i3, i4, inside, classifier) \
|
|
{ \
|
|
inside = false; \
|
|
Orientation(data, tmp, i2, i3, i4, tmp.numerics_h_ori1); \
|
|
if (tmp.numerics_h_ori1 >= 0) { \
|
|
Orientation(data, tmp, i1, i2, i4, tmp.numerics_h_ori1); \
|
|
if (tmp.numerics_h_ori1 > 0) { \
|
|
Orientation(data, tmp, i3, i1, i4, tmp.numerics_h_ori1); \
|
|
if (tmp.numerics_h_ori1 > 0) { \
|
|
inside = true; \
|
|
classifier = 0; \
|
|
} \
|
|
else if (tmp.numerics_h_ori1 == 0) { \
|
|
inside = true; \
|
|
classifier = 1; \
|
|
} \
|
|
} \
|
|
else if (tmp.numerics_h_ori1 == 0) { \
|
|
Orientation(data, tmp, i3, i1, i4, tmp.numerics_h_ori1); \
|
|
if (tmp.numerics_h_ori1 > 0) { \
|
|
inside = true; \
|
|
classifier = 2; \
|
|
} \
|
|
else if (tmp.numerics_h_ori1 == 0) { \
|
|
inside = true; \
|
|
classifier = 3; \
|
|
} \
|
|
} \
|
|
} \
|
|
}
|
|
|
|
|
|
#define ComputeRatio(A, B, C, RATIO) \
|
|
{ \
|
|
if (Min(B, C) <= ZERO) { \
|
|
RATIO = Min(B, C); \
|
|
} \
|
|
else { \
|
|
RATIO = CD_1_0 + (A - B - C) / (CD_2_0 * sqrt(B * C)); \
|
|
} \
|
|
}
|
|
|