SaraP
337 lines
14 KiB
C++
337 lines
14 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"
|
|
|
|
/* #define TRACE */
|
|
|
|
|
|
#define Determine2Sidedness(ai, bi, ci, pnt, dist, eps) \
|
|
{\
|
|
dist = PntLineDist(ai, bi, ci, pnt); \
|
|
if (eq(dist, eps)) { \
|
|
restart = true; \
|
|
} \
|
|
else if (dist < 0.0) { \
|
|
ai = -ai; \
|
|
bi = -bi; \
|
|
ci = -ci; \
|
|
} \
|
|
}
|
|
|
|
|
|
/* */
|
|
/* computes the center cntr and the radius r2 of the circle tangential */
|
|
/* to the segments segs[i], segs[j] and passing through the point pnts[k].*/
|
|
/* it returns false if the radius is too large to be computed numerically, */
|
|
/* or if it is obvious that problems have occurred. */
|
|
/* */
|
|
vr_bool vroniObject::SegSegPntCntr(int i, int j, int k, int e, coord *cntr,
|
|
double *r2, vr_bool *problematic, int *site)
|
|
{
|
|
coord p, q, u, v, w, Q, centers[VRONI_MAXSOL];
|
|
int n1;
|
|
double ai, bi, ci, aj = 0.0, bj = 0.0, cj = 0.0;
|
|
double alpha, beta, dist, x, y, z;
|
|
double roots[2], eps;
|
|
#ifdef TRACE
|
|
vr_bool debug = false;
|
|
coord pi, pj;
|
|
#endif
|
|
int i_in, j_in;
|
|
int num_sol = 0, best_sol = NIL, m;
|
|
double radii[VRONI_MAXSOL];
|
|
double lgth_i, lgth_j;
|
|
|
|
|
|
assert(InSegsList(i));
|
|
assert(InSegsList(j));
|
|
assert(InPntsList(k));
|
|
assert(!(eq(GetSegLgth(i), ZERO)));
|
|
assert(!(eq(GetSegLgth(j), ZERO)));
|
|
|
|
eps = ZERO;
|
|
i_in = i;
|
|
j_in = j;
|
|
|
|
*site = NIL;
|
|
*problematic = false;
|
|
|
|
/* */
|
|
/* we sort the two segment indices in order to guarantee that the center */
|
|
/* of these three sites is always computed consistently. */
|
|
/* */
|
|
SortTwoNumbers(i, j, n1);
|
|
|
|
#ifdef TRACE
|
|
if (debug) {
|
|
printf("SegSegPntCntr - seg1=%d seg2=%d\n", i, j);
|
|
p = GetSegStartCoord(i);
|
|
printf("start (%d) - (%20.16f %20.16f)\n", i, p.x, p.y);
|
|
p = GetSegEndCoord(i);
|
|
printf(" end (%d) - (%20.16f %20.16f)\n", i, p.x, p.y);
|
|
p = GetSegStartCoord(j);
|
|
printf("start (%d) - (%20.16f %20.16f)\n", j, p.x, p.y);
|
|
p = GetSegEndCoord(j);
|
|
printf(" end (%d) - (%20.16f %20.16f)\n", j, p.x, p.y);
|
|
printf("SegSegPntCntr - pnt %d\n", k);
|
|
printf("pnt[%d] - (%20.16f %20.16f)\n", k, pnts[k].p.x, pnts[k].p.y);
|
|
}
|
|
#endif
|
|
|
|
/* */
|
|
/* check whether the two segments are identical. */
|
|
/* */
|
|
if (i == j) {
|
|
VD_Dbg_Warning("SegSegPntCntr() - the same segment!");
|
|
return false;
|
|
}
|
|
n1 = Get1stVtx(i, SEG);
|
|
if (IsSegStartPnt(j, n1) || IsSegEndPnt(j, n1)) {
|
|
n1 = Get2ndVtx(i, SEG);
|
|
if (IsSegStartPnt(j, n1) || IsSegEndPnt(j, n1)) {
|
|
SetInvalidSites(i, SEG, j, SEG, ZERO);
|
|
restart = false;
|
|
return false; /* duplicate segs ! */
|
|
}
|
|
}
|
|
|
|
/* */
|
|
/* check whether the point pnts[k] is an endpoint of segs[i] or */
|
|
/* segs[j]. */
|
|
/* */
|
|
p = GetPntCoords(k);
|
|
if (IsSegStartPnt(j, k) || IsSegEndPnt(j, k)) VroniSwap(i, j, n1);
|
|
if (IsSegStartPnt(j, k) || IsSegEndPnt(j, k)) {
|
|
/* */
|
|
/* pnts[k] is an endpoint of both segs[i] and segs[j]. */
|
|
/* */
|
|
*cntr = p;
|
|
*r2 = 0.0;
|
|
*site = k;
|
|
return true;
|
|
}
|
|
|
|
lgth_i = GetSegLgth(i);
|
|
lgth_j = GetSegLgth(j);
|
|
|
|
while (eps <= ZERO_MAX) {
|
|
|
|
num_sol = 0;
|
|
restart = false;
|
|
|
|
if (IsSegStartPnt(i, k) || IsSegEndPnt(i, k)) {
|
|
/* */
|
|
/* compute the center cntr and the radius r2 of the circle */
|
|
/* tangential to the segments segs[i], segs[j] and passing */
|
|
/* through the point pnts[k], where pnts[k] is an endpoint of */
|
|
/* segs[i]. */
|
|
/* */
|
|
/* get the precomputed line data. */
|
|
/* */
|
|
GetSegEqnData(i, &ai, &bi, &ci);
|
|
if (IsLftRgtSite(e, i , SEG)) {
|
|
n1 = GetStartNode(e);
|
|
if (!IsNodeDeleted(n1)) n1 = GetEndNode(e);
|
|
w = GetNodeCoord(n1);
|
|
Determine2Sidedness(ai, bi, ci, w, dist, eps);
|
|
n1 = 1;
|
|
}
|
|
else {
|
|
n1 = 0;
|
|
}
|
|
GetSegEqnData(j, &aj, &bj, &cj);
|
|
w = MakeVec( ai, bi);
|
|
|
|
/* */
|
|
/* we intersect the bisector between segs[i] and p with the */
|
|
/* offset of segs[j] in order to find the appropriate offset */
|
|
/* value. in general, we will get two solutions. */
|
|
/* */
|
|
if (!IntersectRayOffsetSeg(p, w, aj, bj, cj, n1, 1, centers, roots,
|
|
&num_sol)) {
|
|
q = GetSegStartCoord(j);
|
|
if (IsInSegConeZero(q, p, aj, bj, lgth_j, eps) == 0) {
|
|
SetInvalidSites(j, SEG, k, PNT, eps);
|
|
return false;
|
|
}
|
|
}
|
|
}
|
|
else {
|
|
/* */
|
|
/* get the precomputed line data and determine sidedness of point */
|
|
/* p relative to the lines. */
|
|
/* */
|
|
n1 = GetStartNode(e);
|
|
if (IsNodeDeleted(n1)) n1 = GetEndNode(e);
|
|
GetSegEqnData(i, &ai, &bi, &ci);
|
|
Determine2Sidedness(ai, bi, ci, p, dist, eps);
|
|
if (restart) {
|
|
u = GetSegStartCoord(i);
|
|
if (IsInSegConeZero(u, p, ai, bi, lgth_i, eps) == 0) {
|
|
SetInvalidSites(i, SEG, k, PNT, eps);
|
|
VD_Dbg_Warning("SegSegPntCntr() - point on segment!");
|
|
restart = false;
|
|
return false;
|
|
}
|
|
restart = false;
|
|
q = GetNodeCoord(n1);
|
|
Determine2Sidedness(ai, bi, ci, q, dist, eps);
|
|
}
|
|
if (!restart) {
|
|
GetSegEqnData(j, &aj, &bj, &cj);
|
|
Determine2Sidedness(aj, bj, cj, p, dist, eps);
|
|
if (restart) {
|
|
u = GetSegStartCoord(j);
|
|
if (IsInSegConeZero(u, p, aj, bj, lgth_j, eps) == 0) {
|
|
SetInvalidSites(j, SEG, k, PNT, eps);
|
|
VD_Dbg_Warning("SegSegPntCntr() - point on segment!");
|
|
restart = false;
|
|
return false;
|
|
}
|
|
restart = false;
|
|
q = GetNodeCoord(n1);
|
|
Determine2Sidedness(aj, bj, cj, q, dist, eps);
|
|
}
|
|
}
|
|
|
|
if (!restart) {
|
|
/* */
|
|
/* we compute the bisector Q + lambda * w of the two lines. */
|
|
/* */
|
|
if (!SegSegBisector(i, j, ai, bi, aj, bj, cj, &Q, &w)) {
|
|
VD_Dbg_Warning("SegSegPntCntr() - cannot compute bisector!");
|
|
}
|
|
else {
|
|
/* */
|
|
/* intersect the line and the circle by solving a */
|
|
/* second-degree equation. */
|
|
/* */
|
|
u = VecSub(Q, p);
|
|
alpha = PntLineDist(ai, bi, ci, Q);
|
|
beta = w.x * ai + w.y * bi;
|
|
x = VecLenSq(w) - beta * beta;
|
|
y = 2.0 * (VecDotProd(u, w) - alpha * beta);
|
|
z = VecLenSq(u) - alpha * alpha;
|
|
Roots2abc(x, y, z, roots, num_sol);
|
|
|
|
if (num_sol == 2) {
|
|
centers[0] = RayPnt(Q, w, roots[0]);
|
|
centers[1] = RayPnt(Q, w, roots[1]);
|
|
|
|
#ifdef TRACE
|
|
if (debug) {
|
|
printf("SegSegPntCntr - u = (%20.16f, %20.16f)\n", u.x, u.y);
|
|
printf("SegSegPntCntr - v = (%20.16f, %20.16f)\n", v.x, v.y);
|
|
}
|
|
#endif
|
|
}
|
|
else if (num_sol == 1) {
|
|
centers[0] = RayPnt(Q, w, roots[0]);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
restart = false;
|
|
|
|
/* */
|
|
/* 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. */
|
|
/* */
|
|
for (m = 0; m < num_sol; ++m) {
|
|
radii[m] = PntPntDist(centers[m], p);
|
|
}
|
|
best_sol = ClassifySolutions(i_in, j_in, k, e, SEG, SEG, PNT,
|
|
false, true, true, centers, radii,
|
|
&num_sol, eps, problematic);
|
|
assert((0 <= best_sol) && (best_sol < num_sol));
|
|
if (!*problematic) {
|
|
*cntr = centers[best_sol];
|
|
*r2 = radii[best_sol];
|
|
return true;
|
|
}
|
|
else {
|
|
eps *= 10.0;
|
|
}
|
|
}
|
|
|
|
if (eq(lgth_i, ZERO_MAX)) {
|
|
/* */
|
|
/* this is a seg with a very small length; we replace it */
|
|
/* */
|
|
VD_Dbg_Warning("SegSegPntCntr() - seg with small length!");
|
|
ReplaceSeg(i);
|
|
restart = true;
|
|
return false;
|
|
}
|
|
else if (eq(lgth_j, ZERO_MAX)) {
|
|
/* */
|
|
/* this is a seg with a very small length; we replace it */
|
|
/* */
|
|
VD_Dbg_Warning("SegSegPntCntr() - seg with small length!");
|
|
ReplaceSeg(j);
|
|
restart = true;
|
|
return false;
|
|
}
|
|
else if (CheckIntersectionsLocally(i, SEG, j, SEG, k, PNT)) {
|
|
restart = true;
|
|
return false;
|
|
}
|
|
|
|
#ifdef VRONI_DBG_WARN
|
|
printf("\nCenter not computed for seg-seg-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);
|
|
q = GetSegStartCoord(j);
|
|
v = GetSegEndCoord(j);
|
|
printf("%20.16f %20.16f %20.16f %20.16f\n", q.x, q.y, v.x, v.y);
|
|
p = GetPntCoords(k);
|
|
printf("%20.16f %20.16f\n", p.x, p.y);
|
|
#endif
|
|
|
|
*cntr = centers[best_sol];
|
|
*r2 = radii[best_sol];
|
|
restart = false;
|
|
|
|
return false;
|
|
}
|