EGSnrc C++ class library  Report PIRS-898 (2021)
Iwan Kawrakow, Ernesto Mainegra-Hing, Frederic Tessier, Reid Townson and Blake Walters
egs_octree.cpp
Go to the documentation of this file.
1 /*
2 ###############################################################################
3 #
4 # EGSnrc egs++ octree geometry
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 #include "egs_octree.h"
172 #include "egs_input.h"
173 
174 void EGS_Octree::printInfo() const {
176  egsInformation(" bounding box minimum = %g %g %g\n", bbxmin, bbymin, bbzmin);
177  egsInformation(" bounding box maximum = %g %g %g\n", bbxmax, bbymax, bbzmax);
178  egsInformation(" bounding box resolution = %d %d %d\n", nx, ny, nz);
179  egsInformation(" octree leaf size = %g %g %g\n", dx, dy, dz);
180  egsInformation(" octree cells (no medium) = %d\n", nreg-nLeaf);
181  egsInformation(" octree cells (medium) = %d\n", nLeaf);
182  egsInformation(" octree average cell size = %.2f\n", nLeafMax/(float)nLeaf);
183  char percent = '%';
184  egsInformation(" octree cell savings = %.1f%c\n", 100*(1-(float)nLeaf/nLeafMax),percent);
185  egsInformation("=======================================================\n");
186 }
187 
188 string EGS_Octree::type("EGS_Octree");
189 
190 static char EGS_OCTREE_LOCAL eoctree_message1[] = "createGeometry(octree): %s\n";
191 static char EGS_OCTREE_LOCAL eoctree_message2[] = "null input?";
192 //static char EGS_OCTREE_LOCAL eoctree_message3[] = "wrong/missing 'octree size' input?";
193 //static char EGS_OCTREE_LOCAL eoctree_message4[] = "expecting 1 or 3 float inputs for 'octree size'";
194 static char EGS_OCTREE_LOCAL eoctree_message5[] = "wrong/missing 'min' or input?";
195 static char EGS_OCTREE_LOCAL eoctree_message6[] = "expecting 3 float inputs for 'box min' input";
196 static char EGS_OCTREE_LOCAL eoctree_message7[] = "wrong/missing 'max' or input?";
197 static char EGS_OCTREE_LOCAL eoctree_message8[] = "expecting 3 float inputs for 'max max' input";
198 static char EGS_OCTREE_LOCAL eoctree_message9[] = "wrong or missing 'resolution' input?";
199 static char EGS_OCTREE_LOCAL eoctree_message10[] = "expecting 3 integer inputs for 'resolution' input";
200 static char EGS_OCTREE_LOCAL eoctree_message11[] = "wrong or missing 'discard child' input?";
201 static char EGS_OCTREE_LOCAL eoctree_message12[] = "expecting 'yes' or 'no' for 'discard child' input";
202 static char EGS_OCTREE_LOCAL eoctree_message13[] = "missing or wrong 'child geometry' input";
203 static char EGS_OCTREE_LOCAL eoctree_message14[] = "undefined child geometry";
204 static char EGS_OCTREE_LOCAL eoctree_message15[] = "you need to define at least one octree box";
205 static char EGS_OCTREE_LOCAL eoctree_message16[] = "wrong 'prune tree' input?";
206 static char EGS_OCTREE_LOCAL eoctree_message17[] = "expecting 'yes' or 'no' for 'prune tree' input?";
207 static char EGS_OCTREE_LOCAL eoctree_key0[] = "octree box";
208 static char EGS_OCTREE_LOCAL eoctree_key1[] = "box min";
209 static char EGS_OCTREE_LOCAL eoctree_key2[] = "box max";
210 static char EGS_OCTREE_LOCAL eoctree_key3[] = "resolution";
211 static char EGS_OCTREE_LOCAL eoctree_key4[] = "child geometry";
212 static char EGS_OCTREE_LOCAL eoctree_key5[] = "discard child";
213 static char EGS_OCTREE_LOCAL eoctree_key6[] = "prune tree";
214 
215 extern "C" {
216 
217  EGS_OCTREE_EXPORT EGS_BaseGeometry *createGeometry(EGS_Input *input) {
218 
219  EGS_Input *i;
220 
221  // check that we have an input
222  if (!input) {
223  egsWarning(eoctree_message1,eoctree_message2);
224  return 0;
225  }
226 
227  // read bounding boxes
228  vector<EGS_Octree_bbox> vBox;
229  while ((i = input->takeInputItem(eoctree_key0))) {
230 
231  // read the bounding box minimum
232  vector<EGS_Float> v;
233  int err = i->getInput(eoctree_key1, v);
234  if (err) {
235  egsWarning(eoctree_message1, eoctree_message5);
236  return 0;
237  }
238  if (v.size() != 3) {
239  egsWarning(eoctree_message1, eoctree_message6);
240  return 0;
241  }
242  EGS_Vector bboxMin(v[0],v[1],v[2]);
243 
244  // read the bounding box maximum
245  err = i->getInput(eoctree_key2, v);
246  if (err) {
247  egsWarning(eoctree_message1, eoctree_message7);
248  return 0;
249  }
250  if (v.size() != 3) {
251  egsWarning(eoctree_message1, eoctree_message8);
252  return 0;
253  }
254  EGS_Vector bboxMax(v[0],v[1],v[2]);
255 
256  // read the bounding box resolution
257  vector<int> bboxRes;
258  err = i->getInput(eoctree_key3, bboxRes);
259  if (err) {
260  egsWarning(eoctree_message1, eoctree_message9);
261  return 0;
262  }
263  if (bboxRes.size() != 3) {
264  egsWarning(eoctree_message1, eoctree_message10);
265  return 0;
266  }
267 
268  EGS_Octree_bbox box = EGS_Octree_bbox(bboxMin, bboxMax, bboxRes);
269  vBox.push_back(box);
270  }
271  if (vBox.size() < 1) {
272  egsWarning(eoctree_message1, eoctree_message15);
273  return 0;
274  }
275 
276  // read discard child option
277  bool discardChild = true;
278  string discard;
279  if (input->getInputItem(eoctree_key5)) {
280  int err = input->getInput(eoctree_key5, discard);
281  if (err) {
282  egsWarning(eoctree_message1, eoctree_message11);
283  return 0;
284  }
285  if (discard.find("yes")==string::npos && discard.find("no")==string::npos) {
286  egsWarning(eoctree_message1, eoctree_message12);
287  return 0;
288  }
289  if (discard.find("no")!=string::npos) {
290  discardChild = false;
291  }
292  }
293 
294  // read prune tree option
295  bool pruneTree = true;
296  string prune;
297  if (input->getInputItem(eoctree_key6)) {
298  int err = input->getInput(eoctree_key6, prune);
299  if (err) {
300  egsWarning(eoctree_message1, eoctree_message16);
301  return 0;
302  }
303  if (prune.find("yes")==string::npos && prune.find("no")==string::npos) {
304  egsWarning(eoctree_message1, eoctree_message17);
305  return 0;
306  }
307  if (prune.find("no")!=string::npos) {
308  pruneTree = false;
309  }
310  }
311 
312  // read and load the child geometry
313  string gname;
314  {
315  int err = input->getInput(eoctree_key4, gname);
316  if (err) {
317  egsWarning(eoctree_message1, eoctree_message13);
318  return 0;
319  }
320  }
322  if (!g) {
323  egsWarning(eoctree_message1, eoctree_message14);
324  return 0;
325  }
326 
327  // create the octree geometry
328  EGS_Octree *octree = new EGS_Octree(vBox, pruneTree, g);
329  octree->setName(input);
330  octree->setBoundaryTolerance(input);
331  octree->setLabels(input);
332  octree->printInfo();
333 
334  if (discardChild) {
335  delete g;
336  }
337  return octree;
338  }
339 }
An octree geometry: header.
virtual void printInfo() const
Print information about this geometry.
An octree geometry.
Definition: egs_octree.h:518
EGS_Input class header file.
A class representing 3D vectors.
Definition: egs_vector.h:56
int setLabels(EGS_Input *input)
Set the labels from an input block.
EGS_GLIB_EXPORT EGS_BaseGeometry * createGeometry(EGS_Input *input)
Definition: egs_glib.cpp:57
Base geometry class. Every geometry class must be derived from EGS_BaseGeometry.
void setBoundaryTolerance(EGS_Input *inp)
Set the value of the boundary tolerance from the input inp.
static EGS_BaseGeometry * getGeometry(const string &Name)
Get a pointer to the geometry named Name.
EGS_Input * getInputItem(const string &key) const
Same as the previous function but now ownership remains with the EGS_Input object.
Definition: egs_input.cpp:245
EGS_InfoFunction EGS_EXPORT egsInformation
Always use this function for reporting the progress of a simulation and any other type of information...
void setName(EGS_Input *inp)
Set the name of the geometry from the input inp.
A class for storing information in a tree-like structure of key-value pairs. This class is used throu...
Definition: egs_input.h:182
int nreg
Number of local regions in this geometry.
EGS_Input * takeInputItem(const string &key, bool self=true)
Get the property named key.
Definition: egs_input.cpp:226
int getInput(const string &key, vector< string > &values) const
Assign values to an array of strings from an input identified by key.
Definition: egs_input.cpp:338
EGS_InfoFunction EGS_EXPORT egsWarning
Always use this function for reporting warnings.