branches/exact/zeo/v_network.cc

   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>
   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. */
  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. */
  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. */
 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. */
 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. */
 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. */
 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. */
 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. */
 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. */
 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. */
 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>
 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. */
 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>
 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>
 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
 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
 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
 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
 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). */
 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. */
 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);