SaraP
230 lines
9.2 KiB
C++
230 lines
9.2 KiB
C++
/*****************************************************************************/
|
|
/* */
|
|
/* Copyright (C) 1999-2025 M. Held */
|
|
/* */
|
|
/* This code is not in the public domain. All rights reserved! Please make */
|
|
/* sure to read the full copyright statement contained in "README.txt" or in */
|
|
/* the "main" file of this code, such as "main.cc". */
|
|
/* */
|
|
/*****************************************************************************/
|
|
/* */
|
|
/* Written by: Martin Held */
|
|
/* */
|
|
/* E-Mail: held@cs.sbg.ac.at */
|
|
/* Fax Mail: (+43 662) 8044-611 */
|
|
/* Voice Mail: (+43 662) 8044-6304 */
|
|
/* Snail Mail: Martin Held */
|
|
/* FB Informatik */
|
|
/* Universitaet Salzburg */
|
|
/* A-5020 Salzburg, Austria */
|
|
/* */
|
|
/*****************************************************************************/
|
|
|
|
/* */
|
|
/* get standard libraries */
|
|
/* */
|
|
#include <stdlib.h>
|
|
#include <stdio.h>
|
|
#include <assert.h>
|
|
#include <math.h>
|
|
#include <float.h>
|
|
|
|
/* */
|
|
/* get my header files */
|
|
/* */
|
|
#include "fpkernel.h"
|
|
#include "vronivector.h"
|
|
#include "vroni_object.h"
|
|
#include "defs.h"
|
|
#include "numerics.h"
|
|
#include "roots.h"
|
|
|
|
|
|
|
|
/* */
|
|
/* computes the center cntr and the radius (squared) r2 of the circle */
|
|
/* passing through the points pnts[i],pnts[j],pnts[k]. it returns false */
|
|
/* if the radius is too large to be computed numerically. */
|
|
/* */
|
|
vr_bool vroniObject::PntPntPntCntr(int i, int j, int k,
|
|
coord *cntr, double *r2)
|
|
{
|
|
double alpha, beta, gamma, delta, t;
|
|
coord AB, BC, CA, P, Q, U, V, A, B, C;
|
|
double aaa[2][2], bbb[2], xy[2];
|
|
int I0, I1, J0, J1, exists, m;
|
|
double ab, bc, ca;
|
|
|
|
assert(InPntsList(i));
|
|
assert(InPntsList(j));
|
|
assert(InPntsList(k));
|
|
|
|
/* */
|
|
/* we sort the three indices in order to guarantee that the center of */
|
|
/* three points is always computed consistently. */
|
|
/* */
|
|
SortThreeNumbers(i, j, k, m);
|
|
|
|
/* */
|
|
/* check whether two points are identical. */
|
|
/* */
|
|
if ((i == j) || (j == k)) {
|
|
VD_Dbg_Warning("PntPntPntCntr() - two identical points!");
|
|
return false;
|
|
}
|
|
|
|
/* */
|
|
/* check whether two points are nearly identical. */
|
|
/* */
|
|
A = GetPntCoords(i);
|
|
B = GetPntCoords(j);
|
|
C = GetPntCoords(k);
|
|
|
|
AB = VecSub(B, A);
|
|
ab = VecLen(AB);
|
|
assert(ab >= 0.0);
|
|
if (le(ab, ZERO)) {
|
|
VD_Dbg_Warning("PntPntPntCntr() - two points too close");
|
|
SetInvalidSites(i, PNT, j, PNT, ZERO);
|
|
return false;
|
|
}
|
|
BC = VecSub(C, B);
|
|
bc = VecLen(BC);
|
|
assert(bc >= 0.0);
|
|
if (le(bc, ZERO)) {
|
|
VD_Dbg_Warning("PntPntPntCntr() - two points too close");
|
|
SetInvalidSites(j, PNT, k, PNT, ZERO);
|
|
return false;
|
|
}
|
|
CA = VecSub(A, C);
|
|
ca = VecLen(CA);
|
|
assert(ca >= 0.0);
|
|
if (le(ca, ZERO)) {
|
|
VD_Dbg_Warning("PntPntPntCntr() - two points too close");
|
|
SetInvalidSites(i, PNT, k, PNT, ZERO);
|
|
return false;
|
|
}
|
|
|
|
/* */
|
|
/* normalize the vectors */
|
|
/* */
|
|
AB = VecDiv(ab, AB);
|
|
BC = VecDiv(bc, BC);
|
|
CA = VecDiv(ca, CA);
|
|
|
|
/* */
|
|
/* find those two vectors that are the least parallel. */
|
|
/* */
|
|
alpha = VecDotProd(AB, BC);
|
|
alpha = Abs(alpha);
|
|
beta = VecDotProd(BC, CA);
|
|
beta = Abs(beta);
|
|
gamma = VecDotProd(CA, AB);
|
|
gamma = Abs(gamma);
|
|
delta = Min3(alpha, beta, gamma);
|
|
|
|
/* */
|
|
/* assume that the angle between AB and AC is closest to a right */
|
|
/* angle. if ab > ca then we construct the line perpendicular to AB */
|
|
/* through the midpoint of A,B and intersect it with the corresponding */
|
|
/* line w.r.t. A,C. the point of intersection will be expressed in */
|
|
/* terms of the line perpendicular to AB. similar for the other cases. */
|
|
/* */
|
|
if (delta == gamma) {
|
|
P = MidPoint(A, B);
|
|
Q = VecCCW(AB);
|
|
U = MidPoint(A, C);
|
|
V = VecCCW(CA);
|
|
if (ab < ca) {
|
|
LineLineIntersection(aaa, bbb, xy, I0, J0, I1, J1, P, Q, U, V, *cntr,
|
|
exists);
|
|
}
|
|
else {
|
|
LineLineIntersection(aaa, bbb, xy, I0, J0, I1, J1, U, V, P, Q, *cntr,
|
|
exists);
|
|
}
|
|
}
|
|
else if (delta == alpha) {
|
|
P = MidPoint(A, B);
|
|
Q = VecCCW(AB);
|
|
U = MidPoint(B, C);
|
|
V = VecCCW(BC);
|
|
if (ab < bc) {
|
|
LineLineIntersection(aaa, bbb, xy, I0, J0, I1, J1, P, Q, U, V, *cntr,
|
|
exists);
|
|
}
|
|
else {
|
|
LineLineIntersection(aaa, bbb, xy, I0, J0, I1, J1, U, V, P, Q, *cntr,
|
|
exists);
|
|
}
|
|
}
|
|
else {
|
|
P = MidPoint(C, A);
|
|
Q = VecCCW(CA);
|
|
U = MidPoint(B, C);
|
|
V = VecCCW(BC);
|
|
if (ca < bc) {
|
|
LineLineIntersection(aaa, bbb, xy, I0, J0, I1, J1, P, Q, U, V, *cntr,
|
|
exists);
|
|
}
|
|
else {
|
|
LineLineIntersection(aaa, bbb, xy, I0, J0, I1, J1, U, V, P, Q, *cntr,
|
|
exists);
|
|
}
|
|
}
|
|
|
|
if (exists == 1) {
|
|
*r2 = PntPntDistSq(A, *cntr);
|
|
return true;
|
|
}
|
|
else {
|
|
/* */
|
|
/* we'll try to find a center using an alternative approach: we find */
|
|
/* a parameter t such that cntr = P + t * Q. thus, we require */
|
|
/* d(A, cntr) = d(C, cntr). */
|
|
/* */
|
|
P = MidPoint(A, B);
|
|
Q = VecCCW(AB); /* note: Q is a unit vector! */
|
|
V = VecSub(P, C);
|
|
delta = VecDotProd(Q, V) * 2.0;
|
|
|
|
if (ne(delta, TINY)) {
|
|
ab = ab * ab / 4.0;
|
|
t = (ab - VecLenSq(V)) / delta;
|
|
*cntr = RayPnt(P, Q, t);
|
|
*r2 = ab + t * t;
|
|
return true;
|
|
}
|
|
else {
|
|
/* */
|
|
/* We could not find a center. Let's check whether two points are */
|
|
/* very close. */
|
|
/* */
|
|
if (le(ab, ZERO_MAX)) {
|
|
VD_Dbg_Warning("PntPntPntCntr() - two points too close");
|
|
SetInvalidSites(i, PNT, j, PNT, ZERO_MAX);
|
|
return false;
|
|
}
|
|
if (le(bc, ZERO_MAX)) {
|
|
VD_Dbg_Warning("PntPntPntCntr() - two points too close");
|
|
SetInvalidSites(j, PNT, k, PNT, ZERO_MAX);
|
|
return false;
|
|
}
|
|
if (le(ca, ZERO_MAX)) {
|
|
VD_Dbg_Warning("PntPntPntCntr() - two points too close");
|
|
SetInvalidSites(i, PNT, k, PNT, ZERO_MAX);
|
|
return false;
|
|
}
|
|
}
|
|
}
|
|
|
|
VD_Warning("PntPntPntCntr() - cannot find center!");
|
|
#ifdef VRONI_DBG_WARN
|
|
printf("A = (%20.16f, %20.16f)\n", A.x, A.y);
|
|
printf("B = (%20.16f, %20.16f)\n", B.x, B.y);
|
|
printf("C = (%20.16f, %20.16f)\n", C.x, C.y);
|
|
#endif
|
|
|
|
return false;
|
|
}
|