SaraP
324 lines
13 KiB
C++
324 lines
13 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 r2 of the circle tangential */
|
|
/* to the segment segs[i] and passing through the points pnts[j],pnts[k]. */
|
|
/* it returns false if the radius is too large to be computed numerically, */
|
|
/* or if it is obvious that problems have occurred. if two centers exist */
|
|
/* then we choose the center that lies on the Voronoi edge e. the vr_bool */
|
|
/* flag problematic is set if a center was computed but its correctness is */
|
|
/* doubtful. */
|
|
/* */
|
|
vr_bool vroniObject::SegPntPntCntr(int i, int j, int k, int e, coord *cntr,
|
|
double *r2, vr_bool *problematic)
|
|
{
|
|
double dist, a, b, c;
|
|
coord p, q, Q, V;
|
|
double t, delta, lgth, eps, d_p_q, d2;
|
|
double alpha, beta, x, y, z, alpha_p, alpha_q;
|
|
coord m, v, n, u, tempc;
|
|
int n1, ii;
|
|
double roots[2];
|
|
#ifdef TRACE
|
|
vr_bool debug = true;
|
|
#endif
|
|
int num_sol = 0, best_sol = NIL;
|
|
double radii[VRONI_MAXSOL];
|
|
coord centers[VRONI_MAXSOL];
|
|
|
|
assert(InSegsList(i));
|
|
assert(InPntsList(j));
|
|
assert(InPntsList(k));
|
|
assert(!(eq(GetSegLgth(i), ZERO)));
|
|
|
|
#ifdef TRACE
|
|
if ((i == 6) && (j == 10) && (k == 11)) {
|
|
printf("\nSegPntPntCntr() - seg %d\n", i);
|
|
printf("start - (%24.20f %24.20f)\n", GetSegStartCoord(i).x,
|
|
GetSegStartCoord(i).y);
|
|
printf("end - (%24.20f %24.20f)\n", GetSegEndCoord(i).x,
|
|
GetSegEndCoord(i).y);
|
|
printf("pnt[%d] - (%24.20f %24.20f)\n", j, pnts[j].p.x, pnts[j].p.y);
|
|
printf("pnt[%d] - (%24.20f %24.20f)\n", k, pnts[k].p.x, pnts[k].p.y);
|
|
}
|
|
#endif
|
|
|
|
*problematic = false;
|
|
|
|
/* */
|
|
/* check whether the two points are identical. */
|
|
/* */
|
|
if (j == k) {
|
|
VD_Dbg_Warning("SegPntPntCntr() - two identical points!");
|
|
return false;
|
|
}
|
|
|
|
/* */
|
|
/* we sort the two point indices in order to guarantee that the center of */
|
|
/* these three sites is always computed consistently. */
|
|
/* */
|
|
SortTwoNumbers(j, k, n1);
|
|
|
|
/* */
|
|
/* check whether one of the points is an endpoint of segs[i]. */
|
|
/* */
|
|
if (IsSegStartPnt(i, k) || IsSegEndPnt(i, k)) VroniSwap(j, k, n1);
|
|
|
|
p = GetPntCoords(j);
|
|
q = GetPntCoords(k);
|
|
d2 = PntPntDistSq(p, q);
|
|
d_p_q = sqrt(d2);
|
|
lgth = GetSegLgth(i);
|
|
eps = ZERO;
|
|
|
|
while (eps <= ZERO_MAX) {
|
|
|
|
num_sol = 0;
|
|
|
|
if (eq(d_p_q, eps)) {
|
|
/* */
|
|
/* both points are very close to each other. */
|
|
/* */
|
|
VD_Dbg_Warning("SegPntPntCntr() - points are too close!");
|
|
SetInvalidSites(j, PNT, k, PNT, eps);
|
|
restart = false;
|
|
return false;
|
|
}
|
|
|
|
if (IsSegStartPnt(i, j) || IsSegEndPnt(i, j)) {
|
|
if (IsSegStartPnt(i, k) || IsSegEndPnt(i, k)) {
|
|
VD_Dbg_Warning("SegPntPntCntr() - both points are endpoints!");
|
|
*problematic = true;
|
|
}
|
|
else {
|
|
GetSegEqnData(i, &a, &b, &c);
|
|
dist = PntLineDist(a, b, c, q);
|
|
|
|
if (eq(dist, eps)) {
|
|
/* */
|
|
/* both points lie on the line! this doesn't make much sense! */
|
|
/* either the input data is corrupted, or we've run into a */
|
|
/* problem... */
|
|
/* */
|
|
Q = GetSegStartCoord(i);
|
|
if (IsInSegConeZero(Q, q, a, b, lgth, eps) == 0) {
|
|
SetInvalidSites(i, SEG, k, PNT, eps);
|
|
restart = false;
|
|
return false;
|
|
}
|
|
else {
|
|
VD_Dbg_Warning("SegPntPntCntr() - both points lie on seg!");
|
|
*problematic = true;
|
|
}
|
|
}
|
|
else {
|
|
V = VecSub(p, q);
|
|
if (dist > 0.0) {
|
|
Q = MakeVec( -a, -b);
|
|
}
|
|
else {
|
|
Q = MakeVec( a, b);
|
|
}
|
|
delta = VecDotProd(Q, V);
|
|
|
|
if (!eq(delta, eps)) {
|
|
t = (- 0.5) * d2 / delta;
|
|
centers[0] = RayPnt(p, Q, t);
|
|
radii[0] = Abs(t);
|
|
num_sol = 1;
|
|
}
|
|
else {
|
|
/* */
|
|
/* the bisector of p,q is parallel to the line!? */
|
|
/* */
|
|
VD_Dbg_Warning("SegPntPntCntr() - cannot find center!");
|
|
*problematic = true;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
else {
|
|
assert(!(eq(GetSegLgth(i), ZERO)));
|
|
p = GetPntCoords(j);
|
|
q = GetPntCoords(k);
|
|
v = VecSub(q, p);
|
|
dist = VecLenSq(v) / 4.0;
|
|
v = VecMult(0.5, v);
|
|
n = VecCCW(v);
|
|
|
|
GetSegEqnData(i, &a, &b, &c);
|
|
|
|
alpha_p = PntLineDist(a, b, c, p);
|
|
if (eq(alpha_p, eps)) {
|
|
/* */
|
|
/* point p lies on the line!? */
|
|
/* */
|
|
u = GetSegStartCoord(i);
|
|
if (IsInSegConeZero(u, p, a, b, lgth, eps) == 0) {
|
|
SetInvalidSites(i, SEG, j, PNT, eps);
|
|
restart = false;
|
|
return false;
|
|
}
|
|
}
|
|
alpha_q = PntLineDist(a, b, c, q);
|
|
if (eq(alpha_q, eps)) {
|
|
/* */
|
|
/* point q lies on the line!? */
|
|
/* */
|
|
u = GetSegStartCoord(i);
|
|
if (IsInSegConeZero(u, q, a, b, lgth, eps) == 0) {
|
|
SetInvalidSites(i, SEG, k, PNT, eps);
|
|
restart = false;
|
|
return false;
|
|
}
|
|
}
|
|
|
|
alpha = (alpha_p + alpha_q) / 2.0;
|
|
beta = n.x * a + n.y * b;
|
|
y = 2 * alpha * beta;
|
|
beta = beta * beta;
|
|
alpha = alpha * alpha;
|
|
x = beta - dist;
|
|
z = alpha - dist;
|
|
m = MidPoint(p, q);
|
|
|
|
/* */
|
|
/* solve the quadratic equation. */
|
|
/* */
|
|
Roots2abc(x, y, z, roots, num_sol);
|
|
|
|
for (ii = 0; ii < num_sol; ++ii) {
|
|
centers[ii] = RayPnt(m, n, roots[ii]);
|
|
radii[ii] = PntPntDist(centers[ii], p);
|
|
#ifdef TRACE
|
|
printf("c[%d]: (%20.16f %20.16f)\n", ii,
|
|
centers[ii].x, centers[ii].y);
|
|
printf("r[%d] = %20.16f\n", ii, radii[ii]);
|
|
#endif
|
|
}
|
|
}
|
|
|
|
#ifdef TRACE
|
|
if (debug) {
|
|
printf("Eqn x = %24.20f\n", x);
|
|
printf("Eqn y = %24.20f\n", y);
|
|
printf("Eqn z = %24.20f\n", z);
|
|
printf("num_sol: %d\n", num_sol);
|
|
printf("t1 = %20.16f\n", roots[0]);
|
|
printf("t2 = %20.16f\n", roots[1]);
|
|
printf("p: (%20.16f, %20.16f)\n", p.x, p.y);
|
|
printf("q: (%20.16f, %20.16f)\n", q.x, q.y);
|
|
printf("start: %d = (%20.16f, %20.16f); rad = %20.16f\n",
|
|
GetStartNode(e),nodes[GetStartNode(e)].p.x,
|
|
nodes[GetStartNode(e)].p.y, nodes[GetStartNode(e)].r2);
|
|
printf("end: %d = (%20.16f, %20.16f); rad = %20.16f\n", GetEndNode(e),
|
|
nodes[GetEndNode(e)].p.x, nodes[GetEndNode(e)].p.y,
|
|
nodes[GetEndNode(e)].r2);
|
|
if (IsNodeDeleted(GetStartNode(e))) printf("start node deleted\n");
|
|
if (IsNodeDeleted(GetEndNode(e))) printf("end node deleted\n");
|
|
}
|
|
#endif
|
|
|
|
/* */
|
|
/* we have up to two possible centers. we select the center that */
|
|
/* (ideally) lies within the cones of influence of the segments and */
|
|
/* on the Voronoi edge e. */
|
|
/* */
|
|
best_sol = ClassifySolutions(i, j, k, e, SEG, PNT, PNT,
|
|
false, true, true, centers, radii,
|
|
&num_sol, eps, problematic);
|
|
assert((0 <= best_sol) && (best_sol < num_sol));
|
|
#ifdef TRACE
|
|
printf("best_sol = %d, problematic = %d\n", best_sol, *problematic);
|
|
#endif
|
|
|
|
if (!*problematic) {
|
|
*cntr = centers[best_sol];
|
|
*r2 = radii[best_sol];
|
|
return true;
|
|
}
|
|
else {
|
|
eps *= 10.0;
|
|
}
|
|
}
|
|
|
|
if (eq(lgth, ZERO_MAX)) {
|
|
/* */
|
|
/* this is a seg with a very small length; we replace it */
|
|
/* */
|
|
VD_Dbg_Warning("SegPntPntCntr() - seg with small length!");
|
|
ReplaceSeg(i);
|
|
restart = true;
|
|
return false;
|
|
}
|
|
else if (eq(d_p_q, ZERO_MAX)) {
|
|
/* */
|
|
/* both points are very close to each other. */
|
|
/* */
|
|
VD_Dbg_Warning("SegPntPntCntr() - points are too close!");
|
|
SetInvalidSites(j, PNT, k, PNT, ZERO_MAX);
|
|
restart = false;
|
|
|
|
return false;
|
|
}
|
|
|
|
#ifdef VRONI_DBG_WARN
|
|
printf("\nCenter not computed for seg-pnt-pnt:\n");
|
|
Q = GetSegStartCoord(i);
|
|
V = GetSegEndCoord(i);
|
|
printf("%20.16f %20.16f %20.16f %20.16f\n", Q.x, Q.y, V.x, V.y);
|
|
p = GetPntCoords(j);
|
|
q = GetPntCoords(k);
|
|
printf("%20.16f %20.16f\n", p.x, p.y);
|
|
printf("%20.16f %20.16f\n", q.x, q.y);
|
|
#endif
|
|
|
|
*cntr = centers[best_sol];
|
|
*r2 = radii[best_sol];
|
|
restart = false;
|
|
|
|
return false;
|
|
}
|