Egalware/fist
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
cbec90699f385a00b9672896b7b4eff5c3ea1414
fist/numerics.cpp
T
SaraP cbec90699f FIST 6.8 : 
- creato il progetto in Visual Studio per compilare come libreria statica
- modifiche al codice originale per integrarlo nelle nostre librerie.
2025-03-04 16:19:35 +01:00

590 lines
20 KiB
C++
Raw Blame History

/*****************************************************************************/
/* */
/* 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 <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <float.h>
#include <assert.h>
/* */
/* 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;
}
}
}
}
Reference in New Issue View Git Blame Copy Permalink