1 #include "v_network.hh"
2
3 /** Initializes the Voronoi network object. The geometry is set up to match a
4 * corresponding container class, and memory is allocated for the network.
5 * \param[in] c a reference to a container or container_poly class. */
6 template<class c_class>
voronoi_network(c_class & c,double net_tol_)7 voronoi_network::voronoi_network(c_class &c,double net_tol_) :
8 bx(c.bx), bxy(c.bxy), by(c.by), bxz(c.bxz), byz(c.byz), bz(c.bz),
9 nx(c.nx), ny(c.ny), nz(c.nz), nxyz(nx*ny*nz),
10 xsp(nx/bx), ysp(ny/by), zsp(nz/bz), net_tol(net_tol_) {
11 int l;
12
13 // Allocate memory for vertex structure
14 pts=new double*[nxyz];
15 idmem=new int*[nxyz];
16 ptsc=new int[nxyz];
17 ptsmem=new int[nxyz];
18 for(l=0;l<nxyz;l++) {
19 pts[l]=new double[4*init_network_vertex_memory];
20 idmem[l]=new int[init_network_vertex_memory];
21 ptsc[l]=0;ptsmem[l]=init_network_vertex_memory;
22 }
23
24 // Allocate memory for network edges and related statistics
25 edc=0;edmem=init_network_vertex_memory*nxyz;
26 ed=new int*[edmem];
27 ne=new int*[edmem];
28 pered=new unsigned int*[edmem];
29 raded=new block*[edmem];
30 nu=new int[edmem];
31 nec=new int[edmem];
32 numem=new int[edmem];
33
34 // Allocate memory for back pointers
35 reg=new int[edmem];
36 regp=new int[edmem];
37
38 // Allocate edge memory
39 for(l=0;l<edmem;l++) {
40 ed[l]=new int[2*init_network_edge_memory];
41 ne[l]=ed[l]+init_network_edge_memory;
42 }
43 for(l=0;l<edmem;l++) raded[l]=new block[init_network_edge_memory];
44 for(l=0;l<edmem;l++) pered[l]=new unsigned int[init_network_edge_memory];
45 for(l=0;l<edmem;l++) {nu[l]=nec[l]=0;numem[l]=init_network_edge_memory;}
46
47 // vertices
48 vmap=new int[4*init_vertices];
49 map_mem=init_vertices;
50 }
51
52 /** The voronoi_network destructor removes the dynamically allocated memory. */
~voronoi_network()53 voronoi_network::~voronoi_network() {
54 int l;
55
56 // Remove Voronoi mapping array
57 delete [] vmap;
58
59 // Remove individual edge arrays
60 for(l=0;l<edmem;l++) delete [] pered[l];
61 for(l=0;l<edmem;l++) delete [] raded[l];
62 for(l=0;l<edmem;l++) delete [] ed[l];
63
64 // Remove back pointers
65 delete [] regp;delete [] reg;
66
67 // Remove memory for edges and related statistics
68 delete [] numem;delete [] nec;delete [] nu;
69 delete [] raded;delete [] pered;
70 delete [] ne;delete [] ed;
71
72 // Remove vertex structure arrays
73 for(l=0;l<nxyz;l++) {
74 delete [] idmem[l];
75 delete [] pts[l];
76 }
77 delete [] ptsmem;delete [] ptsc;
78 delete [] idmem;delete [] pts;
79 }
80
81 /** Increase network memory for a particular region. */
add_network_memory(int l)82 void voronoi_network::add_network_memory(int l) {
83 ptsmem[l]<<=1;
84
85 // Check to see that an absolute maximum in memory allocation
86 // has not been reached, to prevent runaway allocation
87 if(ptsmem[l]>max_network_vertex_memory)
88 voro_fatal_error("Container vertex maximum memory allocation exceeded",VOROPP_MEMORY_ERROR);
89
90 // Allocate new arrays
91 double *npts(new double[4*ptsmem[l]]);
92 int *nidmem(new int[ptsmem[l]]);
93
94 // Copy the contents of the old arrays into the new ones
95 for(int i=0;i<4*ptsc[l];i++) npts[i]=pts[l][i];
96 for(int i=0;i<ptsc[l];i++) nidmem[i]=idmem[l][i];
97
98 // Delete old arrays and update pointers to the new ones
99 delete [] pts[l];delete [] idmem[l];
100 pts[l]=npts;idmem[l]=nidmem;
101 }
102
103 /** Increase edge network memory. */
add_edge_network_memory()104 void voronoi_network::add_edge_network_memory() {
105 int i;
106 edmem<<=1;
107
108 // Allocate new arrays
109 int **ned(new int*[edmem]);
110 int **nne(new int*[edmem]);
111 block **nraded(new block*[edmem]);
112 unsigned int **npered(new unsigned int*[edmem]);
113 int *nnu(new int[edmem]);
114 int *nnec(new int[edmem]);
115 int *nnumem(new int[edmem]);
116 int *nreg(new int[edmem]);
117 int *nregp(new int[edmem]);
118
119 // Copy the contents of the old arrays into the new ones
120 for(i=0;i<edc;i++) {
121 ned[i]=ed[i];
122 nne[i]=ne[i];
123 nraded[i]=raded[i];
124 npered[i]=pered[i];
125 nnu[i]=nu[i];
126 nnec[i]=nec[i];
127 nnumem[i]=numem[i];
128 nreg[i]=reg[i];
129 nregp[i]=regp[i];
130 }
131
132 // Carry out new allocation
133 while(i<edmem) {
134 ned[i]=new int[2*init_network_edge_memory];
135 nne[i]=ned[i]+init_network_edge_memory;
136 nnu[i]=nnec[i]=0;nnumem[i]=init_network_edge_memory;
137 nraded[i]=new block[init_network_edge_memory];
138 npered[i++]=new unsigned int[init_network_edge_memory];
139 }
140
141 // Delete old arrays and update pointers to the new ones
142 delete [] ed;ed=ned;
143 delete [] ne;ne=nne;
144 delete [] raded;raded=nraded;
145 delete [] pered;pered=npered;
146 delete [] nu;nu=nnu;
147 delete [] nec;nec=nnec;
148 delete [] numem;numem=nnumem;
149 delete [] reg;reg=nreg;
150 delete [] regp;regp=nregp;
151 }
152
153 /** Increase a particular vertex memory. */
add_particular_vertex_memory(int l)154 void voronoi_network::add_particular_vertex_memory(int l) {
155 numem[l]<<=1;
156
157 // Check that the vertex allocation does not exceed a maximum safe
158 // limit
159 if(numem[l]>max_vertex_order)
160 voro_fatal_error("Particular vertex maximum memory allocation exceeded",VOROPP_MEMORY_ERROR);
161
162 // Allocate new arrays
163 int *ned(new int[2*numem[l]]);
164 int *nne(ned+numem[l]);
165 block *nraded(new block[numem[l]]);
166 unsigned int *npered(new unsigned int[numem[l]]);
167
168 // Copy the contents of the old arrays into the new ones
169 for(int i=0;i<nu[l];i++) {
170 ned[i]=ed[l][i];
171 nraded[i]=raded[l][i];
172 npered[i]=pered[l][i];
173 }
174 for(int i=0;i<nec[l];i++) nne[i]=ne[l][i];
175
176 // Delete old arrays and update pointers to the new ones
177 delete [] ed[l];ed[l]=ned;ne[l]=nne;
178 delete [] raded[l];raded[l]=nraded;
179 delete [] pered[l];pered[l]=npered;
180 }
181
182
183 /** Increases the memory for the vertex mapping.
184 * \param[in] pmem the amount of memory needed. */
add_mapping_memory(int pmem)185 void voronoi_network::add_mapping_memory(int pmem) {
186 do {map_mem<<=1;} while(map_mem<pmem);
187 delete [] vmap;
188 vmap=new int[4*map_mem];
189 }
190
191 /** Clears the class of all vertices and edges. */
clear_network()192 void voronoi_network::clear_network() {
193 int l;
194 edc=0;
195 for(l=0;l<nxyz;l++) ptsc[l]=0;
196 for(l=0;l<edmem;l++) nu[l]=0;
197 }
198
199 /** Outputs the network in a format that can be read by gnuplot.
200 * \param[in] fp a file handle to write to. */
draw_network(FILE * fp)201 void voronoi_network::draw_network(FILE *fp) {
202 voronoicell c;
203 int l,q,ai,aj,ak;
204 double x,y,z,*ptsp;
205 for(l=0;l<edc;l++) {
206 ptsp=pts[reg[l]]+4*regp[l];
207 x=*(ptsp++);y=*(ptsp++);z=*ptsp;
208 for(q=0;q<nu[l];q++) {
209 unpack_periodicity(pered[l][q],ai,aj,ak);
210 if(ed[l][q]<l&&ai==0&&aj==0&&ak==0) continue;
211 ptsp=pts[reg[ed[l][q]]]+4*regp[ed[l][q]];
212 fprintf(fp,"%g %g %g\n%g %g %g\n\n\n",x,y,z,
213 *ptsp+bx*ai+bxy*aj+bxz*ak,
214 ptsp[1]+by*aj+byz*ak,ptsp[2]+bz*ak);
215 }
216 }
217 }
218
219 /** Prints out the network.
220 * \param[in] fp a file handle to write to.
221 * \param[in] reverse_remove a boolean value, setting whether or not to remove
222 * reverse edges. */
print_network(FILE * fp,bool reverse_remove)223 void voronoi_network::print_network(FILE *fp,bool reverse_remove) {
224 int ai,aj,ak,j,l,ll,q;
225 double x,y,z,x2,y2,z2,*ptsp;
226
227 // Print the vertex table
228 fprintf(fp,"Vertex table:\n%d\n",edc);
229 //os << edc << "\n";
230 for(l=0;l<edc;l++) {
231 ptsp=pts[reg[l]];j=4*regp[l];
232 fprintf(fp,"%d %g %g %g %g",l,ptsp[j],ptsp[j+1],ptsp[j+2],ptsp[j+3]);
233 for(ll=0;ll<nec[l];ll++) fprintf(fp," %d",ne[l][ll]);
234 fputs("\n",fp);
235 }
236
237 // Print out the edge table, loop over vertices
238 fputs("\nEdge table:\n",fp);
239 for(l=0;l<edc;l++) {
240
241 // Store the position of this vertex
242 ptsp=pts[reg[l]];j=4*regp[l];
243 x=ptsp[j];y=ptsp[j+1];z=ptsp[j+2];
244
245 // Loop over edges of this vertex
246 for(q=0;q<nu[l];q++) {
247
248 unpack_periodicity(pered[l][q],ai,aj,ak);
249
250 // If this option is enabled, then the code will not
251 // print edges from i to j for j<i.
252 if(reverse_remove) if(ed[l][q]<l&&ai==0&&aj==0&&ak==0) continue;
253
254 fprintf(fp,"%d -> %d",l,ed[l][q]);
255 raded[l][q].print(fp);
256
257 // Compute and print the length of the edge
258 ptsp=pts[reg[ed[l][q]]];j=4*regp[ed[l][q]];
259 x2=ptsp[j]+ai*bx+aj*bxy+ak*bxz-x;
260 y2=ptsp[j+1]+aj*by+ak*byz-y;
261 z2=ptsp[j+2]+ak*bz-z;
262 fprintf(fp," %d %d %d %g\n",ai,aj,ak,sqrt(x2*x2+y2*y2+z2*z2));
263 }
264 }
265 }
266
267 // Converts three periodic image displacements into a single unsigned integer.
268 // \param[in] i the periodic image in the x direction.
269 // \param[in] j the periodic image in the y direction.
270 // \param[in] k the periodic image in the z direction.
271 // \return The packed integer. */
pack_periodicity(int i,int j,int k)272 inline unsigned int voronoi_network::pack_periodicity(int i,int j,int k) {
273 unsigned int pa=((unsigned int) (127+i));
274 pa<<=8;pa+=((unsigned int) (127+j));
275 pa<<=8;pa+=((unsigned int) (127+k));
276 return pa;
277 }
278
279 /** Unpacks an unsigned integer into three periodic image displacements.
280 * \param[in] pa the packed integer.
281 // \param[out] i the periodic image in the x direction.
282 // \param[out] j the periodic image in the y direction.
283 // \param[out] k the periodic image in the z direction. */
unpack_periodicity(unsigned int pa,int & i,int & j,int & k)284 inline void voronoi_network::unpack_periodicity(unsigned int pa,int &i,int &j,int &k) {
285 i=((signed int) (pa>>16))-127;
286 j=((signed int) ((pa>>8)&255))-127;
287 k=((signed int) (pa&255))-127;
288 }
289
290 /** Adds a Voronoi cell to the network structure. The routine first checks all
291 * of the Voronoi cell vertices and merges them with existing ones where
292 * possible. Edges are then added to the structure.
293 * \param[in] c a reference to a Voronoi cell.
294 * \param[in] (x,y,z) the position of the Voronoi cell.
295 * \param[in] idn the ID number of the particle associated with the cell. */
296 template<class v_cell>
add_to_network_internal(v_cell & c,int idn,double x,double y,double z,double rad,int * cmap)297 void voronoi_network::add_to_network_internal(v_cell &c,int idn,double x,double y,double z,double rad,int *cmap) {
298 int i,j,k,ijk,l,q,ai,aj,ak,*vmp(cmap);
299 double gx,gy,vx,vy,vz,crad,*cp(c.pts);
300
301 // Loop over the vertices of the Voronoi cell
302 for(l=0;l<c.p;l++,vmp+=4) {
303
304 // Compute the real position of this vertex, and evaluate its
305 // position along the non-rectangular axes
306 vx=x+cp[3*l]*0.5;vy=y+cp[3*l+1]*0.5;vz=z+cp[3*l+2]*0.5;
307 gx=vx-vy*(bxy/by)+vz*(bxy*byz-by*bxz)/(by*bz);
308 gy=vy-vz*(byz/bz);
309
310 // Compute the adjusted radius, which will be needed either way
311 crad=0.5*sqrt(cp[3*l]*cp[3*l]+cp[3*l+1]*cp[3*l+1]+cp[3*l+2]*cp[3*l+2])-rad;
312
313 // Check to see if a vertex very close to this one already
314 // exists in the network
315 if(search_previous(gx,gy,vx,vy,vz,ijk,q,vmp[1],vmp[2],vmp[3])) {
316
317 // If it does, then just map the Voronoi cell
318 // vertex to it
319 *vmp=idmem[ijk][q];
320
321 // Store this radius if it smaller than the current
322 // value
323 if(pts[ijk][4*q+3]>crad) pts[ijk][4*q+3]=crad;
324 } else {
325 k=step_int(vz*zsp);if(k<0||k>=nz) {ak=step_div(k,nz);vx-=bxz*ak;vy-=byz*ak;vz-=bz*ak;k-=ak*nz;} else ak=0;
326 j=step_int(gy*ysp);if(j<0||j>=ny) {aj=step_div(j,ny);vx-=bxy*aj;vy-=by*aj;j-=aj*ny;} else aj=0;
327 i=step_int(gx*xsp);if(i<0||i>=nx) {ai=step_div(i,nx);vx-=bx*ai;i-=ai*nx;} else ai=0;
328
329 vmp[1]=ai;vmp[2]=aj;vmp[3]=ak;
330 ijk=i+nx*(j+ny*k);
331
332 if(edc==edmem) add_edge_network_memory();
333 if(ptsc[ijk]==ptsmem[ijk]) add_network_memory(ijk);
334
335 reg[edc]=ijk;regp[edc]=ptsc[ijk];
336 pts[ijk][4*ptsc[ijk]]=vx;
337 pts[ijk][4*ptsc[ijk]+1]=vy;
338 pts[ijk][4*ptsc[ijk]+2]=vz;
339 pts[ijk][4*ptsc[ijk]+3]=crad;
340 idmem[ijk][ptsc[ijk]++]=edc;
341 *vmp=edc++;
342 }
343
344 // Add the neighbor information to this vertex
345 add_neighbor(*vmp,idn);
346 }
347
348 add_edges_to_network(c,x,y,z,rad,cmap);
349 }
350
351 /** Adds a neighboring particle ID to a vertex in the Voronoi network, first
352 * checking that the ID is not already recorded.
353 * \param[in] k the Voronoi vertex.
354 * \param[in] idn the particle ID number. */
add_neighbor(int k,int idn)355 inline void voronoi_network::add_neighbor(int k,int idn) {
356 for(int i=0;i<nec[k];i++) if(ne[k][i]==idn) return;
357 if(nec[k]==numem[k]) add_particular_vertex_memory(k);
358 ne[k][nec[k]++]=idn;
359 }
360
361 /** Adds edges to the network structure, after the vertices have been
362 * considered. This routine assumes that the cmap array has a mapping between
363 * Voronoi cell vertices and Voronoi network vertices. */
364 template<class v_cell>
add_edges_to_network(v_cell & c,double x,double y,double z,double rad,int * cmap)365 void voronoi_network::add_edges_to_network(v_cell &c,double x,double y,double z,double rad,int *cmap) {
366 int i,j,ai,aj,ak,bi,bj,bk,k,l,q,*vmp;unsigned int cper;
367 double vx,vy,vz,wx,wy,wz,dx,dy,dz,dis;double *pp;
368 for(l=0;l<c.p;l++) {
369 vmp=cmap+4*l;k=*(vmp++);ai=*(vmp++);aj=*(vmp++);ak=*vmp;
370 pp=pts[reg[k]]+4*regp[k];
371 vx=pp[0]+ai*bx+aj*bxy+ak*bxz;
372 vy=pp[1]+aj*by+ak*byz;
373 vz=pp[2]+ak*bz;
374 for(q=0;q<c.nu[l];q++) {
375 i=c.ed[l][q];
376 vmp=cmap+4*i;
377 j=*(vmp++);bi=*(vmp++);bj=*(vmp++);bk=*vmp;
378
379 // Skip if this is a self-connecting edge
380 if(j==k&&bi==ai&&bj==aj&&bk==ak) continue;
381 cper=pack_periodicity(bi-ai,bj-aj,bk-ak);
382 pp=pts[reg[j]]+(4*regp[j]);
383 wx=pp[0]+bi*bx+bj*bxy+bk*bxz;
384 wy=pp[1]+bj*by+bk*byz;
385 wz=pp[2]+bk*bz;
386 dx=wx-vx;dy=wy-vy;dz=wz-vz;
387 dis=(x-vx)*dx+(y-vy)*dy+(z-vz)*dz;
388 dis/=dx*dx+dy*dy+dz*dz;
389 if(dis<0) dis=0;if(dis>1) dis=1;
390 wx=vx-x+dis*dx;wy=vy-y+dis*dy;wz=vz-z+dis*dz;
391 int nat=not_already_there(k,j,cper);
392 if(nat==nu[k]) {
393 if(nu[k]==numem[k]) add_particular_vertex_memory(k);
394 ed[k][nu[k]]=j;
395 raded[k][nu[k]].first(sqrt(wx*wx+wy*wy+wz*wz)-rad,dis);
396 pered[k][nu[k]++]=cper;
397 } else {
398 raded[k][nat].add(sqrt(wx*wx+wy*wy+wz*wz)-rad,dis);
399 }
400 }
401 }
402 }
403
404 template<class v_cell>
add_to_network_rectangular_internal(v_cell & c,int idn,double x,double y,double z,double rad,int * cmap)405 void voronoi_network::add_to_network_rectangular_internal(v_cell &c,int idn,double x,double y,double z,double rad,int *cmap) {
406 int i,j,k,ijk,l,q,ai,aj,ak,*vmp(cmap);
407 double vx,vy,vz,crad,*cp(c.pts);
408
409 for(l=0;l<c.p;l++,vmp+=4) {
410 vx=x+cp[3*l]*0.5;vy=y+cp[3*l+1]*0.5;vz=z+cp[3*l+2]*0.5;
411 crad=0.5*sqrt(cp[3*l]*cp[3*l]+cp[3*l+1]*cp[3*l+1]+cp[3*l+2]*cp[3*l+2])-rad;
412 if(safe_search_previous_rect(vx,vy,vz,ijk,q,vmp[1],vmp[2],vmp[3])) {
413 *vmp=idmem[ijk][q];
414
415 // Store this radius if it smaller than the current
416 // value
417 if(pts[ijk][4*q+3]>crad) pts[ijk][4*q+3]=crad;
418 } else {
419 k=step_int(vz*zsp);
420 if(k<0||k>=nz) {
421 ak=step_div(k,nz);
422 vz-=ak*bz;vy-=ak*byz;vx-=ak*bxz;k-=ak*nz;
423 } else ak=0;
424 j=step_int(vy*ysp);
425 if(j<0||j>=ny) {
426 aj=step_div(j,ny);
427 vy-=aj*by;vx-=aj*bxy;j-=aj*ny;
428 } else aj=0;
429 i=step_int(vx*xsp);
430 if(i<0||i>=nx) {
431 ai=step_div(i,nx);
432 vx-=ai*bx;i-=ai*nx;
433 } else ai=0;
434 vmp[1]=ai;
435 vmp[2]=aj;
436 vmp[3]=ak;
437 ijk=i+nx*(j+ny*k);
438 if(edc==edmem) add_edge_network_memory();
439 if(ptsc[ijk]==ptsmem[ijk]) add_network_memory(ijk);
440 reg[edc]=ijk;regp[edc]=ptsc[ijk];
441 pts[ijk][4*ptsc[ijk]]=vx;
442 pts[ijk][4*ptsc[ijk]+1]=vy;
443 pts[ijk][4*ptsc[ijk]+2]=vz;
444 pts[ijk][4*ptsc[ijk]+3]=crad;
445 idmem[ijk][ptsc[ijk]++]=edc;
446 *vmp=edc++;
447 }
448
449 add_neighbor(*vmp,idn);
450 }
451
452 add_edges_to_network(c,x,y,z,rad,cmap);
453 }
454
not_already_there(int k,int j,unsigned int cper)455 int voronoi_network::not_already_there(int k,int j,unsigned int cper) {
456 for(int i=0;i<nu[k];i++) if(ed[k][i]==j&&pered[k][i]==cper) return i;
457 return nu[k];
458 }
459
search_previous(double gx,double gy,double x,double y,double z,int & ijk,int & q,int & pi,int & pj,int & pk)460 bool voronoi_network::search_previous(double gx,double gy,double x,double y,double z,int &ijk,int &q,int &pi,int &pj,int &pk) {
461 int ai=step_int((gx-net_tol)*xsp),bi=step_int((gx+net_tol)*xsp);
462 int aj=step_int((gy-net_tol)*ysp),bj=step_int((gy+net_tol)*ysp);
463 int ak=step_int((z-net_tol)*zsp),bk=step_int((z+net_tol)*zsp);
464 int i,j,k,mi,mj,mk;
465 double px,py,pz,px2,py2,px3,*pp;
466
467 for(k=ak;k<=bk;k++) {
468 pk=step_div(k,nz);px=pk*bxz;py=pk*byz;pz=pk*bz;mk=k-nz*pk;
469 for(j=aj;j<=bj;j++) {
470 pj=step_div(j,ny);px2=px+pj*bxy;py2=py+pj*by;mj=j-ny*pj;
471 for(i=ai;i<=bi;i++) {
472 pi=step_div(i,nx);px3=px2+pi*bx;mi=i-nx*pi;
473 ijk=mi+nx*(mj+ny*mk);
474 pp=pts[ijk];
475 for(q=0;q<ptsc[ijk];q++,pp+=4) if(abs(*pp+px3-x)<net_tol&&abs(pp[1]+py2-y)<net_tol&&abs(pp[2]+pz-z)<net_tol) return true;
476 }
477 }
478 }
479 return false;
480 }
481
safe_search_previous_rect(double x,double y,double z,int & ijk,int & q,int & ci,int & cj,int & ck)482 bool voronoi_network::safe_search_previous_rect(double x,double y,double z,int &ijk,int &q,int &ci,int &cj,int &ck) {
483 const double tol(0.5*net_tol);
484 if(search_previous_rect(x+tol,y+tol,z+tol,ijk,q,ci,cj,ck)) return true;
485 if(search_previous_rect(x-tol,y+tol,z+tol,ijk,q,ci,cj,ck)) return true;
486 if(search_previous_rect(x+tol,y-tol,z+tol,ijk,q,ci,cj,ck)) return true;
487 if(search_previous_rect(x-tol,y-tol,z+tol,ijk,q,ci,cj,ck)) return true;
488 if(search_previous_rect(x+tol,y+tol,z-tol,ijk,q,ci,cj,ck)) return true;
489 if(search_previous_rect(x-tol,y+tol,z-tol,ijk,q,ci,cj,ck)) return true;
490 if(search_previous_rect(x+tol,y-tol,z-tol,ijk,q,ci,cj,ck)) return true;
491 return search_previous_rect(x-tol,y-tol,z-tol,ijk,q,ci,cj,ck);
492 }
493
search_previous_rect(double x,double y,double z,int & ijk,int & q,int & ci,int & cj,int & ck)494 bool voronoi_network::search_previous_rect(double x,double y,double z,int &ijk,int &q,int &ci,int &cj,int &ck) {
495 int k=step_int(z*zsp);
496 if(k<0||k>=nz) {
497 ck=step_div(k,nz);
498 z-=ck*bz;y-=ck*byz;x-=ck*bxz;k-=ck*nz;
499 } else ck=0;
500
501 int j=step_int(y*ysp);
502 if(j<0||j>=ny) {
503 cj=step_div(j,ny);
504 y-=cj*by;x-=cj*bxy;j-=cj*ny;
505 } else cj=0;
506
507 ijk=step_int(x*xsp);
508 if(ijk<0||ijk>=nx) {
509 ci=step_div(ijk,nx);
510 x-=ci*bx;ijk-=ci*nx;
511 } else ci=0;
512
513 ijk+=nx*(j+ny*k);double *pp(pts[ijk]);
514 for(q=0;q<ptsc[ijk];q++,pp+=4)
515 if(abs(*pp-x)<net_tol&&abs(pp[1]-y)<net_tol&&abs(pp[2]-z)<net_tol) return true;
516 return false;
517 }
518
519 /** Custom int function, that gives consistent stepping for negative numbers.
520 * With normal int, we have (-1.5,-0.5,0.5,1.5) -> (-1,0,0,1).
521 * With this routine, we have (-1.5,-0.5,0.5,1.5) -> (-2,-1,0,1). */
step_int(double a)522 inline int voronoi_network::step_int(double a) {
523 return a<0?int(a)-1:int(a);
524 }
525
526 /** Custom integer division function, that gives consistent stepping for
527 * negative numbers. */
step_div(int a,int b)528 inline int voronoi_network::step_div(int a,int b) {
529 return a>=0?a/b:-1+(a+1)/b;
530 }
531
532 // Explicit instantiation
533 template voronoi_network::voronoi_network(container_periodic&, double);
534 template voronoi_network::voronoi_network(container_periodic_poly&, double);
535 template void voronoi_network::add_to_network<voronoicell>(voronoicell&, int, double, double, double, double);
536 template void voronoi_network::add_to_network<voronoicell_neighbor>(voronoicell_neighbor&, int, double, double, double, double);
537 template void voronoi_network::add_to_network_rectangular<voronoicell>(voronoicell&, int, double, double, double, double);
538 template void voronoi_network::add_to_network_rectangular<voronoicell_neighbor>(voronoicell_neighbor&, int, double, double, double, double);
539