Egalware/vroni
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
f3e15b8c8d6b7f1efad2a9da1edfa65b5c194c88
vroni/desperate.cc
T
SaraP f3e15b8c8d vroni 7.6 : 
- aggiunti file della libreria e progetto visual studio.
2023-09-06 15:44:02 +02:00

817 lines
23 KiB
C++
Raw Blame History

/*****************************************************************************/
/* */
/* Copyright (C) 1999-2023 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 <stdlib.h>
#include <stdio.h>
#include <assert.h>
#include <math.h>
#include <float.h>
/* */
/* 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;
}
Reference in New Issue View Git Blame Copy Permalink