Skip to content

Instantly share code, notes, and snippets.

@ialhashim

ialhashim/powercrust.h

Created December 26, 2014 03:18
Show Gist options
  • Star 2 You must be signed in to star a gist
  • Fork 1 You must be signed in to fork a gist
  • Save ialhashim/373a08ed160c43f028bb to your computer and use it in GitHub Desktop.
Save ialhashim/373a08ed160c43f028bb to your computer and use it in GitHub Desktop.
Download ZIP
Portable Powercrust in C++
Raw
powercrust.h
// Portable version of powercrust, adapted from https://github.com/timhutton/vtkpowercrust
/*===================================================================================================
vtkPowerCrustSurfaceReconstruction algorithm reconstructs surfaces from unorganized point data.
Copyright (C) 2014 Arash Akbarinia, Tim Hutton, Bruce Lamond Dieter Pfeffer, Oliver Moss
====================================================================================================*/
/*
#include "CellArray.h"
#include "FloatArray.h"
#include "Information.h"
#include "InformationVector.h"
#include "ObjectFactory.h"
#include "PointData.h"
*/
#include <assert.h>
#include <float.h>
#include <vector>
#include <iostream>
#include <fstream>
#include <sstream>
typedef double Coord;
typedef Coord* point;
typedef point site;
#define MAXBLOCKS 10000
#define Nobj 10000
#define STORAGE_GLOBALS(X) \
\
extern size_t X##_size; \
extern X *X##_list; \
extern X *new_block_##X(int); \
void free_##X##_storage(void); \
#define INCP(X,p,k) ((X*) ( (char*)p + (k) * X##_size)) /* portability? */
#define STORAGE(X) \
\
size_t X##_size; \
X *X##_list = 0; \
\
X *new_block_##X(int make_blocks) \
{ int i; \
static X *X##_block_table[MAXBLOCKS]; \
X *xlm, *xbt; \
static int num_##X##_blocks; \
if (make_blocks) { \
assert(num_##X##_blocks<MAXBLOCKS); \
DEB(0, before) DEBEXP(0, Nobj * X##_size) \
\
xbt = X##_block_table[num_##X##_blocks++] = (X*)malloc(Nobj * X##_size); \
memset(xbt,0,Nobj * X##_size); \
if (!xbt) { \
DEBEXP(-10,num_##X##_blocks) \
} \
assert(xbt); \
\
xlm = INCP(X,xbt,Nobj); \
for (i=0;i<Nobj; i++) { \
xlm = INCP(X,xlm,(-1)); \
xlm->next = X##_list; \
X##_list = xlm; \
} \
\
return X##_list; \
} \
\
for (i=0; i<num_##X##_blocks; i++) \
free(X##_block_table[i]); \
num_##X##_blocks = 0; \
X##_list = 0; \
return 0; \
} \
\
void free_##X##_storage(void) {new_block_##X(0);} \
/*end of STORAGE*/
#define NEWL(X,p) \
{ \
p = X##_list ? X##_list : new_block_##X(1); \
assert(p); \
X##_list = p->next; \
} \
#define NEWLRC(X,p) \
{ \
p = X##_list ? X##_list : new_block_##X(1); \
assert(p); \
X##_list = p->next; \
p->ref_count = 1; \
} \
#define FREEL(X,p) \
{ \
memset((p),0,X##_size); \
(p)->next = X##_list; \
X##_list = p; \
} \
#define dec_ref(X,v) {if ((v) && --(v)->ref_count == 0) FREEL(X,(v));}
#define inc_ref(X,v) {if (v) v->ref_count++;}
#define NULLIFY(X,v) {dec_ref(X,v); v = NULL;}
#define mod_refs(op,s) \
{ \
int i; \
neighbor *mrsn; \
\
for (i=-1,mrsn=s->neigh-1;i<cdim;i++,mrsn++) \
op##_ref(basis_s, mrsn->basis); \
}
#define copy_simp(new,s) \
{ NEWL(simplex,new); \
memcpy(new,s,simplex_size); \
mod_refs(inc,s); \
} \
#define DEBS(qq) {if (DEBUG>qq) {
#define EDEBS }}
#define DEBOUT 0
#define DEB(ll,mes) DEBS(ll) if(DEBOUT){fprintf(DEBOUT,#mes "\n");fflush(DEBOUT);} EDEBS
#define DEBEXP(ll,exp) DEBS(ll) if(DEBOUT){fprintf(DEBOUT,#exp "=%G\n", (double) exp); fflush(DEBOUT);} EDEBS
#define MAXDIM 8
#define BLOCKSIZE 100000
#define DEBUG -7
#define EXACT 1 /* sunghee */
#define CNV 0 /* sunghee : status of simplex, if it's on convex hull */
#define VV 1 /* sunghee : if it's regular simplex */
#define SLV -1 /* sunghee : if orient3d=0, sliver simplex */
#define AV 2 /* if av contains the averaged pole vector */
#define POLE_OUTPUT 3 /* VV is pole and it's ouput */
#define SQ(a) ((a)*(a)) /* sunghee */
#define BAD_POLE -1
// next two lines added by TJH to avoid name collision
#undef IN
#undef OUT
#define IN 2
#define OUT 1
#define INIT 0
#define NO 1
#define FIRST_EDGE 0
#define POW 1
#define VISITED 3
#define VALIDEDGE 24
#define ADDAXIS 13
/* for priority queue */
#define LEFT(i) ((i) * 2)
#define RIGHT(i) ((i) * 2 + 1)
#define PARENT(i) ((i) / 2)
typedef struct basis_s
{
struct basis_s *next; /* free list */
int ref_count; /* storage management */
int lscale; /* the log base 2 of total scaling of vector */
Coord sqa, sqb; /* sums of squared norms of a part and b part */
Coord vecs[1]; /* the actual vectors, extended by malloc'ing bigger */
} basis_s;
typedef struct neighbor
{
site vert; /* vertex of simplex */
/* short edgestatus[3]; FIRST_EDGE if not visited
NOT_POW if not dual to powercrust faces
POW if dual to powercrust faces */
struct simplex *simp; /* neighbor sharing all vertices but vert */
basis_s *basis; /* derived vectors */
} neighbor;
typedef struct simplex
{
simplex() : isVvNull(false) {}
simplex *next; /* used in free list */
short mark;
double vv[3];
bool isVvNull;
double sqradius; /* squared radius of Voronoi ball */
short status;/* sunghee : 0(CNV) if on conv hull so vv contains normal vector;
1(VV) if vv points to circumcenter of simplex;
-1(SLV) if cond=0 so vv points to hull
2(AV) if av contains averaged pole */
long poleindex; /* for 1st DT, if status==POLE_OUTPUT, contains poleindex; for 2nd, contains vertex index for powercrust output for OFF file format */
short edgestatus[6]; /* edge status :(01)(02)(03)(12)(13)(23)
FIRST_EDGE if not visited
VISITED
NOT_POW if not dual to powercrust faces
POW if dual to powercrust faces */
/* short tristatus[4]; triangle status :
FIRST if not visited
NO if not a triangle
DEG if degenerate triangle
SURF if surface triangle
NORM if fails normal test
VOR if falis voronoi edge test
VOR_NORM if fails both test */
/* NOTE!!! neighbors has to be the LAST field in the simplex stucture,
since it's length gets altered by some tricky Clarkson-move.
Also peak has to be the one before it.
Don't try to move these babies!! */
long visit; /* number of last site visiting this simplex */
basis_s* normal; /* normal vector pointing inward */
neighbor peak; /* if null, remaining vertices give facet */
neighbor neigh[1]; /* neighbors of simplex */
} simplex;
/* structure for list of opposite poles, opplist. */
typedef struct plist
{
long pid;
double angle;
plist *next;
} plist;
/* regular triangulation edge, between pole pid to center of simp? */
typedef struct edgesimp
{
short kth;
double angle; /* angle between balls */
simplex *simp;
long pid;
edgesimp *next;
} edgesimp;
/* additional info about poles: label for pole, pointer to list of regular
triangulation edges, squared radius of polar ball. adjlist is an
array of polelabels. */
typedef struct polelabel
{
edgesimp *eptr;
short bad;
short label;
double in; /* 12/7/99 Sunghee for priority queue */
double out; /* 12/7/99 Sunghee for priority queue */
int hid; /* 0 if not in the heap, otherwise heap index 1..heap_size*/
double sqradius;
double oppradius; /* minimum squared radius of this or any opposite ball */
double samp_distance;
int grafindex; /* index in thinning graph data structure */
} polelabel;
typedef struct fg_node fg;
typedef struct tree_node Tree;
struct tree_node
{
Tree *left, *right;
site key;
int size; /* maintained to be the number of nodes rooted here */
fg *fgs;
Tree *next; /* freelist */
};
typedef struct fg_node
{
Tree *facets;
double dist, vol; /* of Voronoi face dual to this */
fg *next; /* freelist */
short mark;
int ref_count;
} fg_node;
typedef struct heap_array
{
int pid;
double pri;
} heap_array;
STORAGE_GLOBALS(basis_s)
STORAGE_GLOBALS(fg)
STORAGE_GLOBALS(Tree)
STORAGE_GLOBALS(simplex)
STORAGE(basis_s)
STORAGE(fg)
STORAGE(Tree)
STORAGE(simplex)
//=====================================================================
class PowerCrustSurfaceReconstructionImpl
{
private:
typedef void* (PowerCrustSurfaceReconstructionImpl::*visit_func) (simplex *, void *);
typedef int (PowerCrustSurfaceReconstructionImpl::*test_func) (simplex *, int, void *);
void ASSERT(int b, const char* message = "");
int pcFALSE;
int pcTRUE;
site p;
std::vector<point> site_blocks;
std::vector<point> site_blocks_pointers;
int num_blocks;
int pdim;
Coord infinity[10]; /* point at infinity for Delaunay triang */
int rdim, /* region dimension: (max) number of sites specifying region */
cdim, /* number of sites currently specifying region */
site_size; /* size of malloc needed for a site */
site get_site_offline(long);
double bound[8][3];
double mult_up;
Coord mins[MAXDIM];
Coord maxs[MAXDIM];
double Huge;
void read_bounding_box(long);
void construct_face(simplex *, short);
void compute_distance(simplex**, int, double*);
#define MAXPOINTS 10000
short mi[MAXPOINTS], mo[MAXPOINTS];
int correct_orientation(double*, double*, double*, double*, double*);
#ifndef _RAND48_H_
#define _RAND48_H_
void _dorand48(unsigned short xseed[3]);
#define RAND48_MULT_0 (0xe66d)
#define RAND48_MULT_1 (0xdeec)
#define RAND48_MULT_2 (0x0005)
#define RAND48_ADD (0x000b)
#endif /* _RAND48_H_ */
double omaxs[3], omins[3]; /* 8 vertices for bounding box */
int num_vtxs;
int num_poles;
/* Data structures for poles */
simplex **pole1, **pole2; /* arrays of poles - per sample*/
std::vector<polelabel> adjlist;
std::vector<plist*> opplist;
double* lfs_lb; /* array of lower bounds for lfs of each sample */
double est_r; /* estimated value of r - user input */
double *pole1_distance, *pole2_distance;
/* for priority queue */
int heap_size;
int scount;
int v1[6];
int v2[6];
int v3[6];
int v4[6];
long num_sites;
short vd_new;
short power_diagram; /* 1 if power diagram */
int dim;
long s_num; /* site number */
double theta; /* input argument - angle defining deep intersection */
double deep; /* input argument.. same as theta for labeling unlabled pole */
long site_numm(site p);
site new_site(site p, long j);
// TJH: trying to replace file use
site read_next_site(long j);
site pole_read_next_site(long j);
std::vector<long> shufmat;
long mat_size; // EPRO set to global to reinitialize
void make_shuffle(void);
long shufflef(long i);
long (PowerCrustSurfaceReconstructionImpl::*shuf) (long);
long (PowerCrustSurfaceReconstructionImpl::*site_num) (site);
site(PowerCrustSurfaceReconstructionImpl::*get_site) ();
site(PowerCrustSurfaceReconstructionImpl::*get_site_n) (long);
site get_next_site(void);
void make_output(simplex *root, void * (PowerCrustSurfaceReconstructionImpl::*visit_gen) (simplex*, visit_func), visit_func visit);
neighbor p_neigh; // EPRO added set to global
basis_s *seesB;
simplex **st_search;
simplex **st_visit_triang_gen;
basis_s *check_perps_b;
#define swap_points(a,b) {point t; t=a; a=b; b=t;}
double alpha_test_alpha;
/* variables for tracking infinite loop */
/* (This should never occur, but it did in early versions) */
int loopStart;
int count;
int lastCount;
#define compare(i,j) (( this->*site_num )(i)-( this->*site_num )(j))
#define node_size(x) ((x) ? ((x)->size) : 0 )
heap_array *heap_A;
int heap_length;
long pnum;
#define push(x, st, tms) *(st + tms++) = x;
#define pop(x, st, tms) x = *(st + --tms);
long visit_triang_gen_vnum;
long visit_triang_gen_ss;
simplex *make_facets_ns;
long search_ss;
simplex *ch_root;
#define NEARZERO(d) ((d) < FLT_EPSILON && (d) > -FLT_EPSILON)
#define SWAP(X,a,b) {X t; t = a; a = b; b = t;}
#define DELIFT 0
int basis_vec_size;
int exact_bits;
double b_err_min, b_err_min_sq;
short vd;
basis_s tt_basis;
basis_s *tt_basisp;
basis_s *infinity_basis;
#define VA(x) ((x)->vecs+rdim)
#define VB(x) ((x)->vecs)
#define two_to(x) ( ((x)<20) ? 1<<(x) : ldexp(1.0,(x)) ) // EPRO added
int sc_lscale;
double sc_max_scale, sc_ldetbound, sc_Sb;
Coord Vec_dot(point x, point y);
Coord Vec_dot_pdim(point x, point y);
Coord Norm2(point x);
void Ax_plus_y(Coord a, point x, point y);
void Ax_plus_y_test(Coord a, point x, point y);
void Vec_scale_test(int n, Coord a, Coord *x);
double sc(basis_s *v, simplex *s, int k, int j);
int reduce_inner(basis_s *v, simplex *s, int k);
void trans(point z, point p, point q);
#define lift(z,s) {if (vd) z[2*rdim-1] =z[rdim-1]= ldexp(Vec_dot_pdim(z,z), -DELIFT);}
int reduce(basis_s **v, point p, simplex *s, int k);
void get_basis_sede(simplex *s);
int out_of_flat(simplex *root, point p);
double cosangle_sq(basis_s* v, basis_s* w);
int check_perps(simplex *s);
void get_normal_sede(simplex *s);
void get_normal(simplex *s);
int sees(site p, simplex *s);
double radsq(simplex *s);
void* zero_marks(simplex* s, void* dum);
void* one_marks(simplex* s, void* dum);
void* conv_facetv(simplex* s, void* dum);
void* mark_points(simplex* s, void* dum);
int alph_test(simplex *s, int i, void *alphap);
void* visit_outside_ashape(simplex *root, visit_func visit);
int check_ashape(simplex *root, double alpha);
simplex *build_convex_hull(short dim, short vdd);
void free_hull_storage(void);
void *compute_vv(simplex *s, void *p);
void *compute_pole2(simplex *s, void *p);
void *compute_3d_power_vv(simplex *s, void *p);
void *compute_axis(simplex *s, void *p);
int close_pole(double* v, double* p, double lfs_lb);
int antiLabel(int label);
/* computes angle between two opposite poles */
double computePoleAngle(simplex* pole1, simplex* pole2, double* samp);
/* Adds a new pair of opposite poles to each other's lists */
void newOpposite(int p1index, int p2index, double pole_angle);
/* Outputs a pole, saving it's squared radius in adjlist */
void outputPole(simplex* pole, int poleid, double* samp, int* num_poles, double distance);
/* Splay using the key i (which may or may not be in the tree.) */
/* The starting root is t, and the tree used is defined by rat */
/* size fields are maintained */
Tree * splay(site i, Tree *t);
Tree * insert(site i, Tree * t);
void free_heap();
void init_heap(int num);
void heapify(int hi);
int extract_max();
int insert_heap(int pi, double pr);
/* make the element heap_A[hi].pr = pr ... */
void update(int hi, double pr);
void *visit_triang_gen(simplex *s, visit_func visit, test_func test);
int truet(simplex *s, int i, void *dum);
void *visit_triang(simplex *root, visit_func visit);
int hullt(simplex *s, int i, void *dummy);
void *facet_test(simplex *s, void *dummy);
/* visit all simplices with facets of the current hull */
void *visit_hull(simplex *root, visit_func visit);
#define lookup(a, b, what) \
{ \
int i; \
neighbor *x; \
for (i = 0, x = a->neigh; (x->what != b) && (i < cdim) ; i++, x++); \
if (i < cdim) \
return x; \
else { \
ASSERT(pcFALSE, "adjacency failure!"); \
return 0; \
} \
} \
neighbor *op_simp(simplex *a, simplex *b);
neighbor *op_vert(simplex *a, site b);
void connect(simplex *s);
simplex *make_facets(simplex *seen);
simplex *extend_simplices(simplex *s);
simplex *search(simplex *root);
point get_another_site(void);
void buildhull(simplex *root);
int propagate();
void opp_update(int pi);
void sym_update(int pi);
void update_pri(int hi, int pi);
void label_unlabeled(int num);
double sqdist(double a[3], double b[3]);
void dir_and_dist(double a[3], double b[3], double dir[3], double* dist);
void tetcircumcenter(double a[3], double b[3], double c[3], double d[3], double circumcenter[3], double *cond);
void tetorthocenter(double a[4], double b[4], double c[4], double d[4], double orthocenter[3], double *cnum);
#define INEXACT /* Nothing */
#define REAL double /* float or double */
#define Absolute(a) ((a) >= 0.0 ? (a) : -(a))
#define Fast_Two_Sum_Tail(a, b, x, y) \
bvirt = x - a; \
y = b - bvirt
#define Fast_Two_Sum(a, b, x, y) \
x = (REAL) (a + b); \
Fast_Two_Sum_Tail(a, b, x, y)
#define Two_Sum_Tail(a, b, x, y) \
bvirt = (REAL) (x - a); \
avirt = x - bvirt; \
bround = b - bvirt; \
around = a - avirt; \
y = around + bround
#define Two_Sum(a, b, x, y) \
x = (REAL) (a + b); \
Two_Sum_Tail(a, b, x, y)
#define Two_Diff_Tail(a, b, x, y) \
bvirt = (REAL) (a - x); \
avirt = x + bvirt; \
bround = bvirt - b; \
around = a - avirt; \
y = around + bround
#define Two_Diff(a, b, x, y) \
x = (REAL) (a - b); \
Two_Diff_Tail(a, b, x, y)
#define Split(a, ahi, alo) \
c = (REAL) (splitter * a); \
abig = (REAL) (c - a); \
ahi = c - abig; \
alo = a - ahi
#define Two_Product_Tail(a, b, x, y) \
Split(a, ahi, alo); \
Split(b, bhi, blo); \
err1 = x - (ahi * bhi); \
err2 = err1 - (alo * bhi); \
err3 = err2 - (ahi * blo); \
y = (alo * blo) - err3
#define Two_Product(a, b, x, y) \
x = (REAL) (a * b); \
Two_Product_Tail(a, b, x, y)
/* Two_Product_Presplit() is Two_Product() where one of the inputs has */
/* already been split. Avoids redundant splitting. */
#define Two_Product_Presplit(a, b, bhi, blo, x, y) \
x = (REAL) (a * b); \
Split(a, ahi, alo); \
err1 = x - (ahi * bhi); \
err2 = err1 - (alo * bhi); \
err3 = err2 - (ahi * blo); \
y = (alo * blo) - err3
/* Macros for summing expansions of various fixed lengths. These are all */
/* unrolled versions of Expansion_Sum(). */
#define Two_One_Diff(a1, a0, b, x2, x1, x0) \
Two_Diff(a0, b , _i, x0); \
Two_Sum( a1, _i, x2, x1)
#define Two_Two_Diff(a1, a0, b1, b0, x3, x2, x1, x0) \
Two_One_Diff(a1, a0, b0, _j, _0, x0); \
Two_One_Diff(_j, _0, b1, x3, x2, x1)
#define Two_One_Product(a1, a0, b, x3, x2, x1, x0) \
Split(b, bhi, blo); \
Two_Product_Presplit(a0, b, bhi, blo, _i, x0); \
Two_Product_Presplit(a1, b, bhi, blo, _j, _0); \
Two_Sum(_i, _0, _k, x1); \
Fast_Two_Sum(_j, _k, x3, x2)
REAL splitter; /* = 2^ceiling(p / 2) + 1. Used to split floats in half. */
REAL epsilon; /* = 2^(-p). Used to estimate roundoff errors. */
/* A set of coefficients used to calculate maximum roundoff errors. */
REAL resulterrbound;
REAL ccwerrboundA, ccwerrboundB, ccwerrboundC;
REAL o3derrboundA, o3derrboundB, o3derrboundC;
void exactinit();
int fast_expansion_sum_zeroelim(int elen, REAL *e, int flen, REAL *f, REAL *h);
int scale_expansion_zeroelim(int elen, REAL *e, REAL b, REAL *h);
REAL estimate(int elen, REAL *e);
REAL orient2dadapt(REAL* pa, REAL* pb, REAL* pc, REAL detsum);
REAL orient3dadapt(REAL* pa, REAL* pb, REAL* pc, REAL* pd, REAL permanent);
REAL orient3d(REAL* pa, REAL* pb, REAL* pc, REAL* pd);
unsigned short X[3];
double double_rand(void);
void init_rand(void);
double logb(double x);
double local_erand48(unsigned short xseed[3]) throw();
unsigned short _rand48_mult[3];
unsigned short _rand48_add;
public:
PowerCrustSurfaceReconstructionImpl(void);
void pcInit();
void freeAll(void);
void adapted_main(double m_mult_up = 1000000);
// these globals are here so we can access them from anywhere in the powercrust code
// if you can find a neat way to improve this then please feel free
struct PolyData
{
double * GetPoint(long j)
{
return &data[j * 3];
}
double GetWeight(long j)
{
return weights[j];
}
long GetNumberOfPoints()
{
return (long)data.size() / 3;
}
void InsertNextPoint(double vx, double vy, double vz)
{
data.push_back(vx);
data.push_back(vy);
data.push_back(vz);
}
void InsertNextPoint(float * v)
{
data.push_back(v[0]);
data.push_back(v[1]);
data.push_back(v[2]);
}
void InsertWeight(double weight)
{
weights.push_back(weight);
}
void InsertTriangle(long v1, long v2, long v3)
{
indices.push_back(v1);
indices.push_back(v2);
indices.push_back(v3);
}
void InsertCell(std::vector<int> & indices)
{
cells.push_back(indices);
}
void readXYZ(std::string filename)
{
std::ifstream ifs(filename);
std::string line;
std::vector<std::string> words;
while (std::getline(ifs, line))
{
std::stringstream linestream(line);
copy(std::istream_iterator<std::string>(linestream),
std::istream_iterator<std::string>(),
std::back_inserter(words));
}
for (auto & word : words)
{
data.push_back( std::atof(word.c_str()) );
}
}
void writeOFF(std::string filename)
{
std::ofstream ofs(filename);
std::ostringstream ss;
// Header
ss << "OFF " << GetNumberOfPoints() << " " << cells.size() << " 0\n";
// Points
for (int i = 0; i < GetNumberOfPoints(); i++){
auto p = GetPoint(i);
ss << p[0] << " " << p[1] << " " << p[2] << "\n";
}
// Faces
for (int i = 0; i < (int)cells.size(); i++){
ss << cells[i].size() << " ";
for (int j = 0; j < (int)cells[i].size(); j++)
ss << cells[i][j] << " ";
ss << "\n";
}
ofs << ss.str();
}
std::vector< std::vector<double> > GetPoints()
{
std::vector< std::vector<double> > pts;
for (int i = 0; i < GetNumberOfPoints(); i++){
std::vector<double> pnt;
auto p = GetPoint(i);
pnt.push_back(p[0]);
pnt.push_back(p[1]);
pnt.push_back(p[2]);
pts.push_back(pnt);
}
return pts;
}
std::vector < std::vector<int> > GetFaces(){ return cells; }
std::vector < double > data;
std::vector < double > weights;
std::vector < long > indices;
std::vector < std::vector<int> > cells;
};
PolyData input;
PolyData output;
PolyData medial_surface;
};
//=====================================================================
PowerCrustSurfaceReconstructionImpl::PowerCrustSurfaceReconstructionImpl(void)
{
}
long PowerCrustSurfaceReconstructionImpl::site_numm(site p)
{
if ((vd_new || power_diagram) && p == infinity) return -1;
if (!p) return -2;
for (int i = 0; i < num_blocks; i++)
{
long j;
if ((j = p - site_blocks[i]) >= 0 && j < BLOCKSIZE * dim)
return j / dim + BLOCKSIZE * i;
}
return -3;
}
site PowerCrustSurfaceReconstructionImpl::new_site(site p, long j)
{
assert(num_blocks + 1 < MAXBLOCKS);
if (0 == (j % BLOCKSIZE))
{
num_blocks++;
site_blocks.resize(num_blocks);
site_blocks[num_blocks - 1] = (site)malloc(BLOCKSIZE * site_size);
site_blocks_pointers.push_back(site_blocks[num_blocks - 1]);
return site_blocks[num_blocks - 1];
}
else
return p + dim;
}
void PowerCrustSurfaceReconstructionImpl::read_bounding_box(long j)
{
int i, k;
double center[3], width;
omaxs[0] = maxs[0];
omins[0] = mins[0];
omaxs[1] = maxs[1];
omins[1] = mins[1];
omaxs[2] = maxs[2];
omins[2] = mins[2];
center[0] = (maxs[0] - mins[0]) / 2;
center[1] = (maxs[1] - mins[1]) / 2;
center[2] = (maxs[2] - mins[2]) / 2;
if ((maxs[0] - mins[0]) > (maxs[1] - mins[1]))
{
if ((maxs[2] - mins[2]) > (maxs[0] - mins[0]))
width = maxs[2] - mins[2];
else width = maxs[0] - mins[0];
}
else
{
if ((maxs[1] - mins[1]) > (maxs[2] - mins[2]))
width = maxs[1] - mins[1];
else width = maxs[2] - mins[2];
}
width = width * 4;
bound[0][0] = center[0] + width;
bound[1][0] = bound[0][0];
bound[2][0] = bound[0][0];
bound[3][0] = bound[0][0];
bound[0][1] = center[1] + width;
bound[1][1] = bound[0][1];
bound[4][1] = bound[0][1];
bound[5][1] = bound[0][1];
bound[0][2] = center[2] + width;
bound[2][2] = bound[0][2];
bound[4][2] = bound[0][2];
bound[6][2] = bound[0][2];
bound[4][0] = center[0] - width;
bound[5][0] = bound[4][0];
bound[6][0] = bound[4][0];
bound[7][0] = bound[4][0];
bound[2][1] = center[1] - width;
bound[3][1] = bound[2][1];
bound[6][1] = bound[2][1];
bound[7][1] = bound[2][1];
bound[1][2] = center[2] - width;
bound[3][2] = bound[1][2];
bound[5][2] = bound[1][2];
bound[7][2] = bound[1][2];
for (k = 0; k < 3; k++)
{
p[k] = bound[0][k];
}
for (i = 1; i < 8; i++)
{
p = new_site(p, j + i);
for (k = 0; k < 3; k++)
{
p[k] = bound[i][k];
}
}
maxs[0] = bound[0][0];
mins[0] = bound[4][0];
maxs[1] = bound[0][1];
mins[1] = bound[2][1];
maxs[2] = bound[0][2];
mins[2] = bound[1][2];
}
site PowerCrustSurfaceReconstructionImpl::read_next_site(long j)
{
ASSERT(j >= 0, "read_next_site ");
p = new_site(p, j);
for (int i = 0; i < dim; i++)
{
p[i] = (double)input.GetPoint(j)[i];
p[i] = floor(mult_up*p[i] + 0.5);
mins[i] = (mins[i] < p[i]) ? mins[i] : p[i];
maxs[i] = (maxs[i] > p[i]) ? maxs[i] : p[i];
}
return p;
}
site PowerCrustSurfaceReconstructionImpl::pole_read_next_site(long j)
{
ASSERT(j >= 0, "pole_read_next_site");
p = new_site(p, j);
for (int i = 0; i < dim; i++)
{
if (i < 3)
p[i] = (double)medial_surface.GetPoint(j)[i];
else
p[i] = medial_surface.GetWeight(j); // (double)medial_surface.GetPointData()->GetScalars()->GetTuple1(j);
p[i] = floor(mult_up*p[i] + 0.5);
mins[i] = (mins[i] < p[i]) ? mins[i] : p[i];
maxs[i] = (maxs[i] > p[i]) ? maxs[i] : p[i];
}
return p;
}
site PowerCrustSurfaceReconstructionImpl::get_site_offline(long i)
{
if (i >= num_sites) return NULL;
else
{
return site_blocks[i / BLOCKSIZE] + (i % BLOCKSIZE) * dim;
}
}
void PowerCrustSurfaceReconstructionImpl::make_shuffle(void)
{
if (mat_size <= num_sites)
{
mat_size = num_sites + 1;
shufmat.resize(mat_size);
}
for (long i = 0; i <= num_sites; i++)
{
shufmat[i] = i;
}
for (long i = 0; i < num_sites; i++)
{
long t = shufmat[i];
long j = i + (long)((num_sites - i) *double_rand()); // cast to long added by TJH
shufmat[i] = shufmat[j];
shufmat[j] = t;
}
}
long PowerCrustSurfaceReconstructionImpl::shufflef(long i)
{
return shufmat[i];
}
site PowerCrustSurfaceReconstructionImpl::get_next_site(void)
{
return (this->*get_site_n) ((this->*shuf) (s_num++));
}
void PowerCrustSurfaceReconstructionImpl::make_output(simplex *root, void * (PowerCrustSurfaceReconstructionImpl::*visit_gen) (simplex*, visit_func), visit_func visit)
{
// ( this->*visit ) ( 0, ( void ( * ) ) out_funcp );
(this->*visit_gen) (root, visit);
}
void PowerCrustSurfaceReconstructionImpl::ASSERT(int b, const char*)
{
if (!b)
{
//our_filter->Error ( message );
}
}
void PowerCrustSurfaceReconstructionImpl::freeAll(void)
{
for (unsigned int i = 0; i < adjlist.size(); i++)
{
edgesimp *tmpEdgesimp;
edgesimp *curEdgesimp;
curEdgesimp = adjlist[i].eptr;
while (curEdgesimp)
{
tmpEdgesimp = curEdgesimp;
curEdgesimp = tmpEdgesimp->next;
free(tmpEdgesimp);
}
}
for (unsigned int i = 0; i < opplist.size(); i++)
{
if (opplist[i] != NULL)
{
plist *tmpPlist;
plist *curPlist;
curPlist = opplist[i];
while (curPlist)
{
tmpPlist = curPlist;
curPlist = tmpPlist->next;
free(tmpPlist);
}
}
}
free(pole1_distance);
free(pole2_distance);
free(pole1);
free(pole2);
free(lfs_lb);
free_hull_storage();
free(p_neigh.basis);
free(seesB);
free(st_search);
free(st_visit_triang_gen);
for (unsigned int i = 0; i < site_blocks_pointers.size(); i++)
free(site_blocks_pointers[i]);
free_heap(); // EPRO added
}
void PowerCrustSurfaceReconstructionImpl::adapted_main(double m_mult_up)
{
int num_poles = 0;
simplex *root;
edgesimp *eindex;
visit_func pr;
mult_up = m_mult_up;
//est_r = our_filter->GetEstimateR();
short bad = 0;
long poleid = 0;
double samp[3];
dim = 3;
if (dim > MAXDIM) ASSERT(pcFALSE, "dimension bound MAXDIM exceeded");
site_size = sizeof(Coord) *dim;
// TJH: we've replaced this file-reading loop with the loop below
for (num_sites = 0; num_sites < input.GetNumberOfPoints(); num_sites++)
{
read_next_site(num_sites);
}
num_sites--;
read_bounding_box(num_sites);
num_sites += 8;
init_rand();
make_shuffle();
shuf = &PowerCrustSurfaceReconstructionImpl::shufflef;
get_site_n = &PowerCrustSurfaceReconstructionImpl::get_site_offline;
/* Step 1: compute DT of input point set */
root = build_convex_hull(dim, vd_new);
/* Step 2: Find poles */
pole1 = (simplex **)calloc(num_sites, sizeof(simplex *));
pole2 = (simplex **)calloc(num_sites, sizeof(simplex *));
lfs_lb = (double*)calloc(num_sites, sizeof(double));
exactinit(); /* Shewchuk's exact arithmetic initialization */
pr = &PowerCrustSurfaceReconstructionImpl::compute_vv;
make_output(root, &PowerCrustSurfaceReconstructionImpl::visit_hull, pr);
pr = &PowerCrustSurfaceReconstructionImpl::compute_pole2;
make_output(root, &PowerCrustSurfaceReconstructionImpl::visit_hull, pr);
/* initialize the sample distance info for the poles */
pole1_distance = (double *)malloc(num_sites*sizeof(double));
pole2_distance = (double *)malloc(num_sites*sizeof(double));
compute_distance(pole1, num_sites - 8, pole1_distance);
compute_distance(pole2, num_sites - 8, pole2_distance);
/* intialize list of lists of pointers to opposite poles */
opplist.resize(2 * num_sites);
/* data about poles; adjacencies, labels, radii */
adjlist.resize(2 * num_sites);
/* loop through sites, writing out poles */
for (int i = 0; i < num_sites - 8; i++)
{
/* rescale the sample to real input coordinates */
for (int k = 0; k < 3; k++)
samp[k] = get_site_offline(i)[k] / mult_up;
/* output poles, both to debugging file and for weighted DT */
/* remembers sqaured radius */
if ((pole1[i] != NULL) && (pole1[i]->status != POLE_OUTPUT))
{
/* if second pole is closer than we think it should be... */
if ((pole2[i] != NULL) && bad && close_pole(samp, pole2[i]->vv, lfs_lb[i]))
{
}
else
{
outputPole(pole1[i], poleid++, samp, &num_poles, pole1_distance[i]);
}
}
if ((pole2[i] != NULL) && (pole2[i]->status != POLE_OUTPUT))
{
/* if pole is closer than we think it should be... */
if (close_pole(samp, pole2[i]->vv, lfs_lb[i]))
{
/* remember opposite bad for late labeling */
if (!bad) adjlist[pole1[i]->poleindex].bad = BAD_POLE;
continue;
}
outputPole(pole2[i], poleid++, samp, &num_poles, pole2_distance[i]);
}
/* keep list of opposite poles for later coloring */
if ((pole1[i] != NULL) && (pole2[i] != NULL) && (pole1[i]->status == POLE_OUTPUT) && (pole2[i]->status == POLE_OUTPUT))
{
double pole_angle = computePoleAngle(pole1[i], pole2[i], samp);
newOpposite(pole1[i]->poleindex, pole2[i]->poleindex, pole_angle);
newOpposite(pole2[i]->poleindex, pole1[i]->poleindex, pole_angle);
}
}
free_hull_storage();
power_diagram = 1;
vd_new = 0;
dim = 4;
// TJH: uncommented this line, seemed like we need to free this memory before num_blocks gets reset
for (unsigned int i = 0; i < site_blocks_pointers.size(); i++)
{
free(site_blocks_pointers[i]);
}
site_blocks_pointers.clear();
num_blocks = 0;
s_num = 0;
scount = 0;
site_size = sizeof(Coord) *dim;
/* save points in order read */
for (num_sites = 0; num_sites < medial_surface.GetNumberOfPoints(); num_sites++)
{
pole_read_next_site(num_sites);
}
//num_sites--;
/* set up the shuffle */
init_rand();
make_shuffle();
shuf = &PowerCrustSurfaceReconstructionImpl::shufflef;
get_site_n = &PowerCrustSurfaceReconstructionImpl::get_site_offline; /* returns stored points, unshuffled */
/* Compute weighted DT */
root = build_convex_hull(dim, vd_new);
/* compute adjacencies and find angles of ball intersections */
pr = &PowerCrustSurfaceReconstructionImpl::compute_3d_power_vv;
make_output(root, &PowerCrustSurfaceReconstructionImpl::visit_hull, pr);
/* labeling */
init_heap(num_poles);
for (int i = 0; i < num_poles; i++)
{
if ((get_site_offline(i)[0]>(2 * omaxs[0] - omins[0])) ||
(get_site_offline(i)[0]< (2 * omins[0] - omaxs[0])) ||
(get_site_offline(i)[1]> (2 * omaxs[1] - omins[1])) ||
(get_site_offline(i)[1]< (2 * omins[1] - omaxs[1])) ||
(get_site_offline(i)[2]> (2 * omaxs[2] - omins[2])) ||
(get_site_offline(i)[2] < (2 * omins[2] - omaxs[2])))
{
adjlist[i].hid = insert_heap(i, 1.0);
adjlist[i].out = 1.0;
adjlist[i].label = OUT;
}
}
while (heap_size != 0)
propagate();
label_unlabeled(num_poles);
for (int i = 0; i < num_poles; i++)
{
if ((adjlist[i].label != IN) && (adjlist[i].label != OUT))
{
}
else
{
eindex = adjlist[i].eptr;
while (eindex != NULL)
{
if ((i < eindex->pid) && (antiLabel(adjlist[i].label) == adjlist[eindex->pid].label))
{
construct_face(eindex->simp, eindex->kth);
}
eindex = eindex->next;
}
}
}
/* compute the medial axis */
pr = &PowerCrustSurfaceReconstructionImpl::compute_axis;
make_output(root, &PowerCrustSurfaceReconstructionImpl::visit_hull, pr);
}
void PowerCrustSurfaceReconstructionImpl::compute_distance(simplex** poles, int size, double* distance)
{
double indices[4][3]; /* the coords of the four vertices of the simplex*/
point v[MAXDIM];
simplex* currSimplex;
double maxdistance = 0;
for (int l = 0; l < size; l++) /* for each pole do*/
{
if (poles[l] != NULL)
{
currSimplex = poles[l];
/* get the coordinates of the four endpoints */
for (int j = 0; j < 4; j++)
{
v[j] = currSimplex->neigh[j].vert;
for (int k = 0; k < 3; k++)
indices[j][k] = v[j][k] / mult_up;
}
/* now compute the actual distance */
maxdistance = 0;
for (int i = 0; i < 4; i++)
{
for (int j = i + 1; j < 4; j++)
{
double currdistance = SQ(indices[i][0] - indices[j][0]) + SQ(indices[i][1] - indices[j][1]) + SQ(indices[i][2] - indices[j][2]);
currdistance = sqrt(currdistance);
if (maxdistance < currdistance)
maxdistance = currdistance;
}
}
distance[l] = maxdistance;
}
}
}
void PowerCrustSurfaceReconstructionImpl::pcInit()
{
pcFALSE = (1 == 0);
pcTRUE = (1 == 1);
//input = NULL;
//output = NULL;
//medial_surface = NULL;
//our_filter = NULL;
seesB = NULL;
check_perps_b = NULL;
p = NULL;
pole1_distance = NULL;
pole2_distance = NULL;
st_search = NULL;
st_visit_triang_gen = NULL;
ch_root = NULL;
pole1 = NULL;
pole2 = NULL;
lfs_lb = NULL;
site_num = NULL;
heap_A = NULL;
make_facets_ns = NULL;
site_blocks.clear();
site_blocks_pointers.clear();
adjlist.clear();
opplist.clear();
num_blocks = 0;
mult_up = 1.0;
num_vtxs = 0;
est_r = 0.6;
num_poles = 0;
heap_size = 0;
scount = 0;
vd_new = 1;
power_diagram = 0;
s_num = 0;
theta = 0.0;
deep = 0.0;
mat_size = 0;
loopStart = -1;
count = 0;
lastCount = 0;
visit_triang_gen_vnum = -1;
visit_triang_gen_ss = 2000;
pdim = 0;
rdim = 0;
cdim = 0;
site_size = 0;
pnum = 0;
num_sites = 0;
dim = 0;
basis_vec_size = 0;
exact_bits = 0;
b_err_min = 0;
b_err_min_sq = 0;
sc_lscale = 0;
sc_max_scale = 0;
sc_ldetbound = 0;
sc_Sb = 0;
alpha_test_alpha = 0;
heap_length = 0;
search_ss = MAXDIM;
Huge = DBL_MAX;
vd = 0;
_rand48_add = RAND48_ADD;
infinity[0] = 57.2;
infinity[1] = 0;
infinity[2] = 0;
infinity[3] = 0;
infinity[4] = 0;
v1[0] = 0;
v1[1] = 0;
v1[2] = 0;
v1[3] = 1;
v1[4] = 1;
v1[5] = 2;
v2[0] = 1;
v2[1] = 2;
v2[2] = 3;
v2[3] = 2;
v2[4] = 3;
v2[5] = 3;
v3[0] = 2;
v3[1] = 3;
v3[2] = 1;
v3[3] = 3;
v3[4] = 0;
v3[5] = 0;
v4[0] = 3;
v4[1] = 1;
v4[2] = 2;
v4[3] = 0;
v4[4] = 2;
v4[5] = 1;
for (int i = 0; i < 8; i++)
{
for (int j = 0; j < 3; j++)
{
bound[i][j] = 0;
}
}
for (int i = 0; i < MAXPOINTS; i++)
{
mi[i] = 0;
mo[i] = 0;
}
omins[0] = 0;
omins[1] = 0;
omins[2] = 0;
omaxs[0] = 0;
omaxs[1] = 0;
omaxs[2] = 0;
p_neigh.vert = 0;
p_neigh.basis = 0;
p_neigh.simp = 0;
X[0] = 0;
X[1] = 10000;
X[2] = 0;
for (int i = 0; i < MAXDIM; i++)
{
mins[i] = DBL_MAX;
maxs[i] = -DBL_MAX;
}
tt_basis.next = 0;
tt_basis.ref_count = 1;
tt_basis.lscale = -1;
tt_basis.sqa = 0;
tt_basis.sqb = 0;
tt_basis.vecs[0] = 0;
tt_basisp = &tt_basis;
infinity_basis = NULL;
X[0] = 0;
X[1] = 0;
X[2] = 0;
_rand48_mult[0] = RAND48_MULT_0;
_rand48_mult[1] = RAND48_MULT_1;
_rand48_mult[2] = RAND48_MULT_2;
_rand48_add = RAND48_ADD;
}
void PowerCrustSurfaceReconstructionImpl::_dorand48(unsigned short xseed[3])
{
unsigned long accu;
unsigned short temp[2];
accu = (unsigned long)_rand48_mult[0] * (unsigned long)xseed[0] + (unsigned long)_rand48_add;
temp[0] = (unsigned short)accu; /* lower 16 bits */
accu >>= sizeof(unsigned short) * 8;
accu += (unsigned long)_rand48_mult[0] * (unsigned long)xseed[1] + (unsigned long)_rand48_mult[1] * (unsigned long)xseed[0];
temp[1] = (unsigned short)accu; /* middle 16 bits */
accu >>= sizeof(unsigned short) * 8;
accu += _rand48_mult[0] * xseed[2] + _rand48_mult[1] * xseed[1] + _rand48_mult[2] * xseed[0];
xseed[0] = temp[0];
xseed[1] = temp[1];
xseed[2] = (unsigned short)accu;
}
double PowerCrustSurfaceReconstructionImpl::local_erand48(unsigned short xseed[3]) throw()
{
_dorand48(xseed);
return ldexp((double)xseed[0], -48) + ldexp((double)xseed[1], -32) + ldexp((double)xseed[2], -16);
}
double PowerCrustSurfaceReconstructionImpl::double_rand(void)
{
return local_erand48(X);
}
void PowerCrustSurfaceReconstructionImpl::init_rand(void)
{
X[1] = 10000;
}
double PowerCrustSurfaceReconstructionImpl::logb(double x)
{
return log(x) / log(2.0); // EPRO added
}
REAL PowerCrustSurfaceReconstructionImpl::orient3dadapt(REAL *pa, REAL *pb, REAL *pc, REAL *pd, REAL permanent)
{
INEXACT REAL adx, bdx, cdx, ady, bdy, cdy, adz, bdz, cdz;
REAL det, errbound;
INEXACT REAL bdxcdy1, cdxbdy1, cdxady1, adxcdy1, adxbdy1, bdxady1;
REAL bdxcdy0, cdxbdy0, cdxady0, adxcdy0, adxbdy0, bdxady0;
REAL bc[4], ca[4], ab[4];
INEXACT REAL bc3, ca3, ab3;
REAL adet[8], bdet[8], cdet[8];
int alen, blen, clen;
REAL abdet[16];
int ablen;
REAL *finnow, *finother, *finswap;
REAL fin1[192], fin2[192];
int finlength;
REAL adxtail, bdxtail, cdxtail;
REAL adytail, bdytail, cdytail;
REAL adztail, bdztail, cdztail;
INEXACT REAL at_blarge, at_clarge;
INEXACT REAL bt_clarge, bt_alarge;
INEXACT REAL ct_alarge, ct_blarge;
REAL at_b[4], at_c[4], bt_c[4], bt_a[4], ct_a[4], ct_b[4];
int at_blen, at_clen, bt_clen, bt_alen, ct_alen, ct_blen;
INEXACT REAL bdxt_cdy1, cdxt_bdy1, cdxt_ady1;
INEXACT REAL adxt_cdy1, adxt_bdy1, bdxt_ady1;
REAL bdxt_cdy0, cdxt_bdy0, cdxt_ady0;
REAL adxt_cdy0, adxt_bdy0, bdxt_ady0;
INEXACT REAL bdyt_cdx1, cdyt_bdx1, cdyt_adx1;
INEXACT REAL adyt_cdx1, adyt_bdx1, bdyt_adx1;
REAL bdyt_cdx0, cdyt_bdx0, cdyt_adx0;
REAL adyt_cdx0, adyt_bdx0, bdyt_adx0;
REAL bct[8], cat[8], abt[8];
int bctlen, catlen, abtlen;
INEXACT REAL bdxt_cdyt1, cdxt_bdyt1, cdxt_adyt1;
INEXACT REAL adxt_cdyt1, adxt_bdyt1, bdxt_adyt1;
REAL bdxt_cdyt0, cdxt_bdyt0, cdxt_adyt0;
REAL adxt_cdyt0, adxt_bdyt0, bdxt_adyt0;
REAL u[4], v[12], w[16];
INEXACT REAL u3;
int vlength, wlength;
REAL negate;
INEXACT REAL bvirt;
REAL avirt, bround, around;
INEXACT REAL c;
INEXACT REAL abig;
REAL ahi, alo, bhi, blo;
REAL err1, err2, err3;
INEXACT REAL _i, _j, _k;
REAL _0;
adx = (REAL)(pa[0] - pd[0]);
bdx = (REAL)(pb[0] - pd[0]);
cdx = (REAL)(pc[0] - pd[0]);
ady = (REAL)(pa[1] - pd[1]);
bdy = (REAL)(pb[1] - pd[1]);
cdy = (REAL)(pc[1] - pd[1]);
adz = (REAL)(pa[2] - pd[2]);
bdz = (REAL)(pb[2] - pd[2]);
cdz = (REAL)(pc[2] - pd[2]);
Two_Product(bdx, cdy, bdxcdy1, bdxcdy0);
Two_Product(cdx, bdy, cdxbdy1, cdxbdy0);
Two_Two_Diff(bdxcdy1, bdxcdy0, cdxbdy1, cdxbdy0, bc3, bc[2], bc[1], bc[0]);
bc[3] = bc3;
alen = scale_expansion_zeroelim(4, bc, adz, adet);
Two_Product(cdx, ady, cdxady1, cdxady0);
Two_Product(adx, cdy, adxcdy1, adxcdy0);
Two_Two_Diff(cdxady1, cdxady0, adxcdy1, adxcdy0, ca3, ca[2], ca[1], ca[0]);
ca[3] = ca3;
blen = scale_expansion_zeroelim(4, ca, bdz, bdet);
Two_Product(adx, bdy, adxbdy1, adxbdy0);
Two_Product(bdx, ady, bdxady1, bdxady0);
Two_Two_Diff(adxbdy1, adxbdy0, bdxady1, bdxady0, ab3, ab[2], ab[1], ab[0]);
ab[3] = ab3;
clen = scale_expansion_zeroelim(4, ab, cdz, cdet);
ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet);
finlength = fast_expansion_sum_zeroelim(ablen, abdet, clen, cdet, fin1);
det = estimate(finlength, fin1);
errbound = o3derrboundB * permanent;
if ((det >= errbound) || (-det >= errbound))
{
return det;
}
Two_Diff_Tail(pa[0], pd[0], adx, adxtail);
Two_Diff_Tail(pb[0], pd[0], bdx, bdxtail);
Two_Diff_Tail(pc[0], pd[0], cdx, cdxtail);
Two_Diff_Tail(pa[1], pd[1], ady, adytail);
Two_Diff_Tail(pb[1], pd[1], bdy, bdytail);
Two_Diff_Tail(pc[1], pd[1], cdy, cdytail);
Two_Diff_Tail(pa[2], pd[2], adz, adztail);
Two_Diff_Tail(pb[2], pd[2], bdz, bdztail);
Two_Diff_Tail(pc[2], pd[2], cdz, cdztail);
if ((adxtail == 0.0) && (bdxtail == 0.0) && (cdxtail == 0.0)
&& (adytail == 0.0) && (bdytail == 0.0) && (cdytail == 0.0)
&& (adztail == 0.0) && (bdztail == 0.0) && (cdztail == 0.0))
{
return det;
}
errbound = o3derrboundC * permanent + resulterrbound * Absolute(det);
det += (adz * ((bdx * cdytail + cdy * bdxtail)
- (bdy * cdxtail + cdx * bdytail))
+ adztail * (bdx * cdy - bdy * cdx))
+ (bdz * ((cdx * adytail + ady * cdxtail)
- (cdy * adxtail + adx * cdytail))
+ bdztail * (cdx * ady - cdy * adx))
+ (cdz * ((adx * bdytail + bdy * adxtail)
- (ady * bdxtail + bdx * adytail))
+ cdztail * (adx * bdy - ady * bdx));
if ((det >= errbound) || (-det >= errbound))
{
return det;
}
finnow = fin1;
finother = fin2;
if (adxtail == 0.0)
{
if (adytail == 0.0)
{
at_b[0] = 0.0;
at_blen = 1;
at_c[0] = 0.0;
at_clen = 1;
}
else
{
negate = -adytail;
Two_Product(negate, bdx, at_blarge, at_b[0]);
at_b[1] = at_blarge;
at_blen = 2;
Two_Product(adytail, cdx, at_clarge, at_c[0]);
at_c[1] = at_clarge;
at_clen = 2;
}
}
else
{
if (adytail == 0.0)
{
Two_Product(adxtail, bdy, at_blarge, at_b[0]);
at_b[1] = at_blarge;
at_blen = 2;
negate = -adxtail;
Two_Product(negate, cdy, at_clarge, at_c[0]);
at_c[1] = at_clarge;
at_clen = 2;
}
else
{
Two_Product(adxtail, bdy, adxt_bdy1, adxt_bdy0);
Two_Product(adytail, bdx, adyt_bdx1, adyt_bdx0);
Two_Two_Diff(adxt_bdy1, adxt_bdy0, adyt_bdx1, adyt_bdx0,
at_blarge, at_b[2], at_b[1], at_b[0]);
at_b[3] = at_blarge;
at_blen = 4;
Two_Product(adytail, cdx, adyt_cdx1, adyt_cdx0);
Two_Product(adxtail, cdy, adxt_cdy1, adxt_cdy0);
Two_Two_Diff(adyt_cdx1, adyt_cdx0, adxt_cdy1, adxt_cdy0,
at_clarge, at_c[2], at_c[1], at_c[0]);
at_c[3] = at_clarge;
at_clen = 4;
}
}
if (bdxtail == 0.0)
{
if (bdytail == 0.0)
{
bt_c[0] = 0.0;
bt_clen = 1;
bt_a[0] = 0.0;
bt_alen = 1;
}
else
{
negate = -bdytail;
Two_Product(negate, cdx, bt_clarge, bt_c[0]);
bt_c[1] = bt_clarge;
bt_clen = 2;
Two_Product(bdytail, adx, bt_alarge, bt_a[0]);
bt_a[1] = bt_alarge;
bt_alen = 2;
}
}
else
{
if (bdytail == 0.0)
{
Two_Product(bdxtail, cdy, bt_clarge, bt_c[0]);
bt_c[1] = bt_clarge;
bt_clen = 2;
negate = -bdxtail;
Two_Product(negate, ady, bt_alarge, bt_a[0]);
bt_a[1] = bt_alarge;
bt_alen = 2;
}
else
{
Two_Product(bdxtail, cdy, bdxt_cdy1, bdxt_cdy0);
Two_Product(bdytail, cdx, bdyt_cdx1, bdyt_cdx0);
Two_Two_Diff(bdxt_cdy1, bdxt_cdy0, bdyt_cdx1, bdyt_cdx0,
bt_clarge, bt_c[2], bt_c[1], bt_c[0]);
bt_c[3] = bt_clarge;
bt_clen = 4;
Two_Product(bdytail, adx, bdyt_adx1, bdyt_adx0);
Two_Product(bdxtail, ady, bdxt_ady1, bdxt_ady0);
Two_Two_Diff(bdyt_adx1, bdyt_adx0, bdxt_ady1, bdxt_ady0,
bt_alarge, bt_a[2], bt_a[1], bt_a[0]);
bt_a[3] = bt_alarge;
bt_alen = 4;
}
}
if (cdxtail == 0.0)
{
if (cdytail == 0.0)
{
ct_a[0] = 0.0;
ct_alen = 1;
ct_b[0] = 0.0;
ct_blen = 1;
}
else
{
negate = -cdytail;
Two_Product(negate, adx, ct_alarge, ct_a[0]);
ct_a[1] = ct_alarge;
ct_alen = 2;
Two_Product(cdytail, bdx, ct_blarge, ct_b[0]);
ct_b[1] = ct_blarge;
ct_blen = 2;
}
}
else
{
if (cdytail == 0.0)
{
Two_Product(cdxtail, ady, ct_alarge, ct_a[0]);
ct_a[1] = ct_alarge;
ct_alen = 2;
negate = -cdxtail;
Two_Product(negate, bdy, ct_blarge, ct_b[0]);
ct_b[1] = ct_blarge;
ct_blen = 2;
}
else
{
Two_Product(cdxtail, ady, cdxt_ady1, cdxt_ady0);
Two_Product(cdytail, adx, cdyt_adx1, cdyt_adx0);
Two_Two_Diff(cdxt_ady1, cdxt_ady0, cdyt_adx1, cdyt_adx0,
ct_alarge, ct_a[2], ct_a[1], ct_a[0]);
ct_a[3] = ct_alarge;
ct_alen = 4;
Two_Product(cdytail, bdx, cdyt_bdx1, cdyt_bdx0);
Two_Product(cdxtail, bdy, cdxt_bdy1, cdxt_bdy0);
Two_Two_Diff(cdyt_bdx1, cdyt_bdx0, cdxt_bdy1, cdxt_bdy0,
ct_blarge, ct_b[2], ct_b[1], ct_b[0]);
ct_b[3] = ct_blarge;
ct_blen = 4;
}
}
bctlen = fast_expansion_sum_zeroelim(bt_clen, bt_c, ct_blen, ct_b, bct);
wlength = scale_expansion_zeroelim(bctlen, bct, adz, w);
finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
finother);
finswap = finnow;
finnow = finother;
finother = finswap;
catlen = fast_expansion_sum_zeroelim(ct_alen, ct_a, at_clen, at_c, cat);
wlength = scale_expansion_zeroelim(catlen, cat, bdz, w);
finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
finother);
finswap = finnow;
finnow = finother;
finother = finswap;
abtlen = fast_expansion_sum_zeroelim(at_blen, at_b, bt_alen, bt_a, abt);
wlength = scale_expansion_zeroelim(abtlen, abt, cdz, w);
finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
finother);
finswap = finnow;
finnow = finother;
finother = finswap;
if (adztail != 0.0)
{
vlength = scale_expansion_zeroelim(4, bc, adztail, v);
finlength = fast_expansion_sum_zeroelim(finlength, finnow, vlength, v,
finother);
finswap = finnow;
finnow = finother;
finother = finswap;
}
if (bdztail != 0.0)
{
vlength = scale_expansion_zeroelim(4, ca, bdztail, v);
finlength = fast_expansion_sum_zeroelim(finlength, finnow, vlength, v,
finother);
finswap = finnow;
finnow = finother;
finother = finswap;
}
if (cdztail != 0.0)
{
vlength = scale_expansion_zeroelim(4, ab, cdztail, v);
finlength = fast_expansion_sum_zeroelim(finlength, finnow, vlength, v,
finother);
finswap = finnow;
finnow = finother;
finother = finswap;
}
if (adxtail != 0.0)
{
if (bdytail != 0.0)
{
Two_Product(adxtail, bdytail, adxt_bdyt1, adxt_bdyt0);
Two_One_Product(adxt_bdyt1, adxt_bdyt0, cdz, u3, u[2], u[1], u[0]);
u[3] = u3;
finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
finother);
finswap = finnow;
finnow = finother;
finother = finswap;
if (cdztail != 0.0)
{
Two_One_Product(adxt_bdyt1, adxt_bdyt0, cdztail, u3, u[2], u[1], u[0]);
u[3] = u3;
finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
finother);
finswap = finnow;
finnow = finother;
finother = finswap;
}
}
if (cdytail != 0.0)
{
negate = -adxtail;
Two_Product(negate, cdytail, adxt_cdyt1, adxt_cdyt0);
Two_One_Product(adxt_cdyt1, adxt_cdyt0, bdz, u3, u[2], u[1], u[0]);
u[3] = u3;
finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
finother);
finswap = finnow;
finnow = finother;
finother = finswap;
if (bdztail != 0.0)
{
Two_One_Product(adxt_cdyt1, adxt_cdyt0, bdztail, u3, u[2], u[1], u[0]);
u[3] = u3;
finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
finother);
finswap = finnow;
finnow = finother;
finother = finswap;
}
}
}
if (bdxtail != 0.0)
{
if (cdytail != 0.0)
{
Two_Product(bdxtail, cdytail, bdxt_cdyt1, bdxt_cdyt0);
Two_One_Product(bdxt_cdyt1, bdxt_cdyt0, adz, u3, u[2], u[1], u[0]);
u[3] = u3;
finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
finother);
finswap = finnow;
finnow = finother;
finother = finswap;
if (adztail != 0.0)
{
Two_One_Product(bdxt_cdyt1, bdxt_cdyt0, adztail, u3, u[2], u[1], u[0]);
u[3] = u3;
finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
finother);
finswap = finnow;
finnow = finother;
finother = finswap;
}
}
if (adytail != 0.0)
{
negate = -bdxtail;
Two_Product(negate, adytail, bdxt_adyt1, bdxt_adyt0);
Two_One_Product(bdxt_adyt1, bdxt_adyt0, cdz, u3, u[2], u[1], u[0]);
u[3] = u3;
finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
finother);
finswap = finnow;
finnow = finother;
finother = finswap;
if (cdztail != 0.0)
{
Two_One_Product(bdxt_adyt1, bdxt_adyt0, cdztail, u3, u[2], u[1], u[0]);
u[3] = u3;
finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
finother);
finswap = finnow;
finnow = finother;
finother = finswap;
}
}
}
if (cdxtail != 0.0)
{
if (adytail != 0.0)
{
Two_Product(cdxtail, adytail, cdxt_adyt1, cdxt_adyt0);
Two_One_Product(cdxt_adyt1, cdxt_adyt0, bdz, u3, u[2], u[1], u[0]);
u[3] = u3;
finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
finother);
finswap = finnow;
finnow = finother;
finother = finswap;
if (bdztail != 0.0)
{
Two_One_Product(cdxt_adyt1, cdxt_adyt0, bdztail, u3, u[2], u[1], u[0]);
u[3] = u3;
finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
finother);
finswap = finnow;
finnow = finother;
finother = finswap;
}
}
if (bdytail != 0.0)
{
negate = -cdxtail;
Two_Product(negate, bdytail, cdxt_bdyt1, cdxt_bdyt0);
Two_One_Product(cdxt_bdyt1, cdxt_bdyt0, adz, u3, u[2], u[1], u[0]);
u[3] = u3;
finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
finother);
finswap = finnow;
finnow = finother;
finother = finswap;
if (adztail != 0.0)
{
Two_One_Product(cdxt_bdyt1, cdxt_bdyt0, adztail, u3, u[2], u[1], u[0]);
u[3] = u3;
finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
finother);
finswap = finnow;
finnow = finother;
finother = finswap;
}
}
}
if (adztail != 0.0)
{
wlength = scale_expansion_zeroelim(bctlen, bct, adztail, w);
finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
finother);
finswap = finnow;
finnow = finother;
finother = finswap;
}
if (bdztail != 0.0)
{
wlength = scale_expansion_zeroelim(catlen, cat, bdztail, w);
finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
finother);
finswap = finnow;
finnow = finother;
finother = finswap;
}
if (cdztail != 0.0)
{
wlength = scale_expansion_zeroelim(abtlen, abt, cdztail, w);
finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
finother);
finswap = finnow;
finnow = finother;
finother = finswap;
}
return finnow[finlength - 1];
}
REAL PowerCrustSurfaceReconstructionImpl::orient3d(REAL *pa, REAL *pb, REAL *pc, REAL *pd)
{
REAL adx, bdx, cdx, ady, bdy, cdy, adz, bdz, cdz;
REAL bdxcdy, cdxbdy, cdxady, adxcdy, adxbdy, bdxady;
REAL det;
REAL permanent, errbound;
adx = pa[0] - pd[0];
bdx = pb[0] - pd[0];
cdx = pc[0] - pd[0];
ady = pa[1] - pd[1];
bdy = pb[1] - pd[1];
cdy = pc[1] - pd[1];
adz = pa[2] - pd[2];
bdz = pb[2] - pd[2];
cdz = pc[2] - pd[2];
bdxcdy = bdx * cdy;
cdxbdy = cdx * bdy;
cdxady = cdx * ady;
adxcdy = adx * cdy;
adxbdy = adx * bdy;
bdxady = bdx * ady;
det = adz * (bdxcdy - cdxbdy)
+ bdz * (cdxady - adxcdy)
+ cdz * (adxbdy - bdxady);
permanent = (Absolute(bdxcdy) + Absolute(cdxbdy)) * Absolute(adz)
+ (Absolute(cdxady) + Absolute(adxcdy)) * Absolute(bdz)
+ (Absolute(adxbdy) + Absolute(bdxady)) * Absolute(cdz);
errbound = o3derrboundA * permanent;
if ((det > errbound) || (-det > errbound))
{
return det;
}
return orient3dadapt(pa, pb, pc, pd, permanent);
}
REAL PowerCrustSurfaceReconstructionImpl::estimate(int elen, REAL *e)
{
REAL Q;
int eindex;
Q = e[0];
for (eindex = 1; eindex < elen; eindex++)
{
Q += e[eindex];
}
return Q;
}
REAL PowerCrustSurfaceReconstructionImpl::orient2dadapt(REAL *pa, REAL *pb, REAL *pc, REAL detsum)
{
INEXACT REAL acx, acy, bcx, bcy;
REAL acxtail, acytail, bcxtail, bcytail;
INEXACT REAL detleft, detright;
REAL detlefttail, detrighttail;
REAL det, errbound;
REAL B[4], C1[8], C2[12], D[16];
INEXACT REAL B3;
int C1length, C2length, Dlength;
REAL u[4];
INEXACT REAL u3;
INEXACT REAL s1, t1;
REAL s0, t0;
INEXACT REAL bvirt;
REAL avirt, bround, around;
INEXACT REAL c;
INEXACT REAL abig;
REAL ahi, alo, bhi, blo;
REAL err1, err2, err3;
INEXACT REAL _i, _j;
REAL _0;
acx = (REAL)(pa[0] - pc[0]);
bcx = (REAL)(pb[0] - pc[0]);
acy = (REAL)(pa[1] - pc[1]);
bcy = (REAL)(pb[1] - pc[1]);
Two_Product(acx, bcy, detleft, detlefttail);
Two_Product(acy, bcx, detright, detrighttail);
Two_Two_Diff(detleft, detlefttail, detright, detrighttail,
B3, B[2], B[1], B[0]);
B[3] = B3;
det = estimate(4, B);
errbound = ccwerrboundB * detsum;
if ((det >= errbound) || (-det >= errbound))
{
return det;
}
Two_Diff_Tail(pa[0], pc[0], acx, acxtail);
Two_Diff_Tail(pb[0], pc[0], bcx, bcxtail);
Two_Diff_Tail(pa[1], pc[1], acy, acytail);
Two_Diff_Tail(pb[1], pc[1], bcy, bcytail);
if ((acxtail == 0.0) && (acytail == 0.0)
&& (bcxtail == 0.0) && (bcytail == 0.0))
{
return det;
}
errbound = ccwerrboundC * detsum + resulterrbound * Absolute(det);
det += (acx * bcytail + bcy * acxtail)
- (acy * bcxtail + bcx * acytail);
if ((det >= errbound) || (-det >= errbound))
{
return det;
}
Two_Product(acxtail, bcy, s1, s0);
Two_Product(acytail, bcx, t1, t0);
Two_Two_Diff(s1, s0, t1, t0, u3, u[2], u[1], u[0]);
u[3] = u3;
C1length = fast_expansion_sum_zeroelim(4, B, 4, u, C1);
Two_Product(acx, bcytail, s1, s0);
Two_Product(acy, bcxtail, t1, t0);
Two_Two_Diff(s1, s0, t1, t0, u3, u[2], u[1], u[0]);
u[3] = u3;
C2length = fast_expansion_sum_zeroelim(C1length, C1, 4, u, C2);
Two_Product(acxtail, bcytail, s1, s0);
Two_Product(acytail, bcxtail, t1, t0);
Two_Two_Diff(s1, s0, t1, t0, u3, u[2], u[1], u[0]);
u[3] = u3;
Dlength = fast_expansion_sum_zeroelim(C2length, C2, 4, u, D);
return (D[Dlength - 1]);
}
int PowerCrustSurfaceReconstructionImpl::scale_expansion_zeroelim(int elen, REAL *e, REAL b, REAL *h)
{
INEXACT REAL Q, sum;
REAL hh;
INEXACT REAL product1;
REAL product0;
int eindex, hindex;
REAL enow;
INEXACT REAL bvirt;
REAL avirt, bround, around;
INEXACT REAL c;
INEXACT REAL abig;
REAL ahi, alo, bhi, blo;
REAL err1, err2, err3;
Split(b, bhi, blo);
Two_Product_Presplit(e[0], b, bhi, blo, Q, hh);
hindex = 0;
if (hh != 0)
{
h[hindex++] = hh;
}
for (eindex = 1; eindex < elen; eindex++)
{
enow = e[eindex];
Two_Product_Presplit(enow, b, bhi, blo, product1, product0);
Two_Sum(Q, product0, sum, hh);
if (hh != 0)
{
h[hindex++] = hh;
}
Fast_Two_Sum(product1, sum, Q, hh);
if (hh != 0)
{
h[hindex++] = hh;
}
}
if ((Q != 0.0) || (hindex == 0))
{
h[hindex++] = Q;
}
return hindex;
}
int PowerCrustSurfaceReconstructionImpl::fast_expansion_sum_zeroelim(int elen, REAL *e, int flen, REAL *f, REAL *h)
{
REAL Q;
INEXACT REAL Qnew;
INEXACT REAL hh;
INEXACT REAL bvirt;
REAL avirt, bround, around;
int eindex, findex, hindex;
REAL enow, fnow;
enow = e[0];
fnow = f[0];
eindex = findex = 0;
if ((fnow > enow) == (fnow > -enow))
{
Q = enow;
enow = e[++eindex];
}
else
{
Q = fnow;
fnow = f[++findex];
}
hindex = 0;
if ((eindex < elen) && (findex < flen))
{
if ((fnow > enow) == (fnow > -enow))
{
Fast_Two_Sum(enow, Q, Qnew, hh);
enow = e[++eindex];
}
else
{
Fast_Two_Sum(fnow, Q, Qnew, hh);
fnow = f[++findex];
}
Q = Qnew;
if (hh != 0.0)
{
h[hindex++] = hh;
}
while ((eindex < elen) && (findex < flen))
{
if ((fnow > enow) == (fnow > -enow))
{
Two_Sum(Q, enow, Qnew, hh);
enow = e[++eindex];
}
else
{
Two_Sum(Q, fnow, Qnew, hh);
fnow = f[++findex];
}
Q = Qnew;
if (hh != 0.0)
{
h[hindex++] = hh;
}
}
}
while (eindex < elen)
{
Two_Sum(Q, enow, Qnew, hh);
enow = e[++eindex];
Q = Qnew;
if (hh != 0.0)
{
h[hindex++] = hh;
}
}
while (findex < flen)
{
Two_Sum(Q, fnow, Qnew, hh);
fnow = f[++findex];
Q = Qnew;
if (hh != 0.0)
{
h[hindex++] = hh;
}
}
if ((Q != 0.0) || (hindex == 0))
{
h[hindex++] = Q;
}
return hindex;
}
void PowerCrustSurfaceReconstructionImpl::exactinit()
{
REAL half;
REAL check, lastcheck;
int every_other;
every_other = 1;
half = 0.5;
epsilon = 1.0;
splitter = 1.0;
check = 1.0;
/* Repeatedly divide `epsilon' by two until it is too small to add to */
/* one without causing roundoff. (Also check if the sum is equal to */
/* the previous sum, for machines that round up instead of using exact */
/* rounding. Not that this library will work on such machines anyway. */
do
{
lastcheck = check;
epsilon *= half;
if (every_other)
{
splitter *= 2.0;
}
every_other = !every_other;
check = 1.0 + epsilon;
} while ((check != 1.0) && (check != lastcheck));
splitter += 1.0;
/* Error bounds for orientation and incircle tests. */
resulterrbound = (3.0 + 8.0 * epsilon) * epsilon;
ccwerrboundA = (3.0 + 16.0 * epsilon) * epsilon;
ccwerrboundB = (2.0 + 12.0 * epsilon) * epsilon;
ccwerrboundC = (9.0 + 64.0 * epsilon) * epsilon * epsilon;
o3derrboundA = (7.0 + 56.0 * epsilon) * epsilon;
o3derrboundB = (3.0 + 28.0 * epsilon) * epsilon;
o3derrboundC = (26.0 + 288.0 * epsilon) * epsilon * epsilon;
}
void PowerCrustSurfaceReconstructionImpl::tetorthocenter(double a[4], double b[4], double c[4], double d[4], double orthocenter[3], double *cnum)
{
double xba, yba, zba, xca, yca, zca, xda, yda, zda, wba, wca, wda;
double balength, calength, dalength;
double xcrosscd, ycrosscd, zcrosscd;
double xcrossdb, ycrossdb, zcrossdb;
double xcrossbc, ycrossbc, zcrossbc;
double denominator;
double xcirca, ycirca, zcirca;
double wa, wb, wc, wd;
wa = a[0] * a[0] + a[1] * a[1] + a[2] * a[2] - a[3];
wb = b[0] * b[0] + b[1] * b[1] + b[2] * b[2] - b[3];
wc = c[0] * c[0] + c[1] * c[1] + c[2] * c[2] - c[3];
wd = d[0] * d[0] + d[1] * d[1] + d[2] * d[2] - d[3];
/* Use coordinates relative to point `a' of the tetrahedron. */
xba = b[0] - a[0];
yba = b[1] - a[1];
zba = b[2] - a[2];
wba = wb - wa;
xca = c[0] - a[0];
yca = c[1] - a[1];
zca = c[2] - a[2];
wca = wc - wa;
xda = d[0] - a[0];
yda = d[1] - a[1];
zda = d[2] - a[2];
wda = wd - wa;
/* Squares of lengths of the edges incident to `a'. */
balength = xba * xba + yba * yba + zba * zba - wba;
calength = xca * xca + yca * yca + zca * zca - wca;
dalength = xda * xda + yda * yda + zda * zda - wda;
/* Cross products of these edges. */
xcrosscd = yca * zda - yda * zca;
ycrosscd = zca * xda - zda * xca;
zcrosscd = xca * yda - xda * yca;
xcrossdb = yda * zba - yba * zda;
ycrossdb = zda * xba - zba * xda;
zcrossdb = xda * yba - xba * yda;
xcrossbc = yba * zca - yca * zba;
ycrossbc = zba * xca - zca * xba;
zcrossbc = xba * yca - xca * yba;
/* Calculate the denominator of the formulae. */
#ifdef EXACT
/* Use orient3d() from http://www.cs.cmu.edu/~quake/robust.html */
/* to ensure a correctly signed (and reasonably accurate) result, */
/* avoiding any possibility of division by zero. */
*cnum = orient3d(b, c, d, a);
denominator = 0.5 / (*cnum);
#else
/* Take your chances with floating-point roundoff. */
denominator = 0.5 / (xba * xcrosscd + yba * ycrosscd + zba * zcrosscd);
#endif
/* Calculate offset (from `a') of circumcenter. */
xcirca = (balength * xcrosscd + calength * xcrossdb + dalength * xcrossbc) *
denominator;
ycirca = (balength * ycrosscd + calength * ycrossdb + dalength * ycrossbc) *
denominator;
zcirca = (balength * zcrosscd + calength * zcrossdb + dalength * zcrossbc) *
denominator;
orthocenter[0] = xcirca;
orthocenter[1] = ycirca;
orthocenter[2] = zcirca;
}
void PowerCrustSurfaceReconstructionImpl::tetcircumcenter(double a[3], double b[3], double c[3], double d[3], double circumcenter[3], double *cond)
{
double xba, yba, zba, xca, yca, zca, xda, yda, zda;
double balength, calength, dalength;
double xcrosscd, ycrosscd, zcrosscd;
double xcrossdb, ycrossdb, zcrossdb;
double xcrossbc, ycrossbc, zcrossbc;
double denominator;
double xcirca, ycirca, zcirca;
/* Use coordinates relative to point `a' of the tetrahedron. */
xba = b[0] - a[0];
yba = b[1] - a[1];
zba = b[2] - a[2];
xca = c[0] - a[0];
yca = c[1] - a[1];
zca = c[2] - a[2];
xda = d[0] - a[0];
yda = d[1] - a[1];
zda = d[2] - a[2];
/* Squares of lengths of the edges incident to `a'. */
balength = xba * xba + yba * yba + zba * zba;
calength = xca * xca + yca * yca + zca * zca;
dalength = xda * xda + yda * yda + zda * zda;
/* Cross products of these edges. */
xcrosscd = yca * zda - yda * zca;
ycrosscd = zca * xda - zda * xca;
zcrosscd = xca * yda - xda * yca;
xcrossdb = yda * zba - yba * zda;
ycrossdb = zda * xba - zba * xda;
zcrossdb = xda * yba - xba * yda;
xcrossbc = yba * zca - yca * zba;
ycrossbc = zba * xca - zca * xba;
zcrossbc = xba * yca - xca * yba;
/* Calculate the denominator of the formulae. */
#ifdef EXACT
/* Use orient3d() from http://www.cs.cmu.edu/~quake/robust.html */
/* to ensure a correctly signed (and reasonably accurate) result, */
/* avoiding any possibility of division by zero. */
*cond = orient3d(b, c, d, a);
denominator = 0.5 / (*cond);
#else
/* Take your chances with floating-point roundoff. */
denominator = 0.5 / (xba * xcrosscd + yba * ycrosscd + zba * zcrosscd);
#endif
/* Calculate offset (from `a') of circumcenter. */
xcirca = (balength * xcrosscd + calength * xcrossdb + dalength * xcrossbc) *
denominator;
ycirca = (balength * ycrosscd + calength * ycrossdb + dalength * ycrossbc) *
denominator;
zcirca = (balength * zcrosscd + calength * zcrossdb + dalength * zcrossbc) *
denominator;
circumcenter[0] = xcirca;
circumcenter[1] = ycirca;
circumcenter[2] = zcirca;
}
double PowerCrustSurfaceReconstructionImpl::sqdist(double a[3], double b[3])
{
/* returns the squared distance between a and b */
return SQ(a[0] - b[0]) + SQ(a[1] - b[1]) + SQ(a[2] - b[2]);
}
void PowerCrustSurfaceReconstructionImpl::dir_and_dist(double a[3], double b[3], double dir[3], double* dist)
{
int k;
for (k = 0; k < 3; k++) dir[k] = b[k] - a[k];
*dist = sqrt(SQ(dir[0]) + SQ(dir[1]) + SQ(dir[2]));
for (k = 0; k<3; k++) dir[k] = dir[k] / (*dist);
}
int PowerCrustSurfaceReconstructionImpl::propagate()
{
int pid = extract_max();
if (adjlist[pid].in > adjlist[pid].out)
adjlist[pid].label = IN;
else adjlist[pid].label = OUT;
if (pid != -1)
{
opp_update(pid);
sym_update(pid);
}
return pid;
}
void PowerCrustSurfaceReconstructionImpl::opp_update(int pi)
{
plist *pindex;
pindex = opplist[pi];
while (pindex != NULL)
{
int npi = pindex->pid;
if (adjlist[npi].label == INIT) /* not yet labeled */
{
if (adjlist[npi].hid == 0) /* not in the heap */
{
if (adjlist[pi].in > adjlist[pi].out)
{
/* propagate in*cos to out */
adjlist[npi].out = (-1.0) * adjlist[pi].in * pindex->angle;
insert_heap(npi, adjlist[npi].out);
}
else if (adjlist[pi].in < adjlist[pi].out)
{
/* propagate out*cos to in */
adjlist[npi].in = (-1.0) * adjlist[pi].out * pindex->angle;
insert_heap(npi, adjlist[npi].in);
}
}
else /* in the heap */
{
int nhi = adjlist[npi].hid;
if (adjlist[pi].in > adjlist[pi].out)
{
/* propagate in*cos to out */
double temp = (-1.0) * adjlist[pi].in * pindex->angle;
if (temp > adjlist[npi].out)
{
adjlist[npi].out = temp;
update_pri(nhi, npi);
}
}
else if (adjlist[pi].in < adjlist[pi].out)
{
/* propagate out*cos to in */
double temp = (-1.0) * adjlist[pi].out * pindex->angle;
if (temp > adjlist[npi].in)
{
adjlist[npi].in = temp;
update_pri(nhi, npi);
}
}
}
}
pindex = pindex->next;
}
}
void PowerCrustSurfaceReconstructionImpl::sym_update(int pi)
{
edgesimp *eindex;
eindex = adjlist[pi].eptr;
while (eindex != NULL)
{
int npi = eindex->pid;
/* try to label deeply intersecting unlabeled neighbors */
if ((adjlist[npi].label == INIT) && (eindex->angle > theta))
{
/* not yet labeled */
if (adjlist[npi].hid == 0) /* not in the heap */
{
if (adjlist[pi].in > adjlist[pi].out)
{
/* propagate in*cos to in */
adjlist[npi].in = adjlist[pi].in * eindex->angle;
insert_heap(npi, adjlist[npi].in);
}
else if (adjlist[pi].in < adjlist[pi].out)
{
/* propagate out*cos to out */
adjlist[npi].out = adjlist[pi].out * eindex->angle;
insert_heap(npi, adjlist[npi].out);
}
}
else /* in the heap */
{
int nhi = adjlist[npi].hid;
if (adjlist[pi].in > adjlist[pi].out)
{
/* propagate in*cos to in */
double temp = adjlist[pi].in * eindex->angle;
if (temp > adjlist[npi].in)
{
adjlist[npi].in = temp;
update_pri(nhi, npi);
}
}
else if (adjlist[pi].in < adjlist[pi].out)
{
/* propagate out*cos to out */
double temp = adjlist[pi].out * eindex->angle;
if (temp > adjlist[npi].out)
{
adjlist[npi].out = temp;
update_pri(nhi, npi);
}
}
}
}
eindex = eindex->next;
}
}
void PowerCrustSurfaceReconstructionImpl::update_pri(int hi, int pi)
{
double pr;
if ((heap_A[hi].pid != pi) || (adjlist[pi].hid != hi))
{
return;
}
if (adjlist[pi].in == 0.0)
{
pr = adjlist[pi].out;
}
else if (adjlist[pi].out == 0.0)
{
pr = adjlist[pi].in;
}
else /* both in/out nonzero */
{
if (adjlist[pi].in > adjlist[pi].out)
{
pr = adjlist[pi].in - adjlist[pi].out - 1;
}
else
{
pr = adjlist[pi].out - adjlist[pi].in - 1;
}
}
update(hi, pr);
}
void PowerCrustSurfaceReconstructionImpl::label_unlabeled(int num)
{
plist *pindex;
edgesimp *eindex;
int opplabel;
for (int i = 0; i < num; i++)
{
if (adjlist[i].label == INIT) /* pole i is unlabeled.. try to label now */
{
opplabel = INIT;
pindex = opplist[i];
if ((pindex == NULL) && (adjlist[i].eptr == NULL))
{
continue;
}
/* check whether there is opp pole */
while (pindex != NULL) /* opp pole */
{
int npi = pindex->pid;
if (adjlist[npi].label != INIT)
{
if (opplabel == INIT) opplabel = adjlist[npi].label;
else if (opplabel != adjlist[npi].label)
{
opplabel = INIT; /* ignore the label of opposite poles */
}
}
pindex = pindex->next;
}
double tangle = -3.0;
double tangle1 = -3.0;
eindex = adjlist[i].eptr;
while (eindex != NULL)
{
int npi = eindex->pid;
if (adjlist[npi].label == IN)
{
if (tangle < eindex->angle)
{
tangle = eindex->angle;
}
}
else if (adjlist[npi].label == OUT)
{
if (tangle1 < eindex->angle)
{
tangle1 = eindex->angle;
}
}
eindex = eindex->next;
}
/* now tangle, tangle 1 are angles of most deeply interesecting in, out poles */
if (tangle == -3.0) /* there was no in poles */
{
if (tangle1 == -3.0) /* there was no out poles */
{
if (opplabel == INIT) /* cannot trust opp pole or no labeled opp pole */
{
}
else if (opplabel == IN)
{
adjlist[i].label = OUT;
}
else adjlist[i].label = IN;
}
else if (tangle1 > deep) /* interesecting deeply only out poles */
{
adjlist[i].label = OUT;
}
else /* no deeply intersecting poles . use opp pole */
{
if (opplabel == INIT) /* cannot trust opp pole or no labeled opp pole */
{
}
else if (opplabel == IN)
{
adjlist[i].label = OUT;
}
else adjlist[i].label = IN;
}
}
else if (tangle1 == -3.0) /* there are in pole but no out pole */
{
if (tangle > deep) /* interesecting deeply only in poles */
{
adjlist[i].label = IN;
}
else /* no deeply intersecting poles . use opp pole */
{
if (opplabel == INIT) /* cannot trust opp pole or no labeled opp pole */
{
}
else if (opplabel == IN)
{
adjlist[i].label = OUT;
}
else adjlist[i].label = IN;
}
}
else /* there are both in/out poles */
{
if (tangle > deep)
{
if (tangle1 > deep) /* intersecting both deeply */
{
/* use opp */
if (opplabel == INIT) /* cannot trust opp pole or no labeled opp pole */
{
/* then give label with bigger angle */
if (tangle > tangle1)
{
adjlist[i].label = IN;
}
else adjlist[i].label = OUT;
}
else if (opplabel == IN)
{
adjlist[i].label = OUT;
}
else adjlist[i].label = IN;
}
else /* intersecting only in deeply */
{
adjlist[i].label = IN;
}
}
else if (tangle1 > deep) /* intersecting only out deeply */
{
adjlist[i].label = OUT;
}
else /* no deeply intersecting poles . use opp pole */
{
if (opplabel == INIT) /* cannot trust opp pole or no labeled opp pole */
{
}
else if (opplabel == IN)
{
adjlist[i].label = OUT;
}
else adjlist[i].label = IN;
}
}
}
}
}
neighbor * PowerCrustSurfaceReconstructionImpl::op_simp(simplex *a, simplex *b)
{
lookup(a, b, simp)
}
neighbor * PowerCrustSurfaceReconstructionImpl::op_vert(simplex *a, site b)
{
lookup(a, b, vert)
}
void PowerCrustSurfaceReconstructionImpl::connect(simplex *s)
{
/* make neighbor connections between newly created simplices incident to p */
site xf, xb, xfi;
simplex *sb, *sf, *seen;
int i;
neighbor *sn;
if (!s)
{
return;
}
assert(!s->peak.vert && s->peak.simp->peak.vert == p && !op_vert(s, p)->simp->peak.vert);
if (s->visit == pnum)
{
return;
}
s->visit = pnum;
seen = s->peak.simp;
xfi = op_simp(seen, s)->vert;
for (i = 0, sn = s->neigh; i < cdim; i++, sn++)
{
xb = sn->vert;
if (p == xb)
{
continue;
}
sb = seen;
sf = sn->simp;
xf = xfi;
if (!sf->peak.vert) /* are we done already? */
{
sf = op_vert(seen, xb)->simp;
if (sf->peak.vert)
{
continue;
}
}
else
{
int LoopCounter = 0;
do
{
if (LoopCounter == 100)
{
ASSERT(pcFALSE, "new adjacency failure!");
}
xb = xf;
xf = op_simp(sf, sb)->vert;
sb = sf;
sf = op_vert(sb, xb)->simp;
LoopCounter++;
} while (sf->peak.vert);
}
sn->simp = sf;
op_vert(sf, xf)->simp = s;
connect(sf);
}
}
simplex * PowerCrustSurfaceReconstructionImpl::make_facets(simplex *seen)
{
/* visit simplices s with sees(p,s), and make a facet for every neighbor
* of s not seen by p */
simplex *n;
neighbor *bn;
int i;
if (!seen) return NULL;
seen->peak.vert = p;
for (i = 0, bn = seen->neigh; i < cdim; i++, bn++)
{
n = bn->simp;
if (pnum != n->visit)
{
n->visit = pnum;
if (sees(p, n)) make_facets(n);
}
if (n->peak.vert) continue;
copy_simp(make_facets_ns, seen);
make_facets_ns->visit = 0;
make_facets_ns->peak.vert = 0;
make_facets_ns->normal = 0;
make_facets_ns->peak.simp = seen;
/* ns->Sb -= ns->neigh[i].basis->sqb; */
NULLIFY(basis_s, make_facets_ns->neigh[i].basis);
make_facets_ns->neigh[i].vert = p;
bn->simp = op_simp(n, seen)->simp = make_facets_ns;
}
return make_facets_ns;
}
simplex * PowerCrustSurfaceReconstructionImpl::extend_simplices(simplex *s)
{
/* p lies outside flat containing previous sites;
* make p a vertex of every current simplex, and create some new simplices */
int i;
int ocdim = cdim - 1;
simplex *ns;
neighbor *nsn;
if (s->visit == pnum) return s->peak.vert ? s->neigh[ocdim].simp : s;
s->visit = pnum;
s->neigh[ocdim].vert = p;
NULLIFY(basis_s, s->normal);
NULLIFY(basis_s, s->neigh[0].basis);
if (!s->peak.vert)
{
s->neigh[ocdim].simp = extend_simplices(s->peak.simp);
return s;
}
else
{
copy_simp(ns, s);
s->neigh[ocdim].simp = ns;
ns->peak.vert = NULL;
ns->peak.simp = s;
ns->neigh[ocdim] = s->peak;
inc_ref(basis_s, s->peak.basis);
for (i = 0, nsn = ns->neigh; i < cdim; i++, nsn++)
nsn->simp = extend_simplices(nsn->simp);
}
return ns;
}
simplex * PowerCrustSurfaceReconstructionImpl::search(simplex *root)
{
/* return a simplex s that corresponds to a facet of the
* current hull, and sees(p, s) */
simplex *s;
neighbor *sn;
int i;
long tms = 0;
if (!st_search)
{
st_search = (simplex **)malloc((search_ss + MAXDIM + 1) *sizeof(simplex*));
}
push(root->peak.simp, st_search, tms);
root->visit = pnum;
if (!sees(p, root))
for (i = 0, sn = root->neigh; i < cdim; i++, sn++)
push(sn->simp, st_search, tms);
while (tms)
{
if (tms > search_ss)
{
st_search = (simplex**)realloc(st_search, ((search_ss += search_ss) + MAXDIM + 1) *sizeof(simplex*));
assert(st_search);
}
pop(s, st_search, tms);
if (s->visit == pnum)
continue;
s->visit = pnum;
if (!sees(p, s))
continue;
if (!s->peak.vert)
{
return s;
}
for (i = 0, sn = s->neigh; i < cdim; i++, sn++)
push(sn->simp, st_search, tms);
}
return NULL;
}
point PowerCrustSurfaceReconstructionImpl::get_another_site(void)
{
point pnext;
if (!(++scount % 1000))
{
}
pnext = (this->*get_site) ();
if (!pnext)
{
return NULL;
}
pnum = (this->*site_num) (pnext)+2;
return pnext;
}
void PowerCrustSurfaceReconstructionImpl::buildhull(simplex *root)
{
while (cdim < rdim)
{
p = get_another_site();
if (!p)
{
return;
}
if (out_of_flat(root, p))
{
extend_simplices(root);
}
else
{
connect(make_facets(search(root)));
}
}
while ((p = get_another_site()) != NULL)
{
connect(make_facets(search(root)));
}
}
int PowerCrustSurfaceReconstructionImpl::reduce(basis_s **v, point p, simplex *s, int k)
{
point z;
point tt = s->neigh[0].vert;
if (!*v) NEWLRC(basis_s, (*v))
else (*v)->lscale = 0;
z = VB(*v);
if (vd || power_diagram)
{
if (p == infinity) memcpy(*v, infinity_basis, basis_s_size);
else
{
trans(z, p, tt);
lift(z, s);
}
}
else trans(z, p, tt);
return reduce_inner(*v, s, k);
}
void PowerCrustSurfaceReconstructionImpl::get_basis_sede(simplex *s)
{
int k = 1;
neighbor *sn = s->neigh + 1,
*sn0 = s->neigh;
if ((vd || power_diagram) && sn0->vert == infinity && cdim > 1)
{
SWAP(neighbor, *sn0, *sn);
NULLIFY(basis_s, sn0->basis);
sn0->basis = tt_basisp;
tt_basisp->ref_count++;
}
else
{
if (!sn0->basis)
{
sn0->basis = tt_basisp;
tt_basisp->ref_count++;
}
else while (k < cdim && sn->basis)
{
k++;
sn++;
}
}
while (k < cdim)
{
NULLIFY(basis_s, sn->basis);
reduce(&sn->basis, sn->vert, s, k);
k++;
sn++;
}
}
int PowerCrustSurfaceReconstructionImpl::out_of_flat(simplex *root, point p)
{
if (!p_neigh.basis)
p_neigh.basis = (basis_s*)malloc(basis_s_size);
p_neigh.vert = p;
cdim++;
root->neigh[cdim - 1].vert = root->peak.vert;
NULLIFY(basis_s, root->neigh[cdim - 1].basis);
get_basis_sede(root);
if ((vd || power_diagram) && root->neigh[0].vert == infinity) return 1;
reduce(&p_neigh.basis, p, root, cdim);
if (p_neigh.basis->sqa != 0) return 1;
cdim--;
return 0;
}
double PowerCrustSurfaceReconstructionImpl::cosangle_sq(basis_s* v, basis_s* w)
{
double dd;
point vv = v->vecs, wv = w->vecs;
dd = Vec_dot(vv, wv);
return dd*dd / Norm2(vv) / Norm2(wv);
}
int PowerCrustSurfaceReconstructionImpl::check_perps(simplex *s)
{
point z, y;
point tt;
for (int i = 1; i < cdim; i++) if (NEARZERO(s->neigh[i].basis->sqb)) return 0;
if (!check_perps_b)
{
check_perps_b = (basis_s*)malloc(basis_s_size);
assert(check_perps_b);
}
else check_perps_b->lscale = 0;
z = VB(check_perps_b);
tt = s->neigh[0].vert;
for (int i = 1; i < cdim; i++)
{
y = s->neigh[i].vert;
if ((vd || power_diagram) && y == infinity) memcpy(check_perps_b, infinity_basis, basis_s_size);
else
{
trans(z, y, tt);
lift(z, s);
}
if (s->normal && cosangle_sq(check_perps_b, s->normal) > b_err_min_sq)
{
return 0;
}
for (int j = i + 1; j < cdim; j++)
{
if (cosangle_sq(check_perps_b, s->neigh[j].basis) > b_err_min_sq)
{
return 0;
}
}
}
return 1;
}
void PowerCrustSurfaceReconstructionImpl::get_normal_sede(simplex *s)
{
neighbor *rn;
int i, j;
get_basis_sede(s);
if (rdim == 3 && cdim == 3)
{
point c;
point a = VB(s->neigh[1].basis);
point b = VB(s->neigh[2].basis);
NEWLRC(basis_s, s->normal);
c = VB(s->normal);
c[0] = a[1] * b[2] - a[2] * b[1];
c[1] = a[2] * b[0] - a[0] * b[2];
c[2] = a[0] * b[1] - a[1] * b[0];
s->normal->sqb = Norm2(c);
for (i = cdim + 1, rn = ch_root->neigh + cdim - 1; i; i--, rn--)
{
for (j = 0; j < cdim && rn->vert != s->neigh[j].vert; j++);
if (j < cdim) continue;
if (rn->vert == infinity)
{
if (c[2] > -b_err_min) continue;
}
else if (!sees(rn->vert, s)) continue;
c[0] = -c[0];
c[1] = -c[1];
c[2] = -c[2];
break;
}
return;
}
for (i = cdim + 1, rn = ch_root->neigh + cdim - 1; i; i--, rn--)
{
for (j = 0; j < cdim && rn->vert != s->neigh[j].vert; j++);
if (j < cdim) continue;
reduce(&s->normal, rn->vert, s, cdim);
if (s->normal->sqb != 0) break;
}
}
void PowerCrustSurfaceReconstructionImpl::get_normal(simplex *s)
{
get_normal_sede(s);
return;
}
int PowerCrustSurfaceReconstructionImpl::sees(site p, simplex *s)
{
point tt, zz;
double dd, dds;
int i;
if (!seesB)
{
seesB = (basis_s*)malloc(basis_s_size);
assert(seesB);
}
else
seesB->lscale = 0;
zz = VB(seesB);
if (cdim == 0)
return 0;
if (!s->normal)
{
get_normal_sede(s);
for (i = 0; i < cdim; i++)
NULLIFY(basis_s, s->neigh[i].basis);
}
tt = s->neigh[0].vert;
if (vd || power_diagram)
{
if (p == infinity)
memcpy(seesB, infinity_basis, basis_s_size);
else
{
trans(zz, p, tt);
lift(zz, s);
}
}
else
trans(zz, p, tt);
for (i = 0; i < 3; i++)
{
dd = Vec_dot(zz, s->normal->vecs);
if (dd == 0.0)
{
return 0;
}
dds = dd*dd / s->normal->sqb / Norm2(zz);
if (dds > b_err_min_sq)
return (dd < 0);
get_basis_sede(s);
reduce_inner(seesB, s, cdim);
}
return 0;
}
double PowerCrustSurfaceReconstructionImpl::radsq(simplex *s)
{
point n;
neighbor *sn;
int i;
/* square of ratio of circumcircle radius to max edge length for Delaunay tetrahedra */
for (i = 0, sn = s->neigh; i < cdim; i++, sn++)
if (sn->vert == infinity) return Huge;
if (!s->normal) get_normal_sede(s);
/* compute circumradius */
n = s->normal->vecs;
if (NEARZERO(n[rdim - 1]))
{
return Huge;
}
return Vec_dot_pdim(n, n) / 4 / n[rdim - 1] / n[rdim - 1];
}
void * PowerCrustSurfaceReconstructionImpl::zero_marks(simplex * s, void *)
{
s->mark = 0;
return NULL;
}
void * PowerCrustSurfaceReconstructionImpl::one_marks(simplex * s, void *)
{
s->mark = 1;
return NULL;
}
int PowerCrustSurfaceReconstructionImpl::alph_test(simplex *s, int i, void *alphap)
{
/*returns 1 if not an alpha-facet */
simplex *si;
neighbor *scn, *sin;
int k;
if (alphap)
{
alpha_test_alpha = *(double*)alphap;
if (!s) return 1;
}
if (i == -1) return 0;
si = s->neigh[i].simp;
scn = s->neigh + cdim - 1;
sin = s->neigh + i;
int nsees = 0;
for (k = 0; k < cdim; k++) if (s->neigh[k].vert == infinity && k != i) return 1;
double rs = radsq(s);
double rsi = radsq(si);
if (rs < alpha_test_alpha && rsi < alpha_test_alpha) return 1;
swap_points(scn->vert, sin->vert);
NULLIFY(basis_s, s->neigh[i].basis);
cdim--;
get_basis_sede(s);
reduce(&s->normal, infinity, s, cdim);
double rsfi = radsq(s);
for (k = 0; k < cdim; k++) if (si->neigh[k].simp == s) break;
int ssees = sees(scn->vert, s);
if (!ssees) nsees = sees(si->neigh[k].vert, s);
swap_points(scn->vert, sin->vert);
cdim++;
NULLIFY(basis_s, s->normal);
NULLIFY(basis_s, s->neigh[i].basis);
if (ssees) return alpha_test_alpha < rs;
if (nsees) return alpha_test_alpha < rsi;
assert(rsfi <= rs + FLT_EPSILON && rsfi <= rsi + FLT_EPSILON);
return alpha_test_alpha <= rsfi;
}
void * PowerCrustSurfaceReconstructionImpl::conv_facetv(simplex *s, void *)
{
for (int i = 0; i < cdim; i++)
if (s->neigh[i].vert == infinity)
{
return s;
}
return NULL;
}
void * PowerCrustSurfaceReconstructionImpl::mark_points(simplex *s, void *)
{
int i, snum;
neighbor *sn;
for (i = 0, sn = s->neigh; i < cdim; i++, sn++)
{
if (sn->vert == infinity) continue;
snum = (this->*site_num) (sn->vert);
if (s->mark) mo[snum] = 1;
else mi[snum] = 1;
}
return NULL;
}
void * PowerCrustSurfaceReconstructionImpl::visit_outside_ashape(simplex *root, visit_func visit)
{
return visit_triang_gen((simplex*)visit_hull(root, &PowerCrustSurfaceReconstructionImpl::conv_facetv), visit, &PowerCrustSurfaceReconstructionImpl::alph_test);
}
int PowerCrustSurfaceReconstructionImpl::check_ashape(simplex *root, double alpha)
{
for (int i = 0; i < MAXPOINTS; i++)
{
mi[i] = mo[i] = 0;
}
visit_hull(root, &PowerCrustSurfaceReconstructionImpl::zero_marks);
alph_test(0, 0, &alpha);
visit_outside_ashape(root, &PowerCrustSurfaceReconstructionImpl::one_marks);
visit_hull(root, &PowerCrustSurfaceReconstructionImpl::mark_points);
for (int i = 0; i < MAXPOINTS; i++) if (mo[i] && !mi[i])
{
return 0;
}
return 1;
}
simplex * PowerCrustSurfaceReconstructionImpl::build_convex_hull(short dim, short vdd)
{
/*
get_s returns next site each call;
hull construction stops when NULL returned;
site_numm returns number of site when given site;
dim dimension of point set;
vdd if (vdd) then return Delaunay triangulation
*/
simplex *s, *root;
cdim = 0;
get_site = &PowerCrustSurfaceReconstructionImpl::get_next_site;
site_num = &PowerCrustSurfaceReconstructionImpl::site_numm;
pdim = dim;
vd = vdd;
exact_bits = (int)(DBL_MANT_DIG*log(FLT_RADIX) / log(2.0)); // cast to int added by TJH // EPRO added
b_err_min = DBL_EPSILON*MAXDIM* (1 << MAXDIM) * MAXDIM * 3.01;
b_err_min_sq = b_err_min * b_err_min;
assert(get_site != NULL);
assert(site_num != NULL);
rdim = vd ? pdim + 1 : pdim;
if (rdim > MAXDIM)
ASSERT(pcFALSE, "dimension bound MAXDIM exceeded");
site_size = sizeof(Coord) *pdim;
basis_vec_size = sizeof(Coord) *rdim;
basis_s_size = sizeof(basis_s) + (2 * rdim - 1) *sizeof(Coord);
simplex_size = sizeof(simplex) + (rdim - 1) *sizeof(neighbor);
Tree_size = sizeof(Tree);
fg_size = sizeof(fg);
root = NULL;
if (vd || power_diagram)
{
p = infinity;
NEWLRC(basis_s, infinity_basis);
infinity_basis->vecs[2 * rdim - 1] = infinity_basis->vecs[rdim - 1] = 1;
infinity_basis->sqa = infinity_basis->sqb = 1;
}
else if (!(p = (this->*get_site) ())) return 0;
NEWL(simplex, root);
ch_root = root;
copy_simp(s, root);
root->peak.vert = p;
root->peak.simp = s;
s->peak.simp = root;
buildhull(root);
return root;
}
void PowerCrustSurfaceReconstructionImpl::free_hull_storage(void)
{
free_basis_s_storage();
free_simplex_storage();
free_Tree_storage();
free_fg_storage();
}
void * PowerCrustSurfaceReconstructionImpl::compute_vv(simplex *s, void *)
{
/* computes Voronoi vertices */
point v[MAXDIM];
int inf = 0;
double cc[3], cond, ta[4][3];
if (!s) return NULL;
for (int j = 0; j < cdim; j++)
{
v[j] = s->neigh[j].vert;
/* v[j] stores coordinates of j'th vertex of simplex s; j=0..3 */
if (v[j] == infinity) /* means simplex s is on the convex hull */
{
inf = 1;
break; /* skip the rest of the for loop, ignore convex hull faces (= bounding box ) */
}
for (int k = 0; k < cdim - 1; k++)
{
ta[j][k] = v[j][k] / mult_up; /* restore original coords */
}
}
if (!inf) /* if not faces on convex hull, compute circumcenter*/
{
tetcircumcenter(ta[0], ta[1], ta[2], ta[3], cc, &cond);
if (cond != 0) /* ignore them if cond = 0 */
{
s->isVvNull = false;
for (int k = 0; k < cdim - 1; k++)
{
s->vv[k] = ta[0][k] + cc[k];
}
s->status = VV;
}
else
{
s->isVvNull = true;
s->status = SLV;
}
}
else /* if on conv hull */
{
s->status = CNV;
}
/* computing poles */
for (int j = 0; j < cdim; j++) /* compute 1st pole for vertex j */
{
int i = (this->*site_num) (s->neigh[j].vert);
if (i == -1) continue;
/* Ignore poles that are too far away to matter - a relic of the
original California-style crust. Probably no longer needed */
if ((s->neigh[j].vert[0] > omaxs[0]) ||
(s->neigh[j].vert[0] < omins[0]) ||
(s->neigh[j].vert[1] > omaxs[1]) ||
(s->neigh[j].vert[1] < omins[1]) ||
(s->neigh[j].vert[2] > omaxs[2]) ||
(s->neigh[j].vert[2] < omins[2]))
{
pole1[i] = NULL;
continue;
}
else
{
if (pole1[i] == NULL)
{
/* the vertex i is encountered for the 1st time */
if (s->status == VV) /* we don't store infinite poles */
{
pole1[i] = s;
continue;
}
}
if ((s->status == VV) && (pole1[i]->status == VV))
{
if (sqdist(pole1[i]->vv, ta[j]) < sqdist(s->vv, ta[j]))
{
pole1[i] = s; /* update 1st pole */
}
}
}
}
return NULL;
}
void * PowerCrustSurfaceReconstructionImpl::compute_pole2(simplex *s, void *p)
{
point v[MAXDIM];
int inf = 0;
double a[3] = { 0, 0, 0 };
site t;
double dir_s[3], dir_p[3], dist_s, dist_p;
if (p)
{
if (!s) return NULL;
}
for (int j = 0; j < cdim; j++)
{
v[j] = s->neigh[j].vert;
int i = (this->*site_num) (v[j]);
if (i == -1) inf = 1;
}
double cos_2r = cos(2 * est_r);
for (int j = 0; j < cdim; j++) /* compute 2nd poles */
{
t = s->neigh[j].vert;
int i = (this->*site_num) (t);
if (i < 0) continue; /* not a vertex */
if (inf) /* on conv hull */
{
if (s->status == CNV)
{
continue;
}
}
if (!pole1[i])
{
continue;
}
if (pole1[i]->isVvNull) //pole1[i]->vv==NULL
{
continue;
}
if (s->isVvNull) //!s->vv
{
if (s->status != SLV)
continue;
}
for (int k = 0; k < cdim - 1; k++) /* a stores orig vertex coord */
{
a[k] = t[k] / mult_up;
}
/* compute direction and length of vector from sample to first pole */
dir_and_dist(a, pole1[i]->vv, dir_p, &dist_p);
/* We have a vertex, and there is a good first pole. */
if ((s->status == VV) && (pole1[i]->status == VV))
{
/* make direction vector from sample to this Voronoi vertex */
dir_and_dist(a, s->vv, dir_s, &dist_s);
/* cosine of angle between angle to vertex and angle to pole */
double cos_sp = dir_s[0] * dir_p[0] + dir_s[1] * dir_p[1] + dir_s[2] * dir_p[2];
/* if there is an estimate for r, use it to estimate lfs */
if (est_r < 1.0)
{
/* near vertices - should be close to sample (a) */
if ((cos_sp < cos_2r) && (cos_sp > -cos_2r))
{
/* use to get lower bound on lfs */
double est_lfs = dist_s / est_r * ((sqrt(1 - cos_sp*cos_sp)) - est_r);
if (est_lfs > lfs_lb[i]) lfs_lb[i] = est_lfs;
}
}
else
{
lfs_lb[i] = 0;
}
if (cos_sp > 0)
{
/* s->vv is in the same side of pole1 */
continue;
}
/* s->vv is a candidate for pole2 */
if (!pole2[i])
{
/* 1st pole2 candidate for vertex i */
pole2[i] = s;
continue;
}
else if (!pole2[i]->vv) /* 2nd pole points null */
{
continue;
}
else if ((pole2[i]->status == VV) && (sqdist(a, pole2[i]->vv) < sqdist(a, s->vv)))
pole2[i] = s; /* update 2nd pole */
}
}
return NULL;
}
int PowerCrustSurfaceReconstructionImpl::close_pole(double* v, double* p, double lfs_lb)
{
return (sqdist(v, p) < lfs_lb * lfs_lb);
}
int PowerCrustSurfaceReconstructionImpl::antiLabel(int label)
{
if (label == IN) return (OUT);
if (label == OUT) return (IN);
return (label);
}
double PowerCrustSurfaceReconstructionImpl::computePoleAngle(simplex* pole1, simplex* pole2, double* samp)
{
return (((pole1->vv[0] - samp[0]) * (pole2->vv[0] - samp[0]) +
(pole1->vv[1] - samp[1]) * (pole2->vv[1] - samp[1]) +
(pole1->vv[2] - samp[2]) * (pole2->vv[2] - samp[2])) /
(sqrt(SQ(pole1->vv[0] - samp[0]) + SQ(pole1->vv[1] - samp[1]) + SQ(pole1->vv[2] - samp[2])) *
sqrt(SQ(pole2->vv[0] - samp[0]) + SQ(pole2->vv[1] - samp[1]) + SQ(pole2->vv[2] - samp[2]))));
}
void PowerCrustSurfaceReconstructionImpl::newOpposite(int p1index, int p2index, double pole_angle)
{
plist* newplist;
newplist = (plist *)malloc(sizeof(plist));
newplist->pid = p2index;
newplist->angle = pole_angle;
newplist->next = opplist[p1index];
opplist[p1index] = newplist;
if (adjlist[p1index].oppradius > adjlist[p2index].sqradius)
{
assert(adjlist[p2index].sqradius > 0.0);
adjlist[p1index].oppradius = adjlist[p2index].sqradius;
}
}
void PowerCrustSurfaceReconstructionImpl::outputPole(simplex* pole, int poleid, double* samp, int* num_poles, double distance)
{
double r2 = SQ(pole->vv[0] - samp[0]) + SQ(pole->vv[1] - samp[1]) + SQ(pole->vv[2] - samp[2]);
double weight = SQ(pole->vv[0]) + SQ(pole->vv[1]) + SQ(pole->vv[2]) - r2;
pole->status = POLE_OUTPUT;
pole->poleindex = poleid;
medial_surface.InsertNextPoint(pole->vv[0], pole->vv[1], pole->vv[2]);
medial_surface.InsertWeight(weight);
/* remember squared radius */
adjlist[poleid].sqradius = r2;
adjlist[poleid].samp_distance = distance;
/* initialize perp dist to MA */
adjlist[poleid].oppradius = r2;
/* initialize */
adjlist[poleid].grafindex = -1;
/* keep count! */
(*num_poles)++;
}
Tree * PowerCrustSurfaceReconstructionImpl::splay(site i, Tree *t)
{
Tree N, *l, *r, *y;
if (!t) return t;
N.left = N.right = NULL;
l = r = &N;
int l_size = 0;
int r_size = 0;
for (;;)
{
int comp = compare(i, t->key);
if (comp < 0)
{
if (!t->left) break;
if (compare(i, t->left->key) < 0)
{
y = t->left; /* rotate right */
t->left = y->right;
y->right = t;
t->size = node_size(t->left) + node_size(t->right) + 1;
t = y;
if (!t->left) break;
}
r->left = t; /* link right */
r = t;
t = t->left;
r_size += 1 + node_size(r->right);
}
else if (comp > 0)
{
if (!t->right) break;
if (compare(i, t->right->key) > 0)
{
y = t->right; /* rotate left */
t->right = y->left;
y->left = t;
t->size = node_size(t->left) + node_size(t->right) + 1;
t = y;
if (!t->right) break;
}
l->right = t; /* link left */
l = t;
t = t->right;
l_size += 1 + node_size(l->left);
}
else break;
}
l_size += node_size(t->left); /* Now l_size and r_size are the sizes of */
r_size += node_size(t->right); /* the left and right trees we just built.*/
t->size = l_size + r_size + 1;
l->right = r->left = NULL;
/* The following two loops correct the size fields of the right path */
/* from the left child of the root and the right path from the left */
/* child of the root. */
for (y = N.right; y != NULL; y = y->right)
{
y->size = l_size;
l_size -= 1 + node_size(y->left);
}
for (y = N.left; y != NULL; y = y->left)
{
y->size = r_size;
r_size -= 1 + node_size(y->right);
}
l->right = t->left; /* assemble */
r->left = t->right;
t->left = N.right;
t->right = N.left;
return t;
}
Tree * PowerCrustSurfaceReconstructionImpl::insert(site i, Tree * t)
{
/* Insert key i into the tree t, if it is not already there. */
/* Return a pointer to the resulting tree. */
Tree *new_tree;
if (t != NULL)
{
t = splay(i, t);
if (compare(i, t->key) == 0)
{
return t; /* it's already there */
}
}
NEWL(Tree, new_tree)
if (!t)
{
new_tree->left = new_tree->right = NULL;
}
else if (compare(i, t->key) < 0)
{
new_tree->left = t->left;
new_tree->right = t;
t->left = NULL;
t->size = 1 + node_size(t->right);
}
else
{
new_tree->right = t->right;
new_tree->left = t;
t->right = NULL;
t->size = 1 + node_size(t->left);
}
new_tree->key = i;
new_tree->size = 1 + node_size(new_tree->left) + node_size(new_tree->right);
return new_tree;
}
void PowerCrustSurfaceReconstructionImpl::free_heap()
{
free(heap_A);
}
void PowerCrustSurfaceReconstructionImpl::init_heap(int num)
{
heap_A = (heap_array *)calloc(num + 1, sizeof(heap_array)); // EPRO added
heap_size = 0;
heap_length = num;
}
void PowerCrustSurfaceReconstructionImpl::heapify(int hi)
{
int largest;
if ((LEFT(hi) <= heap_size) && (heap_A[LEFT(hi)].pri > heap_A[hi].pri))
largest = LEFT(hi);
else largest = hi;
if ((RIGHT(hi) <= heap_size) && (heap_A[RIGHT(hi)].pri > heap_A[largest].pri))
largest = RIGHT(hi);
if (largest != hi)
{
int temp = heap_A[hi].pid;
heap_A[hi].pid = heap_A[largest].pid;
adjlist[heap_A[hi].pid].hid = hi;
heap_A[largest].pid = temp;
adjlist[heap_A[largest].pid].hid = largest;
double td = heap_A[hi].pri;
heap_A[hi].pri = heap_A[largest].pri;
heap_A[largest].pri = td;
heapify(largest);
}
}
int PowerCrustSurfaceReconstructionImpl::extract_max()
{
if (heap_size < 1) return -1;
int max = heap_A[1].pid;
heap_A[1].pid = heap_A[heap_size].pid;
heap_A[1].pri = heap_A[heap_size].pri;
adjlist[heap_A[1].pid].hid = 1;
heap_size--;
heapify(1);
return max;
}
int PowerCrustSurfaceReconstructionImpl::insert_heap(int pi, double pr)
{
heap_size++;
int i = heap_size;
while ((i > 1) && (heap_A[PARENT(i)].pri < pr))
{
heap_A[i].pid = heap_A[PARENT(i)].pid;
heap_A[i].pri = heap_A[PARENT(i)].pri;
adjlist[heap_A[i].pid].hid = i;
i = PARENT(i);
};
heap_A[i].pri = pr;
heap_A[i].pid = pi;
adjlist[pi].hid = i;
return i;
}
void PowerCrustSurfaceReconstructionImpl::update(int hi, double pr)
{
heap_A[hi].pri = pr;
int pi = heap_A[hi].pid;
if (pr > heap_A[PARENT(hi)].pri)
{
int i = hi;
while ((i > 1) && (heap_A[PARENT(i)].pri < pr))
{
heap_A[i].pid = heap_A[PARENT(i)].pid;
heap_A[i].pri = heap_A[PARENT(i)].pri;
adjlist[heap_A[i].pid].hid = i;
i = PARENT(i);
};
heap_A[i].pri = pr;
heap_A[i].pid = pi;
adjlist[pi].hid = i;
}
else heapify(hi);
}
void * PowerCrustSurfaceReconstructionImpl::visit_triang_gen(simplex *s, visit_func visit, test_func test)
{
/* starting at s, visit simplices t such that test(s,i,0) is true,
* and t is the i'th neighbor of s;
* apply visit function to all visited simplices;
* when visit returns nonNULL, exit and return its value */
neighbor *sn;
void *v;
simplex *t;
int i;
long tms = 0;
visit_triang_gen_vnum--;
if (!st_visit_triang_gen)
{
st_visit_triang_gen = (simplex**)malloc((visit_triang_gen_ss + MAXDIM + 1) * sizeof(simplex*));
assert(st_visit_triang_gen);
}
if (s) push(s, st_visit_triang_gen, tms);
while (tms)
{
if (tms > visit_triang_gen_ss)
{
st_visit_triang_gen = (simplex**)realloc(st_visit_triang_gen, ((visit_triang_gen_ss += visit_triang_gen_ss) + MAXDIM + 1) * sizeof(simplex*));
assert(st_visit_triang_gen);
}
pop(t, st_visit_triang_gen, tms);
if (!t || t->visit == visit_triang_gen_vnum) continue;
t->visit = visit_triang_gen_vnum;
if ((v = (this->*visit) (t, 0)) != NULL)
{
return v;
}
for (i = -1, sn = t->neigh - 1; i < cdim; i++, sn++)
if ((sn->simp->visit != visit_triang_gen_vnum) && sn->simp && (this->*test) (t, i, 0))
push(sn->simp, st_visit_triang_gen, tms);
}
return NULL;
}
int PowerCrustSurfaceReconstructionImpl::truet(simplex *, int , void *)
{
return 1;
}
void * PowerCrustSurfaceReconstructionImpl::visit_triang(simplex *root, visit_func visit)
{
return visit_triang_gen(root, visit, &PowerCrustSurfaceReconstructionImpl::truet);
}
int PowerCrustSurfaceReconstructionImpl::hullt(simplex *, int i, void *)
{
return i > -1;
}
void * PowerCrustSurfaceReconstructionImpl::facet_test(simplex *s, void *)
{
return (!s->peak.vert) ? s : NULL;
}
void * PowerCrustSurfaceReconstructionImpl::visit_hull(simplex *root, visit_func visit)
{
return visit_triang_gen((simplex*)visit_triang(root, &PowerCrustSurfaceReconstructionImpl::facet_test), visit, &PowerCrustSurfaceReconstructionImpl::hullt);
}
Coord PowerCrustSurfaceReconstructionImpl::Vec_dot(point x, point y)
{
Coord sum = 0;
for (int i = 0; i < rdim; i++) sum += x[i] * y[i];
return sum;
}
Coord PowerCrustSurfaceReconstructionImpl::Vec_dot_pdim(point x, point y)
{
Coord sum = 0;
for (int i = 0; i < pdim; i++) sum += x[i] * y[i];
return sum;
}
Coord PowerCrustSurfaceReconstructionImpl::Norm2(point x)
{
Coord sum = 0;
for (int i = 0; i < rdim; i++) sum += x[i] * x[i];
return sum;
}
void PowerCrustSurfaceReconstructionImpl::Ax_plus_y(Coord a, point x, point y)
{
for (int i = 0; i < rdim; i++)
{
*y++ += a * *x++;
}
}
void PowerCrustSurfaceReconstructionImpl::Ax_plus_y_test(Coord a, point x, point y)
{
for (int i = 0; i < rdim; i++)
{
*y++ += a * *x++;
}
}
void PowerCrustSurfaceReconstructionImpl::Vec_scale_test(int n, Coord a, Coord *x)
{
register Coord *xx = x, *xend = xx + n;
while (xx != xend)
{
*xx *= a;
xx++;
}
}
double PowerCrustSurfaceReconstructionImpl::sc(basis_s *v, simplex *s, int k, int j)
{
/* amount by which to scale up vector, for reduce_inner */
if (j < 10)
{
double labound = logb(v->sqa) / 2;
sc_max_scale = exact_bits - labound - 0.66* (k - 2) - 1 - DELIFT;
if (sc_max_scale < 1)
{
sc_max_scale = 1;
}
if (j == 0)
{
int i;
neighbor *sni;
basis_s *snib;
sc_ldetbound = DELIFT;
sc_Sb = 0;
for (i = k - 1, sni = s->neigh + k - 1; i > 0; i--, sni--)
{
snib = sni->basis;
sc_Sb += snib->sqb;
sc_ldetbound += logb(snib->sqb) / 2 + 1;
sc_ldetbound -= snib->lscale;
}
}
}
if (sc_ldetbound - v->lscale + logb(v->sqb) / 2 + 1 < 0)
{
return 0;
}
else
{
sc_lscale = (int)(logb(2 * sc_Sb / (v->sqb + v->sqa*b_err_min)) / 2); // cast to int added by TJH
if (sc_lscale > sc_max_scale)
{
sc_lscale = (int)sc_max_scale; // cast added by TJH (is lscale really meant to be int, not double?)
}
else if (sc_lscale < 0) sc_lscale = 0;
v->lscale += sc_lscale;
return two_to(sc_lscale);
}
}
int PowerCrustSurfaceReconstructionImpl::reduce_inner(basis_s *v, simplex *s, int k)
{
point va = VA(v), vb = VB(v);
int i;
basis_s *snibv;
neighbor *sni;
v->sqa = v->sqb = Norm2(vb);
if (k <= 1)
{
memcpy(vb, va, basis_vec_size);
return 1;
}
for (int j = 0; j < 250; j++)
{
memcpy(vb, va, basis_vec_size);
for (i = k - 1, sni = s->neigh + k - 1; i > 0; i--, sni--)
{
snibv = sni->basis;
double dd = -Vec_dot(VB(snibv), vb) / snibv->sqb;
Ax_plus_y(dd, VA(snibv), vb);
}
v->sqb = Norm2(vb);
v->sqa = Norm2(va);
if (2 * v->sqb >= v->sqa)
{
return 1;
}
Vec_scale_test(rdim, sc(v, s, k, j), va);
for (i = k - 1, sni = s->neigh + k - 1; i > 0; i--, sni--)
{
snibv = sni->basis;
double dd = -Vec_dot(VB(snibv), va) / snibv->sqb;
dd = floor(dd + 0.5);
Ax_plus_y_test(dd, VA(snibv), va);
}
}
return 0;
}
void PowerCrustSurfaceReconstructionImpl::trans(point z, point p, point q)
{
for (int i = 0; i < pdim; i++)
z[i + rdim] = z[i] = p[i] - q[i];
}
void * PowerCrustSurfaceReconstructionImpl::compute_3d_power_vv(simplex *s, void *p)
{
point v[MAXDIM];
int inf = 0;
int index = 0;
double cc[3], cond, ta[4][4];
edgesimp *pindex;
if (p)
{
if (!s) return NULL;
}
for (int j = 0; j < cdim; j++)
{
v[j] = s->neigh[j].vert;
/* v[j] stores coordinates of j'th vertex of simplex s; j=0..3 */
if (v[j] == infinity) /* means simplex s is on the convex hull */
{
inf = 1;
continue; /* skip the rest of the for loop; process next vertex */
}
for (int k = 0; k < 4; k++)
{
ta[index][k] = v[j][k] / mult_up; /* restore original coords */
}
index++;
}
/* if not faces on convex hull, process */
if (!inf)
{
/* build structure for each edge, including angle of intersection */
for (int k = 0; k < 6; k++)
{
if (s->edgestatus[k] == FIRST_EDGE) /* not visited edge */
{
pindex = adjlist[(this->*site_num) (v[v1[k]])].eptr;
bool visited_edge = false;
while (pindex != NULL)
{
if (pindex->pid == (this->*site_num) (v[v2[k]])) /* already in the list */
{
visited_edge = true;
break;
}
pindex = pindex->next;
}
if (!visited_edge)
{
double d = sqdist(ta[v1[k]], ta[v2[k]]);
double r1 = SQ(ta[v1[k]][0]) + SQ(ta[v1[k]][1]) + SQ(ta[v1[k]][2]) - ta[v1[k]][3];
double r2 = SQ(ta[v2[k]][0]) + SQ(ta[v2[k]][1]) + SQ(ta[v2[k]][2]) - ta[v2[k]][3];
double e = 2 * sqrt(r1) * sqrt(r2);
edgesimp *newplist1;
newplist1 = (edgesimp *)malloc(sizeof(edgesimp));
newplist1->simp = s;
newplist1->kth = k;
newplist1->angle = (r1 + r2 - d) / e;
newplist1->pid = (this->*site_num) (v[v1[k]]);
newplist1->next = adjlist[(this->*site_num) (v[v2[k]])].eptr;
adjlist[(this->*site_num) (v[v2[k]])].eptr = newplist1;
edgesimp *newplist2;
newplist2 = (edgesimp *)malloc(sizeof(edgesimp));
newplist2->simp = s;
newplist2->kth = k;
newplist2->angle = (r1 + r2 - d) / e;
newplist2->pid = (this->*site_num) (v[v2[k]]);
newplist2->next = adjlist[(this->*site_num) (v[v1[k]])].eptr;
adjlist[(this->*site_num) (v[v1[k]])].eptr = newplist2;
s->edgestatus[k] = VISITED;
}
}
}
tetorthocenter(ta[0], ta[1], ta[2], ta[3], cc, &cond);
/* cc is the displacement of orthocenter from ta[0] */
/* cond is the denominator ( orient2d ) value */
if (cond != 0) /* ignore them if cond = 0 */
{
s->isVvNull = false;
for (int k = 0; k < 3; k++)
{
s->vv[k] = ta[0][k] + cc[k];
}
s->status = VV;
}
else /* if cond=0, s is SLIVER */
{
s->isVvNull = true;
s->status = SLV;
}
}
else /* if on conv hull, ignore */
{
s->isVvNull = true;
s->status = CNV;
}
return NULL;
}
void * PowerCrustSurfaceReconstructionImpl::compute_axis(simplex *s, void *p)
{
point v[MAXDIM];
point point1, point2;
int edgedata[6];
int indices[6];
if (p)
{
if (!s) return NULL;
}
if ((s->status == CNV) || (s->status == SLV)) return NULL; /* skip inf faces */
for (int j = 0; j < cdim; j++)
{
v[j] = s->neigh[j].vert;
}
for (int k = 0; k < 6; k++) /* for each edge */
{
edgedata[k] = 0;
if ((s->edgestatus[k] != POW)) /* not dual to a power face */
{
point1 = v[v1[k]];
point2 = v[v2[k]];
int pindex = (this->*site_num) (point1);
int qindex = (this->*site_num) (point2);
if (adjlist[pindex].label == IN && adjlist[qindex].label == IN)
{
if (s->edgestatus[k] != ADDAXIS)
{
}
edgedata[k] = VALIDEDGE;
indices[v1[k]] = pindex;
indices[v2[k]] = qindex;
s->edgestatus[k] = ADDAXIS;
}
/* now start adding triangles if present */
}
}
if ((edgedata[0] == VALIDEDGE) && (edgedata[1] == VALIDEDGE) && (edgedata[3] == VALIDEDGE))
{
{
medial_surface.InsertTriangle(indices[v1[0]], indices[v2[1]], indices[v1[3]]);
}
}
if ((edgedata[1] == VALIDEDGE) && (edgedata[2] == VALIDEDGE) && (edgedata[5] == VALIDEDGE))
{
{
medial_surface.InsertTriangle(indices[v1[1]], indices[v2[2]], indices[v1[5]]);
}
}
if ((edgedata[0] == VALIDEDGE) && (edgedata[2] == VALIDEDGE) && (edgedata[4] == VALIDEDGE))
{
{
medial_surface.InsertTriangle(indices[v1[0]], indices[v2[2]], indices[v1[4]]);
}
}
if ((edgedata[3] == VALIDEDGE) && (edgedata[4] == VALIDEDGE) && (edgedata[5] == VALIDEDGE))
{
{
medial_surface.InsertTriangle(indices[v1[3]], indices[v2[4]], indices[v1[5]]);
}
}
return NULL;
}
void PowerCrustSurfaceReconstructionImpl::construct_face(simplex *s, short k)
{
site edge0, edge1, nextv, remv, prevv, outsite, insite;
simplex *prevs, *nexts;
int j, numedges, l1, l2, nk, l, nedge0 = 0, nedge1 = 0, nremv = 0, nnextv = 0, i;
char indface[1024][32]; /* the indices of the face */
double plane[3][3];
double outpole[3], inpole[3];
edge0 = s->neigh[v1[k]].vert;
edge1 = s->neigh[v2[k]].vert;
if (adjlist[(this->*site_num) (edge0)].label == OUT)
{
outsite = edge0;
insite = edge1;
}
else
{
outsite = edge1;
insite = edge0;
}
for (j = 0; j < 3; j++)
{
outpole[j] = outsite[j] / mult_up;
inpole[j] = insite[j] / mult_up;
}
nextv = s->neigh[v3[k]].vert;
/* nextv is the opposite vtx of the next simplex */
remv = s->neigh[v4[k]].vert;
/* remv is a vtx of the next simplex with edge0, edge1 */
prevv = remv;
/* prevv is the vtx shared by prevs and nexts besides edge0, edge1 */
/* construct its dual power face */
s->edgestatus[k] = POW;
/* visit the next simplex */
prevs = s;
nexts = s->neigh[v3[k]].simp;
numedges = 0;
while (nexts != s)
{
if (nexts->status == CNV)
{
break;
}
else
{
if (prevs->status != POLE_OUTPUT)
{
/* this vertex is not yet output */
prevs->status = POLE_OUTPUT;
prevs->poleindex = num_vtxs++;
// TJH: PC contains the points for the powercrust surface
// so we hijack the data and pipe it to our structure
{
float vp[3];
vp[0] = prevs->vv[0];
vp[1] = prevs->vv[1];
vp[2] = prevs->vv[2];
output.InsertNextPoint(vp);
}
}
if (numedges < 3)
{
plane[numedges][0] = prevs->vv[0];
plane[numedges][1] = prevs->vv[1];
plane[numedges][2] = prevs->vv[2];
}
sprintf(indface[numedges], "%ld ", prevs->poleindex);
numedges++;
/* find edgenumber k of nexts for this edge */
for (l = 0; l < 4; l++)
{
if (nexts->neigh[l].vert == edge0)
{
nedge0 = l;
continue;
}
else if (nexts->neigh[l].vert == edge1)
{
nedge1 = l;
continue;
}
else if (nexts->neigh[l].vert == prevv)
{
nremv = l;
continue;
}
else if (nexts->neigh[l].vert == nextv)
{
nnextv = l;
continue;
}
else
{
nnextv = l;
}
}
if (nedge0 > nedge1)
{
l1 = nedge1;
l2 = nedge0;
}
else
{
l2 = nedge1;
l1 = nedge0;
}
if (l1 == 0)
{
if (l2 == 1) nk = 0;
else if (l2 == 2) nk = 1;
else nk = 2;
}
else if (l1 == 1)
{
if (l2 == 2) nk = 3;
else nk = 4;
}
else nk = 5;
/* found nk for the edge */
nexts->edgestatus[nk] = POW; /* record that it's visited */
/* visit next simplex (opposite vertex ns )*/
prevs = nexts;
prevv = nexts->neigh[nnextv].vert;
nexts = nexts->neigh[nremv].simp;
}
}
if (prevs->status != POLE_OUTPUT)
{
prevs->status = POLE_OUTPUT;
prevs->poleindex = num_vtxs++;
// TJH: PC contains the points for the powercrust surface
// so we hijack the data and pipe it to our structure
{
float vp[3];
vp[0] = prevs->vv[0];
vp[1] = prevs->vv[1];
vp[2] = prevs->vv[2];
output.InsertNextPoint(vp);
}
}
if (numedges < 3)
{
plane[numedges][0] = prevs->vv[0];
plane[numedges][1] = prevs->vv[1];
plane[numedges][2] = prevs->vv[2];
}
sprintf(indface[numedges], "%ld ", prevs->poleindex);
numedges++;
//output.InsertNextCell(numedges);
std::vector<int> indices;
if (!correct_orientation(plane[0], plane[1], plane[2], inpole, outpole))
for (i = numedges - 1; i >= 0; i--)
{
{
//output.GetPolys()->InsertCellPoint(atoi(indface[i]));
indices.push_back(atoi(indface[i]));
}
}
else
{
for (i = 0; i < numedges; i++)
{
{
//output.GetPolys()->InsertCellPoint(atoi(indface[i]));
indices.push_back(atoi(indface[i]));
}
}
}
output.InsertCell(indices);
}
int PowerCrustSurfaceReconstructionImpl::correct_orientation(double *p1, double *p2, double *p3, double *inp, double *outp)
{
double normal[3];
double v1[3], v2[3];
double xcross, ycross, zcross;
int numplus = 0, numminus = 0;
normal[0] = outp[0] - inp[0];
normal[1] = outp[1] - inp[1];
normal[2] = outp[2] - inp[2];
v1[0] = p2[0] - p1[0];
v1[1] = p2[1] - p1[1];
v1[2] = p2[2] - p1[2];
v2[0] = p3[0] - p2[0];
v2[1] = p3[1] - p2[1];
v2[2] = p3[2] - p2[2];
xcross = v1[1] * v2[2] - v1[2] * v2[1];
ycross = v1[2] * v2[0] - v1[0] * v2[2];
zcross = v1[0] * v2[1] - v1[1] * v2[0];
if ((xcross*normal[0]) > 0)
numplus++;
else
numminus++;
if ((ycross*normal[1]) > 0)
numplus++;
else
numminus++;
if ((zcross*normal[2]) > 0)
numplus++;
else
numminus++;
if (numplus > numminus)
return 1;
else
return 0;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment