/*****************************************************************************/ /* */ /* 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. */ /* */ /*****************************************************************************/ /* */ /* get standard libraries */ /* */ #include #include #include #include #include /* */ /* get my header files */ /* */ #include "fpkernel.h" #include "martin.h" #include "defs.h" #include "list.h" #include "header.h" #include "basic.h" #include "numerics.h" /* */ /* prototypes of functions provided in this file */ /* */ boolean SegIntersect(global_struct *all, int i1, int i2, int i3, int i4, int i5); boolean SegIntersection(global_struct *all, int i1, int i2, int i3, int i4, int i5); machine_double GetRatio(datadef *data, int i, int j, int k); machine_double GetDoubleRatio(datadef *data, int i, int j, int k, int m); int SpikeAngle(listdef *list, datadef *data, int i, int j, int k, list_ind ind); int RecSpikeAngle(listdef *list, datadef *data, int i1, int i2, int i3, list_ind ind1, list_ind ind3); double Angle(point p, point p1, point p2); boolean PntInTri(datadef *data, int i1, int i2, int i3, int i4, double *area); #if defined(WITH_MPFRBACKEND) || defined(WITH_EXPR_WRAPPER) double EPS = 1.0e-15; double ZERO = 1.0e-16; double M_PI = 3.14159265358979323846; double M_2PI = 6.28318530717958647693; double M_PI_180 = 0.01745329251994329576; void set_mpfrbackend_prec(mpfr_prec_t mpfr_prec, boolean verbose) { mpfr_set_default_prec(mpfr_prec); SET_PRECISION(M_PI, mpfr_prec); SET_MPFR_KERNEL_PI(M_PI); SET_PRECISION(M_2PI, mpfr_prec); M_2PI = M_PI * 2; SET_PRECISION(M_PI_180, mpfr_prec); M_PI_180 = M_PI / 180; SET_CORE_MPFR_KERNEL(EPS, 0, 1e-15 / (exp2(100 * (sqrt(double(mpfr_prec) / 53) - 1))), "MP_EPS"); SET_CORE_MPFR_KERNEL(ZERO, 0, 1e-16 / (exp2(100 * (sqrt(double(mpfr_prec) / 53) - 1))), "MP_ZERO"); if(verbose) { FP_printf("Setting EPS to %e\n", FP_PRNTARG(EPS)); FP_printf("Setting ZERO to %e\n", FP_PRNTARG(ZERO)); } } #endif /* */ /* checks whether the line segments i1, i2 and i3, i4 intersect. no */ /* intersection is reported if they intersect at a common vertex. */ /* the function assumes that i1 <= i2 and i3 <= i4. if i3 or i4 lies */ /* on i1, i2 then an intersection is reported, but no intersection is */ /* reported if i1 or i2 lies on i3, i4. this function is not symmetric! */ /* */ boolean SegIntersect(global_struct *all, int i1, int i2, int i3, int i4, int i5) { datadef *data = &all->c_data; griddef *grid = &all->c_grid; tmp_data_def tmp_data; int ori1, ori2, ori3, ori4; assert(InPointsList(data, i1)); assert(InPointsList(data, i2)); assert(InPointsList(data, i3)); assert(InPointsList(data, i4)); assert(i1 <= i2); assert(i3 <= i4); if ((i1 == i2) || (i3 == i4)) return false; if ((i1 == i3) && (i2 == i4)) return true; if ((i3 == i5) || (i4 == i5)) { #ifdef PARTITION_FIST #pragma omp atomic #endif ++grid->ident_cntr; } Orientation(data, tmp_data, i1, i2, i3, ori3); Orientation(data, tmp_data, i1, i2, i4, ori4); if (((ori3 == 1) && (ori4 == 1)) || ((ori3 == -1) && (ori4 == -1))) return false; if (ori3 == 0) { if (StrictlyInBetween(i1, i2, i3)) return true; if (ori4 == 0) { if (StrictlyInBetween(i1, i2, i4)) return true; } else return false; } else if (ori4 == 0) { if (StrictlyInBetween(i1, i2, i4)) return true; else return false; } Orientation(data, tmp_data,i3, i4, i1, ori1); Orientation(data, tmp_data,i3, i4, i2, ori2); if (((ori1 <= 0) && (ori2 <= 0)) || ((ori1 >= 0) && (ori2 >= 0))) return false; return true; } /* */ /* checks whether the line segments i1, i2 and i3, i4 intersect. no */ /* intersection is reported if they intersect at a common vertex or if they */ /* are collinear. the function assumes that i1 <= i2 and i3 <= i4. */ /* */ boolean SegIntersection(global_struct *all, int i1, int i2, int i3, int i4, int i5) { datadef *data = &all->c_data; griddef *grid = &all->c_grid; tmp_data_def tmp_data; int ori1, ori2, ori3, ori4; assert(InPointsList(data, i1)); assert(InPointsList(data, i2)); assert(InPointsList(data, i3)); assert(InPointsList(data, i4)); assert(i1 <= i2); assert(i3 <= i4); if ((i1 == i2) || (i3 == i4)) return false; if (i1 == i3) return (i2 == i4); else if (i2 == i4) return false; if ((i3 == i5) || (i4 == i5)) ++grid->ident_cntr; Orientation(data, tmp_data, i1, i2, i3, ori3); Orientation(data, tmp_data, i1, i2, i4, ori4); if (((ori3 == 1) && (ori4 == 1)) || ((ori3 == -1) && (ori4 == -1))) return false; if (ori3 == 0) { if (StrictlyInBetween(i1, i2, i3)) return true; if (ori4 == 0) { if (StrictlyInBetween(i1, i2, i4)) return true; } else return false; } else if (ori4 == 0) { if (StrictlyInBetween(i1, i2, i4)) return true; else return false; } Orientation(data, tmp_data, i3, i4, i1, ori1); Orientation(data, tmp_data, i3, i4, i2, ori2); if (((ori1 <= 0) && (ori2 <= 0)) || ((ori1 >= 0) && (ori2 >= 0))) return false; return true; } /* */ /* this function returns a quality measure ratio of a triangle i, j, k, */ /* where it is assumed that i, j defines the new diagonal. the smaller */ /* ratio is, the earlier the triangle will be clipped. (all ratios are */ /* positive.) */ /* */ machine_double GetRatio(datadef *data, int i, int j, int k) { machine_double a, b, c, ratio = CD_0_0; point p, q, r; md_point side; assert(InPointsList(data, i)); assert(InPointsList(data, j)); assert(InPointsList(data, k)); p = data->points[i]; q = data->points[j]; /* */ /* compute the squared distance a of the diagonal of the new triangle */ /* */ PntPntDistSqd(p, q, side, a); if (a <= ZERO) return (a); r = data->points[k]; PntPntDistSqd(p, r, side, b); PntPntDistSqd(r, q, side, c); /* */ /* compute the ratio */ /* */ #ifdef BOUNDARY /* */ /* ratio equals the two edge lengths of the sides of the ear */ /* */ ratio = sqrt(b) + sqrt(c); #else /* */ /* we consider the angle, gamma, opposite of a and compute a measure, */ /* ratio, for its size. the smaller the angle and the shorter a, the */ /* smaller (i.e., "better") the ratio is: */ /* 0.0 <= ratio := 1.0 - cos[gamma] <= 2.0 */ /* this value is then scaled according to the length of a, in order to */ /* avoid clipping long diagonals in early stages of the triangulation. */ /* */ ComputeRatio(a, b, c, ratio); #endif assert(gt(data->bbox_diagonal_sqd, ZERO)); return (ratio * ratio * (1.0 + sqrt((a / data->bbox_diagonal_sqd)))); } /* */ /* this function computes a quality measure ratio of the two triangles */ /* i, j, k and i, j, m, where it is assumed that i, j defines the new */ /* diagonal, and m is a point opposite of the diagonal. the smaller */ /* ratio is, the earlier the triangle will be clipped. (all ratios are */ /* positive.) */ /* */ machine_double GetDoubleRatio(datadef *data, int i, int j, int k, int m) { machine_double a, b, c, ratio1 = 0.0, ratio2 = 0.0; point p, q, r, u; md_point side; assert(InPointsList(data, i)); assert(InPointsList(data, j)); assert(InPointsList(data, k)); p = data->points[i]; q = data->points[j]; /* */ /* compute the squared distance a of the diagonal of the new triangle */ /* */ PntPntDistSqd(p, q, side, a); if (a <= ZERO) return (a); r = data->points[k]; PntPntDistSqd(p, r, side, b); PntPntDistSqd(r, q, side, c); /* */ /* compute the ratio of i, j, k: */ /* */ ComputeRatio(a, b, c, ratio1); /* */ /* compute the ratio of i, j, m: */ /* */ u = data->points[m]; PntPntDistSqd(p, u, side, b); PntPntDistSqd(u, q, side, c); ComputeRatio(a, b, c, ratio2); ratio1 = Max(ratio1, ratio2); return (ratio1 * ratio1 * (1.0 + sqrt((a / data->bbox_diagonal_sqd)))); } int SpikeAngle(listdef *list, datadef *data, int i, int j, int k, list_ind ind) { list_ind ind1, ind2, ind3; #ifndef NDEBUG int i1, i2, i3; #endif int RecSpikeAngle(listdef *list, datadef *data, int i1, int i2, int i3, list_ind ind1, list_ind ind3); assert(InPointsList(data, i)); assert(InPointsList(data, j)); assert(InPointsList(data, k)); ind2 = ind; #ifndef NDEBUG i2 = GetIndex(list, ind2); assert(i2 == j); #endif ind1 = GetPrevNode(list, ind2); #ifndef NDEBUG i1 = GetIndex(list, ind1); assert(i1 == i); #endif ind3 = GetNextNode(list, ind2); #ifndef NDEBUG i3 = GetIndex(list, ind3); assert(i3 == k); #endif return RecSpikeAngle(list, data, i, j, k, ind1, ind3); } int RecSpikeAngle(listdef *list, datadef *data, int i1, int i2, int i3, list_ind ind1, list_ind ind3) { int ori, ori1, ori2, i0, ii1, ii2; point pq, pr; double dot; tmp_data_def tmp_data; if (GetAngle(list, ind1) == -2) return -2; if (ind1 == ind3) { /* */ /* all points are collinear??? well, then it does not really matter */ /* which angle is returned. perhaps, -2 is the best bet as my code */ /* likely regards this contour as a hole. */ /* */ return -2; } if (i1 != i3) { if (i1 < i2) { ii1 = i1; ii2 = i2; } else { ii1 = i2; ii2 = i1; } if (InBetween(ii1, ii2, i3)) { i2 = i3; ind3 = GetNextNode(list, ind3); i3 = GetIndex(list, ind3); if (ind1 == ind3) return -2; Orientation(data, tmp_data, i1, i2, i3, ori); if (ori > 0) return 2; else if (ori < 0) return -2; else return RecSpikeAngle(list, data, i1, i2, i3, ind1, ind3); } else { i2 = i1; ind1 = GetPrevNode(list, ind1); i1 = GetIndex(list, ind1); if (ind1 == ind3) return -2; Orientation(data, tmp_data, i1, i2, i3, ori); if (ori > 0) return 2; else if (ori < 0) return -2; else return RecSpikeAngle(list, data, i1, i2, i3, ind1, ind3); } } else { i0 = i2; i2 = i1; ind1 = GetPrevNode(list, ind1); i1 = GetIndex(list, ind1); if (ind1 == ind3) return -2; ind3 = GetNextNode(list, ind3); i3 = GetIndex(list, ind3); if (ind1 == ind3) return -2; Orientation(data, tmp_data, i1, i2, i3, ori); if (ori > 0) { Orientation(data, tmp_data, i1, i2, i0, ori1); if (ori1 > 0) { Orientation(data, tmp_data, i2, i3, i0, ori2); if (ori2 > 0) return -2; } return 2; } else if (ori < 0) { Orientation(data, tmp_data, i2, i1, i0, ori1); if (ori1 > 0) { Orientation(data, tmp_data, i3, i2, i0, ori2); if (ori2 > 0) return 2; } return -2; } else { VectorSub2D(data->points[i1], data->points[i2], pq); VectorSub2D(data->points[i3], data->points[i2], pr); dot = DotProduct2D(pq, pr); if (dot < C_0_0) { Orientation(data, tmp_data, i2, i1, i0, ori); if (ori > 0) return 2; else return -2; } else { return RecSpikeAngle(list, data, i1, i2, i3, ind1, ind3); } } } } /* */ /* computes the signed angle between p, p1 and p, p2. */ /* */ /* warning: this function does not handle a 180-degree angle correctly! */ /* (this is no issue in our application, as we will always compute */ /* the angle centered at the mid-point of a valid diagonal.) */ /* */ double Angle(point p, point p1, point p2) { int sign; double macro_var; double angle1, angle2, angle; point v1, v2; sign = SignEps(macro_var,Det2D(p2, p, p1), EPS); if (sign == 0) return C_0_0; VectorSub2D(p1, p, v1); VectorSub2D(p2, p, v2); angle1 = atan2(TO_MDOUBLE(v1.y), TO_MDOUBLE(v1.x)); angle2 = atan2(TO_MDOUBLE(v2.y), TO_MDOUBLE(v2.x)); #ifdef WITH_COREBACKEND if (angle1 < C_0_0) angle1 += C_2_0 * acos(CD_m1_0); if (angle2 < C_0_0) angle2 += C_2_0 * acos(CD_m1_0); #else if (angle1 < C_0_0) angle1 += M_2PI; if (angle2 < C_0_0) angle2 += M_2PI; #endif angle = angle1 - angle2; #ifdef WITH_COREBACKEND if (angle > acos(CD_m1_0)) angle = C_2_0 * acos(CD_m1_0) - angle; else if (-angle > acos(CD_m1_0)) angle = C_2_0 * acos(CD_m1_0) + angle; #else if (angle > M_PI) angle = M_2PI - angle; else if (-angle > M_PI) angle = M_2PI + angle; #endif if (sign == 1) { if (angle < C_0_0) return -angle; else return angle; } else { if (angle > C_0_0) return -angle; else return angle; } } boolean PntInTri(datadef *data, int i1, int i2, int i3, int i4, double *area) { tmp_data_def tmp_data; StableDet2D(data, tmp_data, i2, i3, i4, *area); if (ge(*area, EPS)) { StableDet2D(data, tmp_data, i1, i2, i4, *area); if (ge(*area, EPS)) { StableDet2D(data, tmp_data, i3, i1, i4, *area); if (ge(*area, EPS)) return true; } *area = -fpkernel_DBL_MAX; } return false; } boolean PntInTriClass(datadef *data, int i1, int i2, int i3, int i4, double *area, int *pnt_on_edge) { tmp_data_def tmp_data; *pnt_on_edge = 0; StableDet2D(data, tmp_data, i2, i3, i4, *area); if (eq(*area, EPS)) { if (i2 <= i3) { if (InBetween(i2, i3, i4)) { *pnt_on_edge = 1; return true; } else { *area = -fpkernel_DBL_MAX; return false; } } else { if (InBetween(i3, i2, i4)) { *pnt_on_edge = 1; return true; } else { *area = -fpkernel_DBL_MAX; return false; } } } else if (*area < C_0_0) { return false; } else { StableDet2D(data, tmp_data, i1, i2, i4, *area); if (eq(*area, EPS)) { if (i1 <= i2) { if (InBetween(i1, i2, i4)) { *pnt_on_edge = 2; return true; } else { *area = -fpkernel_DBL_MAX; return false; } } else { if (InBetween(i2, i1, i4)) { *pnt_on_edge = 2; return true; } else { *area = -fpkernel_DBL_MAX; return false; } } } else if (*area < C_0_0) { *area = -fpkernel_DBL_MAX; return false; } else { StableDet2D(data, tmp_data, i3, i1, i4, *area); if (eq(*area, EPS)) { if (i3 <= i1) { if (InBetween(i3, i1, i4)) { *pnt_on_edge = 2; return true; } else { *area = TO_REAL(-fpkernel_DBL_MAX); return false; } } else { if (InBetween(i1, i3, i4)) { *pnt_on_edge = 2; return true; } else { *area = -fpkernel_DBL_MAX; return false; } } } else if (*area < C_0_0) { *area = -fpkernel_DBL_MAX; return false; } else { return true; } } } }