/*****************************************************************************/ /* */ /* 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-172 */ /* Voice Mail: (+43 662) 8044-6304 */ /* Snail Mail: Martin Held */ /* FB Informatik */ /* Universitaet Salzburg */ /* A-5020 Salzburg, Austria */ /* */ /*****************************************************************************/ /* */ /* get standard libraries */ /* */ #include #include #include #include #include /* */ /* get my header files */ /* */ #include "fpkernel.h" #include "vronivector.h" #include "vroni_object.h" #include "defs.h" #include "numerics.h" #include "offset.h" #include "stack.h" #include "wmat.h" /* */ /* local stack */ /* */ void vroniObject::StoreActiveEdge(int e) { assert(num_active_edges >= 0); assert(num_active_edges < max_num_active_edges); assert(InEdgesList(e)); active_edges[num_active_edges] = e; ++num_active_edges; } vr_bool vroniObject::InEdgeDataList(int i) { return ((i >= 0) && (i < num_edge_data)); } vr_bool vroniObject::InActiveEdgeList(int i) { return ((i >= 0) && (i < num_active_edges)); } vr_bool vroniObject::InEdgeFlagList(int i) { return ((i >= 0) && (i < num_edge_flags)); } /* void MarkUnboundedRegion(int e, int n) { double t; int e1; while (GetEdgeFlagNew(e)) { SetEdgeFlagNew(e, false); n = GetOtherNode(e, n); t = GetNodeParam(n); if (gt(t, ZERO)) { e1 = GetCCWEdge(e, n); if (GetEdgeFlagNew(e1)) { e = GetCWEdge(e, n); if (GetEdgeFlagNew(e) && (e1 != e)) { Push(e); Push(n); } e = e1; } else { e = GetCWEdge(e, n); while (!GetEdgeFlagNew(e) && StackIsNotEmpty) { Pop(n); Pop(e); } } } else { while (!GetEdgeFlagNew(e) && StackIsNotEmpty) { Pop(n); Pop(e); } } } return; } */ /* */ /* for t1 <= t3 <= t2 this function computes */ /* */ /* lambda = (t1 - t3) / (t2 - t1) */ /* */ /* iteratively by means of bisection. hence it works also if t1 and t2 are */ /* nearly equal. */ /* */ double DegenerateInterpolation(double_arg t1, double_arg t2, double_arg t3) { double lambda = 0.5, lambda_min = 0.0, lambda_max = 1.0, t; int cntr = 0; assert(t1 <= t2); if (t3 <= t1) return 0.0; else if (t3 >= t2) return 1.0; while (cntr < 60) { assert((lambda_min <= lambda) && (lambda <= lambda_max)); lambda = (lambda_min + lambda_max) / 2.0; t = (1.0 - lambda) * t1 + lambda * t2; if (t < t3) { if (t <= t1) return lambda_max; else { if ((lambda < lambda_max) && (lambda > lambda_min)) lambda_min = lambda; else return lambda_max; } } else if (t > t3) { if (t >= t2) return lambda_min; else { if ((lambda < lambda_max) && (lambda > lambda_min)) lambda_max = lambda; else return lambda_min; } } else { return lambda; } ++cntr; } return lambda; } vr_bool vroniObject::IsCorrectSide(int i, int type, int n, vr_bool left_offset) { coord p, q; double a, b, c, r, dist; assert(type==ARC || InSegsList(i)); assert(type==SEG || InArcsList(i)); assert(InNodesList(n)); /* */ /* let a user influence this function by replacing my sidedness test by */ /* a user's sidedness test. */ /* */ ExtApplIsCorrectSide; /* */ /* get position of node */ /* */ p = GetNodeCoord(n); if( type == SEG ) { /* */ /* It's a segment... */ /* */ /* my sidedness test is based on the relative position of node n */ /* with respect to the line segment i. */ /* */ /* */ GetSegEqnData(i, &a, &b, &c); if (left_offset) return (PntLineDist(a, b, c, p) > 0.0); else return (PntLineDist(a, b, c, p) < 0.0); } else { /* */ /* It's a circular arc... */ /* */ q = GetArcCenter(i); r = GetArcRadius(i); dist = PntPntDist(q,p); //printf("arc a%d, node n%d, left_offset=%d. dist=%e, r=%e.\n", i, n, left_offset, dist, r); /* */ /* depending on the orientation of the arc, determine the sidedness of */ /* the point... */ /* */ if(GetArcOrientation(i)) { //printf("\tCCW orientation of arc, correct side: %d\n", left_offset ? distr); return (left_offset ? (distr)); } else { //printf("\tCW orientation of arc, correct side: %d\n", left_offset ? dist>r : distr) : (dist *max_rad) { *max_rad = t; *max_edge = e; } } else { StoreActiveEdge(e); /* */ /* user-specific code, e.g., in order to output the boundary of */ /* the face of the arrangement which is traversed by this scan */ /* */ ExtApplScanVDEdges; } } if (gt(t, ZERO)) { e1 = GetCCWEdge(e, n); if (GetEdgeFlagNew(e1)) { e = GetCWEdge(e, n); if (GetEdgeFlagNew(e) && (e1 != e)) { Push(e); Push(n); } e = e1; } else { e = GetCWEdge(e, n); while (!GetEdgeFlagNew(e) && StackIsNotEmpty) { Pop(n); Pop(e); } } } else { while (!GetEdgeFlagNew(e) && StackIsNotEmpty) { Pop(n); Pop(e); } } } return; } /* */ /* this function determines the correct bisectors if one-sided offsetting or */ /* or one-sided WMAT computation is required, and allocates and initializes */ /* the appropriate arrays. if requested, it also determines the center of */ /* the maximum inscribed/empty circle. */ /* */ void vroniObject::InitializeEdgeData(vr_bool unrestricted, vr_bool left_offset, vr_bool compute_mic, int *mic_edge) { int e, e1, e2, i, n, num_edges; t_site i_type; double t1, t2; vr_bool offset_side; double max_rad = 0.0; int max_edge = NIL; *mic_edge = NIL; num_edges = GetNumberOfEdges(); if (num_edges >= max_num_edge_data) { num_edge_data = num_edges; max_num_edge_data = num_edge_data + 1; gentlyResizeSTLVector(edge_data, max_num_edge_data, "data_off:edge_data"); for (e = 0; e < num_edges; ++e) { SetEdgeDataInit(e, true); #ifdef MAT SetWmatEdge(e, false); #endif } data_off_VD_ID = GetVDIdentifier(); } else { num_edge_data = num_edges; i = GetVDIdentifier(); if (data_off_VD_ID != i) { data_off_VD_ID = i; for (e = 0; e < num_edges; ++e) SetEdgeDataInit(e, true); #ifdef MAT if (!GetWMATStatus()) { for (e = 0; e < num_edges; ++e) SetWmatEdge(e, false); } #endif } } if (num_edges >= max_num_edge_flags) { num_edge_flags = num_edges; max_num_edge_flags = num_edge_flags + 1; edge_flags = (e_flag*) ReallocateArray(edge_flags, max_num_edge_flags, sizeof(e_flag), "data_off:edge_flags"); } else { num_edge_flags = num_edges; } if (num_edges >= max_num_active_edges) { max_num_active_edges = num_edges + 1; active_edges = (int*) ReallocateArray(active_edges, max_num_active_edges, sizeof(int), "data_off:active_edges"); } /* */ /* all Voronoi edges that might be intersected by an offset curve are put */ /* into a list of "active" edges. note that we deal only with those edges */ /* that are not defined by the four dummy points; see IsEdgeGenuine(). */ /* */ num_active_edges = 0; min_pnt_ind = 1; max_pnt_ind = num_pnts - 2; if (unrestricted) { /* */ /* offsetting to be performed at both sides of the input line sites */ /* */ for (e = 4; e < num_edges; ++e) { SetEdgeFlagNew(e, true); if (IsEdgeGenuine(e, min_pnt_ind, max_pnt_ind)) StoreActiveEdge(e); } SetEdgeFlagNew(0, false); SetEdgeFlagNew(1, false); SetEdgeFlagNew(2, false); SetEdgeFlagNew(3, false); } else { /* */ /* offsetting to be performed at only one side of the input line sites */ /* */ for (e = 0; e < num_edges; ++e) { SetEdgeFlagNew(e, true); } /* * The following code _s_h_o_u_l_d_ be good enough to discard all VD * edges of the (one) unbounded region. By making IsCorrecSide() return * true, this modification can be used to computed offsets within all * bounded regions of the partition induced by a (possibly * self-intersecting) input. * e = 0; GetEdgeParam(e, &t1, &t2); if (t1 < t2) n = GetStartNode(e); else n = GetEndNode(e); if (n != NIL) { MarkUnboundedRegion(e, n); } else { VD_Warning("cannot mark unbounded edges!\n"); } */ ExtApplInitializeEdgeData; SetEdgeFlagNew(0, false); SetEdgeFlagNew(1, false); SetEdgeFlagNew(2, false); SetEdgeFlagNew(3, false); for (e = 4; e < num_edges; ++e) { /* */ /* check for every VD edge on which side of the input it lies. once */ /* the relative position of one VD edge is known, the positions of */ /* all neighboring edges are known, too. we do not cross input */ /* boundaries, though. */ /* */ if (GetEdgeFlagNew(e) && IsEdgeGenuine(e, min_pnt_ind, max_pnt_ind)) { GetLftSiteData(e, &i, &i_type); if (i_type == PNT) GetRgtSiteData(e, &i, &i_type); if (i_type == SEG || i_type == ARC ) { //Let n be the node of e with the higher clearance. //and NIL if both clearances are zero. GetEdgeParam(e, &t1, &t2); if (t1 >= t2) { if (gt(t1, ZERO)) n = GetStartNode(e); else n = NIL; } else { if (gt(t2, ZERO)) n = GetEndNode(e); else n = NIL; } if (n != NIL) { /* */ /* we have not yet checked this VD edge. determine its */ /* relative position, and mark all neighboring edges. */ /* */ offset_side = IsCorrectSide(i, i_type, n, left_offset); if (compute_mic && offset_side) { if (t1 > max_rad) { max_rad = t1; max_edge = e; } if (t2 > max_rad) { max_rad = t2; max_edge = e; } } ScanVDEdges(e, n, offset_side, compute_mic, &max_edge, &max_rad); e1 = GetCCWEdge(e, n); if (e1 >= 4) ScanVDEdges(e1, n, offset_side, compute_mic, &max_edge, &max_rad); e2 = GetCWEdge(e, n); if ((e2 >= 4) && (e1 != e2)) ScanVDEdges(e2, n, offset_side, compute_mic, &max_edge, &max_rad); } } } } for (e = 4; e < num_edges; ++e) { SetEdgeFlagNew(e, true); } } e_data_init = false; *mic_edge = max_edge; return; } void vroniObject::ResetEdgeData(void) { e_data_init = true; num_edge_flags = 0; num_edge_data = 0; num_active_edges = 0; num_stack = 0; return; } void vroniObject::FreeEdgeData(void) { edge_data.clear(); FreeMemory((void**) &edge_flags, "data_off:edge_flags"); FreeMemory((void**) &active_edges, "data_off:active_edges"); FreeMemory((void**) &stack, STACK_STRING); num_stack = 0; max_num_stack = 0; e_data_init = true; num_edge_flags = 0; num_edge_data = 0; num_active_edges = 0; max_num_edge_flags = 0; max_num_edge_data = 0; max_num_active_edges = 0; return; } vr_bool vroniObject::ComputeParabolaData(int i, e_formula *coeff) { int i1, i2; t_site type, arc_type; coord p, u, v, w, u_max, u_min; double a, b, c, A, B, t1, t2, dist, du, dv, dw, r, t_max, t_min, dx, dy; int n1, n2, k1, k2; n1 = GetStartNode(i); n2 = GetEndNode(i); GetNodeData(n1, &u, &t1); GetNodeData(n2, &v, &t2); dist = t1 - t2; if (eq(dist, ZERO)) { /* */ /* degenerate Voronoi edge: most likely, this is a very short edge */ /* */ return false; } GetLftSiteData(i, &i1, &type); if (type == SEG) { GetRgtSiteData(i, &i2, &arc_type); assert((arc_type == PNT) || (arc_type == ARC)); } else { arc_type = type; assert((arc_type == PNT) || (arc_type == ARC)); i2 = i1; GetRgtSiteData(i, &i1, &type); assert(type == SEG); } if (t1 > t2) { u_max = u; u_min = v; t_max = t1; t_min = t2; } else { u_max = v; u_min = u; t_max = t2; t_min = t1; } GetSegEqnData(i1, &a, &b, &c); if (arc_type == PNT) { p = GetPntCoords(i2); r = 0.0; k1 = -1; dist = du = PntLineDist(a, b, c, p); if (eq(dist, ZERO)) du = PntLineDist(a, b, c, u_max); } else { r = GetArcRadius(i2); p = GetArcCenter(i2); dist = PntLineDist(a, b, c, p); du = PntLineDist(a, b, c, u_max); dw = PntPntDist(u_max, p); if (dw > r) k1 = 1; else k1 = -1; } assert(!eq(du, ZERO)); if (du > 0.0) k2 = -1; else k2 = 1; /* */ /* parabola formula w.r.t. clearance distance t: */ /* x(t) = p.x - a * dist - k2 * a * t + */ /* + sign * b * sqrt(r^2 + 2 * t * (r * k1 - dist * k2) - dist^2) */ /* y(t) = p.y - b * dist - k2 * b * t - */ /* + sign * a * sqrt(r^2 + 2 * t * (r * k1 - dist * k2) - dist^2) */ /* i.e., */ /* x(t) = A - k2 * a * t + sign * b * sqrt(discriminant) */ /* y(t) = B - k2 * b * t - sign * a * sqrt(discriminant) */ /* where A := p.x - a * dist, B := p.y - b * dist, and */ /* discriminant := r^2 + 2 * t * (r * k1 - dist * k2) - dist^2 */ /* */ A = p.x - a * dist; B = p.y - b * dist; xParabola(A, a, b, 1.0, dist, t_max, k1, k2, r, v.x); yParabola(B, a, b, 1.0, dist, t_max, k1, k2, r, v.y); xParabola(A, a, b, -1.0, dist, t_max, k1, k2, r, w.x); yParabola(B, a, b, -1.0, dist, t_max, k1, k2, r, w.y); dv = PntPntDistSq(u_max, v); dw = PntPntDistSq(u_max, w); if (eq(dv, ZERO_MAX) && eq(dw, ZERO_MAX)) { xParabola(A, a, b, 1.0, dist, t_min, k1, k2, r, v.x); yParabola(B, a, b, 1.0, dist, t_min, k1, k2, r, v.y); xParabola(A, a, b, -1.0, dist, t_min, k1, k2, r, w.x); yParabola(B, a, b, -1.0, dist, t_min, k1, k2, r, w.y); dx = PntPntDistSq(u_min, v); dy = PntPntDistSq(u_min, w); if (eq(dx, ZERO_MAX) && eq(dy, ZERO_MAX)) { if (VroniMax(dv, dx) > VroniMax(dw, dy)) coeff->sign = false; else coeff->sign = true; } else { if (dx < dy) coeff->sign = true; else coeff->sign = false; } } else { if (dv < dw) coeff->sign = true; else coeff->sign = false; } coeff->A = A; coeff->B = B; coeff->d = dist; if (k1 == 1) coeff->k1 = true; else coeff->k1 = false; if (k2 == 1) coeff->k2 = true; else coeff->k2 = false; return true; } /* */ /* compute (and store) the coefficients of my parameterization formula for a */ /* parabolic VD edge defined by a line segment and a point. note that the */ /* parameterization uses the clearance radius as parameter. */ /* */ void vroniObject::CreateParabolaData(int i) { assert(!e_data_init); assert(InEdgeFlagList(i)); assert(InEdgeDataList(i)); if (ComputeParabolaData(i, &(edge_data[i]))) { SetEdgeDataInit(i, false); SetEdgeFlagDeg(i, false); } else { SetEdgeFlagDeg(i, true); SetEdgeDataInit(i, false); } return; } /* */ /* this is an attempt to intersect a VD edge defined by two nearly parallel */ /* line segments with an offset segment at offset distance t. */ /* */ void vroniObject::EvaluateDegenerateEdge(int i, double_arg t, coord *w) { int n1, n2; coord u, v, p; double delta, t1, t2; assert(InEdgesList(i)); assert(!GetEdgeDataInit(i)); n1 = GetStartNode(i); n2 = GetEndNode(i); GetNodeData(n1, &u, &t1); GetNodeData(n2, &v, &t2); if (t1 == t) { *w = u; return; } else if (t2 == t) { *w = v; return; } if (t1 > t2) { VroniSwap(u, v, p); VroniSwap(t1, t2, delta); } assert(t1 <= t2); assert((t1 <= t) && (t <= t2)); delta = DegenerateInterpolation(t1, t2, t); *w = LinearComb(u, v, delta); return; } /* */ /* evaluate the parameterization formula for the parabolic VD edge i that */ /* has the input segment i1 as one of its defining sites. */ /* */ void vroniObject::EvaluateParabolaData(int i, int i1, double_arg t, coord *w) { int i2; double sign, a, b, c, A, B, dist, k1, k2, r; t_site arc_type; #ifndef NDEBUG int n1, n2; double t1, t2; coord u, v; #endif #ifndef NDEBUG assert(InEdgesList(i)); assert(InSegsList(i1)); assert(!GetEdgeDataInit(i)); n1 = GetStartNode(i); n2 = GetEndNode(i); GetNodeData(n1, &u, &t1); GetNodeData(n2, &v, &t2); if (t1 < t2) assert((t1 <= t) && (t <= t2)); else assert((t2 <= t) && (t <= t1)); #endif if (GetEdgeFlagDeg(i)) { /* */ /* degenerate Voronoi edge: */ /* */ EvaluateDegenerateEdge(i, t, w); } else { GetSegEqnData(i1, &a, &b, &c); A = edge_data[i].A; B = edge_data[i].B; dist = edge_data[i].d; if (edge_data[i].sign) sign = 1.0; else sign = -1.0; if (edge_data[i].k1) k1 = 1.0; else k1 = -1.0; if (edge_data[i].k2) k2 = 1.0; else k2 = -1.0; if (IsLftSite(i, i1, SEG)) { GetRgtSiteData(i, &i2, &arc_type); } else { GetLftSiteData(i, &i2, &arc_type); } if (arc_type == PNT) r = 0.0; else r = GetArcRadius(i2); xParabola(A, a, b, sign, dist, t, k1, k2, r, w->x); yParabola(B, a, b, sign, dist, t, k1, k2, r, w->y); } return; } /* */ /* evaluate the parameterization formula for the hyperbolic VD edge i */ /* defined by two points. */ /* */ void vroniObject::EvaluateHyperbolaData(int i, double t, coord *w) { double dist, delta; int n1, n2; coord u, v; double t1, t2; assert(InEdgesList(i)); assert(!GetEdgeDataInit(i)); n1 = GetStartNode(i); n2 = GetEndNode(i); GetNodeData(n1, &u, &t1); GetNodeData(n2, &v, &t2); #ifndef NDEBUG if (t1 < t2) assert((t1 <= t) && (t <= t2)); else assert((t2 <= t) && (t <= t1)); #endif if (GetEdgeFlagDeg(i)) { /* */ /* degenerate Voronoi edge */ /* */ EvaluateDegenerateEdge(i, t, w); } else { if (t1 > t2) { t1 = t2; u = v; } v.x = edge_data[i].A; v.y = edge_data[i].B; dist = edge_data[i].d; delta = t1 * t1 - dist; if (delta < 0.0) delta = 0.0; dist = t * t - dist; if (dist < 0.0) { #ifdef VRONI_WARN if (lt(dist, ZERO)) printf("warning in EvaluateHyperbolaBisector() - dist = %f\n", dist); #endif dist = 0.0; } t = sqrt(dist) - sqrt(delta); *w = RayPnt(u, v, t); } return; } vr_bool vroniObject::ComputeHyperbolaData(int i, e_formula *coeff) { int i1, i2; t_site ltype, rtype; coord p1, p2, u, v, w; double lgth, t1, t2, dist; int n1, n2; /* */ /* this is a special hyperbolic bisector that degenerates to a line */ /* */ n1 = GetStartNode(i); n2 = GetEndNode(i); GetNodeData(n1, &u, &t1); GetNodeData(n2, &v, &t2); dist = t1 - t2; GetLftSiteData(i, &i1, <ype); assert(ltype == PNT); GetRgtSiteData(i, &i2, &rtype); assert(rtype == PNT); p1 = GetPntCoords(i1); p2 = GetPntCoords(i2); dist = PntPntDistSq(p1, p2) / 4.0; if (t1 > t2) w = VecSub(u, v); else w = VecSub(v, u); lgth = VecLen(w); if (eq(lgth, ZERO)) { return false; } SetEdgeDataInit(i, false); SetEdgeFlagDeg(i, false); coeff->A = w.x / lgth; coeff->B = w.y / lgth; coeff->d = dist; return true; } /* */ /* compute (and store) the coefficients of my parameterization formula for a */ /* hyperbolic VD edge defined by two point sites. note that the */ /* parameterization uses the clearance radius as parameter. */ /* */ void vroniObject::CreateHyperbolaData(int i) { assert(!e_data_init); assert(InEdgeFlagList(i)); assert(InEdgeDataList(i)); if (ComputeHyperbolaData(i, &(edge_data[i]))) { SetEdgeDataInit(i, false); SetEdgeFlagDeg(i, false); } else { SetEdgeFlagDeg(i, true); SetEdgeDataInit(i, false); } return; } /* */ /* evaluate my parameterization formula for an elliptic or hyperbolic arc */ /* defined by two arcs or a point and an arc. */ /* */ /* */ void vroniObject::EvaluateDegenerateHyperEll(int i, double t, coord *w) { int i1, i2, n1, n2; t_site ltype, rtype; coord u, v, p1, p2, p; double r1, r2, dx, dy, dist, delta, radius, alpha, angle_s, angle_e; double t1, t2; GetLftSiteData(i, &i1, <ype); GetRgtSiteData(i, &i2, &rtype); assert((ltype == PNT) || (ltype == ARC)); assert((rtype == PNT) || (rtype == ARC)); n1 = GetStartNode(i); n2 = GetEndNode(i); GetNodeData(n1, &u, &t1); GetNodeData(n2, &v, &t2); if (t1 == t) { *w = u; return; } else if (t2 == t) { *w = v; return; } if (ltype == PNT) { r1 = 0.0; p1 = GetPntCoords(i1); } else { r1 = GetArcRadius(i1); p1 = GetArcCenter(i1); } if (rtype == PNT) { r2 = 0.0; p2 = GetPntCoords(i2); } else { r2 = GetArcRadius(i2); p2 = GetArcCenter(i2); } dx = p2.x - p1.x; dy = p2.y - p1.y; dist = sqrt(dx * dx + dy * dy); if (eq(dist, ZERO_MAX)) { dist = r1 - r2; if (eq(dist, ZERO_MAX)) { /* */ /* the bisector is a line segment! */ /* */ dist = t1 - t2; if (eq(dist, ZERO)) { EvaluateDegenerateEdge(i, t, w); } else { if (t1 > t2) { t = (t - t2) / (t1 - t2); u = v; } else { t = (t - t1) / (t2 - t1); } v.x = edge_data[i].A; v.y = edge_data[i].B; *w = RayPnt(u, v, t); } } else { /* */ /* degenerate elliptic edge */ /* */ p = MidPoint(p1, p2); if (VecDet(p, u, v) < 0.0) { VroniSwap(u, v, *w); VroniSwap(t1, t2, delta); } angle_s = atan2(u.y - p.y, u.x - p.x); angle_e = atan2(v.y - p.y, v.x - p.x); if (angle_s < 0.0) angle_s += M_2PI; if (angle_e < 0.0) angle_e += M_2PI; if (angle_e < angle_s) angle_e += M_2PI; if (t1 <= t2) delta = DegenerateInterpolation(t1, t2, t); else delta = 1.0 - DegenerateInterpolation(t2, t1, t); assert((0.0 <= delta) && (delta <= 1.0)); radius = (1.0 - delta) * PntPntDist(p, u) + delta * PntPntDist(p, v); alpha = (1.0 - delta) * angle_s + delta * angle_e; w->x = p.x + radius * cos(alpha); w->y = p.y + radius * sin(alpha); } } else { /* */ /* degenerate hyperbolic edge */ /* */ EvaluateDegenerateEdge(i, t, w); } return; } /* */ /* evaluate the parameterization formula for the hyperbolic or elliptic VD */ /* edge i defined by two arcs or a point and an arc. */ /* */ #ifdef GENUINE_ARCS void vroniObject::EvaluateHyperbolaEllipseData(int i, double_arg t, coord *w) { int i1, i2, i3; t_site ltype, rtype, itype; double r1, r2; coord p1, p2; double A, B, C, D, dist; double r1t, r2t, dx, dy, ht; double k1, k2, sign; #ifndef NDEBUG int n1, n2; coord u, v; double t1, t2; #endif #ifndef NDEBUG assert(InEdgesList(i)); assert(!GetEdgeDataInit(i)); n1 = GetStartNode(i); n2 = GetEndNode(i); GetNodeData(n1, &u, &t1); GetNodeData(n2, &v, &t2); if (t1 < t2) assert((t1 <= t) && (t <= t2)); else assert((t2 <= t) && (t <= t1)); #endif if (GetEdgeFlagDeg(i)) { /* */ /* degenerate Voronoi edge: most likely, this is a very short edge */ /* */ EvaluateDegenerateHyperEll(i, t, w); } else { GetLftSiteData(i, &i1, <ype); GetRgtSiteData(i, &i2, &rtype); assert((ltype == PNT) || (ltype == ARC)); assert((rtype == PNT) || (rtype == ARC)); if (i2 > i1) { VroniSwap(i1, i2, i3); VroniSwap(ltype, rtype, itype); } else if (i1 == i2) { if (rtype == PNT) { VroniSwap(i1, i2, i3); VroniSwap(ltype, rtype, itype); } } if (ltype == PNT) { r1 = 0.0; p1 = GetPntCoords(i1); } else { r1 = GetArcRadius(i1); p1 = GetArcCenter(i1); } if (rtype == PNT) { r2 = 0.0; p2 = GetPntCoords(i2); } else { r2 = GetArcRadius(i2); p2 = GetArcCenter(i2); } dx = (p2.x - p1.x); dy = (p2.y - p1.y); dist = edge_data[i].d; if (eq(dist, ZERO)) { EvaluateDegenerateHyperEll(i, t, w); return; } dx /= dist; dy /= dist; A = edge_data[i].A; B = edge_data[i].B; C = edge_data[i].C; D = edge_data[i].D; if (edge_data[i].k1) k1 = 1.0; else k1 = -1.0; if (edge_data[i].k2) k2 = 1.0; else k2 = -1.0; if (edge_data[i].sign) sign = 1.0; else sign = -1.0; xHyperEll(A, C, sign, dist, dy, t, k1, k2, r1, r2, r1t, r2t, ht, w->x); yHyperEll(B, D, sign, dist, dx, t, k1, k2, r1, r2, r1t, r2t, ht, w->y); } return; } vr_bool vroniObject::ComputeHyperbolaEllipseData(int i, e_formula *coeff) { int i1, i2, i3; t_site ltype, rtype, itype; double r1, r2; coord p1, p2, u, v, w, u_max, u_min; double A, B, C, D, t1, t2, dist, dv, dw, t_max, t_min; double r1t, r2t, dx, dy, h, ht, delta; int n1, n2, k1, k2; #ifdef TRACE double sign; #endif n1 = GetStartNode(i); n2 = GetEndNode(i); GetNodeData(n1, &u, &t1); GetNodeData(n2, &v, &t2); GetLftSiteData(i, &i1, <ype); GetRgtSiteData(i, &i2, &rtype); assert((ltype == PNT) || (ltype == ARC)); assert((rtype == PNT) || (rtype == ARC)); if (i2 > i1) { VroniSwap(i1, i2, i3); VroniSwap(ltype, rtype, itype); } else if (i1 == i2) { if (rtype == PNT) { VroniSwap(i1, i2, i3); VroniSwap(ltype, rtype, itype); } } if (t1 > t2) { u_max = u; u_min = v; t_max = t1; t_min = t2; } else { u_max = v; u_min = u; t_max = t2; t_min = t1; } if (ltype == PNT) { r1 = 0.0; k1 = 1; p1 = GetPntCoords(i1); } else { r1 = GetArcRadius(i1); p1 = GetArcCenter(i1); dw = PntPntDist(u_max, p1); if (dw > r1) k1 = 1; else k1 = -1; } if (rtype == PNT) { r2 = 0.0; k2 = 1; p2 = GetPntCoords(i2); } else { r2 = GetArcRadius(i2); p2 = GetArcCenter(i2); dw = PntPntDist(u_max, p2); if (dw > r2) k2 = 1; else k2 = -1; } /* */ /* ellipse/hyperbola formula w.r.t. clearance distance t: */ /* x(t) = A - C * t + sign * dy * sqrt(r1(t)^2 - h(t)^2) */ /* y(t) = B - D * t - sign * dx * sqrt(r1(t)^2 - h(t)^2) */ /* where dist := PntPntDist(p1, p2) */ /* dx := (p2.x - p1.x) / dist */ /* dy := (p2.y - p1.y) / dist */ /* delta := (k2 * r2 - k1 * r1) / dist */ /* h := (r2 * r2 - r1 * r1 - dist * dist) / (2 * dist) */ /* A := p1.x - dx * h */ /* B := p1.y - dy * h */ /* C := dx * delta */ /* D := dy * delta */ /* r1t := r1 + k1 * t */ /* r2t := r2 + k2 * t */ /* ht := (r2t^2 - r1t^2 - dist^2) / (2 * dist) */ /* */ /* */ dx = p2.x - p1.x; dy = p2.y - p1.y; dist = sqrt(dx * dx + dy * dy); if (eq(dist, ZERO_MAX)) { dist = r1 - r2; if (eq(dist, ZERO_MAX)) { if (t1 > t2) w = VecSub(u, v); else w = VecSub(v, u); coeff->A = w.x; coeff->B = w.y; } return false; } dx /= dist; dy /= dist; delta = (k2 * r2 - k1 * r1) / dist; h = (r2 * r2 - r1 * r1 - dist * dist) / (2.0 * dist); A = p1.x - dx * h; B = p1.y - dy * h; C = dx * delta; D = dy * delta; xHyperEll(A, C, 1.0, dist, dy, t_max, k1, k2, r1, r2, r1t, r2t, ht, v.x); yHyperEll(B, D, 1.0, dist, dx, t_max, k1, k2, r1, r2, r1t, r2t, ht, v.y); xHyperEll(A, C, -1.0, dist, dy, t_max, k1, k2, r1, r2, r1t, r2t, ht, w.x); yHyperEll(B, D, -1.0, dist, dx, t_max, k1, k2, r1, r2, r1t, r2t, ht, w.y); dv = PntPntDistSq(u_max, v); dw = PntPntDistSq(u_max, w); if (eq(dv, ZERO_MAX) && eq(dw, ZERO_MAX)) { xHyperEll(A, C, 1.0, dist, dy, t_min, k1, k2, r1, r2, r1t, r2t, ht, v.x); yHyperEll(B, D, 1.0, dist, dx, t_min, k1, k2, r1, r2, r1t, r2t, ht, v.y); xHyperEll(A, C, -1.0, dist, dy, t_min, k1, k2, r1, r2, r1t, r2t, ht, w.x); yHyperEll(B, D, -1.0, dist, dx, t_min, k1, k2, r1, r2, r1t, r2t, ht, w.y); dx = PntPntDistSq(u_min, v); dy = PntPntDistSq(u_min, w); if (eq(dx, ZERO_MAX) && eq(dy, ZERO_MAX)) { /* if ((((dx == Max(dx, dy)) && (dw == Max(dv, dw)) && (dx > ZERO * 10.0) && (dw > ZERO * 10.0)) || ((dy == Max(dx, dy)) && (dv == Max(dv, dw)) && (dy > ZERO * 10.0) && (dv > ZERO * 10.0))) && (eq(t1 - t2, ZERO * 10.0))) { if (PntPntDist(u, v) > 1e-4) { if (dist > 1e-6) { FP_printf("\ndv = %20.16f\ndw = %20.16f\n", FP_PRNTARG(sqrt(dv)), FP_PRNTARG(sqrt(dw))); FP_printf("dx = %20.16f\ndy = %20.16f\n", FP_PRNTARG(sqrt(dx)), FP_PRNTARG(sqrt(dy))); FP_printf("dt = %20.16f\n", FP_PRNTARG(Abs(t1-t2))); FP_printf("t1 = %20.16f\n", FP_PRNTARG(t1)); FP_printf("t2 = %20.16f\n", FP_PRNTARG(t2)); FP_printf("uv = %20.16f\n", FP_PRNTARG(PntPntDist(u, v))); FP_printf("di = %20.16f\n", FP_PRNTARG(dist)); FP_printf("det1 = %20.16f\n", FP_PRNTARG(VecDet(p1, p2, u))); FP_printf("det2 = %20.16f\n", FP_PRNTARG(VecDet(p1, p2, v))); } } } */ if (VroniMax(dv, dx) > VroniMax(dw, dy)) coeff->sign = false; else coeff->sign = true; } else { if (dx < dy) coeff->sign = true; else coeff->sign = false; } } else { if (dv < dw) coeff->sign = true; else coeff->sign = false; } coeff->A = A; coeff->B = B; coeff->C = C; coeff->D = D; coeff->d = dist; if (k1 == 1) coeff->k1 = true; else coeff->k1 = false; if (k2 == 1) coeff->k2 = true; else coeff->k2 = false; #ifdef TRACE if (coeff->sign) sign = 1.0; else sign = -1.0; xHyperEll(A, C, sign, dist, dy, t_max, k1, k2, r1, r2, r1t, r2t, ht, v.x); yHyperEll(B, D, sign, dist, dx, t_max, k1, k2, r1, r2, r1t, r2t, ht, v.y); xHyperEll(A, C, sign, dist, dy, t_min, k1, k2, r1, r2, r1t, r2t, ht, w.x); yHyperEll(B, D, sign, dist, dx, t_min, k1, k2, r1, r2, r1t, r2t, ht, w.y); #endif return true; } /* */ /* compute (and store) the coefficients of my parameterization formula for a */ /* hyperbolic or elliptic VD edge defined by two arcs/points. note that the */ /* parameterization uses the clearance radius as parameter. */ /* */ void vroniObject::CreateHyperbolaEllipseData(int i) { if (ComputeHyperbolaEllipseData(i, &(edge_data[i]))) { SetEdgeFlagDeg(i, false); } else { SetEdgeFlagDeg(i, true); } SetEdgeDataInit(i, false); return; } #endif vr_bool vroniObject::ComputeLineData(int i, e_formula *coeff) { coord u, v, w; double t1, t2, dist; int n1, n2; n1 = GetStartNode(i); n2 = GetEndNode(i); GetNodeData(n1, &u, &t1); GetNodeData(n2, &v, &t2); dist = t1 - t2; if (eq(dist, ZERO)) { /* */ /* degenerate Voronoi edge; most likely, this is a very short edge */ /* */ return false; } if (t1 > t2) w = VecSub(u, v); else w = VecSub(v, u); SetEdgeDataInit(i, false); SetEdgeFlagDeg(i, false); coeff->A = w.x; coeff->B = w.y; return true; } /* */ /* compute (and store) the coefficients of my parameterization formula for a */ /* linear VD edge defined by two straight-line segments. note that the */ /* parameterization uses the clearance radius as parameter. */ /* */ void vroniObject::CreateLineData(int i) { assert(InEdgeFlagList(i)); assert(InEdgeDataList(i)); if (ComputeLineData(i, &(edge_data[i]))) { SetEdgeDataInit(i, false); SetEdgeFlagDeg(i, false); } else { SetEdgeFlagDeg(i, true); SetEdgeDataInit(i, false); } return; } /* */ /* evaluate the parameterization formula for the linear VD edge i. */ /* */ void vroniObject::EvaluateLineData(int i, double t, coord *w) { int n1, n2; coord u, v; double t1, t2; assert(InEdgesList(i)); assert(!GetEdgeDataInit(i)); n1 = GetStartNode(i); n2 = GetEndNode(i); GetNodeData(n1, &u, &t1); GetNodeData(n2, &v, &t2); #ifndef NDEBUG if (t1 < t2) assert((t1 <= t) && (t <= t2)); else assert((t2 <= t) && (t <= t1)); #endif if (GetEdgeFlagDeg(i)) { /* */ /* degenerate Voronoi edge */ /* */ EvaluateDegenerateEdge(i, t, w); } else { assert(!eq(t1-t2, ZERO)); if (t1 > t2) { t = (t - t2) / (t1 - t2); u = v; } else { t = (t - t1) / (t2 - t1); } v.x = edge_data[i].A; v.y = edge_data[i].B; *w = RayPnt(u, v, t); } return; } #ifdef MAT void vroniObject::OutputMAEdges(FILE *output, int e, int n) { int e1, edge_class, parent_id; t_site l_type, r_type; double t; coord p; while (GetEdgeFlagNew(e)) { parent_id = n; SetEdgeFlagNew(e, false); n = GetOtherNode(e, n); t = GetNodeParam(n); if (IsWmatEdge(e)) { p = GetNodeCoord(n); if (IsStartNode(e, n)) { l_type = GetLftSiteType(e); r_type = GetRgtSiteType(e); } else { l_type = GetRgtSiteType(e); r_type = GetLftSiteType(e); } if (l_type == SEG) { if (r_type == SEG) edge_class = 2; else if (r_type == PNT) edge_class = 4; else { edge_class = 0; VD_Warning("OutputMAEdges() - this is not a polygonal input!"); } } else if (l_type == PNT) { if (r_type == SEG) edge_class = 3; else if (r_type == PNT) edge_class = 1; else { edge_class = 0; VD_Warning("OutputMAEdges() - this is not a polygonal input!"); } } else { edge_class = 0; VD_Warning("OutputMAEdges() - this is not a polygonal input!"); } FP_fprintf(output, "%7d %20.16f %20.16f %20.16f %7d %1d\n", n, FP_PRNTARG(UnscaleX(p.x)), FP_PRNTARG(UnscaleY(p.y)), FP_PRNTARG(UnscaleV(t)), parent_id, edge_class); if (gt(t, ZERO)) { e1 = GetCCWEdge(e, n); if (GetEdgeFlagNew(e1) && IsWmatEdge(e1)) { e = GetCWEdge(e, n); if (GetEdgeFlagNew(e) && (e1 != e) && IsWmatEdge(e)) { Push(e); Push(n); } e = e1; } else { e = GetCWEdge(e, n); while (!GetEdgeFlagNew(e) && !IsWmatEdge(e) && StackIsNotEmpty) { Pop(n); Pop(e); } } } else if (StackIsNotEmpty) { do { Pop(n); Pop(e); } while (!GetEdgeFlagNew(e) && !IsWmatEdge(e) && StackIsNotEmpty); } } else if (StackIsNotEmpty) { do { Pop(n); Pop(e); } while (!GetEdgeFlagNew(e) && !IsWmatEdge(e) && StackIsNotEmpty); } } return; } void vroniObject::OutputMA(char output_file[]) { int n = 0, e, e1, e2, num_edges; vr_bool at_boundary = true; coord p; double t = 0.0; FILE *output; assert(scale_factor > 0.0); if (!(GetVDStatus() || GetWMATStatus())) { VD_Warning("OutputMA() - medial axis not known!\n"); return; } ResetStack; num_edges = GetNumberOfEdges(); for (e = 0; e < num_edges; ++e) { SetEdgeFlagNew(e, true); } /* */ /* find a MA node that is not on the boundary */ /* */ e = 0; while (at_boundary && (e < num_edges)) { if (IsWmatEdge(e)) { n = GetStartNode(e); t = GetNodeParam(n); if (gt(t, ZERO)) { at_boundary = false; } else { n = GetOtherNode(e, n); t = GetNodeParam(n); if (gt(t, ZERO)) { at_boundary = false; } else { ++e; } } } else { ++e; } } if (at_boundary) { VD_Warning("OutputMA() - medial axis not known or corrupted data!\n"); return; } output = OpenFileVD(output_file, "w"); p = GetNodeCoord(n); FP_fprintf(output, "%7d %20.16f %20.16f %20.16f %7d %1d\n", n, FP_PRNTARG(UnscaleX(p.x)), FP_PRNTARG(UnscaleY(p.y)), FP_PRNTARG(UnscaleV(t)), -1, 9); OutputMAEdges(output, e, n); e1 = GetCCWEdge(e, n); if ((e1 >= 4) && IsWmatEdge(e1)) OutputMAEdges(output, e1, n); e2 = GetCWEdge(e, n); if ((e1 != e2) && (e2 >= 4) && IsWmatEdge(e2)) OutputMAEdges(output, e2, n); fclose(output); return; } #endif #ifdef WRITE_VD #define INCR 0.02 inline void vroniObject::StoreActivePnt(coord3D **active_pnts, int *num_active_pnts, int *max_num_active_pnts, int i, vr_bool isPnt) { coord p; vr_bool isNew; isNew = false; if (isPnt) { // p = GetPntCoords(i); //printf("pnt i = %d with x = %f and y = %f\n", i, // UnscaleX(p.x), UnscaleY(p.y)); //printf("idx = %d\n", GetPntIndex(i)); if (GetPntIndex(i) == NIL) { ++(*num_active_pnts); SetPntIndex(i, *num_active_pnts); p = GetPntCoords(i); isNew = true; } } else { //p = GetNodeCoord(i); //printf("nde i = %d with x = %f and y = %f\n", i, // UnscaleX(p.x), UnscaleY(p.y)); //printf("idx = %d\n", GetNodeIndex(i)); if (GetNodeIndex(i) == NIL) { ++(*num_active_pnts); SetNodeIndex(i, *num_active_pnts); p = GetNodeCoord(i); isNew = true; } } if (isNew) { if (*num_active_pnts >= *max_num_active_pnts) { *max_num_active_pnts *= 2; (*active_pnts) = (coord3D*) ReallocateArray(*active_pnts, *max_num_active_pnts, sizeof(coord3D), "data_off:active_pnts"); } (*active_pnts)[*num_active_pnts].x = p.x; (*active_pnts)[*num_active_pnts].y = p.y; if (isPnt) { (*active_pnts)[*num_active_pnts].z = 0.0; } else { (*active_pnts)[*num_active_pnts].z = GetNodeParam(i); } } return; } void vroniObject::RecApxVDEdge(int e, double delta, coord3D **active_pnts, int *num_active_pnts, int *max_num_active_pnts, coord p1, coord p2, double t1, double t2) { double t0; coord p0; int n0; t0 = (t1 + t2) / 2.0; EvaluateBisectorData(e, t0, &p0); n0 = StoreNode(p0.x, p0.y, t0); if (PntPntDist(p0, p2) > delta) RecApxVDEdge(e, delta, active_pnts, num_active_pnts, max_num_active_pnts, p0, p2, t0, t2); InsertDummyNode(e, n0, false); StoreActivePnt(active_pnts, num_active_pnts, max_num_active_pnts, n0, false); if (PntPntDist(p1, p0) > delta) RecApxVDEdge(e, delta, active_pnts, num_active_pnts, max_num_active_pnts, p1, p0, t1, t0); return; } void vroniObject::ApproximateVDEdge(int e, double delta, coord3D **active_pnts, int *num_active_pnts, int *max_num_active_pnts) { int n1, n2, i1, i2; coord p1, p2; double t1, t2; t_site ltype, rtype; assert(delta > 0.0); n1 = GetStartNode(e); n2 = GetEndNode(e); GetNodeData(n1, &p1, &t1); GetNodeData(n2, &p2, &t2); if (t1 > 0.0) StoreActivePnt(active_pnts, num_active_pnts, max_num_active_pnts, n1, false); if (t2 > 0.0) StoreActivePnt(active_pnts, num_active_pnts, max_num_active_pnts, n2, false); GetLftSiteData(e, &i1, <ype); GetRgtSiteData(e, &i2, &rtype); if (((ltype == SEG) && (rtype == PNT)) || ((rtype == SEG) && (ltype == PNT))) { if ((t1 > 0.0) && (t2 > 0.0)) { /* */ /* parabolic bisector which may need to be approximated */ /* */ if (PntPntDist(p1, p2) > delta) RecApxVDEdge(e, delta, active_pnts, num_active_pnts, max_num_active_pnts, p1, p2, t1, t2); } } else if ((ltype == PNT) && (rtype == PNT)) { if ((t1 > 0.0) && (t2 > 0.0)) { /* */ /* hyperbolic bisector which may need to be approximated */ /* */ if (PntPntDist(p1, p2) > delta) RecApxVDEdge(e, delta, active_pnts, num_active_pnts, max_num_active_pnts, p1, p2, t1, t2); } } return; } void vroniObject::ScanVDCell(int e, int i, t_site type, FILE *fp) { int i1, n1, n2; t_site t; /* set fan_tri to TRUE in order to triangulate the Voronoi region of a */ /* reflex vertex */ vr_bool fan_tri = false; GetLftSiteData(e, &i1, &t); if ((i == i1) && (type == t)) { n1 = GetStartNode(e); n2 = GetEndNode(e); } else { GetRgtSiteData(e, &i1, &t); assert((i == i1) && (type == t)); n2 = GetStartNode(e); n1 = GetEndNode(e); } if (type == PNT) { if (fan_tri) { i1 = GetPntIndex(i); n1 = n2; e = GetCWEdge(e, n1); n2 = GetOtherNode(e, n1); while (GetNodeParam(n2) > 0.0) { if (GetNodeIndex(n1) != GetNodeIndex(n2)) fprintf(fp, "f %d %d %d", i1, GetNodeIndex(n1), GetNodeIndex(n2)); n1 = n2; e = GetCWEdge(e, n1); n2 = GetOtherNode(e, n1); } } else { fprintf(fp, "f %d ", GetPntIndex(i)); } } else { assert(type == SEG); if (GetNodeSite(n1) == Get1stVtx(i, SEG)) { fprintf(fp, "f %d ", GetPntIndex(Get2ndVtx(i, SEG))); fprintf(fp, "%d ", GetPntIndex(Get1stVtx(i, SEG))); } else { assert(GetNodeSite(n1) == Get2ndVtx(i, SEG)); fprintf(fp, "f %d ", GetPntIndex(Get1stVtx(i, SEG))); fprintf(fp, "%d ", GetPntIndex(Get2ndVtx(i, SEG))); } } while (GetNodeParam(n2) > 0.0) { if (GetNodeIndex(n1) != GetNodeIndex(n2)) fprintf(fp, "%d ", GetNodeIndex(n2)); n1 = n2; e = GetCWEdge(e, n1); n2 = GetOtherNode(e, n1); } fprintf(fp, "\n"); return; } /* */ /* This function writes a polygonal approximation of all "interior" Voronoi */ /* regions of a consistently oriented multiply-connected (polygonal) area in */ /* OBJ format to a file. Each Voronoi region is represented as one OBJ face, */ /* where every vertex of a face has a clearance radius as its z-coordinate. */ /* Hence, effectively this function outputs a "Voronoi mountain". */ /* */ /* Note: My indices start at 0, whereas OBJ's indices start at 1! */ /* */ /* Note: Currently only polygonal input and parabolic arcs are handled! */ /* */ void vroniObject::OutputVD(char output_file[], double vd_apx_dist, vr_bool left_vd, vr_bool right_vd) { FILE *fp; int n1, n2, i, j, e, i1, not_needed; int num_active_pnts, max_num_active_pnts; t_site type; coord3D *active_pnts = (coord3D*)NULL; double delta; coord p1, p2; if (!left_vd && !right_vd) { VD_Warning("OutputVD() - no side requested for VD output\n"); return; } else if (left_vd && right_vd) { VD_Warning("OutputVD() - only one side can be requested for VD output\n"); return; } if (num_arcs > 0) { VD_Warning("OutputVD() - no arcs supported yet\n"); return; } VD_Info("...computing a polygonal approximation of VD -- "); if (!GetVDStatus()) { VD_Warning("OutputVD() - Voronoi diagram not known!\n"); return; } if (GetNodesSqdFlag()) TakeSquareOfNodes(); if (!GetDeg2Flag()) InsertDegreeTwoNodes(); if (le(vd_apx_dist, ZERO)) { delta = INCR / sqrt(sqrt(sqrt(((double) num_pnts)))); } else { delta = ScaleV(vd_apx_dist); } /* */ /* allocate memory and initialize the flags for the edge data. */ /* */ InitializeEdgeData(false, left_vd, false, ¬_needed); /* */ /* allocate memory for the active vertices and nodes */ /* */ num_active_pnts = 0; max_num_active_pnts = num_pnts + num_active_edges + 1; active_pnts = (coord3D*) ReallocateArray(active_pnts, max_num_active_pnts, sizeof(coord3D), "data_off:active_pnts"); /* */ /* discard all zero-length VD edges */ /* */ j = 0; while (j < num_active_edges) { e = active_edges[j]; ++j; n1 = GetStartNode(e); n2 = GetEndNode(e); p1 = GetNodeCoord(n1); p2 = GetNodeCoord(n2); if ((GetNodeParam(n1) > 0.0) || (GetNodeParam(n2) > 0.0)) { if (PntPntDist(p1, p2) <= ZERO_MAX) { if (GetNodeIndex(n1) != NIL) { SetNodeIndex(n2, GetNodeIndex(n1)); } else { StoreActivePnt(&active_pnts, &num_active_pnts, &max_num_active_pnts, n2, false); SetNodeIndex(n1, GetNodeIndex(n2)); } } } } /* */ /* approximate all conic VD arcs */ /* */ j = 0; while (j < num_active_edges) { e = active_edges[j]; ++j; n1 = GetStartNode(e); n2 = GetEndNode(e); GetLftSiteData(e, &i, &type); assert(type != ARC); if ((GetNodeParam(n1) == 0.0) && (GetNodeParam(n2) > 0.0)) SetVDEdge(i, type, e); if (type == PNT) { // printf("pnt %d\n", i); StoreActivePnt(&active_pnts, &num_active_pnts, &max_num_active_pnts, i, true); } else { assert(type == SEG); // printf("seg %d\n", i); i1 = Get1stVtx(i, type); StoreActivePnt(&active_pnts, &num_active_pnts, &max_num_active_pnts, i1, true); i1 = Get2ndVtx(i, type); StoreActivePnt(&active_pnts, &num_active_pnts, &max_num_active_pnts, i1, true); } GetRgtSiteData(e, &i, &type); assert(type != ARC); if ((GetNodeParam(n2) == 0.0) && (GetNodeParam(n1) > 0.0)) SetVDEdge(i, type, e); if (type == PNT) { // printf("\nVDEdge(%d) = %d\n", i, GetVDEdge(i, PNT)); // printf("pnt %d for edge %d\n", i, e); // printf("VDEdge(%d) = %d\n", i, GetVDEdge(i, PNT)); StoreActivePnt(&active_pnts, &num_active_pnts, &max_num_active_pnts, i, true); } else { assert(type == SEG); // printf("seg %d\n", i); i1 = Get1stVtx(i, type); StoreActivePnt(&active_pnts, &num_active_pnts, &max_num_active_pnts, i1, true); i1 = Get2ndVtx(i, type); StoreActivePnt(&active_pnts, &num_active_pnts, &max_num_active_pnts, i1, true); } ApproximateVDEdge(e, delta, &active_pnts, &num_active_pnts, &max_num_active_pnts); } fp = OpenFileVD(output_file, "w"); fprintf(fp, "# .obj file generated by approximating all Voronoi regions "); fprintf(fp, "computed by\n"); fprintf(fp, "# Martin Held's VRONI. "); fprintf(fp, "email address: held@cs.sbg.ac.at.\n"); fprintf(fp, "# every z-coordinate corresponds to a clearance radius.\n"); /* */ /* output all active points */ /* */ fprintf(fp, "\n# %d vertices\n", num_active_pnts); for (i = 1; i <= num_active_pnts; ++i) { fprintf(fp, "v %20.16f %20.16f %20.16f\n", UnscaleX(active_pnts[i].x), UnscaleY(active_pnts[i].y), UnscaleV(active_pnts[i].z)); } /* */ /* output all VD faces */ /* */ j = 0; while (j < num_active_edges) { e = active_edges[j]; ++j; GetLftSiteData(e, &i, &type); if (GetVDEdge(i, type) == e) ScanVDCell(e, i, type, fp); GetRgtSiteData(e, &i, &type); if (GetVDEdge(i, type) == e) ScanVDCell(e, i, type, fp); } fclose(fp); free(active_pnts); VD_Info(" done!\n"); return; } #endif #include "ext_appl_data_off.cc"