00001 // Voro++, a 3D cell-based Voronoi library 00002 // 00003 // Author : Chris H. Rycroft (LBL / UC Berkeley) 00004 // Email : [email protected] 00005 // Date : July 1st 2008 00006 00007 /** \file cell.hh 00008 * \brief Header file for the voronoicell_base template and related classes. */ 00009 00010 #ifndef VOROPP_CELL_HH 00011 #define VOROPP_CELL_HH 00012 00013 #include "config.hh" 00014 #include <cstdlib> 00015 #include <iostream> 00016 #include <fstream> 00017 #include <cmath> 00018 using namespace std; 00019 00020 /** \brief Function for printing fatal error messages and exiting. 00021 * 00022 * Function for printing fatal error messages and exiting. 00023 * \param[in] p a pointer to the message to print. 00024 * \param[in] status the status code to return with. */ 00025 void voropp_fatal_error(const char *p,int status) { 00026 cerr << "voro++: " << p << endl; 00027 exit(status); 00028 } 00029 00030 /** \brief A class to reliably carry out floating point comparisons, storing 00031 * marginal cases for future reference. 00032 * 00033 * Floating point comparisons can be unreliable on some processor 00034 * architectures, and can produce unpredictable results. On a number of popular 00035 * Intel processors, floating point numbers are held to higher precision when 00036 * in registers than when in memory. When a register is swapped from a register 00037 * to memory, a truncation error, and in some situations this can create 00038 * circumstances where for two numbers c and d, the program finds c>d first, 00039 * but later c<d. The programmer has no control over when the swaps between 00040 * memory and registers occur, and recompiling with slightly different code can 00041 * give different results. One solution to avoid this is to force the compiler 00042 * to evaluate everything in memory (e.g. by using the -ffloat-store option in 00043 * the GNU C++ compiler) but this could be viewed overkill, since it slows the 00044 * code down, and the extra register precision is useful. 00045 * 00046 * In the plane cutting routine of the voronoicell class, we need to reliably 00047 * know whether a vertex lies inside, outside, or on the cutting plane, since 00048 * if it changed during the tracing process there would be confusion. This 00049 * class makes these tests reliable, by storing the results of marginal cases, 00050 * where the vertex lies within tolerance2 of the cutting plane. If that vertex 00051 * is tested again, then code looks up the value of the table in a buffer, 00052 * rather than doing the floating point comparison again. Only vertices which 00053 * are close to the plane are stored and tested, so this routine should create 00054 * minimal computational overhead. 00055 */ 00056 class suretest { 00057 public: 00058 /** This is a pointer to the array in the voronoicell class 00059 * which holds the vertex coordinates.*/ 00060 fpoint *p; 00061 suretest(); 00062 ~suretest(); 00063 inline void init(fpoint x,fpoint y,fpoint z,fpoint rsq); 00064 inline int test(int n,fpoint &ans); 00065 private: 00066 int check_marginal(int n,fpoint &ans); 00067 /** This stores the current memory allocation for the marginal 00068 * cases. */ 00069 int current_marginal; 00070 /** This stores the total number of marginal points which are 00071 * currently in the buffer. */ 00072 int sc; 00073 /** This array contains a list of the marginal points, and also 00074 * the outcomes of the marginal tests. */ 00075 int *sn; 00076 /** The x coordinate of the normal vector to the test plane. */ 00077 fpoint px; 00078 /** The y coordinate of the normal vector to the test plane. */ 00079 fpoint py; 00080 /** The z coordinate of the normal vector to the test plane. */ 00081 fpoint pz; 00082 /** The magnitude of the normal vector to the test plane. */ 00083 fpoint prsq; 00084 }; 00085 00086 class neighbor_track; 00087 00088 /** \brief A class encapsulating all the routines for storing and calculating 00089 * a single Voronoi cell. 00090 * 00091 * This class encapsulates all the routines for storing and calculating a 00092 * single Voronoi cell. The cell can first be initialized by the init() function 00093 * to be a rectangular box. The box can then be successively cut by planes 00094 * using the plane function. Other routines exist for outputting the cell, 00095 * computing its volume, or finding the largest distance of a vertex from the 00096 * cell center. The cell is described by two arrays. pts[] is a floating point 00097 * array which holds the vertex positions. ed[] holds the table of edges, and 00098 * also a relation table that determines how two vertices are connected to one 00099 * another. The relation table is redundant, but helps speed up the 00100 * computation. The function check_relations() checks that the relational table 00101 * is valid. */ 00102 template <class n_option> 00103 class voronoicell_base { 00104 public: 00105 /** This holds the current size of the arrays ed and nu, which 00106 * hold the vertex information. If more vertices are created 00107 * than can fit in this array, then it is dynamically extended 00108 * using the add_memory_vertices routine. */ 00109 int current_vertices; 00110 /** This holds the current maximum allowed order of a vertex, 00111 * which sets the size of the mem, mep, and mec arrays. If a 00112 * vertex is created with more vertices than this, the arrays 00113 * are dynamically extended using the add_memory_vorder routine. 00114 */ 00115 int current_vertex_order; 00116 /** This sets the size of the main delete stack. */ 00117 int current_delete_size; 00118 /** This sets the size of the auxiliary delete stack. */ 00119 int current_delete2_size; 00120 /** This sets the total number of vertices in the current cell. 00121 */ 00122 int p; 00123 /** This is the index of particular point in the cell, which is 00124 * used to start the tracing routines for plane intersection 00125 * and cutting. These routines will work starting from any 00126 * point, but it's often most efficient to start from the last 00127 * point considered, since in many cases, the cell construction 00128 * algorithm may consider many planes with similar vectors 00129 * concurrently. */ 00130 int up; 00131 /** This is a two dimensional array that holds information 00132 * about the edge connections of the vertices that make up the 00133 * cell. The two dimensional array is not allocated in the 00134 * usual method. To account for the fact the different vertices 00135 * have different orders, and thus require different amounts of 00136 * storage, the elements of ed[i] point to one-dimensional 00137 * arrays in the mep[] array of different sizes. 00138 * 00139 * More specifically, if vertex i has order m, then ed[i] 00140 * points to a one-dimensional array in mep[m] that has 2*m+1 00141 * entries. The first m elements hold the neighboring edges, so 00142 * that the jth edge of vertex i is held in ed[i][j]. The next 00143 * m elements hold a table of relations which is redundant but 00144 * helps speed up the computation. It satisfies the relation 00145 * ed[ed[i][j]][ed[i][m+j]]=i. The final entry holds a back 00146 * pointer, so that ed[i+2*m]=i. These are used when 00147 * rearranging the memory. */ 00148 int **ed; 00149 /** This array holds the order of the vertices in the Voronoi 00150 * cell. This array is dynamically allocated, with its current 00151 * size held by current_vertices. */ 00152 int *nu; 00153 /** This in an array with size 3*current_vertices for holding 00154 * the positions of the vertices. */ 00155 fpoint *pts; 00156 /** This is a class used in the plane routine for carrying out 00157 * reliable comparisons of whether points in the cell are 00158 * inside, outside, or on the current cutting plane. */ 00159 suretest sure; 00160 voronoicell_base(); 00161 ~voronoicell_base(); 00162 void init(fpoint xmin,fpoint xmax,fpoint ymin,fpoint ymax,fpoint zmin,fpoint zmax); 00163 inline void init_octahedron(fpoint l); 00164 inline void init_tetrahedron(fpoint x0,fpoint y0,fpoint z0,fpoint x1,fpoint y1,fpoint z1,fpoint x2,fpoint y2,fpoint z2,fpoint x3,fpoint y3,fpoint z3); 00165 void draw_pov(ostream &os,fpoint x,fpoint y,fpoint z); 00166 inline void draw_pov(const char *filename,fpoint x,fpoint y,fpoint z); 00167 inline void draw_pov(fpoint x,fpoint y,fpoint z); 00168 void draw_pov_mesh(ostream &os,fpoint x,fpoint y,fpoint z); 00169 inline void draw_pov_mesh(const char *filename,fpoint x,fpoint y,fpoint z); 00170 inline void draw_pov_mesh(fpoint x,fpoint y,fpoint z); 00171 void draw_gnuplot(ostream &os,fpoint x,fpoint y,fpoint z); 00172 inline void draw_gnuplot(const char *filename,fpoint x,fpoint y,fpoint z); 00173 inline void draw_gnuplot(fpoint x,fpoint y,fpoint z); 00174 fpoint volume(); 00175 fpoint max_radius_squared(); 00176 fpoint total_edge_distance(); 00177 fpoint surface_area(); 00178 void centroid(fpoint &cx,fpoint &cy,fpoint &cz); 00179 int number_of_faces(); 00180 int number_of_edges(); 00181 void output_vertex_orders(ostream &os); 00182 void output_vertices(ostream &os); 00183 void output_vertices(ostream &os,fpoint x,fpoint y,fpoint z); 00184 void output_face_areas(ostream &os); 00185 void output_face_orders(ostream &os); 00186 void output_face_freq_table(ostream &os); 00187 void output_face_vertices(ostream &os); 00188 void output_face_perimeters(ostream &os); 00189 void output_normals(ostream &os); 00190 void output_neighbors(ostream &os,bool later=false); 00191 bool nplane(fpoint x,fpoint y,fpoint z,fpoint rs,int p_id); 00192 inline bool nplane(fpoint x,fpoint y,fpoint z,int p_id); 00193 inline bool plane(fpoint x,fpoint y,fpoint z,fpoint rs); 00194 inline bool plane(fpoint x,fpoint y,fpoint z); 00195 bool plane_intersects(fpoint x,fpoint y,fpoint z,fpoint rs); 00196 bool plane_intersects_guess(fpoint x,fpoint y,fpoint z,fpoint rs); 00197 inline void init_test(int n); 00198 void add_vertex(fpoint x,fpoint y,fpoint z,int a); 00199 void add_vertex(fpoint x,fpoint y,fpoint z,int a,int b); 00200 void add_vertex(fpoint x,fpoint y,fpoint z,int a,int b,int c); 00201 void add_vertex(fpoint x,fpoint y,fpoint z,int a,int b,int c,int d); 00202 void add_vertex(fpoint x,fpoint y,fpoint z,int a,int b,int c,int d,int e); 00203 void construct_relations(); 00204 void check_relations(); 00205 void check_duplicates(); 00206 void print_edges(); 00207 void label_facets(); 00208 void check_facets(); 00209 inline void perturb(fpoint r); 00210 private: 00211 /** This a one dimensional array that holds the current sizes 00212 * of the memory allocations for them mep array.*/ 00213 int *mem; 00214 /** This is a one dimensional array that holds the current 00215 * number of vertices of order p that are stored in the mep[p] 00216 * array. */ 00217 int *mec; 00218 /** This is a two dimensional array for holding the information 00219 * about the edges of the Voronoi cell. mep[p] is a 00220 * one-dimensional array for holding the edge information about 00221 * all vertices of order p, with each vertex holding 2*p+1 00222 * integers of information. The total number of vertices held 00223 * on mep[p] is stored in mem[p]. If the space runs out, the 00224 * code allocates more using the add_memory() routine. */ 00225 int **mep; 00226 /** This is the delete stack, used to store the vertices which 00227 * are going to be deleted during the plane cutting procedure. 00228 */ 00229 int *ds; 00230 /** This is the auxiliary delete stack, which has size set by 00231 * current_delete2_size. */ 00232 int *ds2; 00233 /** This holds the number of points currently on the auxiliary 00234 * delete stack. */ 00235 int stack2; 00236 /** This object contains all the functions required to carry 00237 * out the neighbor computation. If the neighbor_none class is 00238 * used for n_option, then all these functions are blank. If 00239 * the neighbor_track class is used, then the neighbor tracking 00240 * is enabled. All the functions for the n_option classes are 00241 * declared inline, so that they should all be completely 00242 * integrated into the routine during compilation. */ 00243 n_option neighbor; 00244 inline int cycle_up(int a,int p); 00245 inline int cycle_down(int a,int p); 00246 void add_memory(int i); 00247 void add_memory_vertices(); 00248 void add_memory_vorder(); 00249 void add_memory_ds(); 00250 void add_memory_ds2(); 00251 inline bool collapse_order1(); 00252 inline bool collapse_order2(); 00253 inline bool delete_connection(int j,int k,bool hand); 00254 inline bool plane_intersects_track(fpoint x,fpoint y,fpoint z,fpoint rs,fpoint g); 00255 inline void reset_edges(); 00256 inline void output_normals_search(ostream &os,int i,int j,int k); 00257 friend class neighbor_track; 00258 }; 00259 00260 /** \brief A class passed to the voronoicell_base template to switch off 00261 * neighbor computation. 00262 * 00263 * This is a class full of empty routines for neighbor computation. If the 00264 * voronoicell_base template is instantiated with this class, then it 00265 * has the effect of switching off all neighbor computation. Since all these 00266 * routines are declared inline, it should have the effect of a zero speed 00267 * overhead in the resulting code. */ 00268 class neighbor_none { 00269 public: 00270 /** This is a blank constructor. */ 00271 neighbor_none(voronoicell_base<neighbor_none> *ivc) {}; 00272 /** This is a blank placeholder function that does nothing. */ 00273 inline void allocate(int i,int m) {}; 00274 /** This is a blank placeholder function that does nothing. */ 00275 inline void add_memory_vertices(int i) {}; 00276 /** This is a blank placeholder function that does nothing. */ 00277 inline void add_memory_vorder(int i) {}; 00278 /** This is a blank placeholder function that does nothing. */ 00279 inline void init() {}; 00280 /** This is a blank placeholder function that does nothing. */ 00281 inline void init_octahedron() {}; 00282 /** This is a blank placeholder function that does nothing. */ 00283 inline void init_tetrahedron() {}; 00284 /** This is a blank placeholder function that does nothing. */ 00285 inline void set_pointer(int p,int n) {}; 00286 /** This is a blank placeholder function that does nothing. */ 00287 inline void copy(int a,int b,int c,int d) {}; 00288 /** This is a blank placeholder function that does nothing. */ 00289 inline void set(int a,int b,int c) {}; 00290 /** This is a blank placeholder function that does nothing. */ 00291 inline void set_aux1(int k) {}; 00292 /** This is a blank placeholder function that does nothing. */ 00293 inline void copy_aux1(int a,int b) {}; 00294 /** This is a blank placeholder function that does nothing. */ 00295 inline void copy_aux1_shift(int a,int b) {}; 00296 /** This is a blank placeholder function that does nothing. */ 00297 inline void set_aux2_copy(int a,int b) {}; 00298 /** This is a blank placeholder function that does nothing. */ 00299 inline void copy_pointer(int a,int b) {}; 00300 /** This is a blank placeholder function that does nothing. */ 00301 inline void set_to_aux1(int j) {}; 00302 /** This is a blank placeholder function that does nothing. */ 00303 inline void set_to_aux2(int j) {}; 00304 /** This is a blank placeholder function that does nothing. */ 00305 inline void print_edges(int i) {}; 00306 /** This is a blank placeholder function that does nothing. */ 00307 inline void allocate_aux1(int i) {}; 00308 /** This is a blank placeholder function that does nothing. */ 00309 inline void switch_to_aux1(int i) {}; 00310 /** This is a blank placeholder function that does nothing. */ 00311 inline void copy_to_aux1(int i,int m) {}; 00312 /** This is a blank placeholder function that does nothing. */ 00313 inline void set_to_aux1_offset(int k,int m) {}; 00314 /** This is a blank placeholder function that does nothing. */ 00315 inline void neighbors(ostream &os,bool later) {}; 00316 /** This is a blank placeholder function that does nothing. */ 00317 inline void label_facets() {}; 00318 /** This is a blank placeholder function that does nothing. */ 00319 inline void check_facets() {}; 00320 }; 00321 00322 /** \brief A class passed to the voronoicell_base template to switch on the 00323 * neighbor computation. 00324 * 00325 * This class encapsulates all the routines which are required to carry out the 00326 * neighbor tracking. If the voronoicell_base template is instantiated with 00327 * this class, then the neighbor computation is enabled. All these routines are 00328 * simple and declared inline, so they should be directly integrated into the 00329 * functions in the voronoicell class during compilation, without zero function 00330 * call overhead. */ 00331 class neighbor_track { 00332 public: 00333 /** This two dimensional array holds the neighbor information 00334 * associated with each vertex. mne[p] is a one dimensional 00335 * array which holds all of the neighbor information for 00336 * vertices of order p. */ 00337 int **mne; 00338 /** This is a two dimensional array that holds the neighbor 00339 * information associated with each vertex. ne[i] points to a 00340 * one-dimensional array in mne[nu[i]]. ne[i][j] holds the 00341 * neighbor information associated with the jth edge of vertex 00342 * i. It is set to the ID number of the plane that made the 00343 * face that is clockwise from the jth edge. */ 00344 int **ne; 00345 neighbor_track(voronoicell_base<neighbor_track> *ivc); 00346 ~neighbor_track(); 00347 /** This is a pointer back to the voronoicell class which 00348 * created this class. It is used to reference the members of 00349 * that class in computations. */ 00350 voronoicell_base<neighbor_track> *vc; 00351 inline void allocate(int i,int m); 00352 inline void add_memory_vertices(int i); 00353 inline void add_memory_vorder(int i); 00354 inline void init(); 00355 inline void init_octahedron(); 00356 inline void init_tetrahedron(); 00357 inline void set_pointer(int p,int n); 00358 inline void copy(int a,int b,int c,int d); 00359 inline void set(int a,int b,int c); 00360 inline void set_aux1(int k); 00361 inline void copy_aux1(int a,int b); 00362 inline void copy_aux1_shift(int a,int b); 00363 inline void set_aux2_copy(int a,int b); 00364 inline void copy_pointer(int a,int b); 00365 inline void set_to_aux1(int j); 00366 inline void set_to_aux2(int j); 00367 inline void print_edges(int i); 00368 inline void allocate_aux1(int i); 00369 inline void switch_to_aux1(int i); 00370 inline void copy_to_aux1(int i,int m); 00371 inline void set_to_aux1_offset(int k,int m); 00372 inline void neighbors(ostream &os,bool later); 00373 inline void label_facets(); 00374 inline void check_facets(); 00375 private: 00376 /** This is an auxiliary pointer which is used in some of the 00377 * low level neighbor operations. */ 00378 int *paux1; 00379 /** This is a second auxiliary pointer which is used in some 00380 * of the low level neighbor operations. */ 00381 int *paux2; 00382 }; 00383 00384 /** The basic voronoicell class. */ 00385 typedef voronoicell_base<neighbor_none> voronoicell; 00386 00387 /** A neighbor-tracking version of the voronoicell. */ 00388 typedef voronoicell_base<neighbor_track> voronoicell_neighbor; 00389 #endif
1.6.1