Files
SaraP 739088af9f Vroni 7.8 :
- aggiornamento versione.
2025-01-29 16:24:30 +01:00

776 lines
25 KiB
C++

/*****************************************************************************/
/* */
/* Copyright (C) 2001-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 */
/* */
/*****************************************************************************/
/* #define TRACE */
#include "vroni_object.h"
void vroniObject::InsertBisectorClose(int i, t_site t, int e, int k,
int e1r, int e2l, int n1r, int n2l)
{
int i0, e0;
t_site t0;
/* */
/* insert a new bisector and close a Voronoi cell. */
/* */
if (GetEndNode(e1r) == n1r) {
GetRgtSiteData(e1r, &i0, &t0);
} else {
assert(GetStartNode(e1r) == n1r);
GetLftSiteData(e1r, &i0, &t0);
}
e0 = StoreEdge(n1r, n2l, i, i0, t, t0);
#ifdef TRACE
if (dbg_center) {
printf("new edge %2d between sites (%2d/0, %2d/%1d) ",
e0, i, i0, t0);
printf("computed\n");
printf("edge %d - start node %d, end node %d\n", e0, n1r, n2l);
}
#endif
CloseVoronoiCell(e1r, e0, e2l, n1r, n2l);
/* */
/* make sure the VD edge pointer is up-to-date */
/* */
SetVDPtr(i0, t0, e0);
return;
}
void vroniObject::InsertReplacement(int i, t_site t, int n, int e, int k,
int* el, int* er, int* nl, int* nr)
{
t_site t1, t2;
int i1, i2;
int n0, site = NIL;
coord q;
double r;
vr_bool problematic;
coord centers[4];
double radii[4];
int num_sol = 0, best_sol = NIL;
/* */
/* the Voronoi node nodes[k] is to be kept. insert a replacement */
/* node for nodes[n] (which will be deleted). */
/* */
if (!IsNodeDummy(k)) MarkNodeUnchecked(k);
#ifdef TRACE
if (dbg_center) {
printf("new center for edge %3d and site %3d of type %d\n", e, i, t);
GetLftSiteData(e, &i1, &t1);
GetRgtSiteData(e, &i2, &t2);
printf(" defining sites of edge: %2d/%1d, %2d/%1d)\n",
i1, t1, i2, t2);
}
#endif
if (!ComputeCenter(e, i, t, &q, &r, &problematic, &site)) {
if (restart) return;
if (IsInvalidData(false))
return;
else {
VD_Dbg_Warning("UpdateRecursive() - desperate mode (center)!");
centers[0] = q;
radii[0] = r;
DesperateComputeCenter(e, i, t, &q, &r);
centers[1] = q;
radii[1] = r;
num_sol = 2;
GetLftSiteData(e, &i1, &t1);
GetRgtSiteData(e, &i2, &t2);
best_sol = ClassifySolutions(i, i2, i1, e, t, t2, t1,
false, true, true, centers, radii,
&num_sol, ZERO_MAX, &problematic);
q = centers[best_sol];
r = radii[best_sol];
}
}
/* */
/* store the new node, and locally prepare the VD for the update. */
/* */
n0 = StoreNode(q.x, q.y, r);
if (site != NIL) SetNodeSite(n0, site);
#ifdef TRACE
if (dbg_center) {
GetLftSiteData(e, &i1, &t1);
GetRgtSiteData(e, &i2, &t2);
printf("new node %2d between sites (%2d/%1d, %2d/%1d, %2d/%1d) ",
n0, i, t, i1, t1, i2, t2);
printf("computed\n");
printf(" (%f, %f), radius2 = %f\n", q.x, q.y, r);
}
#endif
UpdateNodeEdgeData(e, n, n0);
*el = *er = e;
*nl = *nr = n0;
return;
}
void vroniObject::FreeVDConstructionData(void)
{
FreeMemory((void**) &pnt_visited, "voronoi:pnt_visited");
FreeMemory((void**) &sites_visited, "voronoi:sites_visited");
FreeMemory((void**) &nodes_checked, "voronoi:nodes_checked");
FreeMemory((void**) &preserve, "voronoi:preserve");
FreeMemory((void**) &loop_stack, "voronoi:loop_stack");
num_visited = 0;
max_num_visited = 0;
num_pnt_visited = 0;
max_num_pnt_visited = 0;
num_checked = 0;
max_num_checked = 0;
num_preserve = 0;
max_num_preserve = 0;
num_loop_stack = 0;
max_num_loop_stack = 0;
return;
}
void vroniObject::ResetVDConstructionData(void)
{
num_visited = 0;
num_pnt_visited = 0;
num_preserve = 0;
num_loop_stack = 0;
num_checked = 0;
return;
}
void vroniObject::ResetNodeStatus(int e, int n)
{
int k, e1, e2;
k = GetOtherNode(e, n);
while (IsDeg2Node(k) && (IsNodeChecked(k) || IsNodeVisited(k))) {
MarkNodeUnchecked(k);
e1 = GetCCWEdge(e, k);
assert(GetCCWEdge(e1, k) == e);
e = e1;
k = GetOtherNode(e, k);
}
if (!(IsNodeChecked(k) || IsNodeVisited(k))) return;
MarkNodeUnchecked(k);
e1 = GetCCWEdge(e, k);
e2 = GetCWEdge(e, k);
ResetNodeStatus(e1, k);
ResetNodeStatus(e2, k);
return;
}
void vroniObject::ResetTreeDeleteStatus(int e, int n)
{
int k, e1, e2;
k = GetOtherNode(e, n);
while (IsDeg2Node(k) && IsNodeDeleted(k)) {
MarkNodeVisited(k);
e1 = GetCCWEdge(e, k);
assert(GetCCWEdge(e1, k) == e);
e = e1;
k = GetOtherNode(e, k);
}
if (!IsNodeDeleted(k)) return;
MarkNodeVisited(k);
e1 = GetCCWEdge(e, k);
e2 = GetCWEdge(e, k);
ResetTreeDeleteStatus(e1, k);
ResetTreeDeleteStatus(e2, k);
return;
}
/* */
/* this routine inserts the closest point between segs[i] and the point */
/* or arc i1 as a degree-2 node if it lies on their common VD edge e. */
/* */
void vroniObject::InsertApexDegreeTwoNodeSeg(int e, int i, int i1, t_site t1)
{
int n, n1, n2;
coord c1, ss, se, p1, p2, v, w, q1, q2, p, u1, u2, apex;
double t, d, r1, r2, a, b, sa, sb, sc;
assert(InEdgesList(e));
assert(t1==PNT || t1==ARC);
assert(t1==ARC || InPntsList(i1));
assert(t1==PNT || InArcsList(i1));
assert(InSegsList(i));
#ifndef GENUINE_ARCS
assert(t1!=ARC);
#endif
if ((t1 == PNT) && (IsSegStartPnt(i, i1) || IsSegEndPnt(i, i1))) {
/* */
/* pnts[i1] is the start point or end point of segs[i]. create a */
/* dummy node at pnts[i1] if the VD edge extends to both sides of */
/* segs[i]. */
/* */
if (IsSegStartPnt(i, i1)) p1 = GetSegStartCoord(i);
else p1 = GetSegEndCoord(i);
if (!HasIncidentSite(i1)) {
/* */
/* segs[i] is the first site incident at pnts[i1] that has been */
/* inserted. */
/* */
InsertSiteNode(e, p1, n, i1);
}
else {
n1 = GetStartNode(e);
n2 = GetEndNode(e);
GetNodeData(n1, &q1, &r1);
GetNodeData(n2, &q2, &r2);
if (gt(r1, ZERO) && gt(r2, ZERO)) {
p = GetPntCoords(i1);
u1 = VecSub(q1, p1);
u2 = VecSub(q2, p1);
if (VecDotProd(u1, u2) < 0.0) {
InsertSiteNode(e, p1, n, i1);
}
}
}
}
else {
/* */
/* get seg and pnt/arc data */
/* */
GetSegEqnData(i, &sa, &sb, &sc);
ss = GetSegStartCoord(i);
se = GetSegEndCoord(i);
c1 = (t1 == PNT) ? GetPntCoords(i1) : GetArcCenter(i1);
/* */
/* project c1 onto seg */
/* */
v = VecSub(se, ss);
d = VecLen(v);
v = VecDiv(d, v);
w = VecSub(c1, ss);
t = VecDotProd(v, w);
if ((t <= 0.0) || (t >= d)) return; /* projection not in CoI */
p2 = RayPnt(ss, v, t);
if (t1 == PNT) {
p1 = c1;
}
else {
r1 = GetArcRadius(i1);
GetChordEqnData(i1, &a, &b, &d);
w = MakeVec(sa, sb);
p1 = RayPnt(c1, w, r1);
if (!IsOnArcStrict(a, b, d, p1)) {
p1 = RayPnt(c1, w, -r1);
if (!IsOnArcStrict(a, b, d, p1)) return;
}
}
apex = MidPoint(p1, p2);
if (IsBetweenVoronoiNodes(e, apex, ZERO)) {
/* */
/* we've found an apex! */
/* */
d = PntPntDist(p1, apex);
n = StoreNode(apex.x, apex.y, d);
InsertDummyNode(e, n, false);
}
}
return;
}
/* */
/* this routine inserts the closest point between arcs[i] and the point/arc*/
/* pnts/arcs[i1] as a degree-2 node if it lies on their common VD edge e. */
/* t1 determines whether i1 is pnt or arc. */
/* */
#ifdef GENUINE_ARCS
void vroniObject::InsertApexDegreeTwoNodeArc(int e, int i, int i1, t_site t1)
{
int n, n1, n2;
coord w, c1, v, c, p1, p2, q1, q2, p, u1, u2;
double d, r, r1, r2;
double a, b;
coord apex;
assert(InEdgesList(e));
assert(InArcsList(i));
assert(t1==PNT || t1==ARC);
assert(t1==PNT || InArcsList(i1));
assert(t1==ARC || InPntsList(i1));
if ((t1 == PNT) && (IsArcStartPnt(i, i1) || IsArcEndPnt(i, i1))) {
/* */
/* pnts[i1] is the start point or end point of segs[i]. create a */
/* dummy node at pnts[i1] if the VD edge extends to both sides of */
/* segs[i]. */
/* */
p = GetPntCoords(i1);
if (!HasIncidentSite(i1)) {
/* */
/* segs[i] is the first site incident at pnts[i1] that has been */
/* inserted. */
/* */
InsertSiteNode(e, p, n, i1);
}
else {
n1 = GetStartNode(e);
n2 = GetEndNode(e);
GetNodeData(n1, &q1, &r1);
GetNodeData(n2, &q2, &r2);
if (gt(r1, ZERO) && gt(r2, ZERO)) {
u1 = VecSub(q1, p);
u2 = VecSub(q2, p);
if (VecDotProd(u1, u2) < 0.0) {
/* */
/* this is a VD edge that extends to both sides of the arc i */
/* */
InsertSiteNode(e, p, n, i1);
}
}
}
}
else {
/* */
/* we check whether i and i1 contain a pair of points p1, p2 */
/* whose midpoint defines an apex of the bisector e. */
/* */
c = GetArcCenter(i);
r = GetArcRadius(i);
c1 = (t1 == PNT) ? GetPntCoords(i1) : GetArcCenter(i1);
v = VecSub(c1, c);
d = VecLen(v);
if (le(d, ZERO)) {
/* */
/* the point is at the center of the arc, or the arcs are */
/* concentric. in either case no apex exists. */
/* */
return;
}
v = VecDiv(d, v);
GetChordEqnData(i, &a, &b, &d);
p2 = RayPnt(c, v, r);
if (!IsOnArcStrict(a, b, d, p2)) {
p2 = RayPnt(c, v, -r);
if (!IsOnArcStrict(a, b, d, p2)) return;
}
if (t1 == PNT) {
p1 = c1;
}
else {
r1 = GetArcRadius(i1);
GetChordEqnData(i1, &a, &b, &d);
p1 = RayPnt(c1, v, r1);
if (!IsOnArcStrict(a, b, d, p1)) {
p1 = RayPnt(c1, v, -r1);
if (!IsOnArcStrict(a, b, d, p1)) return;
}
}
apex = MidPoint(p1, p2);
if (IsBetweenVoronoiNodes(e, apex, ZERO)) {
/* */
/* we've found an apex! */
/* */
d = PntPntDist(p1, apex);
n = StoreNode(apex.x, apex.y, d);
InsertDummyNode(e, n, false);
}
}
return;
}
#endif
/* */
/* the loop stack contains edges leading from n1 to n3, plus the actual */
/* loop. we'll delete the prefix from n1 to n3. */
/* */
void vroniObject::PurifyLoop(int n1, int n3)
{
int number, k, j, e = 0, k0;
GetLoopStackSize(number);
for (j = 0; j < number; ++j) {
GetLoopStackElement(j, e);
k = GetStartNode(e);
k0 = GetEndNode(e);
if ((k == n3) || (k0 == n3)) {
if (n1 != n3) ++j;
ResetLoopStackMin(j);
return;
}
}
assert(1 == 0);
return;
}
/* */
/* check whether a loop of VD nodes marked for deletion exists. we start at */
/* VD edge e at node n and visit neighboring nodes that have also been */
/* flagged as "visited" in the previous step of the insertion of a new line */
/* segment. (those nodes will be deleted afterwards.) during the scan we */
/* convert the status of those nodes from "visited" to "deleted". if we */
/* encounter during this (recursive) scan a node that carries already the */
/* status "deleted" then we know that a loop exists. all edges traversed */
/* are stored in a stack, which is reset whenever it is clear that some set */
/* of edges traversed so far does not form a loop. */
/* */
vr_bool vroniObject::LoopExists(int e, int n, int *m)
{
int k, e1, e2, size;
k = GetOtherNode(e, n);
GetLoopStackSize(size);
PushOntoLoopStack(e);
while (IsDeg2Node(k) && IsNodeVisited(k)) {
MarkNodeDeleted(k);
e1 = GetCCWEdge(e, k);
assert(GetCCWEdge(e1, k) == e);
e = e1;
k = GetOtherNode(e, k);
PushOntoLoopStack(e);
}
if (IsNodeDeleted(k)) {
*m = k;
return true;
}
else if (!IsNodeVisited(k)) {
ResetLoopStackSize(size);
return false;
}
MarkNodeDeleted(k);
e1 = GetCCWEdge(e, k);
if (LoopExists(e1, k, m))
return true;
e2 = GetCWEdge(e, k);
if (LoopExists(e2, k, m))
return true;
else {
ResetLoopStackSize(size);
return false;
}
}
vr_bool vroniObject::ComputeCenter(int e, int i3, t_site t3,
coord *w, double *r2,
vr_bool *problematic, int *site)
{
int i1, i2;
t_site t1, t2;
vr_bool computed = false;
coord centers[VRONI_MAXSOL];
double radii[VRONI_MAXSOL];
int num_sol = 0, best_sol = NIL;
*problematic = false;
*site = NIL;
GetLftSiteData(e, &i1, &t1);
GetRgtSiteData(e, &i2, &t2);
#ifdef TRACE
if ((i3 == 13) && (t3 == ARC)) {
printf("ComputeCenter(): e = %d\n", e);
printf("ComputeCenter(): i1 = %d, t1 = %d\n", i1, t1);
printf("ComputeCenter(): i2 = %d, t2 = %d\n", i2, t2);
printf("ComputeCenter(): i3 = %d, t3 = %d\n", i3, t3);
}
#endif
if (t1 == PNT) {
if (t2 == PNT) {
if (t3 == PNT) {
computed = PntPntPntCntr(i1, i2, i3, w, r2);
}
else if (t3 == SEG) {
computed = SegPntPntCntr(i3, i1, i2, e, w, r2, problematic);
}
#ifdef GENUINE_ARCS
else if (t3 == ARC) {
computed = ArcPntPntCntr(i3, i1, i2, e, w, r2, problematic);
}
#endif
else {
VD_Warning("ComputeCenter() - site type not available");
}
}
else if (t2 == SEG) {
if (t3 == PNT) {
computed = SegPntPntCntr(i2, i1, i3, e, w, r2, problematic);
}
else if (t3 == SEG) {
computed = SegSegPntCntr(i3, i2, i1, e, w, r2, problematic, site);
}
#ifdef GENUINE_ARCS
else if (t3 == ARC) {
computed = ArcSegPntCntr(i3, i2, i1, e, w, r2, problematic, site);
}
#endif
else {
VD_Warning("ComputeCenter() - site type not available");
}
}
#ifdef GENUINE_ARCS
else if (t2 == ARC) {
if (t3 == ARC) {
computed = ArcArcPntCntr(i3, i2, i1, e, w, r2, problematic, site);
}
else {
VD_Warning("ComputeCenter() - site type not available");
}
}
#endif
else {
VD_Warning("ComputeCenter() - site type not available");
}
}
else if (t1 == SEG) {
if (t2 == PNT) {
if (t3 == SEG) {
computed = SegSegPntCntr(i3, i1, i2, e, w, r2, problematic, site);
}
#ifdef GENUINE_ARCS
else if (t3 == ARC) {
computed = ArcSegPntCntr(i3, i1, i2, e, w, r2, problematic, site);
}
#endif
else {
VD_Warning("ComputeCenter() - site type not available");
}
}
else if (t2 == SEG) {
if (t3 == SEG) {
computed = SegSegSegCntr(i1, i2, i3, e, w, r2, problematic, site);
}
#ifdef GENUINE_ARCS
else if (t3 == ARC) {
computed = ArcSegSegCntr(i3, i2, i1, e, w, r2, problematic, site);
}
#endif
else {
VD_Warning("ComputeCenter() - site type not available");
}
}
#ifdef GENUINE_ARCS
else if (t2 == ARC) {
if (t3 == ARC) {
computed = ArcArcSegCntr(i3, i2, i1, e, w, r2, problematic, site);
}
else {
VD_Warning("ComputeCenter() - site type not available");
}
}
#endif
else {
VD_Warning("ComputeCenter() - site type not available");
}
}
#ifdef GENUINE_ARCS
else if (t1 == ARC) {
if (t2 == PNT) {
if (t3 == ARC) {
computed = ArcArcPntCntr(i3, i1, i2, e, w, r2, problematic, site);
}
else {
VD_Warning("ComputeCenter() - site type not available");
}
}
else if (t2 == SEG) {
if (t3 == ARC) {
computed = ArcArcSegCntr(i3, i1, i2, e, w, r2, problematic, site);
}
else {
VD_Warning("ComputeCenter() - site type not available");
}
}
else if (t2 == ARC) {
if (t3 == ARC) {
computed = ArcArcArcCntr(i3, i2, i1, e, w, r2, problematic, site);
}
else {
VD_Warning("ComputeCenter() - site type not available");
}
}
else {
VD_Warning("ComputeCenter() - site type not available");
}
}
#endif
else {
VD_Warning("ComputeCenter() - site type not available");
}
if (computed) {
if (*problematic) {
#ifdef VRONI_DBG_WARN
printf("Using DesperateComputeCenter for e = %d and i3 = %d\n",
e, i3);
printf("i1 = %d of type %d, i2 = %d of type %d\n", i1, t1, i2, t2);
#endif
centers[0] = *w;
radii[0] = *r2;
DesperateComputeCenter(e, i3, t3, w, r2);
centers[1] = *w;
radii[1] = *r2;
num_sol = 2;
best_sol = ClassifySolutions(i3, i2, i1, e, t3, t2, t1,
false, true, true, centers, radii,
&num_sol, ZERO_MAX, problematic);
*w = centers[best_sol];
*r2 = radii[best_sol];
if (best_sol == 1) VD_Dbg_Warning("ComputeCenter() - DesperateComputeCenter() improved solution!\n");
}
return true;
}
return false;
}
void vroniObject::UpdateRecursive(int i, t_site t, int n, int e,
int *el, int *er, int *nl, int *nr)
{
#ifdef TRACE
vr_bool dbg_center = false;
#endif
int k, e1, e2;
int e1r, n1r, e2l, n2l;
coord q;
assert((t == SEG) || (t == ARC));
k = GetOtherNode(e, n);
#ifdef TRACE
if (i == 2) dbg_center = true;
else dbg_center = false;
if (dbg_center) {
printf("UpdateRecursive - edge e=%d, node n=%d\n", e, n);
printf("UpdateRecursive - checking node k=%d of degree %d",
k, 3 - IsDeg2Node(k));
if (IsNodeDeleted(k)) printf(" --- to be deleted\n");
else printf(" --- to be kept\n");
}
#endif
if (IsNodeDeleted(k)) {
/* */
/* the Voronoi edge edges[e] is to be deleted; proceed recursively */
/* */
if (IsDeg2Node(k)) {
e1 = GetCCWEdge(e, k);
UpdateRecursive(i, t, k, e1, el, er, nl, nr);
if (restart) return;
}
else {
e1 = GetCCWEdge(e, k);
e2 = GetCWEdge(e, k);
UpdateRecursive(i, t, k, e1, el, &e1r, nl, &n1r);
if (restart) return;
UpdateRecursive(i, t, k, e2, &e2l, er, &n2l, nr);
if (restart) return;
/* */
/* insert a new bisector and close a Voronoi cell. */
/* */
InsertBisectorClose(i, t, e, k, e1r, e2l, n1r, n2l);
}
/* */
/* delete edges[e] and nodes[k]. */
/* */
DeleteEdge(e);
DeleteNode(k);
#ifdef TRACE
if (dbg_center) {
printf("node %2d has been deleted!\n", k);
printf("edge %2d has been deleted!\n", e);
}
#endif
}
else {
/* */
/* the Voronoi node nodes[k] is to be kept. insert a replacement */
/* node for nodes[n] (which will be deleted). */
/* */
InsertReplacement(i, t, n, e, k, el, er, nl, nr);
}
return;
}