/*****************************************************************************/ /* */ /* Copyright (C) 1999-2025 M. Held */ /* */ /* This code is not in the public domain. All rights reserved! Please make */ /* sure to read the full copyright statement contained in "README.txt" or in */ /* the "main" file of this code, such as "main.cc". */ /* */ /*****************************************************************************/ /* */ /* Written by: Martin Held */ /* */ /* E-Mail: held@cs.sbg.ac.at */ /* Fax Mail: (+43 662) 8044-611 */ /* Voice Mail: (+43 662) 8044-6304 */ /* Snail Mail: Martin Held */ /* FB Informatik */ /* Universitaet Salzburg */ /* A-5020 Salzburg, Austria */ /* */ /*****************************************************************************/ /* */ /* get standard libraries */ /* */ #include #include #include #include #include /* */ /* get my header files */ /* */ #include "fpkernel.h" #include "vroni_object.h" #include "roots.h" #include "numerics.h" #include "geom.h" #define SAMPLE 100 vr_bool vroniObject::DesperatePntPntPntCntr(const coord & A, const coord & B, const coord & C, coord *cntr, double *r2) { int i, j = NIL; coord AB, BC, CA; coord u, v, w, p, dir, mid; double ab, bc, ca; double s = 0.0, t_min, t_max, delta; double t[SAMPLE+1]; double dist, new_dist = LARGE, old_dist = LARGE; AB = VecSub(B, A); ab = VecLen(AB); assert(ab >= 0.0); BC = VecSub(C, B); bc = VecLen(BC); assert(bc >= 0.0); CA = VecSub(A, C); ca = VecLen(CA); assert(ca >= 0.0); if (ab <= ca) { if (ca <= bc) { /* ab, ca <= bc */ u = B; v = C; w = A; dir = VecDiv(bc, BC); } else { /* ab, bc <= ca */ u = A; v = C; w = B; dir = VecDiv(ca, CA); } } else { if (ab <= bc) { /* ab, ca <= bc */ u = B; v = C; w = A; dir = VecDiv(bc, BC); } else { /* ca, bc <= ab */ u = A; v = B; w = C; dir = VecDiv(ab, AB); } } dir = VecCCW(dir); mid = MidPoint(u, v); t_min = - too_large; t_max = too_large; #ifdef VRONI_DBG_WARNINGS printf("DesperatePntPntPntCntr:\n"); printf("t_min = %20.16f\n", t_min); printf("t_max = %20.16f\n", t_max); printf("A = (%20.16f, %20.16f)\n", A.x, A.y); printf("B = (%20.16f, %20.16f)\n", B.x, B.y); printf("C = (%20.16f, %20.16f)\n", C.x, C.y); printf("dir = (%20.16f, %20.16f)\n", dir.x, dir.y); printf("mid = (%20.16f, %20.16f)\n", mid.x, mid.y); #endif do { old_dist = new_dist; delta = (t_max - t_min) / SAMPLE; t[0] = t_min; t[SAMPLE] = t_max; for (i = 1; i < SAMPLE; ++i) t[i] = t_min + i * delta; for (i = 0; i <= SAMPLE; ++i) { p.x = mid.x + t[i] * dir.x; p.y = mid.y + t[i] * dir.y; dist = PntPntDist(u, p) - PntPntDist(w, p); if (dist < 0.0) dist = -dist; if (dist < new_dist) { new_dist = dist; j = i; s = t[i]; } } if (new_dist < old_dist) { if (j == 0) { t_min = t[0]; t_max = t[2]; } else if (j == SAMPLE) { t_min = t[SAMPLE-2]; t_max = t[SAMPLE]; } else { t_min = t[j-1]; t_max = t[j+1]; } } #ifdef VRONI_DBG_WARNINGS printf("t_min = %20.16f\n", t_min); printf("t_max = %20.16f\n", t_max); printf("s = %20.16f\n", s); printf("dist = %20.16f\n", new_dist); #endif } while (new_dist < old_dist); p.x = mid.x + s * dir.x; p.y = mid.y + s * dir.y; *r2 = (PntPntDist(u, p) + PntPntDist(v, p) + PntPntDist(w, p)) / 3.0; *cntr = p; #ifdef VRONI_DBG_WARNINGS printf("cntr = (%20.16f, %20.16f)\n", cntr->x, cntr->y); printf("rad = %20.16f\n", *r2); #endif *r2 = *r2 * *r2; return (j != NIL); } vr_bool vroniObject::NumericalPntPntPntCntr(const coord & A, const coord & B, const coord & C, coord *cntr, double *r2) { double alpha, beta, gamma, delta, t; coord AB, BC, CA, P, Q, U, V; double aaa[2][2], bbb[2], xy[2]; int I0, I1, J0, J1, exists; double ab, bc, ca; vr_bool success; /* */ /* check whether two points are nearly identical. we return a circle */ /* whose diameter is given by those points which do not coincide. */ /* */ AB = VecSub(B, A); ab = VecLen(AB); assert(ab >= 0.0); if (le(ab, ZERO_MAX)) { AB = MidPoint(A, B); *cntr = MidPoint(AB, C); *r2 = PntPntDist(AB, C) / 2.0; return true; } BC = VecSub(C, B); bc = VecLen(BC); assert(bc >= 0.0); if (le(bc, ZERO)) { BC = MidPoint(B, C); *cntr = MidPoint(BC, A); *r2 = PntPntDist(BC, A) / 2.0; return true; } CA = VecSub(A, C); ca = VecLen(CA); assert(ca >= 0.0); if (le(ca, ZERO)) { CA = MidPoint(C, A); *cntr = MidPoint(CA, B); *r2 = PntPntDist(CA, B) / 2.0; return true; } /* */ /* normalize the vectors */ /* */ AB = VecDiv(ab, AB); BC = VecDiv(bc, BC); CA = VecDiv(ca, CA); /* */ /* find those two vectors that are the least parallel. */ /* */ alpha = VecDotProd(AB, BC); alpha = Abs(alpha); beta = VecDotProd(BC, CA); beta = Abs(beta); gamma = VecDotProd(CA, AB); gamma = Abs(gamma); delta = Min3(alpha, beta, gamma); /* */ /* assume that the angle between AB and AC is closest to a right */ /* angle. if ab > ca then we construct the line perpendicular to AB */ /* through the midpoint of A,B and intersect it with the corresponding */ /* line w.r.t. A,C. the point of intersection will be expressed in */ /* terms of the line perpendicular to AB. similar for the other cases. */ /* */ if (delta == gamma) { P = MidPoint(A, B); Q = VecCCW(AB); U = MidPoint(A, C); V = VecCCW(CA); if (ab < ca) { LineLineIntersection(aaa, bbb, xy, I0, J0, I1, J1, P, Q, U, V, *cntr, exists); } else { LineLineIntersection(aaa, bbb, xy, I0, J0, I1, J1, U, V, P, Q, *cntr, exists); } } else if (delta == alpha) { P = MidPoint(A, B); Q = VecCCW(AB); U = MidPoint(B, C); V = VecCCW(BC); if (ab < bc) { LineLineIntersection(aaa, bbb, xy, I0, J0, I1, J1, P, Q, U, V, *cntr, exists); } else { LineLineIntersection(aaa, bbb, xy, I0, J0, I1, J1, U, V, P, Q, *cntr, exists); } } else { P = MidPoint(C, A); Q = VecCCW(CA); U = MidPoint(B, C); V = VecCCW(BC); if (ca < bc) { LineLineIntersection(aaa, bbb, xy, I0, J0, I1, J1, P, Q, U, V, *cntr, exists); } else { LineLineIntersection(aaa, bbb, xy, I0, J0, I1, J1, U, V, P, Q, *cntr, exists); } } if (exists == 1) { *r2 = PntPntDist(A, *cntr); return true; } else { /* */ /* we'll try to find a center using an alternative approach: we find */ /* a parameter t such that cntr = P + t * Q. thus, we require */ /* d(A, cntr) = d(C, cntr). */ /* */ P = MidPoint(A, B); Q = VecCCW(AB); /* note: Q is a unit vector! */ V = VecSub(P, C); delta = VecDotProd(Q, V) * 2.0; if (ne(delta, ZERO)) { ab = ab * ab / 4.0; t = (ab - VecLenSq(V)) / delta; *cntr = RayPnt(P, Q, t); *r2 = ab + t * t; *r2 = sqrt(*r2); return true; } else { success = DesperatePntPntPntCntr(A, B, C, cntr, r2); if (success) *r2 = sqrt(*r2); return success; } } } #define NS 10 void vroniObject::DesperateComputeCenter(int e, int i, t_site ti, coord *w, double *r2) { int n, j, k, num_samples[3]; int i_sites[3], min_ind[3], new_min, new_max; t_site s_types[3]; coord p, u, v, cntr; coord samples[3][NS+1]; coord start[3], end[3]; double t, angle_s, angle_e, angle, du, dv, radius; double dist, dist1, dist2, dist3, dist_diff = LARGE, old_diff; #ifdef VR_TRACE int counter = 1; #endif #ifdef VR_TRACE #endif #ifdef VR_TRACE printf("in DesperateComputeCenter()\n"); #endif GetLftSiteData(e, &i_sites[0], &s_types[0]); GetRgtSiteData(e, &i_sites[1], &s_types[1]); i_sites[2] = i; s_types[2] = ti; *r2 = LARGE; #ifdef VR_TRACE printf("site %d/%d, site %d/%d, site %d/%d\n", i, ti, i_sites[0], s_types[0], i_sites[1], s_types[1]); int i_sites[3], min_ind[3], new_min, new_max; t_site s_types[3]; #endif for (j = 0; j < 3; j += 1) { if (s_types[j] == PNT) { start[j] = end[j] = GetPntCoords(i_sites[j]); } else if (s_types[j] == SEG) { start[j] = GetSegStartCoord(i_sites[j]); end[j] = GetSegEndCoord(i_sites[j]); } else if (s_types[j] == ARC) { start[j] = GetArcStartCoord(i_sites[j]); end[j] = GetArcEndCoord(i_sites[j]); } } do { #ifdef VR_TRACE printf("s1: (%20.16f,%20.16f)\ne1: (%20.16f,%20.16f)\n", start[0].x, start[0].y, end[0].x, end[0].y); printf("s2: (%20.16f,%20.16f)\ne2: (%20.16f,%20.16f)\n", start[1].x, start[1].y, end[1].x, end[1].y); printf("s3: (%20.16f,%20.16f)\ne3: (%20.16f,%20.16f)\n\n", start[2].x, start[2].y, end[2].x, end[2].y); #endif old_diff = dist_diff; for (j = 0; j < 3; j += 1) { if (s_types[j] == PNT) { samples[j][0] = start[j]; num_samples[j] = 1; } else if (s_types[j] == SEG) { num_samples[j] = NS + 1; samples[j][0] = u = start[j]; samples[j][NS] = v = end[j]; for (k = 1; k <= NS-1; ++k) { t = ((double) k) / NS; samples[j][k] = LinearComb(u, v, t); } } else if (s_types[j] == ARC) { num_samples[j] = NS + 1; samples[j][0] = u = start[j]; samples[j][NS] = v = end[j]; p = GetArcCenter(i_sites[j]); 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; du = PntPntDist(u, p); dv = PntPntDist(v, p); for (k = 1; k <= NS-1; ++k) { t = ((double) k) / NS; angle = (1.0 - t) * angle_s + t * angle_e; radius = (1.0 - t) * du + t * dv; samples[j][k].x = p.x + radius * cos(angle); samples[j][k].y = p.y + radius * sin(angle); } } } for (n = 0; n < num_samples[0]; ++n) { for (j = 0; j < num_samples[1]; ++j) { for (k = 0; k < num_samples[2]; ++k) { if (NumericalPntPntPntCntr(samples[0][n], samples[1][j], samples[2][k], &cntr, &radius)) { dist1 = ComputeSiteDistance(i_sites[0], s_types[0], cntr); dist2 = ComputeSiteDistance(i_sites[1], s_types[1], cntr); dist3 = ComputeSiteDistance(i_sites[2], s_types[2], cntr); dist = Abs(dist1 - radius) + Abs(dist2 - radius) + Abs(dist3 - radius); if (dist < dist_diff) { dist_diff = dist; min_ind[0] = n; min_ind[1] = j; min_ind[2] = k; *r2 = (dist1 + dist2 + dist3) / 3.0; *w = cntr; } } } } } #ifdef VR_TRACE printf("new radius after %2d round: %20.16f for (%d/%d, %d/%d, %d/%d)\n", counter, *r2, i_sites[0], s_types[0], i_sites[1], s_types[1], i_sites[2], s_types[2]); printf(" dist_diff = %20.16f\n", dist_diff); printf(" center = (%20.16f %20.16f)\n", w->x, w->y); ++counter; #endif if (dist_diff < old_diff) { for (j = 0; j < 3; j += 1) { if (s_types[j] == PNT) { start[j] = end[j] = samples[j][0]; } else { if (min_ind[j] == 0) { new_min = 0; new_max = 2; } else if (min_ind[j] == NS) { new_min = NS - 2; new_max = NS; } else { new_min = min_ind[j] - 1; new_max = min_ind[j] + 1; } u = start[j]; v = end[j]; if (s_types[j] == SEG) { t = new_min / ((double) NS); start[j] = LinearComb(u, v, t); t = new_max / ((double) NS); end[j] = LinearComb(u, v, t); } else { p = GetArcCenter(i_sites[j]); 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; du = PntPntDist(u, p); dv = PntPntDist(v, p); t = ((double) new_min) / NS; angle = (1.0 - t) * angle_s + t * angle_e; radius = (1.0 - t) * du + t * dv; start[j].x = p.x + radius * cos(angle); start[j].y = p.y + radius * sin(angle); t = ((double) new_max) / NS; angle = (1.0 - t) * angle_s + t * angle_e; radius = (1.0 - t) * du + t * dv; end[j].x = p.x + radius * cos(angle); end[j].y = p.y + radius * sin(angle); } } } } } while (dist_diff < old_diff); /* */ /* compare to and choose among old nodes */ /* */ n = GetStartNode(e); GetNodeData(n, &cntr, &radius); dist1 = ComputeSiteDistance(i_sites[0], s_types[0], cntr); dist2 = ComputeSiteDistance(i_sites[1], s_types[1], cntr); dist3 = ComputeSiteDistance(i_sites[2], s_types[2], cntr); dist = Abs(dist1 - radius) + Abs(dist2 - radius) + Abs(dist3 - radius); if (dist < old_diff) { *r2 = (dist1 + dist2 + dist3) / 3.0; *w = cntr; old_diff = dist; } n = GetEndNode(e); GetNodeData(n, &cntr, &radius); dist1 = ComputeSiteDistance(i_sites[0], s_types[0], cntr); dist2 = ComputeSiteDistance(i_sites[1], s_types[1], cntr); dist3 = ComputeSiteDistance(i_sites[2], s_types[2], cntr); dist = Abs(dist1 - radius) + Abs(dist2 - radius) + Abs(dist3 - radius); if (dist < old_diff) { *r2 = (dist1 + dist2 + dist3) / 3.0; *w = cntr; } return; } int vroniObject::DesperateFindMostViolatedNodeSeg(int i1, int j) { int e, k1, k2, n; coord p, q1, q2; double dist1, dist2, radius1, radius2; double a, b, c, lgth; e = GetVDPtr(i1, PNT); k1 = GetEndNode(e); k2 = GetStartNode(e); GetNodeData(k1, &q1, &radius1); GetNodeData(k2, &q2, &radius2); GetSegEqnData(j, &a, &b, &c); lgth = GetSegLgth(j); p = GetSegStartCoord(j); if (GetNodeSite(k1) != i1) dist1 = AbsPntSegDist(p, q1, a, b, c, lgth, TINY) - radius1; else dist1 = 0.0; if (GetNodeSite(k2) != i1) dist2 = AbsPntSegDist(p, q2, a, b, c, lgth, TINY) - radius2; else dist2 = 0.0; if (dist1 < dist2) n = StoreNode(q1.x, q1.y, radius1); else n = StoreNode(q2.x, q2.y, radius2); InsertDummyNode(e, n, false); desperate = true; return n; } void vroniObject::DesperatePreserveVoronoiCellSeg(int i, t_site t, int i0) { int e, n, e0, k, k0, e_max; double a, b, c, r_2, d, r_max, d_max; coord q, w_max; troubles = true; e = GetVDPtr(i, t); assert(InEdgesList(e)); assert(IsLftSite(e, i, t) || IsRgtSite(e, i, t)); GetSegEqnData(i0, &a, &b, &c); if (IsRgtSite(e, i, t)) { k = GetEndNode(e); k0 = GetStartNode(e); } else { k = GetStartNode(e); k0 = GetEndNode(e); } GetNodeData(k0, &q, &r_2); d_max = AbsPntLineDist(a, b, c, q); e_max = e; w_max = q; r_max = r_2; e0 = e; do { GetNodeData(k, &q, &r_2); d = AbsPntLineDist(a, b, c, q); if (d > d_max) { d_max = d; w_max = q; e_max = e; r_max = r_2; } e = GetCCWEdge(e, k); k = GetOtherNode(e, k); } while (e != e0); n = StoreNode(w_max.x, w_max.y, r_max); InsertDummyNode(e_max, n, true); desperate = true; return; } void vroniObject::DesperateBreakLoopSeg(int i0) { int start, number, n, k, j, e, k0, e_max; double a, b, c, d_max, r_max, d1, d0, r1_2, r0_2; coord q1, q0, w_max; troubles = true; number = GetLoopStackSizeFunc(); start = GetLoopStackMinFunc(); GetSegEqnData(i0, &a, &b, &c); d_max = -1.0; e_max = NIL; r_max = 0.0; w_max = GetSegStartCoord(i0); for (j = start; j < number; ++j) { e = GetLoopStackElementFunc(j); k = GetStartNode(e); k0 = GetEndNode(e); GetNodeData(k0, &q0, &r0_2); d0 = AbsPntLineDist(a, b, c, q0); GetNodeData(k, &q1, &r1_2); d1 = AbsPntLineDist(a, b, c, q1); if (d0 > d_max) { d_max = d0; e_max = e; w_max = q0; r_max = r0_2; } if (d1 > d_max) { d_max = d1; e_max = e; w_max = q1; r_max = r1_2; } } assert(e_max != NIL); n = StoreNode(w_max.x, w_max.y, r_max); InsertDummyNode(e_max, n, true); desperate = true; return; } int vroniObject::DesperateFindMostViolatedNodeArc(int i1, int j) { int e, k1, k2, n; coord p, q1, q2, ns, ne; double dist1, dist2, r, radius1, radius2; e = GetVDPtr(i1, PNT); k1 = GetEndNode(e); k2 = GetStartNode(e); GetNodeData(k1, &q1, &radius1); GetNodeData(k2, &q2, &radius2); ns = GetArcStartNormal(j); ne = GetArcEndNormal(j); p = GetArcCenter(j); r = GetArcRadius(j); if (GetNodeSite(k1) != i1) { q1 = VecSub(q1, p); dist1 = AbsPntArcDist(ns, ne, q1, r, TINY) - radius1; } else dist1 = 0.0; if (GetNodeSite(k2) != i1) { q2 = VecSub(q2, p); dist2 = AbsPntArcDist(ns, ne, q2, r, TINY) - radius2; } else dist2 = 0.0; if (dist1 < dist2) n = StoreNode(q1.x, q1.y, radius1); else n = StoreNode(q2.x, q2.y, radius2); InsertDummyNode(e, n, false); desperate = true; return n; } void vroniObject::DesperatePreserveVoronoiCellArc(int i, t_site t, int i0) { int e, n, e0, k, k0, e_max; double r, r_2, d, r_max, d_max; coord p, q, w_max; troubles = true; e = GetVDPtr(i, t); assert(InEdgesList(e)); assert(IsLftSite(e, i, t) || IsRgtSite(e, i, t)); if (IsRgtSite(e, i, t)) { k = GetEndNode(e); k0 = GetStartNode(e); } else { k = GetStartNode(e); k0 = GetEndNode(e); } GetNodeData(k0, &q, &r_2); p = GetArcCenter(i0); r = GetArcRadius(i0); d_max = AbsPntCircleDist(p, r, q); e_max = e; w_max = q; r_max = r_2; e0 = e; do { GetNodeData(k, &q, &r_2); d = AbsPntCircleDist(p, r, q); if (d > d_max) { d_max = d; w_max = q; e_max = e; r_max = r_2; } e = GetCCWEdge(e, k); k = GetOtherNode(e, k); } while (e != e0); n = StoreNode(w_max.x, w_max.y, r_max); InsertDummyNode(e_max, n, true); desperate = true; return; } void vroniObject::DesperateBreakLoopArc(int i0) { int start, number, n, k, j, e, k0, e_max; double d_max, r_max, d1, d0, r1_2, r0_2, r; coord q1, q0, w_max, c; troubles = true; number = GetLoopStackSizeFunc(); start = GetLoopStackMinFunc(); c = GetArcCenter(i0); r = GetArcRadius(i0); d_max = -1.0; e_max = NIL; r_max = 0.0; w_max = GetArcStartCoord(i0); for (j = start; j < number; ++j) { e = GetLoopStackElementFunc(j); k = GetStartNode(e); k0 = GetEndNode(e); GetNodeData(k0, &q0, &r0_2); GetNodeData(k, &q1, &r1_2); d0 = PntPntDist(c,q0); d0 = Abs(d0-r); d1 = PntPntDist(c,q1); d1 = Abs(d1-r); if (d0 > d_max) { d_max = d0; e_max = e; w_max = q0; r_max = r0_2; } if (d1 > d_max) { d_max = d1; e_max = e; w_max = q1; r_max = r1_2; } } assert(e_max != NIL); n = StoreNode(w_max.x, w_max.y, r_max); InsertDummyNode(e_max, n, true); desperate = true; return; }