/*****************************************************************************/ /* */ /* Copyright (C) 2007-2025 S. Huber, 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: Stefan Huber, Martin Held */ /* */ /* E-Mail: held@cs.sbg.ac.at */ /* Fax Mail: (+43 662) 8044-611 */ /* Voice Mail: (+43 662) 8044-6304 */ /* Snail Mail: Martin Held */ /* FB Informatik */ /* Universitaet Salzburg */ /* A-5020 Salzburg, Austria */ /* */ /*****************************************************************************/ /* */ /* get standard libraries */ /* */ #include #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 "roots.h" /* */ /* Calculate the Voronoi center of arc and two segs */ /* */ vr_bool vroniObject::ArcSegSegCntr(int i, int j, int k, int e, coord* cntr, double *r2, vr_bool *problematic, int *site) { int m; coord c1; double rr1; double s1a, s1b, s1c; double s2a, s2b, s2c; coord centers[VRONI_MAXSOL]; double radii[VRONI_MAXSOL]; int num_sol = 0, best_sol = NIL; double eps = ZERO; coord v, w, sp, ep, spi, epi, spj, epj; vr_bool counter_tangential; #ifdef VRONI_DBG_WARN coord spk, epk; #endif /* */ /* Check whether the three cites are all in lists */ /* */ assert(InArcsList(i)); assert(InSegsList(j)); assert(InSegsList(k)); *problematic = false; /* */ /* we sort the two seg indices in order to guarantee that the center of */ /* these three sites is always computed consistently. */ /* */ SortTwoNumbers(j, k, m); /* */ /* check whether the three sites share a common end point */ /* */ m = Get1stVtx(k, SEG); if (IsArcStartPnt(i, m) || IsArcEndPnt(i, m)) { if (IsSegStartPnt(j, m) || IsSegEndPnt(j, m)) { *cntr = GetPntCoords(m); *r2 = 0.0; *site = m; return true; /* common end point */ } } m = Get2ndVtx(k, SEG); if (IsArcStartPnt(i, m) || IsArcEndPnt(i, m)) { if (IsSegStartPnt(j, m) || IsSegEndPnt(j, m)) { *cntr = GetPntCoords(m); *r2 = 0.0; *site = m; return true; /* common end point */ } } c1 = GetArcCenter(i); rr1 = GetArcRadius(i); GetSegEqnData(j, &s1a, &s1b, &s1c); s1c*= -1.0; GetSegEqnData(k, &s2a, &s2b, &s2c); s2c*= -1.0; if (((s1a < -ZERO) && (s2a > ZERO)) || ((s1a > ZERO) && (s2a < ZERO))) { s2a *= -1.0; s2b *= -1.0; s2c *= -1.0; } else if (((s1b < -ZERO) && (s2b > ZERO)) || ((s1b > ZERO) && (s2b < ZERO))) { s2a *= -1.0; s2b *= -1.0; s2c *= -1.0; } while (eps <= ZERO_MAX) { num_sol = 0; //If second segment is tangential to arc --> let it be segment j if ((Abs(Abs(s2a*c1.x+s2b*c1.y-s2c)-rr1) <= eps) && ((PntPntDist(GetArcStartCoord(i), GetSegStartCoord(k)) <= eps) || (PntPntDist(GetArcStartCoord(i), GetSegEndCoord(k)) <= eps) || (PntPntDist(GetArcEndCoord(i), GetSegStartCoord(k)) <= eps) || (PntPntDist(GetArcEndCoord(i), GetSegEndCoord(k)) <= eps)) && !(IsPntInArcCone(i,GetSegStartCoord(k)) && IsPntInArcCone(i,GetSegEndCoord(k)))) { int tmpi; double tmpd; VroniSwap(j, k, tmpi); VroniSwap(s1a, s2a, tmpd); VroniSwap(s1b, s2b, tmpd); VroniSwap(s1c, s2c, tmpd); } spi = GetArcStartCoord(i); epi = GetArcEndCoord(i); spj = GetSegStartCoord(j); epj = GetSegEndCoord(j); #ifdef TRACE if ((i == 1305) && (((j == 1605) && (k == 1604)) || ((j == 1604) && (k == 1605)))) { printf("arc%d-seg%d-seg%d: eps = %20.16f\n", i, j, k, eps); printf("start arc: (%20.16f %20.16f)\n", spi.x, spi.y); printf("end arc: (%20.16f %20.16f)\n", epi.x, epi.y); printf("center arc: (%20.16f %20.16f)\n", c1.x, c1.y); printf("start seg 1: (%20.16f %20.16f)\n", spj.x, spj.y); printf("end seg 1: (%20.16f %20.16f)\n", epj.x, epj.y); printf("start seg 2: (%20.16f %20.16f)\n", GetSegStartCoord(k).x, GetSegStartCoord(k).y); printf("end seg 2: (%20.16f %20.16f)\n", GetSegEndCoord(k).x, GetSegEndCoord(k).y); } #endif //Segment tangential to arc! //Do not change the eps to eps_MAX here! See arcs_347 if (Abs(Abs(s1a*c1.x+s1b*c1.y-s1c)-rr1) <= eps) { counter_tangential = true; if (PntPntDist(spi,spj) <= eps) { ep = spj; sp = epj; v = VecSub(sp, ep); w = VecSub(epi, ep); counter_tangential = (VecDotProd(v, w) < 0.0) ? false : true; } else if (PntPntDist(epi,spj) <= eps) { ep = spj; sp = epj; v = VecSub(sp, ep); w = VecSub(spi, ep); counter_tangential = (VecDotProd(v, w) < 0.0) ? false : true; } else if (PntPntDist(spi,epj) <= eps) { ep = epj; sp = spj; v = VecSub(sp, ep); w = VecSub(epi, ep); counter_tangential = (VecDotProd(v, w) < 0.0) ? false : true; } else if (PntPntDist(epi,epj) <= eps) { ep = epj; sp = spj; v = VecSub(sp, ep); w = VecSub(spi, ep); counter_tangential = (VecDotProd(v, w) < 0.0) ? false : true; } if (!counter_tangential) { //Intersect the ray and get all centers v.x = s1a; v.y = s1b; num_sol = IntersectRaySegment(ep, v, s2a, s2b, s2c, centers, radii, eps); } else { num_sol = CircSegSegCenters(c1, rr1, s1a,s1b,s1c, s2a,s2b,s2c, centers, radii, eps); } } //parallel segments with common endpoint else if ((Abs(s1a-s2a) <= eps) && (Abs(s1b-s2b) <= eps) && (Abs(s1c-s2c) <= eps) && ((PntPntDist(spj,GetSegStartCoord(k)) <= eps) || (PntPntDist(spj,GetSegEndCoord(k)) <= eps) || (PntPntDist(epj,GetSegStartCoord(k)) <= eps) || (PntPntDist(epj,GetSegEndCoord(k)) <= eps))) { //Get the common endpoint if( PntPntDist(spj, GetSegStartCoord(k))<= eps || PntPntDist(spj, GetSegEndCoord(k))<= eps ) ep = spj; else ep = epj; v = MakeVec(s1a, s1b); num_sol = IntersectRayCircle(ep, v, c1, rr1, centers, radii, eps); } else { //Calculate all equidistant points num_sol = CircSegSegCenters(c1, rr1, s1a,s1b,s1c, s2a,s2b,s2c, centers, radii, eps); } /* */ /* classify all solutions and pick the best solution */ /* */ best_sol = ClassifySolutions(i, j, k, e, ARC, SEG, SEG, false, true, true, centers, radii, &num_sol, eps, problematic); assert((0 <= best_sol) && (best_sol < num_sol)); #ifdef TRACE if ((i == 40) && (((j == 2) && (k == 42)) || ((j == 42) && (k == 2)))) { printf("arc %d, seg %d, seg %d\n", i, j, k); for (m = 0; m < num_sol; ++m) { printf("centers[%d]: (%20.16f %20.16f)\n", m, centers[m].x, centers[m].y); printf("radii[%d]: %20.16f\n", m, radii[m]); } printf("problematic = %d, eps = %20.16f\n", *problematic, eps); } #endif if (!*problematic) { *cntr = centers[best_sol]; *r2 = radii[best_sol]; return true; } else { eps *= 10.0; } } if (eq(rr1, ZERO_MAX)) { /* */ /* this is an arc with a very small radius; we replace it */ /* by a seg */ /* */ VD_Dbg_Warning("ArcSegSegCntr() - arc with small radius!"); ReplaceArc(i); restart = true; return false; } else if (eq(PntPntDist(GetArcStartCoord(i), GetArcEndCoord(i)), ZERO_MAX)) { /* */ /* this is an arc with a very short cord; we replace it by a seg */ /* */ VD_Dbg_Warning("ArcArcArcCntr() - arc with short cord!"); ReplaceArc(i); restart = true; return false; } else if (eq(GetSegLgth(j), ZERO_MAX)) { /* */ /* this is a seg with a very small length; we replace it */ /* */ VD_Dbg_Warning("ArcSegSegCntr() - seg with small length!"); ReplaceSeg(j); restart = true; return false; } else if (eq(GetSegLgth(k), ZERO_MAX)) { /* */ /* this is a seg with a very small length; we replace it */ /* */ VD_Dbg_Warning("ArcSegSegCntr() - seg with small length!"); ReplaceSeg(k); restart = true; return false; } else if (CheckIntersectionsLocally(i, ARC, j, SEG, k, SEG)) { restart = true; return false; } #ifdef VRONI_DBG_WARN spi = GetArcStartCoord(i); epi = GetArcEndCoord(i); c1 = GetArcCenter(i); printf("\nCenter not computed for arc-seg-seg: %d-%d-%d\n", i, j, k); printf("%20.16f %20.16f %20.16f %20.16f %20.16f %20.16f\n", spi.x, spi.y, epi.x, epi.y, c1.x, c1.y); epj = GetSegStartCoord(j); spj = GetSegEndCoord(j); epk = GetSegStartCoord(k); spk = GetSegEndCoord(k); printf("%20.16f %20.16f %20.16f %20.16f\n", spj.x, spj.y, epj.x, epj.y); printf("%20.16f %20.16f %20.16f %20.16f\n", spk.x, spk.y, epk.x, epk.y); #endif *cntr = centers[best_sol]; *r2 = radii[best_sol]; restart = false; return false; }