/*****************************************************************************/ /* */ /* 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 */ /* */ /*****************************************************************************/ /* */ /* get my header files */ /* */ #include "vroni_object.h" vr_bool vroniObject::IsSameSide(int i, t_site t, int j, int i1, coord p_node) { coord p, q1, q2, q3; double r, det1, det2, side1, side2; double a, b, c; //printf("checking sidedness relative to site %d/%d\n", i, t); p = GetArcCenter(j); if (IsArcStartPnt(j, i1)) { q1 = GetArcStartCoord(j); q2 = GetArcEndCoord(j); } else { q1 = GetArcEndCoord(j); q2 = GetArcStartCoord(j); } if (t == SEG) { if (IsSegStartPnt(i, i1)) q3 = GetSegEndCoord(i); else q3 = GetSegStartCoord(i); } else { assert(t == ARC); if (IsArcStartPnt(i, i1)) q3 = GetArcEndCoord(i); else q3 = GetArcStartCoord(i); } det1 = Det2D(p, q1, q2); det2 = Det2D(p, q1, q3); //printf("det1 = %20.16f, det2 = %20.16f\n", det1, det2); if (SignEps(det1, ZERO) * SignEps(det2, ZERO) < 0) return true; if (t == SEG) { GetSegEqnData(i, &a, &b, &c); side1 = PntLineDist(a, b, c, q2); side2 = PntLineDist(a, b, c, p_node); } else { p = GetArcCenter(i); r = GetArcRadius(i); side1 = PntCircleDist(p, r, q2); side2 = PntCircleDist(p, r, p_node); } //printf("side1 = %20.16f, side2 = %20.16f\n", side1, side2); if (SignEps(side1, ZERO) * SignEps(side2, ZERO) < 0) return false; else return true; } /** * Recursively search for the most violated point in the VD-cell of * point i1 when inserting arc j (where i1 is an endpoint). The * VD-Node k0 is the start of rec. search, k is the current one * and *n and *dist hold the best candidate. * * Tactics: * - If we reached k0 again, return (n,dist) = (NIL,DBL_MAX) * - If we are still at a node associated to i1, then go on recursively CW and CCW * - Otherwise calc *n and *dist. However, beware of special case of e being on * boundary of CI of arc j. */ void vroniObject::FindMostViolatedRecursiveArc(int e, int k0, int i1, int j, int k, int *n, double *dist) { // MODIF questa funzione non dovrebbe mai scorrere tutti i nodi, quindi se si entra nella ricorsione troppe volte // viene lanciata eccezione m_nMostViolatedRecursiveArcCount ++ ; if ( m_nMostViolatedRecursiveArcCount > 2 * num_nodes) throw std::runtime_error("VRONI error: ComputeVD() - FindMostViolatedRecursveArc infinite loop"); coord p, q1, ns, ne, pend, in1center; double dist2, radius1, r; int e1, e2, k1, n2, in1; t_site int1; assert(InEdgesList(e)); assert(InPntsList(i1)); assert(InNodesList(k0)); assert(InNodesList(k)); assert(InArcsList(j)); assert( IsArcStartPnt(j, i1) || IsArcEndPnt(j, i1) ); k1 = GetOtherNode(e, k); assert(InNodesList(k1)); //printf("FMRV-arc: e=%d, k0=%d, i1=%d, j=%d, k=%d, k1=%d\n", // e, k0, i1, j, k, k1); if (k1 == k0) { /* */ /* returned to start or reached a dummy node: stop */ /* */ *dist = DBL_MAX; *n = NIL; return; } q1 = GetNodeCoord(k1); if ( PntPntDist(q1, GetPntCoords(i1)) < ZERO ) { /* */ /* node k1 coincides with point i1: keep on searching */ /* */ e1 = GetCCWEdge(e, k1); e2 = GetCWEdge(e, k1); FindMostViolatedRecursiveArc(e1, k0, i1, j, k1, n, dist); if( e1 != e2) { /* */ /* follow also second path away from point i1 */ /* */ FindMostViolatedRecursiveArc(e2, k0, i1, j, k1, &n2, &dist2); if (dist2 < *dist) { *n = n2; *dist = dist2; } } } else { /* */ /* the Voronoi edge e from k to k1 points away from point i1 */ /* */ GetNodeData(k1, &q1, &radius1); ns = GetArcStartNormal(j); ne = GetArcEndNormal(j); p = GetArcCenter(j); r = GetArcRadius(j); //Get distance to arc *dist = AbsPntArcDist(ns, ne, VecSub(q1,p), r, TINY) - radius1; *n = k1; //printf(" FMRV-arc: node %d has *dist=%e\n", *n, *dist); if (eq(*dist, zero)) { /* */ /* special case: k1 lies on common boundary of CI of two sites */ /* that share the same tangent line */ /* */ //printf(" FMVR-arc: Consider node %d, which is on CI-boundary (*dist=%e)\n", k1, *dist); pend = IsArcStartPnt(j,i1) ? GetArcEndCoord(j) : GetArcStartCoord(j); if (IsLftSite(e, i1, PNT)) GetRgtSiteData(e, &in1, &int1); else GetLftSiteData(e, &in1, &int1); assert( int1==SEG || int1==ARC); //printf(" other site 1: %d, %d\n", in1, int1); if (!IsSameSide(in1, int1, j, i1, q1)) { *dist = DBL_MAX; *n = NIL; return; } if (IsLftSite(e, in1, int1)) GetRgtSiteData(e, &in1, &int1); else GetLftSiteData(e, &in1, &int1); /* */ /* This site could again be the point, but maybe another edge e2 */ /* points back to e1 and its opposite site has to be checked. */ /* */ if (int1 == PNT ) { e2 = GetCCWEdge(e, k1); if (!IsLftSite(e2, i1, PNT) && !IsRgtSite(e2, i1, PNT)) /* */ /* CCW-edge does not belong to point i1 --> set e2 to CW-edge */ /* */ e2 = GetCWEdge(e, k1); if (IsLftSite(e2, i1, PNT) || IsRgtSite(e2, i1, PNT)) { /* */ /* edge e2 belongs to point i1; select the opposite site */ /* */ if (IsLftSite(e2, i1, PNT) ) GetRgtSiteData(e, &in1, &int1); else GetLftSiteData(e, &in1, &int1); } } if (int1 != PNT ) { //printf(" other site 2: %d, %d\n", in1, int1); if (!IsSameSide(in1, int1, j, i1, q1)) { *dist = DBL_MAX; *n = NIL; return; } } } } return; } /* */ /* find the Voronoi node on the Voronoi boundary of the point pnts[i1] */ /* whose Voronoi disk is violated the most by the insertion of arcs[j]. */ /* note that pnts[i1] and pnts[i2] are the endpoints of arcs[j]. if */ /* the Voronoi cells of pnts[i1] and pnts[i2] share nodes then we take */ /* one of those nodes, and set node_shared = true. */ /* */ int vroniObject::FindMostViolatedNodeArc(int i1, int i2, int j, vr_bool *node_shared) { int e, e1, e0, k, k_min = NIL, n, n1; coord p, q, ns, ne; double r, dist, dist1, dist_min, rad; vr_bool shared, shared1; #ifdef TRACE vr_bool trace_mvna = true; #endif assert(InPntsList(i1)); assert(InPntsList(i2)); assert(InArcsList(j)); assert(IsArcStartPnt(j, i1) || IsArcEndPnt(j, i1)); assert(IsArcStartPnt(j, i2) || IsArcEndPnt(j, i2)); *node_shared = false; /* */ /* we'll scan the Voronoi region of the point in order to find the */ /* node whose Voronoi disk is most violated by the insertion of the */ /* arc arcs[j]. */ /* */ e = GetVDPtr(i1, PNT); assert(InEdgesList(e)); assert(IsLftSite(e, i1, PNT) || IsRgtSite(e, i1, PNT)); ns = GetArcStartNormal(j); ne = GetArcEndNormal(j); p = GetArcCenter(j); r = GetArcRadius(j); #ifdef TRACE if (trace_mvna) { printf("FindMostViolatedNodeArc() -\n"); printf("center p: (%20.16f, %20.16f)\n", p.x, p.y); printf("radius r: %20.16f\n", r); printf("ns: (%20.16f, %20.16f)\n", ns.x, ns.y); printf("ne: (%20.16f, %20.16f)\n", ne.x, ne.y); } #endif if (IsRgtSite(e, i1, PNT)) k = GetEndNode(e); else k = GetStartNode(e); assert(InNodesList(k)); e0 = e; dist_min = LARGE; *node_shared = false; shared = IsLftRgtSite(e, i1, PNT) && IsLftRgtSite(e, i2, PNT); #ifdef TRACE if (trace_mvna) printf("i1 = %d, i2 = %d, j = %d\n", i1, i2, j); #endif do { e1 = GetCCWEdge(e, k); shared1 = IsLftRgtSite(e1, i1, PNT) && IsLftRgtSite(e1, i2, PNT); shared = shared || shared1; #ifdef TRACE if (trace_mvna) PrintNodeData(k, "checked"); #endif /* */ /* check whether q is in the cone of influence of arcs[j]. */ /* */ if (GetNodeSite(k) != i1) { GetNodeData(k, &q, &rad); q = VecSub(q, p); dist = AbsPntArcDist(ns, ne, q, r, TINY) - rad; if (*node_shared && shared && (dist < dist_min)) { #ifdef TRACE if (trace_mvna) { printf("node %d accepted (shared)\n", k); printf("old distance: %f, new distance: %f\n", dist_min, dist); PrintNodeData(k, "FindMostViolatedArc"); } #endif dist_min = dist; k_min = k; } else if (dist < dist_min) { #ifdef TRACE if (trace_mvna) { if (shared) printf("node %d accepted (shared)\n", k); else printf("node %d accepted\n", k); printf("old distance: %20.16f\nnew distance: %20.16f\n", dist_min, dist); PrintNodeData(k, "FindMostViolatedArc"); } #endif dist_min = dist; k_min = k; if (shared) *node_shared = true; } } e = e1; k = GetOtherNode(e, k); shared = shared1; } while (e != e0); if (((k_min == NIL) || (ge(dist_min, ZERO))) && HasIncidentSite(i1)) { /* */ /* at least one seg or arc is already incident at this point. */ /* */ *node_shared = false; k = GetPntNode(i1); assert(InNodesList(k)); e = GetIncidentEdge(k); assert(InEdgesList(e)); //Restore data for calls coming now... k = GetPntNode(i1); e = GetIncidentEdge(k); e1 = GetCCWEdge(e, k); m_nMostViolatedRecursiveArcCount = 0 ; // MODIF azzero il contatore per la nuova chiamata FindMostViolatedRecursiveArc(e, k, i1, j, k, &n, &dist); m_nMostViolatedRecursiveArcCount = 0 ; // MODIF azzero il contatore per la nuova chiamata FindMostViolatedRecursiveArc(e1, k, i1, j, k, &n1, &dist1); if (dist < dist1) k_min = n; else k_min = n1; #ifdef TRACE if (trace_mvna) { printf("node %d accepted\n", k_min); printf("old distance: %20.16f, new distance: %20.16f\n", dist_min, VroniMin(dist, dist1)); PrintNodeData(k_min, "FindMostViolatedArc"); } #endif } return k_min; } /* */ /* computes a unit normal vector w relative to point p and the site */ /* (seg, arc) that points away from the site. it also computes the distance */ /* of p from the supporting element (line, circle) of the site. */ /* */ void vroniObject::GetPntSiteNormal(int i, t_site site, coord p, coord* v, double* dist) { double len; if (site == SEG) { double a, b, c; GetSegEqnData(i, &a, &b, &c); *dist = a*p.x + b*p.y + c ; if (*dist < 0.0) { v->x = -a; v->y = -b; } else { v->x = a; v->y = b; } } else if (site == ARC) { coord c = GetArcCenter(i); double r = GetArcRadius(i); *v = VecSub(p, c); len = VecLen(*v); if (len > zero) { /* */ /* point is not in the center of the arc */ /* */ if (len < r) { *v = VecDiv(-len, *v); *dist = r - len; } else { *v = VecDiv(len, *v); *dist = len - r; } } else { /* */ /* point is in the center of the arc */ /* */ coord s, e; s = GetArcStartCoord(i); e = GetArcEndCoord(i); s = VecSub(c, s); e = VecSub(c, e); //Add start and end vector *v = VecAdd(s, e); len = VecLen(*v); *v = VecDiv(len, *v); *dist = Abs(r - len); } } else assert(0); return; } void vroniObject::ProjectPntOnSite( int i, t_site site, coord p, coord* w) { coord v; double len; //Projection on point is always the point itself if( site==PNT ) { *w = GetPntCoords(i); } //Projecting on segment else if( site==SEG ) { double a, b, c; GetSegEqnData(i, &a, &b, &c); v.x = a; v.y = b; if( IsInSegConeZero( GetSegStartCoord(i), p, a, b, GetSegLgth(i), zero) == 0 ) { const double dist = a*p.x + b*p.y + c ; v = VecMult(-dist, v); *w = VecAdd(v, p); assert( Abs(w->x*a + w->y*b + c) < ZERO_MAX ); } else { w->x = LARGE; w->y = LARGE; } } else if( site==ARC ) { if( IsPntInArcConeEps(i, p, zero) ) { coord c = GetArcCenter(i); v = VecSub( p, c); len = VecLen(v); //Point not in center if( len > zero ) { len = GetArcRadius(i) / len; v = VecMult(len, v); *w = VecAdd(v, c); } //Point is in center of arc else { //Get start and end vectors coord s, e; s = GetArcStartCoord(i); e = GetArcEndCoord(i); s = VecSub(s, c); e = VecSub(e, c); //Add start and end vector and scale //it to radius-length v = VecAdd(s, e); len = VecLen(v); v = VecMult(GetArcRadius(i)/len, v); //This is our projection *w = VecAdd(v, c); } assert( IsPntInArcCone(i, *w)); } else { w->x = LARGE; w->y = LARGE; } } else assert(0); return; } /* */ /* this routine first checks whether all VD nodes on the boundary of the */ /* Voronoi cell of site i of type t have been marked for deletion due to */ /* the insertion of the circular arc i0. (as a matter of fact, t == PNT !) */ /* it then constructs a point on the boundary of the Voronoi cell of i */ /* which is guaranteed to be closer to pnt i than to arc i0. */ /* */ vr_bool vroniObject::PreserveVoronoiCellArc(int i, t_site t, int i0) { int e, n, e0, k, k0, s, s0; double dist, r0_2, r_2, r2, ap, bp, cp, d, d0; coord p, q, q0, v, w, u; #ifdef TRACE if (debug_flag_VD) { #ifdef GRAPHICS if (graphics) AddToCurBuffer(i, t, AlertColor); #endif printf("in PreserveVoronoiCellArc for site %d of type %d and arc %d\n", i, t, i0); printf("start [%4d]: (%20.16f, %20.16f)\n", arcs[i0].i1, pnts[arcs[i0].i1].p.x, pnts[arcs[i0].i1].p.y); printf("end [%4d]: (%20.16f, %20.16f)\n", arcs[i0].i2, pnts[arcs[i0].i2].p.x, pnts[arcs[i0].i2].p.y); printf("center: (%20.16f, %20.16f)\n", arcs[i0].c.x, arcs[i0].c.y); printf("pnt [%4d]: (%20.16f, %20.16f)\n", i, pnts[i].p.x, pnts[i].p.y); } #endif assert((t == PNT) && InPntsList(i)); /* */ /* we'll find a bisector on the boundary of the Voronoi cell of site i */ /* for which we can guarantee that some part of it will survive the */ /* insertion of the circular arc arcs[i0]. we will insert a dummy node */ /* on this bisector in order to prevent its removal. after the insertion */ /* of arcs[i0] we will remove the dummy node. */ /* */ /* scan the Voronoi cell in CCW order and check whether all nodes */ /* really are to be deleted. */ /* */ e = GetVDPtr(i, t); assert(InEdgesList(e)); assert(IsLftSite(e, i, t) || IsRgtSite(e, i, t)); if (IsRgtSite(e, i, t)) k = GetEndNode(e); else k = GetStartNode(e); e0 = e; do { if (!IsNodeVisited(k)) { /* */ /* hm... this Voronoi cell has undeleted nodes! */ /* */ #ifdef TRACE if (debug_flag_VD) { printf("the node %d of this Voronoi cell is undeleted\n", k); } #endif return true; } e = GetCCWEdge(e, k); k = GetOtherNode(e, k); } while (e != e0); /* */ /* indeed, all Voronoi nodes of this cell are to be deleted. */ /* */ if (IsArcStartPnt(i0, i) || IsArcEndPnt(i0, i)) { /* */ /* check whether there exists a node with radius 0.0. */ /* */ e = GetVDPtr(i, t); if (IsRgtSite(e, i, t)) k = GetEndNode(e); else k = GetStartNode(e); e0 = e; do { if (i == GetNodeSite(k)) { GetNodeData(k, &q, &r2); assert(r2 == 0.0); InsertSiteNode(e, q, n, i); MarkNodeChecked(n); #ifdef TRACE if (debug_flag_VD) { printf("dummy node with radius 0 inserted\n"); #ifdef GRAPHICS if (graphics) AddToCurBuffer(n, VDN, AlertColor); #endif } #endif return true; } e = GetCCWEdge(e, k); k = GetOtherNode(e, k); } while (e != e0); return false; } /* */ /* compute the coefficients of the equation of the line through p and */ /* the direction vector v of a ray that starts at p, points away from */ /* arcs[i0], and is perpendicular to arcs[i0]. */ /* */ p = GetPntCoords(i); GetPntSiteNormal(i0, ARC, p, &v, &dist); if (eq(dist, zero)) { if (IsPntInArcConeEps(i0, p, zero)) { SetInvalidSites(i0, ARC, i, PNT, zero); VD_Dbg_Warning("PreserveVoronoiCellArc() - point lies on arc!"); if (IsInvalidData(false)) { restart = true; return false; } } } ap = - v.y; bp = v.x; cp = - ap * p.x - bp * p.y; /* */ /* scan the boundary of the VD cell of p and find a VD edge whose */ /* end points lie on different sides of the ray through p */ /* */ e = GetVDPtr(i, t); if (IsRgtSite(e, i, t)) { k = GetEndNode(e); k0 = GetStartNode(e); } else { k = GetStartNode(e); k0 = GetEndNode(e); } GetNodeData(k0, &q0, &r0_2); d0 = PntLineDist(ap, bp, cp, q0); s0 = SignEps(d0, zero); e0 = e; do { GetNodeData(k, &q, &r_2); d = PntLineDist(ap, bp, cp, q); s = SignEps(d, ZERO); if ((s * s0 < 0) && gt(r_2, zero) && gt(r0_2, zero)) { /* */ /* the nodes q0 and q lie on different sides of the line */ /* through p. */ /* */ if (IntersectRayBisector(i, v, e, &w, &r2)) { if (IsBetweenVoronoiNodes(e, w, zero)) { n = StoreNode(w.x, w.y, r2); InsertDummyNode(e, n, true); MarkNodeDummy(n); #ifdef TRACE if (debug_flag_VD) { printf("dummy node inserted in PreserveVoronoiCellArc\n"); #ifdef GRAPHICS if (graphics) AddToCurBuffer(n, VDN, AlertColor); #endif } #endif return true; } } } s0 = s; q0 = q; r0_2 = r_2; e = GetCCWEdge(e, k); k = GetOtherNode(e, k); } while (e != e0); return false; } vr_bool vroniObject::BreakLoopArc(int i0) { int start, number, i,j, e; coord v; coord ca; double ra; int sn; int en; double rs, re; coord cs, ce; t_site lsite, rsite; int li, ri; coord pl, pr, p; t_site osite; int oi; double distl; double distr; double cap; coord centers[VRONI_MAXSOL]; double radii[VRONI_MAXSOL]; int sol; //printf("BreakLoopArc, i0=%d\n", i0); assert( InArcsList(i0) ); ca = GetArcCenter(i0); ra = GetArcRadius(i0); GetLoopStackSize(number); GetLoopStackMin(start); //No loop for us, return true. if (start >= number) return true; //Check every edge in the loop for breaking the edge up for (j = start; j < number; ++j) { //The edge e is the current edge... GetLoopStackElement(j, e); //printf("Check edge %d\n", e); //Get data of start and end node sn = GetStartNode(e); en = GetStartNode(e); GetNodeData( sn, &cs, &rs ); GetNodeData( en, &ce, &re ); //Aha, do not split those edges... if( rs LARGE/2.0 ) continue; //Get normalized direction and length of center of arc to v v = VecSub(p, ca); cap = VecLen(v); v = VecDiv(cap, v); //Get the intersections if( osite==PNT || osite==ARC) { coord c = osite==PNT ? GetPntCoords(oi) : GetArcCenter(oi); double r = osite==PNT ? 0 : GetArcRadius(oi); sol = IntersectRayCircle(p, v, c, r, centers, radii, ZERO ); } else { double a,b,c; GetSegEqnData(oi, &a, &b, &c); //The procedures in arc_common.c need the form a.x+b.y = c and not //the form a.x + b.y + c = 0 c *= -1.0; sol = IntersectRaySegment(p, v, a,b,c, centers, radii, ZERO ); } //Check those solutions for( i=0; i=r2 && IsBetweenVoronoiNodes(e,w,zero) ) { int n = StoreNode(w.x, w.y, r2); InsertDummyNode(e, n, true); MarkNodeDummy(n); #ifdef TRACE #ifdef GRAPHICS if (graphics) AddToCurBuffer(n, VDN, AlertColor); #endif #endif //printf("splitted\n"); return true; } } } //printf("not splitted\n"); return false; } void vroniObject::MarkNonRecursiveArc(int e, int n, int i, std::stack& e_stack, std::stack& n_stack, std::stack& i_stack) { int i0, i1, i2, k, e1, e2, k1, k2; t_site t1, t2; coord p, q, pq, ns, ne; double dist, rad, r; #ifdef TRACE vr_bool mark_rec = true; #endif coord qn; double radn; #ifdef TRACE if (mark_rec) { printf("MarkNonRecursiveArc edge e = %d, site i = %d\n", e, i); PrintNodeData(n, "search started here"); } #endif k = GetOtherNode(e, n); do { #ifdef TRACE if (mark_rec) { printf("MarkNonRecursiveArc - new node k = %d\n", k); PrintNodeData(k, "new"); } #endif if (!IsNodeUnchecked(k)) /* */ /* this node has already been visited; make sure to avoid a loop. */ /* */ return; if (IsDeg2Node(k)) { e1 = GetCCWEdge(e, k); assert(GetCWEdge(e, k) == e1); k1 = GetOtherNode(e1, k); i0 = GetNodeSite(k); if (IsArcStartPnt(i, i0) || IsArcEndPnt(i, i0)) { MarkNodeChecked(k); return; } GetNodeData(k, &q, &rad); ns = GetArcStartNormal(i); ne = GetArcEndNormal(i); p = GetArcCenter(i); r = GetArcRadius(i); pq = VecSub(q, p); dist = AbsPntArcDist(ns, ne, pq, r, TINY) - rad; //printf("\tdist: %e, strictly in cone: %d rad: %e\n", dist, IsPntInArcCone(i, q), rad); //printf("\tp=(%e, %e), %e), q=(%e, %e)\n", p.x, p.y, r, q.x, q.y); //We want, that the point is in CI if ( le(dist, zero) ) { MarkNodeVisited(k); #ifdef TRACE #ifdef GRAPHICS if (graphics) AddToCurBuffer(k, VDN, CurrColor); #endif if (mark_rec) { printf("node %2d to be deleted\n", k); } #endif /* */ /* make sure that we don't delete an entire Voronoi cell. */ /* */ #ifdef TRACE if (mark_rec) { PrintNodeData(k1, "might be checked next"); PrintEdgeData(e1, "incident edge"); } #endif if (IsNodeVisited(k1)) { GetLftSiteData(e1, &i1, &t1); GetRgtSiteData(e1, &i2, &t2); if (t1 == PNT) if (!HasIncidentSite(i1)) AddToPreserveBuffer(i1, t1); if (t2 == PNT) if (!HasIncidentSite(i2)) AddToPreserveBuffer(i2, t2); return; } else { if (IsNodeUnchecked(k1)) { e = e1; n = k; k = k1; #ifdef TRACE if (mark_rec) { printf("MarkNonRecursiveArc edge e = %d, ", e); printf("node n = %d of status %d, site i = %d\n", n, GetNodeStatus(n), i); printf("MarkNonRecursiveArc new node k =%d\n", k); PrintNodeData(k, "new"); } #endif } else { return; } } } else { MarkNodeChecked(k); #ifdef TRACE if (mark_rec) PrintNodeData(k, "checked"); #endif return; } } } while (IsDeg2Node(k)); #ifdef TRACE if (mark_rec) printf("moved on to deg-3 node k = %d\n", k); #endif i0 = GetNodeSite(k); if (IsArcStartPnt(i, i0) || IsArcEndPnt(i, i0)) { MarkNodeChecked(k); return; } GetNodeData(k, &q, &rad); ns = GetArcStartNormal(i); ne = GetArcEndNormal(i); p = GetArcCenter(i); r = GetArcRadius(i); pq = VecSub(q, p); dist = AbsPntArcDist(ns, ne, pq, r, TINY) - rad; #ifdef TRACE if (mark_rec) { printf("\tdist: %e, in cone: %d rad: %e\n", dist, IsPntInArcCone(i, q), rad); printf("\tp=(%e, %e), %e), q=(%e, %e)\n", p.x, p.y, r, q.x, q.y); } #endif GetNodeData(n, &qn, &radn); //We want that the point is in CI //Please note: When we calc VD of 4 arcs building a circle //then there are several nodes in the center. We do not want to //delete them all, especiall we stop marking, if we would delete //a edge with length zero. //printf(" dist = %d, dist=%e\n", le(dist, zero), dist ); if( le(dist, zero) ) { MarkNodeVisited(k); #ifdef TRACE #ifdef GRAPHICS if (graphics) AddToCurBuffer(k, VDN, CurrColor); #endif if (mark_rec) { printf("node %2d to be deleted\n", k); PrintNodeData(k, "deleted"); } #endif e1 = GetCCWEdge(e, k); e2 = GetCWEdge(e, k); k1 = GetEndNode(e1); if (k1 == k) { GetRgtSiteData(e1, &i1, &t1); } else { assert(k == GetStartNode(e1)); GetLftSiteData(e1, &i1, &t1); } assert((IsLftSite(e2, i1, t1)) || (IsRgtSite(e2, i1, t1))); k1 = GetOtherNode(e1, k); k2 = GetOtherNode(e2, k); /* */ /* make sure that we don't delete an entire Voronoi cell. */ /* */ if (IsNodeVisited(k1)) { GetOtherSiteData(e1, i1, t1, &i2, &t2); if (t2 == PNT) if (!HasIncidentSite(i2)) AddToPreserveBuffer(i2, t2); if (IsNodeVisited(k2)) { if (t1 == PNT) if (!HasIncidentSite(i1)) AddToPreserveBuffer(i1, t1); GetOtherSiteData(e2, i1, t1, &i2, &t2); if (t2 == PNT) if(!HasIncidentSite(i2)) AddToPreserveBuffer(i2, t2); } } else if (IsNodeVisited(k2)) { GetOtherSiteData(e2, i1, t1, &i2, &t2); if (t2 == PNT) if (!HasIncidentSite(i2)) AddToPreserveBuffer(i2, t2); } #ifdef TRACE if (mark_rec) { PrintNodeData(k1, "might be checked next"); PrintEdgeData(e1, "incident edge"); PrintNodeData(k2, "might be checked next"); PrintEdgeData(e2, "incident edge"); } #endif if (IsNodeUnchecked(k2)) { e_stack.push(e2); n_stack.push(k); i_stack.push(i); #ifdef TRACE ++stack_size; #endif } if (IsNodeUnchecked(k1)) { e_stack.push(e1); n_stack.push(k); i_stack.push(i); #ifdef TRACE ++stack_size; #endif } #ifdef TRACE if (stack_size > max_stack_size) max_stack_size = stack_size; #endif } else { MarkNodeChecked(k); #ifdef TRACE if (mark_rec) PrintNodeData(k, "checked"); #endif } return; } void vroniObject::MarkRecursiveArc(int e, int n, int i) { /* */ /* since deep recursions turned out to be a problem on some platforms, we */ /* resort to a stack to avoid explicit recursive calls. */ /* */ std::stack e_stack, n_stack, i_stack; #ifdef TRACE stack_size = 0; max_stack_size = 0; #endif MarkNonRecursiveArc(e, n, i, e_stack, n_stack, i_stack); while (!e_stack.empty()) { assert((e_stack.size() == n_stack.size()) && (n_stack.size() == i_stack.size())); e = e_stack.top(); e_stack.pop(); n = n_stack.top(); n_stack.pop(); i = i_stack.top(); i_stack.pop(); #ifdef TRACE --stack_size; #endif MarkNonRecursiveArc(e, n, i, e_stack, n_stack, i_stack); } #ifdef TRACE printf("max_stack_size = %16d\n", max_stack_size); #endif return; } void vroniObject::InsertArcIntoVD(int i) { int j, k, j1, j2, n1, n2, n3; int e, e0, e1, e2, e3; t_site t0, t1, t2, t3; int i0, i1, i2, i3; int e1l, e1r, n1l, n1r, e2l, e2r, n2l, n2r, e3l, e3r, n3l, n3r; vr_bool node_shared, broken, preserved; #ifdef TRACE vr_bool trace_arc = true, trace_mark = true; vr_bool center = false; #endif #ifndef NDEBUG coord p1, p2, p3, p4, ns, ne; double a, b, c, radius1, radius2, radius, orient; #endif #ifndef NDEBUG p1 = GetArcStartCoord(i); p2 = GetArcEndCoord(i); p3 = GetArcCenter(i); radius = GetArcRadius(i); orient = VecDet(p3, p1, p2); assert(orient > 0.0); radius1 = PntPntDist(p1, p3); radius2 = PntPntDist(p2, p3); assert(eq(radius1 - radius2, ZERO_MAX)); assert(eq(radius1 - radius, ZERO_MAX)); assert(eq(radius2 - radius, ZERO_MAX)); GetChordEqnData(i, &a, &b, &c); p4 = MidPoint(p1, p2); assert(IsOnArc(a, b, c, p4)); ns = GetArcStartNormal(i); ne = GetArcEndNormal(i); p4 = VecSub(p4, p3); assert(IsInArcCone(ns, ne, p4) == 0); #endif #ifdef GRAPHICS if (graphics) { ResetCurBuffer(); ResetVDBuffer(); MarkDrawSite(i, ARC); AddToCurBuffer(i, ARC, SiteCurrColor); } #endif ResetPreserveBuffer; /* */ /* get the indices of the start point and of the end point of the arc. */ /* */ j1 = Get1stVtx(i, ARC); assert(InPntsList(j1)); j2 = Get2ndVtx(i, ARC); assert(InPntsList(j2)); #ifdef TRACE if (trace_arc) { printf("\ninserting arc %d (%d->%d):\n", i, j1, j2); printf("start: %20.16f %20.16f\nend: %20.16f %20.16f\n", pnts[j1].p.x, pnts[j1].p.y, pnts[j2].p.x, pnts[j2].p.y); printf("center: %20.16f %20.16f\n", (GetArcCenter(i)).x, (GetArcCenter(i)).y); } #endif SetThreshold(zero, TINY); /* */ /* find the nearest Voronoi nodes in the Voronoi cells of j1 and j2. */ /* */ do { n1 = FindMostViolatedNodeArc(j1, j2, i, &node_shared); if (n1 == NIL) { IncreaseThreshold(zero); if (gt(zero, ZERO_MAX)) { if (IsInvalidData(true)) { restart = true; return; } else { VD_Warning("InsertArcIntoVD() - desperate mode (MVN)"); n1 = DesperateFindMostViolatedNodeArc(j1, i); node_shared = false; assert(n1 != NIL); } } } #ifdef VRONI_DBG_WARN else if (gt(zero, ZERO) && !quiet) { printf("InsertArcIntoVD() - "); printf("most violated node not found till zero = %20.16f for arc %d", zero, i); printf("\n"); } #endif } while (n1 == NIL); #ifdef VRONI_INFO UpdateRelaxationMemory(zero); #endif SetThreshold(zero, TINY); assert(InNodesList(n1)); MarkNodeUnchecked(n1); if (node_shared) { n2 = n1; } else { do { n2 = FindMostViolatedNodeArc(j2, j1, i, &node_shared); if (n2 == NIL) { IncreaseThreshold(zero); if (gt(zero, ZERO_MAX)) { if (IsInvalidData(true)) { restart = true; return; } else { VD_Warning("InsertArcIntoVD() - desperate mode (MVN)"); n2 = DesperateFindMostViolatedNodeArc(j2, i); node_shared = false; assert(n2 != NIL); } } } #ifdef VRONI_DBG_WARN else if (gt(zero, ZERO) && !quiet) { printf("InsertArcIntoVD() - "); printf("most violated node not found till zero = %20.16f for ", zero); printf("arc %d\n", i); } #endif } while (n2 == NIL); assert(InNodesList(n2)); MarkNodeUnchecked(n2); } #ifdef VRONI_INFO UpdateRelaxationMemory(zero); #endif SetThreshold(zero, TINY); #ifdef TRACE if (trace_arc) printf("InsertArcIntoVD - start node is %d, end node is %d\n", n1, n2); #endif do { /* */ /* mark all Voronoi nodes that are to be deleted. */ /* */ #ifdef GRAPHICS if (graphics) { #ifdef TRACE AddToCurBuffer(n1, VDN, CirColor); AddToCurBuffer(n2, VDN, VDColor); #endif } #endif if (IsDeg2Node(n1)) { e1 = GetIncidentEdge(n1); e2 = GetCCWEdge(e1, n1); e3 = e1; GetLftSiteData(e1, &i1, &t1); GetRgtSiteData(e1, &i2, &t2); MarkNodeVisited(n1); #ifdef TRACE if (trace_mark) { printf("\n+++++ First call of MarkRecursiveArc()\n"); } #endif MarkRecursiveArc(e1, n1, i); #ifdef TRACE if (trace_mark) { printf("\n+++++ Second call of MarkRecursiveArc()\n"); } #endif MarkRecursiveArc(e2, n1, i); } else { e1 = GetIncidentEdge(n1); e2 = GetCCWEdge(e1, n1); e3 = GetCCWEdge(e2, n1); MarkNodeVisited(n1); GetLftSiteData(e1, &i1, &t1); GetRgtSiteData(e1, &i2, &t2); GetLftSiteData(e2, &i3, &t3); if (((i3 == i1) && (t3 == t1)) || ((i3 == i2) && (t3 == t2))) GetRgtSiteData(e2, &i3, &t3); assert(!(((i3 == i1) && (t3 == t1)) && ((i3 == i2) && (t3 == t2)))); #ifdef TRACE if (trace_mark) { printf("\n+++++ First call of MarkRecursiveArc()\n"); } #endif MarkRecursiveArc(e1, n1, i); #ifdef TRACE if (trace_mark) { printf("\n+++++ Second call of MarkRecursiveArc()\n"); } #endif MarkRecursiveArc(e2, n1, i); #ifdef TRACE if (trace_mark) { printf("\n+++++ Third call of MarkRecursiveArc()\n"); } #endif MarkRecursiveArc(e3, n1, i); } /* */ /* check whether we've reached the end node n2. if it has not been */ /* reached then we'll enlarge the tolerance threshold zero and we'll */ /* try again. */ /* */ if (!IsNodeVisited(n2)) { IncreaseThreshold(zero); if (gt(zero, ZERO_MAX)) { if (IsInvalidData(true)) { restart = true; return; } else { VD_Warning("InsertArcIntoVD() - desperate mode (path)"); MarkNodeUnchecked(n2); n2 = n1; desperate = true; troubles = true; } } else { ResetPreserveBuffer; ResetNodeStatus(e1, n1); ResetNodeStatus(e2, n1); if (!IsDeg2Node(n1)) ResetNodeStatus(e3, n1); } } #ifdef VRONI_DBG_WARN else if (gt(zero, ZERO) && !quiet) { printf("InsertArcIntoVD() - "); printf("not reached end node till zero = %20.16f for arc %d!\n", zero, i); } #endif } while (!IsNodeVisited(n2)); #ifdef VRONI_INFO UpdateRelaxationMemory(zero); #endif /* */ /* make sure that no VD cell is completely deleted. */ /* */ for (j = 0; j < num_preserve; ++j) { #ifdef TRACE printf("to be preserved: site %d of type %d\n", preserve[j].site, preserve[j].type); #endif SetThreshold(zero, ZERO); do { preserved = PreserveVoronoiCellArc(preserve[j].site, preserve[j].type, i); if (!preserved) { if (restart) return; IncreaseThreshold(zero); if (gt(zero, ZERO_MAX)) { VD_Dbg_Warning("InsertArcIntoVD() - cell not preserved!"); if (IsInvalidData(true)) { restart = true; return; } else { VD_Warning("InsertArcIntoVD() - desperate mode (preserve)"); DesperatePreserveVoronoiCellArc(preserve[j].site, preserve[j].type, i); preserved = true; } } } #ifdef VRONI_DBG_WARN else if (gt(zero, ZERO) && !quiet) { printf("InsertArcIntoVD() - "); printf("cell not preserved till zero = %20.16f for seg %d!\n", zero, i); } #endif } while (!preserved); #ifdef VRONI_INFO UpdateRelaxationMemory(zero); #endif } ResetPreserveBuffer; /* */ /* make sure that the structure that will be removed forms a tree. */ /* */ e1 = GetIncidentEdge(n1); MarkNodeDeleted(n1); ResetLoopStack(); while (LoopExists(e1, n1, &n3)) { PurifyLoop(n1, n3); SetThreshold(zero, ZERO); do { broken = BreakLoopArc(i); if (!broken) { if (restart) return; IncreaseThreshold(zero); if (gt(zero, ZERO_MAX)) { if (IsInvalidData(true)) { restart = true; return; } else { VD_Warning("InsertArcIntoVD() - desperate mode (loop)!"); DesperateBreakLoopArc(i); broken = true; } } } if (broken) { #ifdef VRONI_DBG_WARN if (gt(zero, ZERO) && le(zero, ZERO_MAX) && !quiet) { printf("InsertArcIntoVD() - "); printf("loop not broken till zero = %20.16f for arc %d!\n", zero, i); } #endif ResetLoopStack(); e1 = GetIncidentEdge(n1); ResetTreeDeleteStatus(e1, n1); if (n3 == n1) { e2 = GetCCWEdge(e1, n1); ResetTreeDeleteStatus(e2, n1); e3 = GetCCWEdge(e2, n1); if (e3 != e1) ResetTreeDeleteStatus(e3, n1); } MarkNodeDeleted(n1); } } while (!broken); #ifdef VRONI_INFO UpdateRelaxationMemory(zero); #endif } e1 = GetIncidentEdge(n1); e2 = GetCCWEdge(e1, n1); ResetLoopStack(); while (LoopExists(e2, n1, &n3)) { PurifyLoop(n1, n3); SetThreshold(zero, ZERO); do { broken = BreakLoopArc(i); if (!broken) { if (restart) return; IncreaseThreshold(zero); if (gt(zero, ZERO_MAX)) { if (IsInvalidData(true)) { restart = true; return; } else { VD_Warning("InsertArcIntoVD() - desperate mode (loop)!"); DesperateBreakLoopArc(i); broken = true; } } } if (broken) { #ifdef VRONI_DBG_WARN if (gt(zero, ZERO) && le(zero, ZERO_MAX) && !quiet) { printf("InsertArcIntoVD() - "); printf("loop not broken till zero = %20.16f for arc %d!\n", zero, i); } #endif ResetLoopStack(); e2 = GetCCWEdge(e1, n1); ResetTreeDeleteStatus(e2, n1); if (n3 == n1) { e3 = GetCCWEdge(e2, n1); if (e3 != e1) ResetTreeDeleteStatus(e3, n1); } } } while (!broken); #ifdef VRONI_INFO UpdateRelaxationMemory(zero); #endif } e2 = GetCCWEdge(e1, n1); e3 = GetCCWEdge(e2, n1); ResetLoopStack(); if (e3 != e1) { while (LoopExists(e3, n1, &n3)) { PurifyLoop(n1, n3); SetThreshold(zero, ZERO); do { broken = BreakLoopArc(i); if (!broken) { if (restart) return; IncreaseThreshold(zero); if (gt(zero, ZERO_MAX)) { if (IsInvalidData(true)) { restart = true; return; } else { VD_Warning("InsertArcIntoVD() - desperate mode (loop)!"); DesperateBreakLoopArc(i); broken = true; } } } if (broken) { #ifdef VRONI_DBG_WARN if (gt(zero, ZERO) && le(zero, ZERO_MAX) && !quiet) { printf("InsertArcIntoVD() - "); printf("loop not broken till zero = %20.16f for arc %d!\n", zero, i); } #endif ResetLoopStack(); e3 = GetCCWEdge(e2, n1); ResetTreeDeleteStatus(e3, n1); } } while (!broken); #ifdef VRONI_INFO UpdateRelaxationMemory(zero); #endif } } /* */ /* recursively discard all nodes and edges that are to be deleted. */ /* also, we insert the new Voronoi nodes and edges. */ /* */ if (IsDeg2Node(n1)) { e1 = GetIncidentEdge(n1); e2 = GetCCWEdge(e1, n1); #ifdef TRACE if (center) printf("*** 1st call of UpdateRecursive - i = %d, n1 = %d, e1 = %d\n", i, n1, e1); #endif UpdateRecursive(i, ARC, n1, e1, &e1l, &e1r, &n1l, &n1r); if (restart) return; #ifdef TRACE if (center) printf("*** 2nd call of UpdateRecursive - i = %d, n1 = %d, e2 = %d\n", i, n1, e2); #endif UpdateRecursive(i, ARC, n1, e2, &e2l, &e2r, &n2l, &n2r); if (restart) return; /* */ /* close the two Voronoi cells around nodes[n], which is to be */ /* deleted. */ /* */ /* cell 1: */ /* */ if (GetEndNode(e1r) == n1r) { GetRgtSiteData(e1r, &i0, &t0); } else { assert(GetStartNode(e1r) == n1r); GetLftSiteData(e1r, &i0, &t0); } e0 = StoreEdge(n1r, n2l, i, i0, ARC, t0); #ifdef TRACE if (center) { printf("new edge %2d between sites (%2d/0, %2d/%1d) computed\n", e0, i, i0, t0); printf("edge %d - start node %d, end node %d\n", e0, n1r, n2l); } #endif SetVDPtr(i0, t0, e0); /* make sure the VD edge pointer is up-to-date! */ CloseVoronoiCell(e1r, e0, e2l, n1r, n2l); /* */ /* cell 2: */ /* */ if (GetEndNode(e2r) == n2r) { GetRgtSiteData(e2r, &i0, &t0); } else { assert(GetStartNode(e2r) == n2r); GetLftSiteData(e2r, &i0, &t0); } e0 = StoreEdge(n2r, n1l, i, i0, ARC, t0); #ifdef TRACE if (center) { printf("new edge %2d between sites (%2d/0, %2d/%1d) computed\n", e0, i, i0, t0); printf("edge %d - start node %d, end node %d\n", e0, n2r, n1l); } #endif SetVDPtr(i0, t0, e0); /* make sure the VD edge pointer is up-to-date! */ CloseVoronoiCell(e2r, e0, e1l, n2r, n1l); } else { e1 = GetIncidentEdge(n1); e2 = GetCCWEdge(e1, n1); e3 = GetCCWEdge(e2, n1); #ifdef TRACE if (center) printf("*** 1st call of UpdateRecursive - i = %d, n1 = %d, e1 = %d\n", i, n1, e1); #endif UpdateRecursive(i, ARC, n1, e1, &e1l, &e1r, &n1l, &n1r); if (restart) return; #ifdef TRACE if (center) printf("*** 2nd call of UpdateRecursive - i = %d, n1 = %d, e2 = %d\n", i, n1, e2); #endif UpdateRecursive(i, ARC, n1, e2, &e2l, &e2r, &n2l, &n2r); if (restart) return; #ifdef TRACE if (center) printf("*** 3rd call of UpdateRecursive - i = %d, n1 = %d, e3 = %d\n", i, n1, e3); #endif UpdateRecursive(i, ARC, n1, e3, &e3l, &e3r, &n3l, &n3r); if (restart) return; /* */ /* close the three Voronoi cells around nodes[n], which is to be */ /* deleted. */ /* */ /* cell 1: */ /* */ if (GetEndNode(e1r) == n1r) { GetRgtSiteData(e1r, &i0, &t0); } else { assert(GetStartNode(e1r) == n1r); GetLftSiteData(e1r, &i0, &t0); } e0 = StoreEdge(n1r, n2l, i, i0, ARC, t0); #ifdef TRACE if (center) { printf("new edge %2d between sites (%2d/0, %2d/%1d) computed\n", e0, i, i0, t0); printf("edge %d - start node %d, end node %d\n", e0, n1r, n2l); } #endif SetVDPtr(i0, t0, e0); /* make sure the VD edge pointer is up-to-date! */ CloseVoronoiCell(e1r, e0, e2l, n1r, n2l); /* */ /* cell 2: */ /* */ if (GetEndNode(e2r) == n2r) { GetRgtSiteData(e2r, &i0, &t0); } else { assert(GetStartNode(e2r) == n2r); GetLftSiteData(e2r, &i0, &t0); } e0 = StoreEdge(n2r, n3l, i, i0, ARC, t0); #ifdef TRACE if (center) { printf("new edge %2d between sites (%2d/0, %2d/%1d) computed\n", e0, i, i0, t0); printf("edge %d - start node %d, end node %d\n", e0, n2r, n3l); } #endif SetVDPtr(i0, t0, e0); /* make sure the VD edge pointer is up-to-date! */ CloseVoronoiCell(e2r, e0, e3l, n2r, n3l); /* */ /* cell 3: */ /* */ if (GetEndNode(e3r) == n3r) { GetRgtSiteData(e3r, &i0, &t0); } else { assert(GetStartNode(e3r) == n3r); GetLftSiteData(e3r, &i0, &t0); } e0 = StoreEdge(n3r, n1l, i, i0, ARC, t0); #ifdef TRACE if (center) { printf("new edge %2d between sites (%2d/0, %2d/%1d) computed\n", e0, i, i0, t0); printf("edge %d - start node %d, end node %d\n", e0, n3r, n1l); } #endif SetVDPtr(i0, t0, e0); /* make sure the VD edge pointer is up-to-date! */ CloseVoronoiCell(e3r, e0, e1l, n3r, n1l); } /* */ /* delete node nodes[n1] */ /* */ DeleteNode(n1); #ifdef TRACE if (center) printf("node %2d has been deleted!\n", n1); #endif /* */ /* set the Voronoi pointer for site arcs[i] */ /* */ SetVDPtr(i, ARC, e0); DeleteDummyNodes(); /* */ /* make sure to have (dummy) degree-2 Voronoi nodes at the apices of the */ /* conic Voronoi arcs around arc[i]. */ /* */ e = GetVDPtr(i, ARC); assert(InEdgesList(e)); assert(IsLftSite(e, i, ARC) || IsRgtSite(e, i, ARC)); if (IsRgtSite(e, i, ARC)) k = GetEndNode(e); else k = GetStartNode(e); e0 = e; GetOtherSiteData(e, i, ARC, &i1, &t1); do { i0 = i1; t0 = t1; if (t1 == PNT || t1 == ARC) InsertApexDegreeTwoNodeArc(e, i, i1, t1); if (t1 == SEG ) InsertApexDegreeTwoNodeSeg(e, i1, i, ARC); e = GetCCWEdge(e, k); k = GetOtherNode(e, k); GetOtherSiteData(e, i, ARC, &i1, &t1); if ((i1 == i0) && (t1 == t0) && (e != e0)) { e = GetCCWEdge(e, k); k = GetOtherNode(e, k); GetOtherSiteData(e, i, ARC, &i1, &t1); } } while (e != e0); #ifndef NDEBUG i1 = Get1stVtx(i, ARC); assert(GetVDPtr(i1, PNT) != NIL); i2 = Get2ndVtx(i, ARC); assert(GetVDPtr(i2, PNT) != NIL); #endif #ifdef GRAPHICS if (graphics) { e = GetVDPtr(i, ARC); assert(InEdgesList(e)); assert(IsLftSite(e, i, ARC) || IsRgtSite(e, i, ARC)); if (IsRgtSite(e, i, ARC)) k = GetEndNode(e); else k = GetStartNode(e); e0 = e; do { GetOtherSiteData(e, i, ARC, &i1, &t1); AddToCurBuffer(e, VDE, VDCurrColor); e = GetCCWEdge(e, k); k = GetOtherNode(e, k); } while (e != e0); } #endif return; }