EGSnrc C++ class library  Report PIRS-898 (2021)
Iwan Kawrakow, Ernesto Mainegra-Hing, Frederic Tessier, Reid Townson and Blake Walters
egs_octree.h
Go to the documentation of this file.
1 /*
2 ###############################################################################
3 #
4 # EGSnrc egs++ octree geometry headers
5 # Copyright (C) 2015 National Research Council Canada
6 #
7 # This file is part of EGSnrc.
8 #
9 # EGSnrc is free software: you can redistribute it and/or modify it under
10 # the terms of the GNU Affero General Public License as published by the
11 # Free Software Foundation, either version 3 of the License, or (at your
12 # option) any later version.
13 #
14 # EGSnrc is distributed in the hope that it will be useful, but WITHOUT ANY
15 # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
16 # FOR A PARTICULAR PURPOSE. See the GNU Affero General Public License for
17 # more details.
18 #
19 # You should have received a copy of the GNU Affero General Public License
20 # along with EGSnrc. If not, see <http://www.gnu.org/licenses/>.
21 #
22 ###############################################################################
23 #
24 # Author: Frederic Tessier, 2008
25 #
26 # Contributors: Iwan Kawrakow
27 # Hubert Ho
28 #
29 ###############################################################################
30 */
31 
32 
33 /*
34 ===============================================================================
35 egs_octree: an octree implementation to the egs++ geometry collection.
36 ===============================================================================
37 
38 An octree is a way to partition space in a tree-like fashion, and it is a direct
39 extension of binary trees to 3D. The idea is to start with a cuboid and divide
40 it in half along each of the three cartesian axes to yield 8 children cuboids
41 (hence the name "octree"), and so on recursively until some desired limit. One
42 can restrict this subdivision scheme for any branch of the octree to obtain
43 leaf cells of varying size: the resolution of the octree can vary spatially.
44 
45 The main motivation for coding an egs++ octree geometry is to increase the
46 performance of large voxelized geometries such as a the human phantom. If we
47 imagine individual voxels in the phantom to be the leaves of the octree, then we
48 can merge voxels of same medium into larger cells and significantly reduce the
49 overall number of regions, thus reducing the number of howfar calls.
50 
51 The geometry can handle multiple bounding boxes in order to model different part
52 of the geometry at different resolutions (see example below). This can be used
53 to model micromatrices, for example.
54 
55 Any egs++ geometry can be "octreefied". One defines a child geometry, a bounding
56 box and the desired resolution, and the geometry will grow and prune the octree
57 accordingly; the child geometry is simply used as an oracle that returns the
58 local medium index for a given point. This approach allows one to use any
59 existing input file and wrap an octree around it, so in effect in can also serve
60 as a generic voxelizer (there is an option to disable pruning).
61 
62 More documentation will follow in the code, but you should be able to try it
63 right away with your own geometries, following the examples I will commit to
64 the egs++/geometry/examples. Please take it out for a spin on your own favorite
65 geometry and let me know if it works well.
66 
67 Enjoy.
68 
69 
70 ===============================================================================
71 EXAMPLE:
72 ===============================================================================
73 
74 Here is a simple example that defines a sphere and then models it with an octree
75 in which the positive octant features a higher resolution:
76 
77 :start geometry definition:
78 
79  :start geometry:
80  library = egs_spheres
81  name = my_sphere
82  :start media input:
83  media = my_medium
84  :stop media input:
85  midpoint = 0 0 0
86  radii = 0.9
87  :stop geometry:
88 
89  :start geometry:
90  library = egs_octree
91  name = my_octree
92  :start octree box:
93  box min = -1 -1 -1
94  box max = 1 1 1
95  resolution = 32 32 32
96  :stop octree box:
97  :start octree box:
98  box min = 0 0 0
99  box max = 1 1 1
100  resolution = 128 128 128
101  :stop octree box:
102  child geometry = my_sphere
103  :stop geometry:
104 
105  simulation geometry = my_octree
106 
107 :stop geometry definition:
108 
109 
110 ===============================================================================
111 PERFORMANCE:
112 ===============================================================================
113 
114 I tested the octree geometry on the fax06 human phantom to check what kind of
115 gains can be achieved. Using the howfar time test from the geometry tester, I
116 found the following:
117 
118 fax06 human phantom:
119 -------------------------------------------------------------------------------
120 model # cells avg cell size avg # steps time / sample (on tux)
121 -------------------------------------------------------------------------------
122 voxels 141716520 0.12 284 14 us
123 octree 9045359 0.30 46 20 us
124 -------------------------------------------------------------------------------
125 GAIN 94 % 84 % 32 %
126 -------------------------------------------------------------------------------
127 
128 There is a major gain in the number of cells and the number of howfar calls,
129 but the performance increase is mitigated by the more intricate logic for
130 navigating in the octree (finding neighbors etc.). So the performance increase
131 is only 32%, as opposed to the 84% that could be expected solely from the
132 decrease of the number of steps. Of course, all these numbers are model
133 specific, and the gain will be better for geometries in which there are large
134 volumes of uniform medium.
135 
136 
137 ===============================================================================
138 MEMORY:
139 ===============================================================================
140 
141 Each node in the final octree (the number of reported regions) requires 46
142 bytes of memory (on 64-bit machines), that is, 38 bytes for each instance of
143 the node EGS_OCtree_node class and 8 bytes for a pointer to the node in the
144 region list. So for example the fax06 phantom requires about 400 MB of memory
145 (but it shows up as using about 600 MB, so maybe there is a leak?).
146 
147 
148 ===============================================================================
149 TODO:
150 ===============================================================================
151 
152 1. Add an option in the VHP geometry to represent the model with an octree, and
153 perhaps implement a bottom-up octree growth algorithm if it can be faster than
154 the current recursive method.
155 
156 2. Add the option to merge voxels when a certain proportion of the children
157 nodes have the same medium.
158 
159 3. Add the option to set the resolution depending on the medium or other
160 criteria than bounding boxes?
161 
162 =============================================================================
163 */
164 
165 
171 #ifndef EGS_OCTREE_
172 #define EGS_OCTREE_
173 
174 #include "egs_functions.h"
175 #include "egs_base_geometry.h"
176 #include <vector>
177 #include <cstring>
178 #include <cmath>
179 
180 
181 #ifdef WIN32
182 
183  #ifdef BUILD_OCTREE_DLL
184  #define EGS_OCTREE_EXPORT __declspec(dllexport)
185  #else
186  #define EGS_OCTREE_EXPORT __declspec(dllimport)
187  #endif
188  #define EGS_OCTREE_LOCAL
189 
190 #else
191 
192  #ifdef HAVE_VISIBILITY
193  #define EGS_OCTREE_EXPORT __attribute__ ((visibility ("default")))
194  #define EGS_OCTREE_LOCAL __attribute__ ((visibility ("hidden")))
195  #else
196  #define EGS_OCTREE_EXPORT
197  #define EGS_OCTREE_LOCAL
198  #endif
199 
200 #endif
201 
202 
204 public:
205  EGS_Vector vmin, vmax;
206  int level, maxlevel;
207  int ixmin, iymin, izmin;
208  int ixmax, iymax, izmax;
209  int nx, ny, nz;
210 
211  // constructor
212  EGS_Octree_bbox(const EGS_Vector &boxMin, const EGS_Vector &boxMax, vector<int> &bboxRes) :
213  vmin(boxMin), vmax(boxMax), nx(bboxRes[0]), ny(bboxRes[1]), nz(bboxRes[2]) {
214  maxlevel = 0;
215  ixmin = iymin = izmin = 0;
216  ixmax = iymax = izmax = 0;
217  }
218 
219  // addition
220  EGS_Octree_bbox operator+ (const EGS_Octree_bbox &b2) const {
221 
222  EGS_Octree_bbox b1(*this);
223 
224  // set the b1 box limits to enclose the two boxes
225  if (b2.vmin.x < b1.vmin.x) {
226  b1.vmin.x = b2.vmin.x;
227  }
228  if (b2.vmin.y < b1.vmin.y) {
229  b1.vmin.y = b2.vmin.y;
230  }
231  if (b2.vmin.z < b1.vmin.z) {
232  b1.vmin.z = b2.vmin.z;
233  }
234  if (b2.vmax.x > b1.vmax.x) {
235  b1.vmax.x = b2.vmax.x;
236  }
237  if (b2.vmax.y > b1.vmax.y) {
238  b1.vmax.y = b2.vmax.y;
239  }
240  if (b2.vmax.z > b1.vmax.z) {
241  b1.vmax.z = b2.vmax.z;
242  }
243 
244  // calculate voxel sizes for the two boxes
245  EGS_Float dx1 = (b1.vmax.x-b1.vmin.x)/b1.nx;
246  EGS_Float dy1 = (b1.vmax.y-b1.vmin.y)/b1.ny;
247  EGS_Float dz1 = (b1.vmax.z-b1.vmin.z)/b1.nz;
248  EGS_Float dx2 = (b2.vmax.x-b2.vmin.x)/b2.nx;
249  EGS_Float dy2 = (b2.vmax.y-b2.vmin.y)/b2.ny;
250  EGS_Float dz2 = (b2.vmax.z-b2.vmin.z)/b2.nz;
251 
252  // calculate the number of voxels in the overall box b1 (using the highest resolution of b1 and b2)
253  if (dx2 < dx1) {
254  b1.nx = (int)((b1.vmax.x - b1.vmin.x + 0.5*dx2) / dx2);
255  }
256  else {
257  b1.nx = (int)((b1.vmax.x - b1.vmin.x +0.5*dx1) / dx1);
258  }
259  if (dy2 < dy1) {
260  b1.ny = (int)((b1.vmax.y - b1.vmin.y + 0.5*dy2) / dy2);
261  }
262  else {
263  b1.ny = (int)((b1.vmax.y - b1.vmin.y +0.5*dy1) / dy1);
264  }
265  if (dz2 < dz1) {
266  b1.nz = (int)((b1.vmax.z - b1.vmin.z + 0.5*dz2) / dz2);
267  }
268  else {
269  b1.nz = (int)((b1.vmax.z - b1.vmin.z +0.5*dz1) / dz1);
270  }
271 
272  // return the overall box b1
273  return b1;
274  }
275 
276  EGS_Octree_bbox &operator+= (const EGS_Octree_bbox &b2) {
277 
278  // set the b1 box limits to enclose the two boxes
279  if (b2.vmin.x < vmin.x) {
280  vmin.x = b2.vmin.x;
281  }
282  if (b2.vmin.y < vmin.y) {
283  vmin.y = b2.vmin.y;
284  }
285  if (b2.vmin.z < vmin.z) {
286  vmin.z = b2.vmin.z;
287  }
288  if (b2.vmax.x > vmax.x) {
289  vmax.x = b2.vmax.x;
290  }
291  if (b2.vmax.y > vmax.y) {
292  vmax.y = b2.vmax.y;
293  }
294  if (b2.vmax.z > vmax.z) {
295  vmax.z = b2.vmax.z;
296  }
297 
298  // calculate voxel sizes for the two boxes
299  EGS_Float dx1 = (vmax.x-vmin.x)/nx;
300  EGS_Float dy1 = (vmax.y-vmin.y)/ny;
301  EGS_Float dz1 = (vmax.z-vmin.z)/nz;
302  EGS_Float dx2 = (b2.vmax.x-b2.vmin.x)/b2.nx;
303  EGS_Float dy2 = (b2.vmax.y-b2.vmin.y)/b2.ny;
304  EGS_Float dz2 = (b2.vmax.z-b2.vmin.z)/b2.nz;
305 
306  // calculate the number of voxels in the overall box b1 (using the highest resolution of b1 and b2)
307  if (dx2 < dx1) {
308  nx = (int)((vmax.x - vmin.x + 0.5*dx2) / dx2);
309  }
310  else {
311  nx = (int)((vmax.x - vmin.x +0.5*dx1) / dx1);
312  }
313  if (dy2 < dy1) {
314  ny = (int)((vmax.y - vmin.y + 0.5*dy2) / dy2);
315  }
316  else {
317  ny = (int)((vmax.y - vmin.y +0.5*dy1) / dy1);
318  }
319  if (dz2 < dz1) {
320  nz = (int)((vmax.z - vmin.z + 0.5*dz2) / dz2);
321  }
322  else {
323  nz = (int)((vmax.z - vmin.z +0.5*dz1) / dz1);
324  }
325 
326  // return the overall box
327  return *this;
328  }
329 };
330 
331 
333 public:
334  int ix, iy, iz;
335  int region;
336  short medium;
337  unsigned short level;
340 
341  // constructor
342  EGS_Octree_node() {
343  medium = -1; // set the medium to -1 by default
344  region = -1; // set region to -1 by default
345  level = 0; // set level to 0 (root) by default
346  ix = iy = iz = 0; // set octree indices to 0 by default
347  child = NULL; // the node has no children initially
348  parent = NULL; // the node has no parent initially
349  }
350 
351  // create children nodes
352  void createChildren() { // create children to the this node
353  if (!child) { // ensure there are no children already
354  child = new EGS_Octree_node [8]; // allocate memory for 8 new children nodes
355  if (!child) {
356  egsFatal("EGS_Octree_node::createChildren(): Memory allocation error");
357  }
358  for (int i=0; i<8; i++) { // loop over all 8 newly created children nodes
359  child[i].level = level+1; // increase level by 1 compared to current level
360  child[i].ix = (ix << 1) | (i>>0 & 0x1); // shift up ix by one, and set new bit to that of the child index's bit 0 (x position)
361  child[i].iy = (iy << 1) | (i>>1 & 0x1); // shift up iy by one, and set new bit to that of the child index's bit 1 (y position)
362  child[i].iz = (iz << 1) | (i>>2 & 0x1); // shift up iz by one, and set new bit to that of the child index's bit 2 (z position)
363  child[i].parent = this; // set the parent pointer of every child to the current node
364  }
365  }
366  }
367 
368  // delete children (recursively)
369  void deleteChildren() { // delete children of this node
370  if (child) { // if this node has children, delete them
371  for (int i=0; i<8; i++) {
372  child[i].deleteChildren(); // recursive calls to delete all children branches
373  }
374  delete [] child; // free the memory allocated for the 8 children
375  child = NULL; // set children pointer to NULL to indicate that there are no children now
376  }
377  }
378 
379  // collapse node
380  int collapseChildren() { // collapse children in this node if they are all the same medium
381  if (child) { // check that we indeed have children
382  for (int i=0; i<8; i++) { // loop over each of the 8 children of this node
383  if (child[i].child) {
384  return 0; // bail out if there are nodes below the children (this could be made recursive)
385  }
386  if (child[i].medium!=child[0].medium) {
387  return 0; // bail out as soon as one children has a different medium
388  }
389  }
390  medium = child[0].medium; // set the medium of this node to that of the children
391  deleteChildren(); // delete the children of this node
392  return 1; // return 1 to indicate that the children were collapsed
393  }
394  return 0;
395  }
396 
397  // insideBBox
398  bool insideBBox(EGS_Octree_bbox &bbox) {
399  int shift = bbox.maxlevel - level;
400  int ii;
401  ii = (ix<<shift);
402  if (ii < bbox.ixmin) {
403  return false;
404  }
405  ii = ~(~ix<<shift);
406  if (ii > bbox.ixmax) {
407  return false;
408  }
409  ii = (iy<<shift);
410  if (ii < bbox.iymin) {
411  return false;
412  }
413  ii = ~(~iy<<shift);
414  if (ii > bbox.iymax) {
415  return false;
416  }
417  ii = iz<<shift;
418  if (ii < bbox.izmin) {
419  return false;
420  }
421  ii = ~(~iz<<shift);
422  if (ii > bbox.izmax) {
423  return false;
424  }
425  return true;
426  }
427 
428  // instersectBBox
429  bool intersectBBox(EGS_Octree_bbox &bbox) {
430 
431  int shift = bbox.maxlevel - level;
432  int iimin, iimax;
433 
434  iimin = (ix<<shift);
435  if (bbox.ixmin > iimin) {
436  iimin = bbox.ixmin;
437  }
438  iimax = ~(~ix<<shift);
439  if (bbox.ixmax < iimax) {
440  iimax = bbox.ixmax;
441  }
442  if (iimax < iimin) {
443  return false;
444  }
445 
446  iimin = (iy<<shift);
447  if (bbox.iymin > iimin) {
448  iimin = bbox.iymin;
449  }
450  iimax = ~(~iy<<shift);
451  if (bbox.iymax < iimax) {
452  iimax = bbox.iymax;
453  }
454  if (iimax < iimin) {
455  return false;
456  }
457 
458  iimin = (iz<<shift);
459  if (bbox.izmin > iimin) {
460  iimin = bbox.izmin;
461  }
462  iimax = ~(~iz<<shift);
463  if (bbox.izmax < iimax) {
464  iimax = bbox.izmax;
465  }
466  if (iimax < iimin) {
467  return false;
468  }
469 
470  return true;
471  }
472 };
473 
474 
518 class EGS_OCTREE_EXPORT EGS_Octree : public EGS_BaseGeometry {
519 
520  EGS_Octree_node *root;
521  EGS_Octree_node **nodeReg;
522  EGS_BaseGeometry *geom;
523  EGS_Float bbxmin, bbymin, bbzmin;
524  EGS_Float bbxmax, bbymax, bbzmax;
525  EGS_Float xmin, ymin, zmin;
526  EGS_Float xmax, ymax, zmax;
527  EGS_Float dx, dy, dz, dxi, dyi, dzi;
528  int ixmin, iymin, izmin;
529  int ixmax, iymax, izmax;
530  int maxlevel, n;
531  int nx, ny, nz;
532  static string type;
533  long int nLeaf, nLeafMax;
534  vector<EGS_Octree_node *> tmp;
535 
536 public:
537 
538  EGS_Octree(vector<EGS_Octree_bbox> &vBox, bool pruneTree, EGS_BaseGeometry *g) : EGS_BaseGeometry(""), geom(g) {
539 
540  // combine bounding boxes to get the overall bounding box
541  EGS_Octree_bbox bbox(vBox[0]);
542  for (int i=1; i<vBox.size(); i++) {
543  bbox += vBox[i];
544  }
545  vBox.push_back(bbox);
546 
547  // save overall resolution
548  nx = bbox.nx;
549  ny = bbox.ny;
550  nz = bbox.nz;
551 
552  // set size of cells at maxlevel
553  dx = (bbox.vmax.x - bbox.vmin.x) / bbox.nx;
554  dy = (bbox.vmax.y - bbox.vmin.y) / bbox.ny;
555  dz = (bbox.vmax.z - bbox.vmin.z) / bbox.nz;
556 
557  // pre-compute cell size inverses to save time
558  dxi = 1.0/dx;
559  dyi = 1.0/dy;
560  dzi = 1.0/dz;
561 
562  // set octree depth (maxlevel)
563  int res = nx;
564  if (ny > res) {
565  res = ny;
566  }
567  if (nz > res) {
568  res = nz;
569  }
570  maxlevel = (int) ceil(log((EGS_Float)res)/0.6931471805599452862);
571 
572  // set octree cell count at maxlevel;
573  n = 1<<maxlevel;
574 
575  // set octree bounds in space
576  xmin = bbox.vmin.x;
577  ymin = bbox.vmin.y;
578  zmin = bbox.vmin.z;
579  xmax = bbox.vmin.x + dx*n;
580  ymax = bbox.vmin.y + dy*n;
581  zmax = bbox.vmin.z + dz*n;
582 
583  // set cell indices range for each box;
584  {
585  for (int i=0; i<vBox.size(); i++) {
586  EGS_Octree_bbox *box = &vBox[i];
587  box->ixmin = (int)((box->vmin.x - xmin + 0.5*dx) * dxi);
588  box->iymin = (int)((box->vmin.y - ymin + 0.5*dy) * dyi);
589  box->izmin = (int)((box->vmin.z - zmin + 0.5*dz) * dzi);
590  box->nx = (int)((box->vmax.x - box->vmin.x + 0.5*dx) * dxi);
591  box->ny = (int)((box->vmax.y - box->vmin.y + 0.5*dy) * dyi);
592  box->nz = (int)((box->vmax.z - box->vmin.z + 0.5*dz) * dzi);
593  box->ixmax = box->ixmin + box->nx - 1;
594  box->iymax = box->iymin + box->ny - 1;
595  box->izmax = box->izmin + box->nz - 1;
596  if (box->ixmin<0) {
597  box->ixmin = 0;
598  }
599  if (box->ixmax>=n) {
600  box->ixmax = n-1;
601  }
602  if (box->iymin<0) {
603  box->iymin = 0;
604  }
605  if (box->iymax>=n) {
606  box->iymax = n-1;
607  }
608  if (box->izmin<0) {
609  box->izmin = 0;
610  }
611  if (box->izmax>=n) {
612  box->izmax = n-1;
613  }
614  box->nx = box->ixmax - box->ixmin + 1;
615  box->ny = box->iymax - box->iymin + 1;
616  box->nz = box->izmax - box->izmin + 1;
617  }
618  }
619 
620 
621  // determine maxlevel in each defined bounding box
622 // {for (int i=0; i<vBox.size(); i++) {
623 // EGS_Octree_bbox *box = &vBox[i];
624 // EGS_Float scale = box->nxtree/(EGS_Float)box->nx;
625 // if (box->nytree/box->ny > scale) scale = box->nytree/(EGS_Float)box->ny;
626 // if (box->nztree/box->nz > scale) scale = box->nztree/(EGS_Float)box->nz;
627 // int difflevel = (int) ceil (log(scale)/0.6931471805599452862);
628 // box->maxlevel = maxlevel;
629 // box->level = maxlevel-difflevel;
630 // }}
631 
632  // avoid indirections for overall bounding box parameters
633  EGS_Octree_bbox *box = &(vBox.back());
634  bbxmin = box->vmin.x;
635  bbxmax = box->vmax.x;
636  bbymin = box->vmin.y;
637  bbymax = box->vmax.y;
638  bbzmin = box->vmin.z;
639  bbzmax = box->vmax.z;
640  ixmin = box->ixmin;
641  ixmax = box->ixmax;
642  iymin = box->iymin;
643  iymax = box->iymax;
644  izmin = box->izmin;
645  izmax = box->izmax;
646 
647  // build octree
648  root = new EGS_Octree_node();
649  growOctree(root, vBox, pruneTree);
650 
651  // allocate memory to hold the region-indexed node pointers in a simple array, and free up the tmp vector
652  nreg = tmp.size();
653  nodeReg = new EGS_Octree_node* [nreg];
654  {
655  for (int i=0; i<nreg; i++) {
656  nodeReg[i] = tmp[i];
657  }
658  }
659  tmp.erase(tmp.begin(),tmp.end());
660 
661  // calculate leaf node statistics
662  nLeaf = 0;
663  nLeafMax = 0;
664  statOctree(root, vBox);
665  }
666 
667 
668  // destructor
669  ~EGS_Octree() {
670  if (root) {
671  root->deleteChildren();
672  delete root;
673  delete [] nodeReg;
674  }
675  }
676 
677 
678  // statOctree
679  void statOctree(EGS_Octree_node *node, vector<EGS_Octree_bbox> &vBox) {
680  if (node->child) {
681  for (int i=0; i<8; i++) {
682  statOctree(node->child+i, vBox);
683  }
684  }
685  else if (node->insideBBox(vBox[vBox.size()-1]) && node->medium>=0) {
686  int shift = maxlevel - node->level;
687  nLeaf++;
688  nLeafMax += 1<<(3*shift);
689  }
690  }
691 
692 
693  // growOctree (recursive)
694  void growOctree(EGS_Octree_node *node, vector<EGS_Octree_bbox> &vBox, bool prune) {
695 
696  // assume node needs no refinement
697  bool refineNode = false;
698 
699  // check if we need to refine this node
700  if (node->level < maxlevel) {
701 
702  // check all bounding boxes to determine local refinement level (priority to boxes defined later)
703  int level = 0;
704  {
705  for (int i=0; i<vBox.size()-1; i++) {
706  int boxlevel = vBox[i].level;
707  if (node->insideBBox(vBox[i])) {
708  level = boxlevel;
709  }
710  }
711  }
712 
713  // node may still need refinement if it intersects a higher resolution bounding box
714  {
715  for (int i=0; i<vBox.size()-1; i++) {
716  int boxlevel = vBox[i].level;
717  if (!node->insideBBox(vBox[i]) && node->intersectBBox(vBox[i])) {
718  if (boxlevel > level) {
719  level = boxlevel;
720  }
721  }
722  }
723  }
724 
725  // if the level of this node is not sufficient, set it for refinement
726  if (node->level < level) {
727  refineNode = true;
728  }
729  }
730 
731  // if this node needs refinement, grow children branches (recursively), and then try to collapse children
732  if (refineNode) {
733  node->createChildren();
734  for (int i=0; i<8; i++) {
735  growOctree(node->child+i, vBox, prune);
736  }
737  if (prune) {
738  if (node->collapseChildren()) {
739  tmp.resize(tmp.size()-8);
740  node->region = tmp.size();
741  tmp.push_back(node);
742  }
743  }
744  }
745  // otherwise this node is a leaf: set the node medium and add it to the region list
746  else {
747  bool insideSomeBox = false;
748  for (int i=0; i<vBox.size()-1; i++) {
749  if (node->insideBBox(vBox[i])) {
750  insideSomeBox = true;
751  break;
752  }
753  }
754  node->medium = -1;
755  if (insideSomeBox) {
756  int shift = maxlevel - node->level;
757  int ix = node->ix << shift;
758  int iy = node->iy << shift;
759  int iz = node->iz << shift;
760  int ireg = geom->isWhere(EGS_Vector(xmin+(ix+0.5)*dx, ymin+(iy+0.5)*dy, zmin+(iz+0.5)*dz));
761  if (ireg>=0) {
762  node->medium = geom->medium(ireg);
763  }
764  }
765  node->region = tmp.size();
766  tmp.push_back(node);
767  }
768  }
769 
770 
771  // getNeighborNodeX
772  EGS_Octree_node *getNeighborNodeX(EGS_Octree_node *node, int ixn, int iyn, int izn) {
773 
774  // check if neighbor index is in range
775  if (ixn < ixmin || ixn > ixmax) {
776  return NULL;
777  }
778 
779  // constrain neighbor indices in y and z
780  int shift = maxlevel - node->level;
781  if (iyn < node->iy<<shift) {
782  iyn = node->iy<<shift;
783  }
784  else if (iyn > ~(~node->iy<<shift)) {
785  iyn = ~(~node->iy<<shift);
786  }
787  if (izn < node->iz<<shift) {
788  izn = node->iz<<shift;
789  }
790  else if (izn > ~(~node->iz<<shift)) {
791  izn = ~(~node->iz<<shift);
792  }
793 
794  // walk up and down the octree to new cell
795  int diff = node->ix ^ (ixn>>shift);
796  shift = 0;
797  while (diff & (1<<shift)) {
798  node = node->parent;
799  shift++;
800  }
801  shift = maxlevel - node->level;
802  while (node->child) {
803  shift--;
804  int childIndex;
805  childIndex = (ixn>>shift & 0x1);
806  childIndex |= (iyn>>shift & 0x1) << 1;
807  childIndex |= (izn>>shift & 0x1) << 2;
808  node = node->child + childIndex;
809  }
810  shift = maxlevel - node->level;
811 
812  return node;
813  }
814 
815 
816  // getNeighborNodeY
817  EGS_Octree_node *getNeighborNodeY(EGS_Octree_node *node, int ixn, int iyn, int izn) {
818 
819  // check if neighbor index is in range
820  if (iyn < iymin || iyn > iymax) {
821  return NULL;
822  }
823 
824  // constrain neighbor indices in x and z
825  int shift = maxlevel - node->level;
826  if (ixn < node->ix<<shift) {
827  ixn = node->ix<<shift;
828  }
829  else if (ixn > ~(~node->ix<<shift)) {
830  ixn = ~(~node->ix<<shift);
831  }
832  if (izn < node->iz<<shift) {
833  izn = node->iz<<shift;
834  }
835  else if (izn > ~(~node->iz<<shift)) {
836  izn = ~(~node->iz<<shift);
837  }
838 
839  // walk up and down the octree to new cell
840  int diff = node->iy ^ (iyn>>shift);
841  shift = 0;
842  while (diff & (1<<shift)) {
843  node = node->parent;
844  shift++;
845  }
846  shift = maxlevel - node->level;
847  while (node->child) {
848  shift--;
849  int childIndex;
850  childIndex = (ixn>>shift & 0x1);
851  childIndex |= (iyn>>shift & 0x1) << 1;
852  childIndex |= (izn>>shift & 0x1) << 2;
853  node = node->child + childIndex;
854  }
855  return node;
856  }
857 
858 
859  // getNeighborNodeZ
860  EGS_Octree_node *getNeighborNodeZ(EGS_Octree_node *node, int ixn, int iyn, int izn) {
861 
862  // check if neighbor index is in range
863  if (izn < izmin || izn > izmax) {
864  return NULL;
865  }
866 
867  // constrain neighbor indices in x and y
868  int shift = maxlevel - node->level;
869  if (ixn < node->ix<<shift) {
870  ixn = node->ix<<shift;
871  }
872  else if (ixn > ~(~node->ix<<shift)) {
873  ixn = ~(~node->ix<<shift);
874  }
875  if (iyn < node->iy<<shift) {
876  iyn = node->iy<<shift;
877  }
878  else if (iyn > ~(~node->iy<<shift)) {
879  iyn = ~(~node->iy<<shift);
880  }
881 
882  // walk up and down the octree to new cell
883  int diff = node->iz ^ (izn>>shift);
884  shift = 0;
885  while (diff & (1<<shift)) {
886  node = node->parent;
887  shift++;
888  }
889  shift = maxlevel - node->level;
890  while (node->child) {
891  shift--;
892  int childIndex;
893  childIndex = (ixn>>shift & 0x1);
894  childIndex |= (iyn>>shift & 0x1) << 1;
895  childIndex |= (izn>>shift & 0x1) << 2;
896  node = node->child + childIndex;
897  }
898  return node;
899  }
900 
901 
902  // getNode
903  EGS_Octree_node *getNode(int ix, int iy, int iz) {
904  if ((ix<ixmin) || (ix>ixmax) || (iy<iymin) || (iy>iymax) || (iz<izmin) || (iz>izmax)) {
905  return NULL;
906  }
907  EGS_Octree_node *node = root;
908  int shift = maxlevel;
909  while (node->child) {
910  shift--;
911  int childIndex;
912  childIndex = (ix>>shift & 0x1);
913  childIndex |= (iy>>shift & 0x1) << 1;
914  childIndex |= (iz>>shift & 0x1) << 2;
915  node = node->child + childIndex;
916  }
917  return node;
918  }
919 
920 
921  // setIndices
922  void setIndices(const EGS_Vector &r, int &ix, int &iy, int &iz) {
923  ix = (int)((r.x-xmin)*dxi);
924  iy = (int)((r.y-ymin)*dyi);
925  iz = (int)((r.z-zmin)*dzi);
926  if (ix<ixmin) {
927  ix=ixmin;
928  }
929  if (ix>ixmax) {
930  ix=ixmax;
931  }
932  if (iy<iymin) {
933  iy=iymin;
934  }
935  if (iy>iymax) {
936  iy=iymax;
937  }
938  if (iz<izmin) {
939  iz=izmin;
940  }
941  if (iz>izmax) {
942  iz=izmax;
943  }
944  }
945 
946 
947  // isInside
948  bool isInside(const EGS_Vector &r) {
949  if (r.x >= bbxmin && r.x <= bbxmax &&
950  r.y >= bbymin && r.y <= bbymax &&
951  r.z >= bbzmin && r.z <= bbzmax) {
952  return true;
953  }
954  return false;
955  }
956 
957 
958  // isWhere
959  int isWhere(const EGS_Vector &r) {
960  if (!isInside(r)) {
961  return -1;
962  }
963  int ix, iy, iz;
964  setIndices(r,ix,iy,iz);
965  return getNode(ix,iy,iz)->region;
966  }
967 
968 
969  // inside (deprecated)
970  int inside(const EGS_Vector &r) {
971  return isWhere(r);
972  }
973 
974 
975  // isWhereFast
976  int isWhereFast(const EGS_Vector &r) {
977  if (!isInside(r)) {
978  return -1;
979  }
980  int ix, iy, iz;
981  setIndices(r, ix, iy, iz);
982  return getNode(ix,iy,iz)->region;
983  }
984 
985 
986  // medium
987  int medium(int ireg) const {
988  return nodeReg[ireg]->medium;
989  }
990 
991 
992  // howfarIn
993  int howfarIn(EGS_Octree_node *node, const EGS_Vector &r, const EGS_Vector &u, EGS_Float &t, EGS_Vector *normal=0) {
994 
995  int ix, iy, iz, tmp;
996  int crossed = -1;
997 
998  if (!node) {
999  return -1;
1000  }
1001 
1002  // set shift and local cell indices
1003  int shift = maxlevel - node->level; // how many levels missing between current level and full depth
1004  ix = node->ix;
1005  iy = node->iy;
1006  iz = node->iz;
1007 
1008  // x direction
1009  if (u.x > 0) {
1010  ix = ~(~ix<<shift); // fill lower level bits with 1's to always consider +x below current node
1011  EGS_Float xBound;
1012  if (ix>=ixmax) {
1013  xBound = bbxmax;
1014  }
1015  else {
1016  xBound = xmin + (ix+1)*dx;
1017  }
1018  EGS_Float d = (xBound-r.x) / u.x;
1019  if (d <= t) {
1020  t = d;
1021  crossed = 0;
1022  ix++;
1023  if (normal) {
1024  *normal = EGS_Vector(-1,0,0);
1025  }
1026  }
1027  }
1028  else if (u.x < 0) {
1029  ix = ix<<shift; // fill lower level bits with 0's to always consider -x below current node
1030  EGS_Float xBound;
1031  if (ix<=ixmin) {
1032  xBound = bbxmin;
1033  }
1034  else {
1035  xBound = xmin + ix*dx;
1036  }
1037  EGS_Float d = (xBound-r.x) / u.x;
1038  if (d <= t) {
1039  t = d;
1040  crossed = 0;
1041  ix--;
1042  if (normal) {
1043  *normal = EGS_Vector(1,0,0);
1044  }
1045  }
1046  }
1047 
1048 
1049  // y direction
1050  if (u.y > 0) {
1051  iy = ~(~iy<<shift); // fill lower level bits with 1's to always consider +y below current node
1052  EGS_Float yBound;
1053  if (iy>=iymax) {
1054  yBound = bbymax;
1055  }
1056  else {
1057  yBound = ymin + (iy+1)*dy;
1058  }
1059  EGS_Float d = (yBound - r.y) / u.y;
1060  if (d <= t) {
1061  t = d;
1062  crossed = 1;
1063  iy++;
1064  if (normal) {
1065  *normal = EGS_Vector(0,-1,0);
1066  }
1067  }
1068  }
1069  else if (u.y < 0) {
1070  iy = iy<<shift; // fill lower level bits with 0's to always consider -y below current node
1071  EGS_Float yBound;
1072  if (iy<=iymin) {
1073  yBound = bbymin;
1074  }
1075  else {
1076  yBound = ymin + iy*dy;
1077  }
1078  EGS_Float d = (yBound-r.y) / u.y;
1079  if (d <= t) {
1080  t = d;
1081  crossed = 1;
1082  iy--;
1083  if (normal) {
1084  *normal = EGS_Vector(0,1,0);
1085  }
1086  }
1087  }
1088 
1089  // z direction
1090  if (u.z > 0) {
1091  iz = ~(~iz<<shift); // fill lower level bits with 1's to always consider +z below current node
1092  EGS_Float zBound;
1093  if (iz>=izmax) {
1094  zBound = bbzmax;
1095  }
1096  else {
1097  zBound = zmin + (iz+1)*dz;
1098  }
1099  EGS_Float d = (zBound-r.z) / u.z;
1100  if (d <= t) {
1101  t = d;
1102  crossed = 2;
1103  iz++;
1104  if (normal) {
1105  *normal = EGS_Vector(0,0,-1);
1106  }
1107  }
1108  }
1109  else if (u.z < 0) {
1110  iz = iz<<shift; // fill lower level bits with 0's to always consider -z below current node
1111  EGS_Float zBound;
1112  if (iz<=izmin) {
1113  zBound = bbzmin;
1114  }
1115  else {
1116  zBound = zmin + iz*dz;
1117  }
1118  EGS_Float d = (zBound-r.z) / u.z;
1119  if (d <= t) {
1120  t = d;
1121  crossed = 2;
1122  iz--;
1123  if (normal) {
1124  *normal = EGS_Vector(0,0,1);
1125  }
1126  }
1127  }
1128 
1129  // get the new region index for the neighbor cell:
1130  // 1) find the new position in the plane perpendicular to the crossing direction
1131  // 2) get the indices for the neighbor cell at maximum depth, corresponding to that position
1132  // 3) call an axis specific function to get the neighbor node
1133  if (crossed==0) {
1134  EGS_Vector ryz(r.x, r.y+t*u.y, r.z+t*u.z);
1135  setIndices(ryz, tmp, iy, iz);
1136  node = getNeighborNodeX(node, ix, iy, iz);
1137  }
1138  else if (crossed==1) {
1139  EGS_Vector rxz(r.x+t*u.x, r.y, r.z+t*u.z);
1140  setIndices(rxz, ix, tmp, iz);
1141  node = getNeighborNodeY(node, ix, iy, iz);
1142  }
1143  else if (crossed==2) {
1144  EGS_Vector rxy(r.x+t*u.x, r.y+t*u.y, r.z);
1145  setIndices(rxy, ix, iy, tmp);
1146  node = getNeighborNodeZ(node, ix, iy, iz);
1147  }
1148 
1149  if (node) {
1150  return node->region;
1151  }
1152  else {
1153  return -1;
1154  }
1155  }
1156 
1157 
1158  // howfarOut
1159  int howfarOut(const EGS_Vector &r, const EGS_Vector &u, EGS_Float &t, EGS_Vector *normal=0) {
1160 
1161  int tmp;
1162  int ix=0, iy=0, iz=0;
1163  EGS_Float d, tlong = 2*t;
1164  EGS_Octree_node *node = NULL;
1165 
1166  // x axis
1167  if (r.x <= bbxmin && u.x > 0) {
1168  ix = ixmin;
1169  d = (bbxmin-r.x) / u.x;
1170  }
1171  else if (r.x >= bbxmax && u.x < 0) {
1172  ix = ixmax;
1173  d = (bbxmax-r.x) / u.x;
1174  }
1175  else {
1176  d = tlong;
1177  }
1178  if (d <= t) {
1179  EGS_Float yy = r.y + u.y*d;
1180  EGS_Float zz = r.z + u.z*d;
1181  if (yy >= bbymin && yy <= bbymax && zz >= bbzmin && zz <= bbzmax) {
1182  EGS_Vector rr(0,yy,zz);
1183  setIndices(rr, tmp, iy, iz);
1184  t = d;
1185  if (normal) {
1186  *normal = (ix == ixmin) ? EGS_Vector(-1,0,0) : EGS_Vector(1,0,0);
1187  }
1188  node = getNode(ix,iy,iz);
1189  return node->region;
1190  }
1191  }
1192 
1193  // y axis
1194  if (r.y <= bbymin && u.y > 0) {
1195  iy = iymin;
1196  d = (bbymin-r.y) / u.y;
1197  }
1198  else if (r.y >= bbymax && u.y < 0) {
1199  iy = iymax;
1200  d = (bbymax-r.y) / u.y;
1201  }
1202  else {
1203  d = tlong;
1204  }
1205  if (d <= t) {
1206  EGS_Float xx = r.x + u.x*d;
1207  EGS_Float zz = r.z + u.z*d;
1208  if (xx >= bbxmin && xx <= bbxmax && zz >= bbzmin && zz <= bbzmax) {
1209  EGS_Vector rr(xx,0,zz);
1210  setIndices(rr, ix, tmp, iz);
1211  t = d;
1212  if (normal) {
1213  *normal = (iy == iymin) ? EGS_Vector(0,-1,0) : EGS_Vector(0,1,0);
1214  }
1215  node = getNode(ix,iy,iz);
1216  return node->region;
1217  }
1218  }
1219 
1220  // z axis
1221  if (r.z <= bbzmin && u.z > 0) {
1222  iz = izmin;
1223  d = (bbzmin-r.z) / u.z;
1224  }
1225  else if (r.z >= bbzmax && u.z < 0) {
1226  iz = izmax;
1227  d = (bbzmax-r.z) / u.z;
1228  }
1229  else {
1230  d = tlong;
1231  }
1232  if (d <= t) {
1233  EGS_Float xx = r.x + u.x*d;
1234  EGS_Float yy = r.y + u.y*d;
1235  if (xx >= bbxmin && xx <= bbxmax && yy >= bbymin && yy <= bbymax) {
1236  EGS_Vector rr(xx,yy,0);
1237  setIndices(rr, ix, iy, tmp);
1238  t = d;
1239  if (normal) {
1240  *normal = (iz == izmin) ? EGS_Vector(0,0,-1) : EGS_Vector(0,0,1);
1241  }
1242  node = getNode(ix,iy,iz);
1243  return node->region;
1244  }
1245  }
1246 
1247  return -1;
1248  }
1249 
1250 
1251  // howfar
1252  int howfar(int ireg, const EGS_Vector &r, const EGS_Vector &u, EGS_Float &t, int *newmed, EGS_Vector *normal=0) {
1253 
1254  int inew = ireg;
1255 
1256  // get new region number
1257  if (ireg==-1) {
1258  inew = howfarOut(r, u, t, normal);
1259  }
1260  else {
1261  inew = howfarIn(nodeReg[ireg], r, u, t, normal);
1262  }
1263 
1264  // set new medium
1265  if (inew>=0 && newmed) {
1266  *newmed = nodeReg[inew]->medium;
1267  }
1268  return inew;
1269  }
1270 
1271 
1272  // hownearIn
1273  EGS_Float hownearIn(int ireg, const EGS_Vector &r) {
1274  EGS_Octree_node *node = nodeReg[ireg];
1275  int shift = maxlevel - node->level;
1276  EGS_Float t1, t2, tx, ty, tz;
1277  int imin, imax;
1278 
1279  // x
1280  imin = node->ix << shift;
1281  imax = ~(~node->ix << shift);
1282  t1 = (r.x-xmin)-dx*imin;
1283  t2 = dx*(imax-imin+1)-t1;
1284  tx = t1 < t2 ? t1 : t2;
1285 
1286  // y
1287  imin = node->iy << shift;
1288  imax = ~(~node->iy << shift);
1289  t1 = (r.y-ymin)-dy*imin;
1290  t2 = dy*(imax-imin+1)-t1;
1291  ty = t1 < t2 ? t1 : t2;
1292 
1293  // z
1294  imin = node->iz << shift;
1295  imax = ~(~node->iz << shift);
1296  t1 = (r.z-zmin)-dz*imin;
1297  t2 = dz*(imax-imin+1)-t1;
1298  tz = t1 < t2 ? t1 : t2;
1299 
1300  return tx<ty && tx<tz ? tx : ty<tz ? ty : tz;
1301  }
1302 
1303 
1304  // hownear
1305  EGS_Float hownear(int ireg, const EGS_Vector &r) {
1306  if (ireg>=0) {
1307  return hownearIn(ireg, r);
1308  }
1309  int nc=0;
1310  EGS_Float s1=0, s2=0;
1311  if (r.x < bbxmin || r.x > bbxmax) {
1312  EGS_Float t = r.x < bbxmin ? bbxmin-r.x : r.x-bbxmax;
1313  nc++;
1314  s1 += t;
1315  s2 += t*t;
1316  }
1317  if (r.y < bbymin || r.y > bbymax) {
1318  EGS_Float t = r.y < bbymin ? bbymin-r.y : r.y-bbymax;
1319  nc++;
1320  s1 += t;
1321  s2 += t*t;
1322  }
1323  if (r.z < bbzmin || r.z > bbzmax) {
1324  EGS_Float t = r.z < bbzmin ? bbzmin-r.z : r.z-bbzmax;
1325  nc++;
1326  s1 += t;
1327  s2 += t*t;
1328  }
1329  return nc == 1 ? s1 : sqrt(s2);
1330  }
1331 
1332 
1333  // getType
1334  const string &getType() const {
1335  return type;
1336  }
1337 
1338  // printInfo
1339  void printInfo() const;
1340 };
1341 
1342 #endif
int iz
octree indices, representing the binary location code of the node in x, y, z;
Definition: egs_octree.h:334
virtual const string & getType() const =0
Get the geometry type.
virtual bool isInside(const EGS_Vector &x)=0
Is the position x inside the geometry?
virtual int medium(int ireg) const
Returns the medium index in region ireg.
virtual void printInfo() const
Print information about this geometry.
EGS_BaseGeometry(const string &Name)
Construct a geometry named Name.
An octree geometry.
Definition: egs_octree.h:518
EGS_Octree_node * child
pointer to children nodes, NULL is there are no children
Definition: egs_octree.h:338
EGS_Float y
y-component
Definition: egs_vector.h:61
A class representing 3D vectors.
Definition: egs_vector.h:56
Global egspp functions header file.
virtual int isWhere(const EGS_Vector &x)=0
In which region is poisition x?
virtual int howfar(int ireg, const EGS_Vector &x, const EGS_Vector &u, EGS_Float &t, int *newmed=0, EGS_Vector *normal=0)=0
Calculate the distance to a boundary from x along the direction u.
short medium
medium index for the node
Definition: egs_octree.h:336
Base geometry class. Every geometry class must be derived from EGS_BaseGeometry.
EGS_Float z
z-component
Definition: egs_vector.h:62
EGS_InfoFunction EGS_EXPORT egsFatal
Always use this function for reporting fatal errors.
virtual int inside(const EGS_Vector &x)=0
Returns the region index, if inside, or -1 if outside (obsolete)
EGS_Octree_node * parent
pointer to the parent node (only root node can have parent set to NULL)
Definition: egs_octree.h:339
virtual EGS_Float hownear(int ireg, const EGS_Vector &x)=0
Calculate the distance to a boundary for position x in any direction.
unsigned short level
depth of the node (root node is level 0)
Definition: egs_octree.h:337
EGS_Float x
x-component
Definition: egs_vector.h:60
int region
region number of the node
Definition: egs_octree.h:335
EGS_BaseGeometry class header file.
int nreg
Number of local regions in this geometry.