Egalware/vroni
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
de65914f437f2635c6c33ed0cc3f69fcb75edd58
vroni/data_off.cc
T
SaraP de65914f43 vroni 7.6 : 
- modifiche per integrare vroni con le nostre librerie
- spostamento header in Extern per evitare duplicati.
2023-11-23 11:26:13 +01:00

1971 lines
63 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-172 */
/* Voice Mail: (+43 662) 8044-6304 */
/* Snail Mail: Martin Held */
/* FB Informatik */
/* Universitaet Salzburg */
/* A-5020 Salzburg, Austria */
/* */
/*****************************************************************************/
/* */
/* get standard libraries */
/* */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <limits.h>
#include <assert.h>
/* */
/* 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 ? dist<r : dist>r);
return (left_offset ? (dist<r) : (dist>r));
}
else {
//printf("\tCW orientation of arc, correct side: %d\n", left_offset ? dist>r : dist<r);
return (left_offset ? (dist>r) : (dist<r));
}
}
}
/* */
/* stack-based pseudo-recursive scan of edges of the VD until no unvisited */
/* VD edge exists. note that the scan does not cross contour boundaries. */
/* depending on the vr_bool offset_side, the edges scanned are stored in a */
/* list of active edges. if compute_mic then the function also updates the */
/* MIC information. */
/* */
void vroniObject::ScanVDEdges(int e, int n, vr_bool offset_side,
vr_bool compute_mic,
int *max_edge, double *max_rad)
{
double t;
int e1;
while (GetEdgeFlagNew(e)) {
SetEdgeFlagNew(e, false);
n = GetOtherNode(e, n);
t = GetNodeParam(n);
if (offset_side && IsEdgeGenuine(e, min_pnt_ind, max_pnt_ind)) {
if (compute_mic) {
/* */
/* determine maximum inscribed circle (MIC) on the fly. */
/* */
if (t > *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 (Max(dv, dx) > Max(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, &ltype);
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, &ltype);
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, &ltype);
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, &ltype);
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 (Max(dv, dx) > Max(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, &ltype);
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, &not_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"
Reference in New Issue View Git Blame Copy Permalink