EGSnrc C++ class library  Report PIRS-898 (2021)
Iwan Kawrakow, Ernesto Mainegra-Hing, Frederic Tessier, Reid Townson and Blake Walters
egs_vhp_geometry.cpp
Go to the documentation of this file.
1 /*
2 ###############################################################################
3 #
4 # EGSnrc egs++ voxelized human phantom 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: Iwan Kawrakow, 2008
25 #
26 # Contributors: Frederic Tessier
27 # Hubert Ho
28 #
29 ###############################################################################
30 #
31 # A "Voxelized Human Phantom" (VHP) geometry.
32 #
33 # Although a VHP geometry is just a XYZ geometry, due to the huge number of
34 # voxels one has to be more careful with memory use for media indices, etc.
35 # That's why this separate geometry implementation. With Richard Kramer's
36 # phantom data the implementation requires ~12 MB to hold all the organ data
37 # in memory. Richard's ASCII phantom data was also transformed into binary
38 # format with the result that loading them is almost instantaneous.
39 #
40 # For now the implementation only provides the macro geometry, i.e., no
41 # spongiosa micro matrices. Adding the micro matrices is the next step.
42 #
43 ###############################################################################
44 */
45 
46 
52 #include "egs_vhp_geometry.h"
53 #include "egs_input.h"
54 #include "egs_functions.h"
55 
56 #include <vector>
57 #include <fstream>
58 #ifndef NO_SSTREAM
59  #include <sstream>
60  #define S_STREAM std::istringstream
61 #else
62  #include <strstream>
63  #define S_STREAM std::istrstream
64 #endif
65 
66 
67 using namespace std;
68 
69 #ifndef SKIP_DOXYGEN
70 EGS_VoxelInfo *EGS_VoxelGeometry::v = 0;
71 int EGS_VoxelGeometry::nv = 0;
72 
73 VHP_OrganData::VHP_OrganData(const char *fname, int slice_min, int slice_max)
74  : nslice(0), ok(false) {
75  ifstream data(fname,ios::binary);
76  if (!data) {
77  egsWarning("VHP_OrganData: failed to open file %s for reading\n",
78  fname);
79  return;
80  }
81  char endian;
82  data.read(&endian,sizeof(endian));
83  char my_endian = egsGetEndian();
84  if (endian != my_endian) {
85  egsWarning("VHP_OrganData: wrong endianess. \n"
86  "This machine has endianess %d, data was written on %d\n"
87  "This requires byte swaping, which has not been implemented yet\n",
88  (int)my_endian,(int)endian);
89  return;
90  }
91  data.read((char *)&dx, sizeof(dx));
92  data.read((char *)&dy, sizeof(dy));
93  data.read((char *)&dz, sizeof(dz));
94  data.read((char *)&nslice,sizeof(nslice));
95  if (slice_min < 0) {
96  slice_min = 0;
97  }
98  if (slice_max > nslice) {
99  slice_max = nslice;
100  }
101  int nwant = slice_max - slice_min;
102  slices = new VHP_SliceInfo [nwant];
103  nx = 0, ny = 0;
104  int is = 0;
105  for (int islice=0; islice<nslice; ++islice) {
106  if (islice >= slice_min) {
107  slices[is].readData(data);
108  int ni, nj;
109  slices[is++].getSliceDimensions(ni,nj);
110  if (ni > nx) {
111  nx = ni;
112  }
113  if (nj > ny) {
114  ny = nj;
115  }
116  }
117  else {
118  VHP_SliceInfo tmp;
119  tmp.readData(data);
120  }
121  if (islice == slice_max - 1) {
122  break;
123  }
124  }
125  nxy = nx*ny;
126  nslice = slice_max - slice_min;
127  if (data.fail()) {
128  egsWarning("VHP_OrganData: I/O error while reading data\n");
129  }
130  else {
131  ok = true;
132  }
133 }
134 
135 VHP_OrganData::~VHP_OrganData() {
136  if (nslice) {
137  delete [] slices;
138  }
139 }
140 
141 void VHP_RowInfo::readData(istream &in) {
142  unsigned short n;
143  in.read((char *)&n,sizeof(n));
144  if (n != nbin) {
145  if (nbin > 0) {
146  delete [] pix;
147  delete [] org;
148  nbin = 0;
149  }
150  }
151  if (n < 1) {
152  return;
153  }
154  if (!nbin) {
155  nbin = n;
156  pix = new unsigned short [nbin+1];
157  org = new unsigned char [nbin];
158  }
159  in.read((char *)pix,(nbin+1)*sizeof(unsigned short));
160  in.read((char *)org,nbin*sizeof(unsigned char));
161 }
162 
163 void VHP_SliceInfo::readData(istream &in) {
164  if (rinfo) {
165  delete [] rinfo;
166  rinfo = 0;
167  }
168  in.read((char *)&firstr,sizeof(firstr));
169  in.read((char *)&lastr,sizeof(lastr));
170  int nrow = lastr-firstr+1;
171  if (nrow > 0) {
172  rinfo = new VHP_RowInfo [nrow];
173  for (int j=0; j<nrow; ++j) {
174  rinfo[j].readData(in);
175  }
176  }
177 }
178 #endif
179 
180 EGS_VHPGeometry::EGS_VHPGeometry(const char *phantom_file,
181  const char *media_file, int slice_min, int slice_max,
182  const string &Name) :
183  EGS_BaseGeometry(Name), vg(0), micros(0), nmicro(0) {
184  organs = new VHP_OrganData(phantom_file,slice_min,slice_max);
185  if (!organs->isOK()) {
186  egsWarning("EGS_VHPGeometry: failed to construct organ data\n");
187  delete organs;
188  organs = 0;
189  return;
190  }
191  ifstream med_data(media_file);
192  if (!med_data) {
193  egsWarning("EGS_VHPGeometry: failed to open media data file %s\n",
194  media_file);
195  delete organs;
196  organs = 0;
197  return;
198  }
199  int norg;
200  med_data >> norg;
201  if (norg < 1) {
202  egsWarning("EGS_VHPGeometry: %d organs in media data file %s?\n",
203  norg,media_file);
204  delete organs;
205  organs = 0;
206  return;
207  }
208  int j;
209  for (j=0; j<256; ++j) {
210  organ_media[j] = -1;
211  organ_names[j] = "Undefined";
212  organ_micro[j] = -1;
213  }
214  char buf[1024];
215  med_data.getline(buf,1023);
216  for (j=0; j<norg; ++j) {
217  int iorg, imed;
218  char c;
219  med_data >> iorg >> imed >> c;
220  med_data.getline(buf,1023);
221  if (iorg < 0 || iorg > 255)
222  egsWarning("EGS_VHPGeometry: bad organ number %d on line %d\n",
223  iorg,j+2);
224  else {
225  organ_media[iorg] = imed;
226  }
227  organ_names[iorg] = c;
228  organ_names[iorg] += buf;
229  }
230  int nmed;
231  med_data >> nmed;
232  med_data.getline(buf,1023);
233  vector<string> media_names;
234  vector<int> media_indeces;
235  for (j=0; j<nmed; ++j) {
236  med_data.getline(buf,1023);
237  string med(buf);
238  media_names.push_back(med);
239  media_indeces.push_back(EGS_BaseGeometry::addMedium(med));
240  }
241  egsInformation("EGS_VHPGeometry: media information\n");
242  for (j=0; j<nmed; ++j) {
243  egsInformation("Medium %s registered as medium number %d\n",
244  media_names[j].c_str(),media_indeces[j]);
245  }
246  for (j=0; j<256; ++j) {
247  if (organ_media[j] >= 0) {
248  if (organ_media[j] >= nmed) {
249  egsWarning("EGS_VHPGeometry: invalid medium index %d\n",
250  organ_media[j]);
251  }
252  else {
253  organ_media[j] = media_indeces[organ_media[j]];
254  }
255  }
256  }
257  EGS_Float dx, dy, dz;
258  int nx, ny, nz;
259  organs->getVoxelSizes(dx,dy,dz);
260  organs->getPhantomDimensions(nx,ny,nz);
261  vg = new EGS_VoxelGeometry(nx,ny,nz,dx,dy,dz);
262  nreg = nx*ny*nz;
263  nmacro = nreg;
264  egsInformation("Phantom size in voxels: x=%d y=%d z=%d\n",nx,ny,nz);
265  egsInformation("Voxel sizes: dx=%g dy=%g dz=%g\n",dx,dy,dz);
266 }
267 
268 void EGS_VHPGeometry::setMicros(EGS_Input *input) {
269  if (!vg || !organs) {
270  egsWarning("EGS_VHPGeometry::setMicros: called for an invalid "
271  "VHP geometry\n");
272  return;
273  }
274  if (!input->getInputItem("micro matrix")) {
275  egsWarning("EGS_VHPGeometry::setMicros: no micro matrix definitions\n");
276  return;
277  }
278  vector<EGS_MicroMatrixCluster *> mv;
279  string dir;
280  input->getInput("micro matrix folder",dir);
281  EGS_Float t_bsc;
282  bool ok = true;
283  const static char *err_msg =
284  "EGS_VHPGeometry::setMicros: wrong/missing '%s' input\n";
285  int err = input->getInput("BSC thickness",t_bsc);
286  if (err || t_bsc < 0) {
287  egsWarning(err_msg,"BSC thickness");
288  ok = false;
289  }
290  string tb_medium, bm_medium;
291  err = input->getInput("TB medium",tb_medium);
292  if (err) {
293  egsWarning(err_msg,"TB medium");
294  ok = false;
295  }
296  err = input->getInput("BM medium",bm_medium);
297  if (err) {
298  egsWarning(err_msg,"BM medium");
299  ok = false;
300  }
301  if (!ok) {
302  egsWarning("EGS_VHPGeometry::setMicros: wrong/missing inputs found\n"
303  " => no micro matrices set for use\n");
304  return;
305  }
306  int tb_med = addMedium(tb_medium), bm_med = addMedium(bm_medium);
307  vector<string> minput;
308  while (!input->getInput("micro matrix",minput)) {
309  delete input->takeInputItem("micro matrix");
310  if (minput.size() < 2) {
311  egsWarning("EGS_VHPGeometry::setMicros: 'micro matrix' input"
312  " requires at least two inputs => ignored\n");
313  continue;
314  }
315  string fname = dir.size() ? egsJoinPath(dir,minput[0]) : minput[0];
316  EGS_MicroMatrixCluster *mcluster = new EGS_MicroMatrixCluster(
317  vg->dx,vg->dy,vg->dz,t_bsc,tb_med,bm_med,fname.c_str());
318  if (!mcluster->isValid()) {
319  egsWarning("failed to construct micro cluster from %s\n",
320  fname.c_str());
321  delete mcluster;
322  }
323  else {
324  bool ok = true;
325  vector<int> orgs;
326  for (int i=1; i<minput.size(); ++i) {
327  S_STREAM s(minput[i].c_str());
328  int org;
329  s >> org;
330  if (s.fail()) {
331  egsWarning("EGS_VHPGeometry::setMicros: invalid input %s "
332  "for micro matrix %s\n",minput[i].c_str(),
333  minput[0].c_str());
334  ok = false;
335  }
336  else {
337  if (org < 0 || org > 255) {
338  egsWarning("EGS_VHPGeometry::setMicros: invalid input"
339  " %d for micro matrix %s\n",org,minput[0].c_str());
340  ok = false;
341  }
342  else {
343  orgs.push_back(org);
344  }
345  }
346  }
347  if (ok) {
348  mv.push_back(mcluster);
349  egsInformation("Using micro matrix data from\n %s\nfor"
350  " organs:",fname.c_str());
351  for (int i=0; i<orgs.size(); ++i) {
352  egsInformation(" %d",orgs[i]);
353  organ_micro[orgs[i]] = mv.size()-1;
354  }
355  egsInformation("\n");
356  }
357  else {
358  delete mcluster;
359  }
360  }
361  }
362  if (mv.size() > 0) {
363  if (nmicro > 0) {
364  for (int i=0; i<nmicro; ++i) {
365  delete micros[i];
366  }
367  delete [] micros;
368  }
369  nmicro = mv.size();
370  micros = new EGS_MicroMatrixCluster* [nmicro];
371  nmax = 0;
372  for (int i=0; i<nmicro; ++i) {
373  micros[i] = mv[i];
374  int nr = micros[i]->nxy*micros[i]->nz;
375  int nm = micros[i]->mxy*micros[i]->mz;
376  if (nm*nr > nmax) {
377  nmax = nm*nr;
378  }
379  }
380  nmacro = vg->nxy*vg->nz;
381  EGS_I64 ntot = nmax;
382  ntot *= nmicro;
383  ntot += nmacro;
384  if (ntot > 2147483647) egsFatal("EGS_VHPGeometry::setMicros: "
385  "too many micro regions\n");
386  nreg = nmacro + nmax*nmicro;
387  egsInformation("Using regions 0...%d to identify macro voxels\n",
388  nmacro-1);
389  egsInformation("Using regions %d...%d to identify micro voxels\n",
390  nmacro,nreg-1);
391  }
392 }
393 
394 EGS_VHPGeometry::~EGS_VHPGeometry() {
395  if (organs) {
396  delete organs;
397  }
398  if (vg) {
399  delete vg;
400  }
401  if (nmicro > 0) {
402  for (int j=0; j<nmicro; ++j) {
403  delete micros[j];
404  }
405  delete [] micros;
406  }
407 
408 }
409 
410 void EGS_VHPGeometry::setMedia(EGS_Input *,int,const int *) {
411  egsWarning("EGS_VHPGeometry::setMedia: don't use this method.\n"
412  " Media are fixed by the organ data defining an VHP geometry\n");
413 }
414 
415 void EGS_VHPGeometry::printInfo() const {
417  if (isOK()) {
418  egsInformation(" voxel sizes: %g %g %g\n",vg->dx,vg->dy,vg->dz);
419  egsInformation(" phantom size: %d %d %d\n",vg->nx,vg->ny,vg->nz);
420  }
421  else {
422  egsInformation(" undefined VHP geometry\n");
423  }
424 }
425 
426 string EGS_VHPGeometry::type = "EGS_VHPGeometry";
427 
428 #ifndef SKIP_DOXYGEN
429 VHPBox EGS_MicroMatrixCluster::bm_boxes[64];
430 bool EGS_MicroMatrixCluster::bm_boxes_initialized = false;
431 
432 EGS_MicroMatrixCluster::EGS_MicroMatrixCluster(
433  EGS_Float Dx, EGS_Float Dy, EGS_Float Dz, EGS_Float bsc_thickness,
434  int tb_med, int bm_med, const char *micro_file) : med_tb(tb_med), med_bm(bm_med),
435  micros(0), vg(0) {
436  //
437  // *** open micro-matrix file
438  //
439  ifstream in(micro_file,ios::binary);
440  if (!in) {
441  egsWarning("EGS_MicroMatrixCluster: failed to open file %s\n",
442  micro_file);
443  return;
444  }
445 
446  //
447  // *** read data
448  //
449  unsigned char Mx, My, Mz;
450  in.read((char *)&Mx,sizeof(Mx));
451  in.read((char *)&My,sizeof(My));
452  in.read((char *)&Mz,sizeof(Mz));
453  mx = Mx;
454  my = My;
455  mz = Mz;
456  mxy = mx*my;
457  nmic = mxy*mz;
458  micros = new unsigned char *[nmic];
459  in.read((char *)&Mx,sizeof(Mx));
460  in.read((char *)&My,sizeof(My));
461  in.read((char *)&Mz,sizeof(Mz));
462  nx = Mx;
463  ny = My;
464  nz = Mz;
465  nxy = nx*ny;
466  int nxyz = nxy*nz, j;
467  nreg = nxyz;
468  {
469  for (int j=0; j<nmic; ++j) {
470  micros[j] = new unsigned char [nxyz];
471  in.read((char *)micros[j],nxyz*sizeof(unsigned char));
472  if (in.fail()) {
473  egsWarning("EGS_MicroMatrixCluster: I/O error while reading "
474  "%d'th micto-matrix from file %s\n",j+1,micro_file);
475  for (int i=0; i<=j; ++i) {
476  delete [] micros[j];
477  }
478  delete [] micros;
479  micros = 0;
480  return;
481  }
482  }
483  }
484 
485  //
486  // *** initialize micro- and macro-voxel sizes and numbers
487  //
488  dx = Dx;
489  dy = Dy;
490  dz = Dz;
491  dxi = 1/dx;
492  dyi = 1/dy;
493  dzi = 1/dz;
494  EGS_Float ddx = Dx/nx, ddy = Dy/ny, ddz = Dz/nz;
495  egsInformation("Macro voxel sizes: %g %g %g\n",dx,dy,dz);
496  egsInformation("Micro voxel sizes: %g %g %g\n",ddx,ddy,ddz);
497  egsInformation("Number of micros: %d %d %d\n",mx,my,mz);
498  egsInformation("Number of micro-voxles: %d %d %d\n",nx,ny,nz);
499 
500  //
501  // *** create an object for geometry computations
502  //
503  vg = new EGS_VoxelGeometry(nx,ny,nz,ddx,ddy,ddz);
504 
505  //
506  // *** initialize the 64 boxes needed for sub-micro-voxel energy deposition
507  // if needed
508  //
509  if (!bm_boxes_initialized) {
510  bm_boxes_initialized = true;
511  for (int ixm=0; ixm<2; ++ixm) {
512  EGS_Float xmin = bsc_thickness*ixm;
513  for (int ixp=0; ixp<2; ++ixp) {
514  EGS_Float xmax = ddx - bsc_thickness*ixp;
515  for (int iym=0; iym<2; ++iym) {
516  EGS_Float ymin = bsc_thickness*iym;
517  for (int iyp=0; iyp<2; ++iyp) {
518  EGS_Float ymax = ddy - bsc_thickness*iyp;
519  for (int izm=0; izm<2; ++izm) {
520  EGS_Float zmin = bsc_thickness*izm;
521  for (int izp=0; izp<2; ++izp) {
522  EGS_Float zmax = ddz - bsc_thickness*izp;
523  int ibox = ixm+2*ixp+4*iym+8*iyp+16*izm+32*izp;
524  bm_boxes[ibox].xmin = xmin;
525  bm_boxes[ibox].xmax = xmax;
526  bm_boxes[ibox].ymin = ymin;
527  bm_boxes[ibox].ymax = ymax;
528  bm_boxes[ibox].zmin = zmin;
529  bm_boxes[ibox].zmax = zmax;
530  }
531  }
532  }
533  }
534  }
535  }
536  }
537 
538  //
539  // *** process micro data. After this loop the mico-matrices will contain
540  // 0 for TB voxels
541  // 1 for BSC-only voxels (possible if BSC thickness > voxel size/2)
542  // 2 for BM-only voxels (i.e., no TB neighbours)
543  // >2 index of VHPBox corresponding to the BSC layers found in the voxel
544  //
545  // we also compute the TB, BM and BSC volumes
546  //
547  v_tb = new EGS_Float [nmic];
548  v_bm = new EGS_Float [nmic];
549  v_bsc = new EGS_Float [nmic];
550  double tot_tb = 0, tot_bm = 0, tot_bsc = 0;
551  for (j=0; j<nmic; ++j) {
552  int iz = j/mxy;
553  int iy = (j-iz*mxy)/mx;
554  int ix = j-iz*mxy-iy*mx;
555  int ixm = ix-1;
556  if (ixm < 0) {
557  ixm = mx-1;
558  }
559  int ixp = ix+1;
560  if (ixp >= mx) {
561  ixp = 0;
562  }
563  int iym = iy-1;
564  if (iym < 0) {
565  iym = my-1;
566  }
567  int iyp = iy+1;
568  if (iyp >= my) {
569  iyp = 0;
570  }
571  int izm = iz-1;
572  if (izm < 0) {
573  izm = mz-1;
574  }
575  int izp = iz+1;
576  if (izp >= mz) {
577  izp = 0;
578  }
579  EGS_Float bsc_vol = 0, bm_vol = 0, tb_vol = 0;
580  for (int jx=0; jx<nx; ++jx) {
581  for (int jy=0; jy<ny; ++jy) {
582  for (int jz=0; jz<nz; ++jz) {
583  int ireg = jx + jy*nx + jz*nxy;
584  if (micros[j][ireg]) {
585  int micxm = !(jx > 0 ? micros[j][ireg-1] :
586  micros[ixm+iy*mx+iz*mxy][nx-1+jy*nx+jz*nxy]);
587  int micxp = !(jx < nx-1 ? micros[j][ireg+1] :
588  micros[ixp+iy*mx+iz*mxy][jy*nx+jz*nxy]);
589  int micym = !(jy > 0 ? micros[j][ireg-nx] :
590  micros[ix+iym*mx+iz*mxy][jx+(ny-1)*nx+jz*nxy]);
591  int micyp = !(jy < ny-1 ? micros[j][ireg+nx] :
592  micros[ix+iyp*mx+iz*mxy][jx+jz*nxy]);
593  int miczm = !(jz > 0 ? micros[j][ireg-nxy] :
594  micros[ix+iy*mx+izm*mxy][jx+jy*nx+(nz-1)*nxy]);
595  int miczp = !(jz < nz-1 ? micros[j][ireg+nxy] :
596  micros[ix+iy*mx+izp*mxy][jx+jy*nx]);
597  if ((micxm && micxp && 2*bsc_thickness>ddx) ||
598  (micym && micyp && 2*bsc_thickness>ddy) ||
599  (miczm && miczp && 2*bsc_thickness>ddz)) {
600  micros[j][ireg] = 1; // i.e. all BSC
601  bsc_vol += 1;
602  }
603  else {
604  int ibox=micxm+2*micxp+4*micym+8*micyp+
605  16*miczm+32*miczp;
606  micros[j][ireg] = ibox + 2;
607  if (!ibox) {
608  bm_vol += 1;
609  }
610  else {
611  double v =
612  (bm_boxes[ibox].xmax-bm_boxes[ibox].xmin)*
613  (bm_boxes[ibox].ymax-bm_boxes[ibox].ymin)*
614  (bm_boxes[ibox].zmax-bm_boxes[ibox].zmin);
615  v /= (ddx*ddy*ddz);
616  bm_vol += v;
617  bsc_vol += 1 - v;
618  }
619  }
620  }
621  else {
622  tb_vol += 1;
623  }
624  }
625  }
626  }
627  egsInformation("Micro-matrix %d from %s:\n",j+1,micro_file);
628  egsInformation(" TB volume fraction: %g\n",tb_vol/nxyz);
629  egsInformation(" BM volume fraction: %g\n",bm_vol/nxyz);
630  egsInformation(" BSC volume fraction: %g\n",bsc_vol/nxyz);
631  egsInformation(" Total: %g\n",(tb_vol+bm_vol+bsc_vol)/nxyz);
632  v_tb[j] = tb_vol*ddx*ddy*ddz;
633  v_bm[j] = bm_vol*ddx*ddy*ddz;
634  v_bsc[j] = bsc_vol*ddx*ddy*ddz;
635  tot_tb += tb_vol;
636  tot_bm += bm_vol;
637  tot_bsc += bsc_vol;
638  }
639  egsInformation("Average volume fractions:\n");
640  egsInformation(" TB volume fraction: %g\n",tot_tb/(nmic*nxyz));
641  egsInformation(" BM volume fraction: %g\n",tot_bm/(nmic*nxyz));
642  egsInformation(" BSC volume fraction: %g\n",tot_bsc/(nmic*nxyz));
643 }
644 
645 EGS_MicroMatrixCluster::~EGS_MicroMatrixCluster() {
646 }
647 #endif
648 
649 EGS_TestMicro::EGS_TestMicro(EGS_Float Dx, EGS_Float Dy, EGS_Float Dz,
650  EGS_Float bsc_t,int tb_med, int bm_med,const char *micro_file,
651  const string &Name) : EGS_BaseGeometry(Name) {
652  micro = new EGS_MicroMatrixCluster(Dx,Dy,Dz,bsc_t,tb_med,bm_med,micro_file);
653  vg = new EGS_VoxelGeometry(micro->mx,micro->my,micro->mz,Dx,Dy,Dz);
654  nr = micro->nx*micro->ny*micro->nz;
655  nreg = nr*micro->mx*micro->my*micro->mz;
656 }
657 
658 void EGS_TestMicro::setMedia(EGS_Input *, int, const int *) {
659  egsFatal("EGS_TestMicro::setMedia(EGS_Input*,int,const int *):\n"
660  " Don't use this method. Media are only set via the micro matrix"
661  " data\n");
662 }
663 
664 EGS_TestMicro::~EGS_TestMicro() {
665 }
666 
667 string EGS_TestMicro::type = "EGS_TestMicro";
668 
669 const static EGS_VHP_LOCAL char *vhp_error_msg1 =
670  "createGeometry(VHP): wrong/missing %s input for micro_test geometry\n";
671 
672 extern "C" {
673 
674  EGS_VHP_EXPORT EGS_BaseGeometry *createGeometry(EGS_Input *input) {
675  string type;
676  int err = input->getInput("type",type);
677  if (!err && type == "micro_test") {
678  bool ok = true;
679  vector<EGS_Float> vs;
680  err = input->getInput("voxel sizes",vs);
681  if (err || vs.size() != 3) {
682  egsWarning(vhp_error_msg1,"'voxel sizes'");
683  ok = false;
684  }
685  EGS_Float bsc_t;
686  err = input->getInput("BSC thickness",bsc_t);
687  if (err) {
688  egsWarning(vhp_error_msg1,"'BSC thickness'");
689  ok = false;
690  }
691  string tb_medium, bm_medium;
692  int err1 = input->getInput("TB medium",tb_medium);
693  int err2 = input->getInput("BM medium",bm_medium);
694  if (err1) {
695  egsWarning(vhp_error_msg1,"'TB medium'");
696  }
697  if (err2) {
698  egsWarning(vhp_error_msg1,"'BM medium'");
699  }
700  if (err1 || err2) {
701  ok = false;
702  }
703  string micro_data;
704  err = input->getInput("micro data file",micro_data);
705  if (err) {
706  egsWarning(vhp_error_msg1,"'micro data file'");
707  ok = false;
708  }
709  if (!ok) {
710  return 0;
711  }
712  int tb_med = EGS_BaseGeometry::addMedium(tb_medium);
713  int bm_med = EGS_BaseGeometry::addMedium(bm_medium);
714  EGS_TestMicro *result = new EGS_TestMicro(vs[0],vs[1],vs[2],bsc_t,
715  tb_med,bm_med,micro_data.c_str());
716  result->setName(input);
717  result->setBoundaryTolerance(input);
718  result->setLabels(input);
719  return result;
720  }
721  string phantom, media;
722  int err1 = input->getInput("phantom data",phantom);
723  int err2 = input->getInput("media data",media);
724  if (err1) egsWarning("createGeometry(EGS_VHPGeometry): "
725  "missing 'phantom data' input\n");
726  if (err2) egsWarning("createGeometry(EGS_VHPGeometry): "
727  "missing 'media data' input\n");
728  if (err1 || err2) {
729  return 0;
730  }
731  vector<int> srange;
732  err1 = input->getInput("slice range",srange);
733  EGS_VHPGeometry *result;
734  if (!err1 && srange.size() == 2 && srange[1] > srange[0]) result =
735  new EGS_VHPGeometry(phantom.c_str(),media.c_str(),srange[0],srange[1]);
736  else {
737  result = new EGS_VHPGeometry(phantom.c_str(),media.c_str());
738  }
739  if (!result->isOK()) {
740  egsWarning("createGeometry(EGS_VHPGeometry): failed to construct "
741  "geometry");
742  delete result;
743  return 0;
744  }
745  result->setName(input);
746  result->setBoundaryTolerance(input);
747  result->setMicros(input);
748  result->setLabels(input);
749  return result;
750  }
751 
752 }
virtual void printInfo() const
Print information about this geometry.
EGS_Input class header file.
string egsJoinPath(const string &first, const string &second)
Join two path variables (or a path and a file name) using the platform specific directory separator a...
int setLabels(EGS_Input *input)
Set the labels from an input block.
Global egspp functions header file.
Voxelized Human Phantom (VHP) geometry: header.
int med
Medium index.
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.
EGS_InfoFunction EGS_EXPORT egsFatal
Always use this function for reporting fatal errors.
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
static int addMedium(const string &medname)
Add a medium or get the index of an existing medium.
A Voxelized Human Phantom (VHP) geometry.
EGS_InfoFunction EGS_EXPORT egsInformation
Always use this function for reporting the progress of a simulation and any other type of information...
int egsGetEndian()
Get the endianess of the machine.
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.