EGSnrc C++ class library  Report PIRS-898 (2021)
Iwan Kawrakow, Ernesto Mainegra-Hing, Frederic Tessier, Reid Townson and Blake Walters
egs_nd_geometry.cpp
Go to the documentation of this file.
1 /*
2 ###############################################################################
3 #
4 # EGSnrc egs++ nd 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, 2005
25 #
26 # Contributors: Frederic Tessier
27 # Reid Townson
28 # Ernesto Mainegra-Hing
29 # Hubert Ho
30 # Randle Taylor
31 #
32 ###############################################################################
33 */
34 
35 
41 #include "egs_nd_geometry.h"
42 #include "egs_input.h"
43 #include "egs_functions.h"
44 #include "egs_transformations.h"
45 
46 #ifdef HAS_GZSTREAM
47  #include "gzstream.h"
48 #endif
49 
50 #include <vector>
51 #include <fstream>
52 
53 using namespace std;
54 
55 EGS_NDGeometry::EGS_NDGeometry(int ng, EGS_BaseGeometry **G,
56  const string &Name, bool O) : EGS_BaseGeometry(Name), N(ng), ortho(O) {
57  nreg = 1;
58  setup();
59  for (int j=0; j<N; j++) {
60  g[j] = G[j];
61  n[j] = nreg;
62  nreg *= g[j]->regions();
63  if (!g[j]->isConvex()) {
64  is_convex = false;
65  }
66  }
67  n[N] = nreg;
68 }
69 
70 EGS_NDGeometry::EGS_NDGeometry(vector<EGS_BaseGeometry *> &G,
71  const string &Name, bool O) : EGS_BaseGeometry(Name), N(G.size()), ortho(O) {
72  nreg = 1;
73  setup();
74  for (int j=0; j<N; j++) {
75  g[j] = G[j];
76  n[j] = nreg;
77  nreg *= g[j]->regions();
78  if (!g[j]->isConvex()) {
79  is_convex = false;
80  }
81  }
82  n[N] = nreg;
83 }
84 
85 EGS_NDGeometry::~EGS_NDGeometry() {
86  for (int j=0; j<N; j++) {
87  int ref = g[j]->deref();
88  if (!ref) {
89  delete g[j];
90  }
91  }
92  delete [] n;
93  delete [] g;
94 }
95 
96 void EGS_NDGeometry::setup() {
97  n = new int [N+1];
98  g = new EGS_BaseGeometry* [N];
99 }
100 
101 void EGS_NDGeometry::printInfo() const {
102  EGS_BaseGeometry::printInfo();
103  egsInformation(" number of dimensions = %d\n",N);
104  for (int j=0; j<N; j++)
105  egsInformation(" dimension %d = %s (type %s)\n",
106  j+1,g[j]->getName().c_str(),g[j]->getType().c_str());
107  egsInformation(
108  "=======================================================\n");
109 }
110 
111 void EGS_NDGeometry::setMedia(EGS_Input *input, int nmed, const int *mind) {
112  EGS_Input *i;
113  med = mind[0];
114  while ((i = input->takeInputItem("set medium"))) {
115  vector<int> inp;
116  int err = i->getInput("set medium",inp);
117  delete i;
118  if (!err) {
119  if (inp.size() == 2) {
120  setMedium(inp[0],inp[0],mind[inp[1]]);
121  }
122  else if (inp.size() == 3) {
123  setMedium(inp[0],inp[1],mind[inp[2]]);
124  }
125  else if (inp.size() == 4) setMedium(inp[0],inp[1],mind[inp[2]],
126  inp[3]);
127  else if (inp.size() == 2*N+1) {
128  setM(0,0,inp,mind[inp[2*N]]);
129  }
130  else egsWarning("EGS_NDGeometry::setMedia(): found %dinputs\n"
131  "in a 'set medium' input. 2, 3, 4, or %d are allowed\n",2*N+1);
132  }
133  else egsWarning("EGS_NDGeometry::setMedia(): wrong 'set medium'"
134  " input\n");
135  }
136 }
137 
138 void EGS_NDGeometry::setM(int ibase, int idim,
139  const vector<int> &ranges, int medium) {
140  int istart = ranges[2*idim], iend = ranges[2*idim+1];
141  if (istart < 0) {
142  istart = 0;
143  }
144  int ndim = n[idim+1] / n[idim];
145  if (iend > ndim) {
146  iend = ndim;
147  }
148  if (idim < N-1) for (int j=istart; j<=iend; j++) {
149  setM(ibase+j*n[idim],idim+1,ranges,medium);
150  }
151  else for (int j=istart; j<=iend; j++) {
152  setMedium(ibase+j*n[idim],ibase+j*n[idim],medium);
153  }
154 }
155 
156 #ifdef EXPLICIT_XYZ
157 
158 string EGS_XYZGeometry::type = "EGS_XYZGeometry";
159 
160 EGS_XYZGeometry::EGS_XYZGeometry(EGS_PlanesX *Xp, EGS_PlanesY *Yp,
161  EGS_PlanesZ *Zp, const string &Name) : EGS_BaseGeometry(Name),
162  xp(Xp), yp(Yp), zp(Zp) {
163  nx = xp->regions();
164  ny = yp->regions();
165  nz = zp->regions();
166  nxy = nx*ny;
167  nreg = nxy*nz;
168  xp->ref();
169  yp->ref();
170  zp->ref();
171  xpos = xp->getPositions();
172  ypos = yp->getPositions();
173  zpos = zp->getPositions();
174 }
175 
176 EGS_XYZGeometry::~EGS_XYZGeometry() {
177  if (!xp->deref()) {
178  delete xp;
179  }
180  if (!yp->deref()) {
181  delete yp;
182  }
183  if (!zp->deref()) {
184  delete zp;
185  }
186 }
187 
188 void EGS_XYZGeometry::printInfo() const {
189  EGS_BaseGeometry::printInfo();
190  xp->printInfo();
191  yp->printInfo();
192  zp->printInfo();
193  egsInformation(
194  "=======================================================\n");
195 }
196 
197 void EGS_XYZGeometry::setMedia(EGS_Input *input, int nmed, const int *mind) {
198  EGS_Input *i;
199  med = mind[0];
200  while ((i = input->takeInputItem("set medium"))) {
201  vector<int> inp;
202  int err = i->getInput("set medium",inp);
203  delete i;
204  if (!err) {
205  if (inp.size() == 2) {
206  setMedium(inp[0],inp[0],mind[inp[1]]);
207  }
208  else if (inp.size() == 3) {
209  setMedium(inp[0],inp[1],mind[inp[2]]);
210  }
211  else if (inp.size() == 4) setMedium(inp[0],inp[1],mind[inp[2]],
212  inp[3]);
213  else if (inp.size() == 7) {
214  int is = inp[0], ie = inp[1];
215  if (is < 0) {
216  is = 0;
217  }
218  if (ie > nx-1) {
219  ie = nx-1;
220  }
221  int js = inp[2], je = inp[3];
222  if (js < 0) {
223  js = 0;
224  }
225  if (je > nxy/nx-1) {
226  je = nxy/nx-1;
227  }
228  int ks = inp[4], ke = inp[5];
229  if (ks < 0) {
230  ks = 0;
231  }
232  if (ke > nreg/nxy-1) {
233  ke = nreg/nxy-1;
234  }
235  //egsInformation("setting medium to %d between %d %d %d %d "
236  // "%d %d\n",mind[inp[6]],is,ie,js,je,ks,ke);
237  for (int i=is; i<=ie; i++) {
238  for (int j=js; j<=je; j++) {
239  for (int k=ks; k<=ke; k++) {
240  int ireg = i + j*nx + k*nxy;
241  //egsWarning("set: %d %d %d %d %d\n",
242  // i,j,k,ireg,mind[inp[6]]);
243  setMedium(ireg,ireg,mind[inp[6]]);
244  }
245  }
246  }
247  }
248  else if (inp.size() == 10) {
249  int is = inp[0], ie = inp[1], idelta = inp[2];
250  if (is < 0) {
251  is = 0;
252  }
253  if (ie > nx-1) {
254  ie = nx-1;
255  }
256  if (idelta <= 0) {
257  idelta = 1;
258  }
259  int js = inp[3], je = inp[4], jdelta = inp[5];
260  if (js < 0) {
261  js = 0;
262  }
263  if (je > nxy/nx-1) {
264  je = nxy/nx-1;
265  }
266  if (jdelta <= 0) {
267  jdelta = 1;
268  }
269  int ks = inp[6], ke = inp[7], kdelta = inp[8];
270  if (ks < 0) {
271  ks = 0;
272  }
273  if (ke > nreg/nxy-1) {
274  ke = nreg/nxy-1;
275  }
276  if (kdelta <= 0) {
277  kdelta = 1;
278  }
279  //egsInformation("setting medium to %d between %d %d %d %d "
280  // "%d %d\n",mind[inp[6]],is,ie,js,je,ks,ke);
281  for (int i=is; i<=ie; i+=idelta) {
282  for (int j=js; j<=je; j+=jdelta) {
283  for (int k=ks; k<=ke; k+=kdelta) {
284  int ireg = i + j*nx + k*nxy;
285  //egsWarning("set: %d %d %d %d %d\n",
286  // i,j,k,ireg,mind[inp[6]]);
287  setMedium(ireg,ireg,mind[inp[9]]);
288  }
289  }
290  }
291  }
292  else egsWarning("EGS_XYZGeometry::setMedia(): found %dinputs\n"
293  "in a 'set medium' input. 2, 3, 4, 7, or 10 are allowed\n");
294  }
295  else egsWarning("EGS_XYZGeometry::setMedia(): wrong 'set medium'"
296  " input\n");
297  }
298 }
299 
300 EGS_XYZGeometry *EGS_XYZGeometry::constructGeometry(const char *dens_file,
301  const char *ramp_file, int dens_or_egsphant_or_interfile) {
302  const static char *func = "EGS_XYZGeometry::constructGeometry";
303  if (!dens_file || !ramp_file) {
304  return 0;
305  }
306  ifstream ramp(ramp_file);
307  if (!ramp) {
308  egsWarning("%s: failed to open CT ramp file %s\n",func,ramp_file);
309  return 0;
310  }
311  vector<string> med_names;
312  vector<EGS_Float> rho_min,rho_max,rho_def;
313 
314  if (dens_or_egsphant_or_interfile == 2) {
315  // interfile format of data
316  int meds;
317  ramp >> meds;
318  for (int i=0; i < meds; i++) {
319  string medname;
320  EGS_Float rmin,rmax,rdef;
321  ramp >> rmin >> rmax >> medname;
322  rdef=(rmin+rmax) / 2.0; // is this correct definition of rdef ??
323  if (ramp.eof() || ramp.fail() || !ramp.good()) {
324  break;
325  }
326  egsInformation("Using medium %s for rho=%g...%g\n",medname.c_str(),
327  rmin,rmax);
328  med_names.push_back(medname);
329  rho_min.push_back(rmin);
330  rho_max.push_back(rmax+1);
331  rho_def.push_back(rdef);
332  }
333  }
334  else {
335  for (EGS_I64 loopCount=0; loopCount<=loopMax; ++loopCount) {
336  if (loopCount == loopMax) {
337  egsFatal("EGS_XYZGeometry::constructGeometry: Too many iterations were required! Input may be invalid, or consider increasing loopMax.");
338  return 0;
339  }
340  // other format of data
341  string medname;
342  EGS_Float rmin,rmax,rdef;
343  ramp >> medname >> rmin >> rmax >> rdef;
344  if (ramp.eof() || ramp.fail() || !ramp.good()) {
345  break;
346  }
347  egsInformation("Using medium %s for rho=%g...%g\n",medname.c_str(),
348  rmin,rmax);
349  med_names.push_back(medname);
350  rho_min.push_back(rmin);
351  rho_max.push_back(rmax);
352  rho_def.push_back(rdef);
353  }
354  }
355  if (med_names.size() < 1) {
356  egsWarning("%s: no media defined in the CT ramp "
357  "file %s!\n",func,ramp_file);
358  return 0;
359  }
360  int Nx, Ny, Nz;
361  EGS_Float *xx, *yy, *zz;
362  float *rho;
363  if (dens_or_egsphant_or_interfile == 0) {
364  int my_endian = egsGetEndian();
365  if (my_endian < 0) egsFatal("%s: machine has an unknown"
366  " endianess!\n",func);
367  ifstream dens(dens_file,ios::binary);
368  if (!dens) {
369  egsWarning("%s: failed to open density matrix file "
370  "%s\n",func,dens_file);
371  return 0;
372  }
373  char their_endian;
374  dens.read(&their_endian,1);
375  if (their_endian != 0 && their_endian != 1)
376  egsFatal("%s: density data created on a machine with"
377  " unknown endianess!\n",func);
378  bool swap = my_endian != their_endian;
379  dens.read((char *) &Nx,sizeof(int));
380  if (swap) {
381  egsSwapBytes(&Nx);
382  }
383  dens.read((char *) &Ny,sizeof(int));
384  if (swap) {
385  egsSwapBytes(&Ny);
386  }
387  dens.read((char *) &Nz,sizeof(int));
388  if (swap) {
389  egsSwapBytes(&Nz);
390  }
391  if (Nx < 1 || Ny < 1 || Nz < 1) {
392  egsWarning("%s: invalid density matrix file: "
393  " Nx=%d Ny=%d Nz=%d\n",func,Nx,Ny,Nz);
394  return 0;
395  }
396  float *x = new float [Nx+1];
397  float *y = new float [Ny+1];
398  float *z = new float [Nz+1];
399  dens.read((char *)x, (Nx+1)*sizeof(float));
400  if (dens.eof() || dens.fail() || !dens.good()) {
401  egsWarning("%s: error while reading %d x-planes from"
402  " file %s\n",func,Nx+1,dens_file);
403  delete [] x;
404  delete [] y;
405  delete [] z;
406  return 0;
407  }
408  dens.read((char *)y, (Ny+1)*sizeof(float));
409  if (dens.eof() || dens.fail() || !dens.good()) {
410  egsWarning("%s: error while reading %d y-planes from"
411  " file %s\n",func,Ny+1,dens_file);
412  delete [] x;
413  delete [] y;
414  delete [] z;
415  return 0;
416  }
417  dens.read((char *)z, (Nz+1)*sizeof(float));
418  if (dens.eof() || dens.fail() || !dens.good()) {
419  egsWarning("%s: error while reading %d z-planes from"
420  " file %s\n",func,Nz+1,dens_file);
421  delete [] x;
422  delete [] y;
423  delete [] z;
424  return 0;
425  }
426  if (swap) {
427  int j;
428  for (j=0; j<Nx; j++) {
429  egsSwapBytes(&x[j]);
430  }
431  for (j=0; j<Ny; j++) {
432  egsSwapBytes(&y[j]);
433  }
434  for (j=0; j<Nz; j++) {
435  egsSwapBytes(&z[j]);
436  }
437  }
438  {
439  int j;
440  xx = new EGS_Float [Nx+1];
441  for (j=0; j<=Nx; j++) {
442  xx[j] = x[j];
443  }
444  yy = new EGS_Float [Ny+1];
445  for (j=0; j<=Ny; j++) {
446  yy[j] = y[j];
447  }
448  zz = new EGS_Float [Nz+1];
449  for (j=0; j<=Nz; j++) {
450  zz[j] = z[j];
451  }
452  delete [] x;
453  delete [] y;
454  delete [] z;
455  }
456  rho = new float [Nx*Ny*Nz];
457  dens.read((char *)rho, Nx*Ny*Nz*sizeof(float));
458  if (dens.fail()) {
459  egsWarning("%s: failed reading mass densities from density matrux file\n",func);
460  delete [] rho;
461  delete [] xx;
462  delete [] yy;
463  delete [] zz;
464  return 0;
465  }
466  if (swap) {
467  for (int j=0; j<Nx*Ny*Nz; j++) {
468  egsSwapBytes(&rho[j]);
469  }
470  }
471  }
472  else if (dens_or_egsphant_or_interfile == 2) {
473  ifstream h_file(dens_file);
474  if (!h_file) {
475  egsWarning("%s: failed to open the interfile header %s\n",func,dens_file);
476  return 0;
477  }
478  string data_file("");
479  int data_type = -1;
480  float scale_x=0.0, scale_y=0.0, scale_z=0.0;
481  for (EGS_I64 loopCount=0; loopCount<=loopMax; ++loopCount) {
482  if (loopCount == loopMax) {
483  egsFatal("EGS_XYZGeometry::constructGeometry: Too many iterations were required! Input may be invalid, or consider increasing loopMax.");
484  return 0;
485  }
486  string line, key, value;
487  size_t pos;
488  getline(h_file, line);
489  if (h_file.eof() || h_file.fail() || !h_file.good()) {
490  break;
491  }
492 
493  pos = line.find(":=");
494  if (pos == string::npos) {
495  continue;
496  }
497  key = line.substr(0, int(pos));
498  value = line.substr(int(pos)+2);
499  while ((key[0] == '!') || (key[0] == ' ')) {
500  key.erase(0,1);
501  }
502  while (key[key.length()-1] == ' ') {
503  key.erase(key.length()-1, 1);
504  }
505  while (value[0] == ' ') {
506  value.erase(0,1);
507  }
508  while (value[value.length()-1] == ' ') {
509  value.erase(value.length()-1, 1);
510  }
511 
512  if (key == "matrix size [1]") {
513  sscanf(value.c_str(), "%u", &Nx);
514  }
515  else if (key == "matrix size [2]") {
516  sscanf(value.c_str(), "%u", &Ny);
517  }
518  else if ((key == "number of slices") || (key == "number of images")) {
519  sscanf(value.c_str(), "%u", &Nz);
520  }
521  else if (key == "scaling factor (mm/pixel) [1]") {
522  sscanf(value.c_str(), "%f", &scale_x);
523  }
524  else if (key == "scaling factor (mm/pixel) [2]") {
525  sscanf(value.c_str(), "%f", &scale_y);
526  }
527  else if (key == "slice thickness (pixels)") {
528  sscanf(value.c_str(), "%f", &scale_z); // ????????????????
529  }
530  else if (key == "name of data file") {
531  data_file = value;
532  }
533  else if (key == "number format") {
534  if ((value == "float") || (value == "FLOAT")) {
535  data_type = 0;
536  }
537  else if ((value == "unsigned integer") || (value == "UNSIGNED INTEGER")) {
538  data_type = 1;
539  }
540  else {
541  egsWarning("%s: unrecognised 'number format' type: %s \n",func,value.c_str());
542  }
543  }
544  }
545  if (Nx < 1 || Ny < 1 || Nz < 1 || scale_x <= 0.0 || scale_y <= 0.0 || scale_z <= 0.0 || data_file == "" || data_type == -1) {
546  egsWarning("%s: invalid interfile header information: "
547  "Nx=%d Ny=%d Nz=%d scale_x=%f scale_y=%f scale_z=%f "
548  "data_file='%s' number_format=%d\n",func,Nx,Ny,Nz,scale_x,scale_y,scale_z,data_file.c_str(),data_type);
549  return 0;
550  }
551  ifstream i_file(data_file.c_str(),ios::binary);
552  if (!i_file) {
553  egsWarning("%s: failed to open interfile data "
554  "%s\n",func,data_file.c_str());
555  return 0;
556  }
557  scale_x /= 10.0;
558  scale_y /= 10.0;
559  scale_z /= 10.0; // convert to cm
560  {
561  int j;
562  xx = new EGS_Float [Nx+1];
563  for (j=0; j<=Nx; j++) {
564  xx[j] = -(float)Nx*scale_x/2.0 + j*scale_x;
565  }
566  yy = new EGS_Float [Ny+1];
567  for (j=0; j<=Ny; j++) {
568  yy[j] = -(float)Ny*scale_y/2.0 + j*scale_y;
569  }
570  zz = new EGS_Float [Nz+1];
571  for (j=0; j<=Nz; j++) {
572  zz[j] = -(float)Nz*scale_z/2.0 + j*scale_z;
573  }
574  }
575  rho = new float [Nx*Ny*Nz];
576  if (data_type == 0) {
577  i_file.read((char *)rho, Nx*Ny*Nz*sizeof(float));
578  }
579  else {
580  unsigned short int *rho_tmp = new unsigned short int [Nx*Ny*Nz];
581  i_file.read((char *)rho_tmp, (Nx*Ny*Nz)*sizeof(unsigned short int));
582  for (int cc = 0; cc<Nx*Ny*Nz; cc++) {
583  rho[cc] = (float)(rho_tmp[cc]);
584  }
585  delete [] rho_tmp;
586  }
587  if (i_file.fail()) {
588  egsWarning("%s: failed reading interfile data\n",func);
589  delete [] rho;
590  delete [] xx;
591  delete [] yy;
592  delete [] zz;
593  return 0;
594  }
595  }
596  else {
597 
598  ifstream tempf(dens_file, ios::binary);
599  istream *data;
600  ifstream textf;
601 #ifdef HAS_GZSTREAM
602  igzstream binf;
603 #endif
604 
605  if (!tempf) {
606  egsWarning("%s: failed to open .egsphant file %s\n",func,dens_file);
607  return 0;
608  }
609 
610  bool is_gzip = (tempf.get() == 0x1f && tempf.get() == 0x8b);
611  tempf.close();
612 
613  if (is_gzip) {
614 #ifdef HAS_GZSTREAM
615  binf.open(dens_file);
616  data = &binf;
617 #else
618  egsWarning("Tried to read gzipped egsphant but egs_ndgeometry was not compiled with gzip support\n");
619  return 0;
620 #endif
621  }
622  else {
623  textf.open(dens_file);
624  data = &textf;
625  }
626 
627 
628  int nmed;
629  (*data) >> nmed;
630  if ((*data).fail()) {
631  egsWarning("%s: failed reading number of media\n",func);
632  return 0;
633  }
634  char buf [1024];
635  int imed;
636  (*data).getline(buf,1023);
637  for (imed=0; imed<nmed; ++imed) {
638  (*data).getline(buf,1023);
639  if ((*data).fail()) {
640  egsWarning("%s: failed reading medium name for %d'th medium\n",func,imed+1);
641  }
642  //egsInformation("Got line <%s>\n",buf);
643  }
644  // ignore estepe
645  EGS_Float dum;
646  for (imed=0; imed<nmed; ++imed) {
647  (*data) >> dum;
648  }
649  (*data) >> Nx >> Ny >> Nz;
650  if ((*data).fail()) {
651  egsWarning("%s: failed reading number of voxels\n",func);
652  return 0;
653  }
654  xx = new EGS_Float [Nx+1];
655  yy = new EGS_Float [Ny+1];
656  zz = new EGS_Float [Nz+1];
657  int j;
658  for (j=0; j<=Nx; ++j) {
659  (*data) >> xx[j];
660  }
661  if ((*data).fail()) {
662  egsWarning("%s: failed reading x-planes\n",func);
663  delete [] xx;
664  delete [] yy;
665  delete [] zz;
666  return 0;
667  }
668  for (j=0; j<=Ny; ++j) {
669  (*data) >> yy[j];
670  }
671  if ((*data).fail()) {
672  egsWarning("%s: failed reading y-planes\n",func);
673  delete [] xx;
674  delete [] yy;
675  delete [] zz;
676  return 0;
677  }
678  for (j=0; j<=Nz; ++j) {
679  (*data) >> zz[j];
680  }
681  if ((*data).fail()) {
682  egsWarning("%s: failed reading z-planes\n",func);
683  delete [] xx;
684  delete [] yy;
685  delete [] zz;
686  return 0;
687  }
688  int nr = Nx*Ny*Nz;
689  //int med;
690  //for(j=0; j<nr; ++j) data >> med;
691  (*data).getline(buf,1023);
692  for (int iz=0; iz<Nz; ++iz) {
693  for (int iy=0; iy<Ny; ++iy) {
694  (*data).getline(buf,1023);
695  }
696  (*data).getline(buf,1023);
697  }
698  if ((*data).fail()) {
699  egsWarning("%s: failed reading media indeces matrix\n",func);
700  delete [] xx;
701  delete [] yy;
702  delete [] zz;
703  return 0;
704  }
705  rho = new float [nr];
706  for (j=0; j<nr; ++j) {
707  (*data) >> rho[j];
708  }
709  if ((*data).fail()) {
710  egsWarning("%s: failed reading mass density matrix\n",func);
711  delete [] xx;
712  delete [] yy;
713  delete [] zz;
714  delete [] rho;
715  return 0;
716  }
717  }
718  EGS_PlanesX *xp = new EGS_PlanesX(Nx+1,xx,"",EGS_XProjector("x-planes"));
719  EGS_PlanesY *yp = new EGS_PlanesY(Ny+1,yy,"",EGS_YProjector("y-planes"));
720  EGS_PlanesZ *zp = new EGS_PlanesZ(Nz+1,zz,"",EGS_ZProjector("z-planes"));
721  EGS_XYZGeometry *result = new EGS_XYZGeometry(xp,yp,zp);
722  EGS_Float rhomin=1e30,rhomax=0;
723  int j;
724  for (j=0; j<Nx*Ny*Nz; j++) {
725  if (rho[j] > rhomax) {
726  rhomax = rho[j];
727  }
728  if (rho[j] < rhomin) {
729  rhomin = rho[j];
730  }
731  }
732  egsInformation("Min. density found: %g\n",rhomin);
733  egsInformation("Max. density found: %g\n",rhomax);
734  int nmed = med_names.size();
735  int *imed = new int [nmed];
736  for (j=0; j<nmed; j++) {
737  imed[j] = -2;
738  }
739  for (j=0; j<Nx*Ny*Nz; j++) {
740  int i;
741  for (i=0; i<nmed; i++) {
742  if (rho[j] >= rho_min[i] && rho[j] < rho_max[i]) {
743  break;
744  }
745  }
746  if (i >= nmed) {
747  egsWarning("%s: no material defined for density %g, leaving voxel %d"
748  " as vacuum\n",func,rho[j],j);
749  result->setMedium(j,j,-1);
750  }
751  else {
752  if (imed[i] < -1) {
753  imed[i] = addMedium(med_names[i]);
754  egsInformation("Using medium %s as mednum %d\n",
755  med_names[i].c_str(),imed[i]);
756  }
757  result->setMedium(j,j,imed[i]);
758  if (rho_def[i] > 0) {
759  EGS_Float rrho = rho[j]/rho_def[i];
760  if (fabs(rrho-1) > epsilon) {
761  result->setRelativeRho(j,j,rrho);
762  }
763  }
764  }
765  }
766  delete [] rho;
767  delete [] imed;
768  return result;
769 }
770 
771 string EGS_DeformedXYZ::def_type = "EGS_DeformedXYZ";
772 
773 char EGS_DeformedXYZ::tetrahedra[] = {0,2,3,6, 0,6,3,7, 0,4,6,7,
774  7,0,1,3, 7,1,0,4, 7,5,1,4
775  };
776 
777 char EGS_DeformedXYZ::plane_order[] = {0,1,2, 0,2,3, 0,3,1, 1,3,2};
778 
779 char EGS_DeformedXYZ::enter_tetra1[] = {0, 3, 4, 0, 0, 2};
780 char EGS_DeformedXYZ::enter_plane1[] = {1, 1, 1, 3, 0, 3};
781 char EGS_DeformedXYZ::enter_tetra2[] = {2, 5, 5, 1, 3, 5};
782 char EGS_DeformedXYZ::enter_plane2[] = {0, 0, 2, 3, 2, 2};
783 
784 
785 char EGS_DeformedXYZ::tetra_data[] = {
786  2,-1, 2, // plane 0 in tetrahedron 0
787  -1, 0, 1, // plane 1 in tetrahedron 0
788  0,-1, 3, // plane 2 in tetrahedron 0
789  1, 1, 4, // plane 3 in tetrahedron 0
790 //------------------------------------------------------------------
791  -1, 0, 0, // plane 0 in tetrahedron 1
792  -1, 0, 3, // plane 1 in tetrahedron 1
793  -1, 0, 2, // plane 2 in tetrahedron 1
794  1, 1, 5, // plane 3 in tetrahedron 1
795 //------------------------------------------------------------------
796  0,-1, 5, // plane 0 in tetrahedron 2
797  -1, 0, 1, // plane 1 in tetrahedron 2
798  -1, 0, 4, // plane 2 in tetrahedron 2
799  2, 1, 0, // plane 3 in tetrahedron 2
800 //------------------------------------------------------------------
801  -1, 0, 4, // plane 0 in tetrahedron 3
802  0, 1, 0, // plane 1 in tetrahedron 3
803  -1, 0, 1, // plane 2 in tetrahedron 3
804  2,-1, 5, // plane 3 in tetrahedron 3
805 //------------------------------------------------------------------
806  -1, 0, 3, // plane 0 in tetrahedron 4
807  -1, 0, 2, // plane 1 in tetrahedron 4
808  -1, 0, 5, // plane 2 in tetrahedron 4
809  1,-1, 0, // plane 3 in tetrahedron 4
810 //------------------------------------------------------------------
811  0, 1, 2, // plane 0 in tetrahedron 5
812  -1, 0, 4, // plane 1 in tetrahedron 5
813  2, 1, 3, // plane 2 in tetrahedron 5
814  1,-1, 1 // plane 3 in tetrahedron 5
815 };
816 
817 
818 EGS_DeformedXYZ::EGS_DeformedXYZ(EGS_PlanesX *Xp, EGS_PlanesY *Yp,
819  EGS_PlanesZ *Zp, const char *defFile, const string &Name) :
820  EGS_XYZGeometry(Xp,Yp,Zp,Name), vectors(0) {
821  setDeformations(defFile);
822  np[0] = nx;
823  np[1] = ny;
824  np[2] = nz;
825  nxyp1 = nx + ny + 1;
826 }
827 
828 int EGS_DeformedXYZ::setDeformations(const char *defFile) {
829  ifstream data(defFile,ios::binary);
830  if (!data) {
831  egsWarning("Failed to open deformations file %s\n",defFile);
832  return 1;
833  }
834  int nvec;
835  data.read((char *)&nvec,sizeof(int));
836  int nneed = (nx+1)*(ny+1)*(nz+1);
837  if (nvec != nneed) {
838  egsWarning("Inconsistent deformations file: %d instead of %d vectors\n",
839  nvec,nneed);
840  return 2;
841  }
842  if (!vectors) {
843  vectors = new EGS_Vector [nvec];
844  }
845  float tmp[3];
846  int i,j,k;
847  for (j=0; j<nvec; ++j) {
848  data.read((char *)tmp,3*sizeof(float));
849  if (data.fail()) {
850  egsWarning("Error while readinf vector %d from file\n",j+1);
851  return 3;
852  }
853  vectors[j] = EGS_Vector(tmp[0],tmp[1],tmp[2]);
854  }
855 
856  // add plane positions to vectors
857  for (k=0; k<=nz; ++k) for (j=0; j<=ny; ++j) for (i=0; i<=nx; ++i) {
858  int index = i + j*(nx+1) + k*(nx+1)*(ny+1);
859  vectors[index] += EGS_Vector(xpos[i],ypos[j],zpos[k]);
860  }
861 
862  for (j=0; j<24; ++j) {
863  int i = tetrahedra[j];
864  switch (i) {
865  case 0:
866  tnodes[j] = 0;
867  break;
868  case 1:
869  tnodes[j] = 1;
870  break;
871  case 2:
872  tnodes[j] = nx+1;
873  break;
874  case 3:
875  tnodes[j] = nx+2;
876  break;
877  case 4:
878  tnodes[j] = (nx+1)*(ny+1);
879  break;
880  case 5:
881  tnodes[j] = (nx+1)*(ny+1)+1;
882  break;
883  case 6:
884  tnodes[j] = (nx+1)*(ny+1)+nx+1;
885  break;
886  case 7:
887  tnodes[j] = (nx+1)*(ny+1)+nx+2;
888  break;
889  }
890  }
891 
892  nreg = 6*nx*ny*nz;
893 
894  return 0;
895 }
896 
897 EGS_DeformedXYZ::~EGS_DeformedXYZ() {
898  if (vectors) {
899  delete [] vectors;
900  }
901 }
902 
903 string EGS_XYZRepeater::type = "EGS_XYZRepeater";
904 
905 EGS_XYZRepeater::EGS_XYZRepeater(EGS_Float xmin, EGS_Float xmax,
906  EGS_Float ymin, EGS_Float ymax,
907  EGS_Float zmin, EGS_Float zmax,
908  int Nx, int Ny, int Nz,
909  EGS_BaseGeometry *G,
910  const string &Name) :
911  EGS_BaseGeometry(Name), g(G), nx(Nx), ny(Ny), nz(Nz) {
912  dx = (xmax - xmin)/nx;
913  dxi = 1/dx;
914  ax = -xmin*dxi;
915  dy = (ymax - ymin)/ny;
916  dyi = 1/dy;
917  ay = -ymin*dyi;
918  dz = (zmax - zmin)/nz;
919  dzi = 1/dz;
920  az = -zmin*dzi;
921  EGS_PlanesX *xp;
922  EGS_PlanesY *yp;
923  EGS_PlanesZ *zp;
924  xp = new EGS_PlanesX(xmin,dx,nx,"",EGS_XProjector("x-planes"));
925  yp = new EGS_PlanesY(ymin,dy,ny,"",EGS_YProjector("y-planes"));
926  zp = new EGS_PlanesZ(zmin,dz,nz,"",EGS_ZProjector("z-planes"));
927  xyz = new EGS_XYZGeometry(xp,yp,zp);
928  xyz->ref();
929  g->ref();
930  nxy = nx*ny;
931  nxyz = nx*ny*nz;
932  ng = g->regions();
933  nreg = nxyz*ng + 1;
934  translation = new EGS_Vector [nxyz];
935  for (int ix=0; ix<nx; ix++)
936  for (int iy=0; iy<ny; iy++)
937  for (int iz=0; iz<nz; iz++)
938  translation[ix+iy*nx+iz*nxy] = EGS_Vector(xmin + dx*(0.5+ix),
939  ymin + dy*(0.5+iy),
940  zmin + dz*(0.5+iz));
941 
942 }
943 
944 EGS_XYZRepeater::~EGS_XYZRepeater() {
945  if (!g->deref()) {
946  delete g;
947  }
948  if (!xyz->deref()) {
949  delete xyz;
950  }
951  delete [] translation;
952 }
953 
954 void EGS_XYZRepeater::printInfo() const {
955  EGS_BaseGeometry::printInfo();
956  egsInformation("%d repetitions between xmin=%g and xmax=%g\n",
957  nx,xyz->getXPositions()[0],xyz->getXPositions()[nx]);
958  egsInformation("%d repetitions between ymin=%g and ymax=%g\n",
959  ny,xyz->getYPositions()[0],xyz->getYPositions()[ny]);
960  egsInformation("%d repetitions between zmin=%g and zmax=%g\n",
961  nz,xyz->getZPositions()[0],xyz->getZPositions()[nz]);
962  const char *med_name = getMediumName(med);
963  if (med_name) {
964  egsInformation("box is filled with %s\n",med_name);
965  }
966  else {
967  egsInformation("box is filled with vacuum\n");
968  }
969  egsInformation("Repeated geometry:\n");
970  g->printInfo();
971 }
972 
973 void EGS_XYZGeometry::voxelizeGeometry(EGS_Input *input) {
974  string gname;
975  int err = input->getInput("voxelize geometry",gname);
976  if (err) {
977  return;
978  }
979  egsInformation("\nEGS_XYZGeometry(%s)\n",name.c_str());
980  egsInformation(" setting media from geometry %s\n",gname.c_str());
981  EGS_BaseGeometry *geometry = EGS_BaseGeometry::getGeometry(gname);
982  if (!geometry) {
983  egsInformation(" this geometry does not exist -> media will be not set.\n");
984  return;
985  }
986  EGS_AffineTransform *T = EGS_AffineTransform::getTransformation(input);
987  if (region_media) {
988  delete [] region_media;
989  }
990  if (rhor) {
991  delete [] rhor;
992  rhor = 0;
993  }
994  if (bfactor) {
995  delete [] bfactor;
996  bfactor = 0;
997  }
998  region_media = new short [nreg];
999 
1000  bool hrs = geometry->hasRhoScaling();
1001  if (hrs) {
1002  has_rho_scaling = true;
1003  rhor = new EGS_Float [nreg];
1004  }
1005  else {
1006  has_rho_scaling = false;
1007  }
1008 
1009  bool hbs = geometry->hasBScaling();
1010  if (hbs) {
1011  has_B_scaling = true;
1012  bfactor = new EGS_Float [nreg];
1013  }
1014  else {
1015  has_B_scaling = false;
1016  }
1017 
1018  EGS_Vector v1(xp->position(0),yp->position(0),zp->position(0));
1019  EGS_Vector v2(xp->position(nx),yp->position(ny),zp->position(nz));
1020  egsInformation(" top/left/front corner : (%g,%g,%g)\n",v1.x,v1.y,v1.z);
1021  egsInformation(" bottom/right/back corner: (%g,%g,%g)\n",v2.x,v2.y,v2.z);
1022  if (T) {
1023  egsInformation(" the above are transformed as follows before looking up media:\n");
1024  T->transform(v1);
1025  T->transform(v2);
1026  egsInformation(" top/left/front corner : (%g,%g,%g)\n",v1.x,v1.y,v1.z);
1027  egsInformation(" bottom/right/back corner: (%g,%g,%g)\n",v2.x,v2.y,v2.z);
1028  }
1029  for (int ix=0; ix<nx; ++ix) for (int iy=0; iy<ny; iy++) for (int iz=0; iz<nz; ++iz) {
1030  int ir = ix + iy*nx + iz*nxy;
1031  EGS_Vector v(0.5*(xp->position(ix)+xp->position(ix+1)),
1032  0.5*(yp->position(iy)+yp->position(iy+1)),
1033  0.5*(zp->position(iz)+zp->position(iz+1)));
1034  if (T) {
1035  T->transform(v);
1036  }
1037  int ireg = geometry->isWhere(v);
1038  if (ireg >= 0) {
1039  region_media[ir] = geometry->medium(ireg);
1040  if (hrs) {
1041  rhor[ir] = geometry->getRelativeRho(ir);
1042  }
1043  if (hbs) {
1044  bfactor[ir] = geometry->getBScaling(ir);
1045  }
1046  }
1047  else {
1048  region_media[ir] = -1;
1049  if (hrs) {
1050  rhor[ir] = 1;
1051  }
1052  if (hbs) {
1053  bfactor[ir] = 1;
1054  }
1055  }
1056  }
1057  vector<string> options;
1058  options.push_back("no");
1059  options.push_back("yes");
1060  bool delete_geometry = input->getInput("delete geometry",options,0);
1061  if (delete_geometry) {
1062  if (geometry->deref() == -1) {
1063  egsInformation(" deleting geometry %s as requested\n",geometry->getName().c_str());
1064  delete geometry;
1065  }
1066  else {
1067  egsInformation(" not deleting geometry %s as requested due to remaining references\n",
1068  geometry->getName().c_str());
1069  geometry->ref();
1070  }
1071  }
1072  egsInformation("\n");
1073 }
1074 
1075 const char *err_msg1 = "createGeometry(EGS_XYZRepeater)";
1076 
1077 #endif
1078 
1079 string EGS_NDGeometry::type = "EGS_NDGeometry";
1080 
1081 extern "C" {
1082 
1083  EGS_NDG_EXPORT EGS_BaseGeometry *createGeometry(EGS_Input *input) {
1084 #ifdef EXPLICIT_XYZ
1085  const static char *func = "createGeometry(XYZ)";
1086  string type;
1087  int is_xyz = input->getInput("type",type);
1088  if (!is_xyz && input->compare("EGS_XYZRepeater",type)) {
1089  string base,medium;
1090  vector<EGS_Float> xr, yr, zr;
1091  int err1 = input->getInput("repeated geometry",base);
1092  int err2 = input->getInput("repeat x",xr);
1093  int err3 = input->getInput("repeat y",yr);
1094  int err4 = input->getInput("repeat z",zr);
1095  int err5 = input->getInput("medium",medium);
1096  EGS_BaseGeometry *g = 0;
1097  if (err1) {
1098  egsWarning("%s: missing 'repeated geometry' input\n",err_msg1);
1099  }
1100  else {
1101  g = EGS_BaseGeometry::getGeometry(base);
1102  if (!g) {
1103  egsWarning("%s: no geometry named %s exists\n",err_msg1,
1104  base.c_str());
1105  err1 = 1;
1106  }
1107  }
1108  if (err2 || xr.size() != 3) {
1109  err1 = 1;
1110  egsWarning("%s: wrong/missing 'repeat x' input\n",err_msg1);
1111  }
1112  if (err3 || yr.size() != 3) {
1113  err1 = 1;
1114  egsWarning("%s: wrong/missing 'repeat y' input\n",err_msg1);
1115  }
1116  if (err4 || zr.size() != 3) {
1117  err1 = 1;
1118  egsWarning("%s: wrong/missing 'repeat z' input\n",err_msg1);
1119  }
1120  if (err1) {
1121  return 0;
1122  }
1123  EGS_Float xmin = xr[0], xmax = xr[1];
1124  int nx = (int)(xr[2]+0.1);
1125  if (xmin >= xmax || nx < 1) {
1126  egsWarning("%s: wrong 'repeat x' input xmin=%g xmax=%g nx=%d\n",
1127  xmin,xmax,nx);
1128  err1 = 1;
1129  }
1130  EGS_Float ymin = yr[0], ymax = yr[1];
1131  int ny = (int)(yr[2]+0.1);
1132  if (ymin >= ymax || ny < 1) {
1133  egsWarning("%s: wrong 'repeat y' input ymin=%g ymax=%g ny=%d\n",
1134  ymin,ymax,ny);
1135  err1 = 1;
1136  }
1137  EGS_Float zmin = zr[0], zmax = zr[1];
1138  int nz = (int)(zr[2]+0.1);
1139  if (zmin >= zmax || nz < 1) {
1140  egsWarning("%s: wrong 'repeat z' input zmin=%g zmax=%g nz=%d\n",
1141  zmin,zmax,nz);
1142  err1 = 1;
1143  }
1144  if (err1) {
1145  return 0;
1146  }
1147 
1148  EGS_XYZRepeater *result =
1149  new EGS_XYZRepeater(xmin,xmax,ymin,ymax,zmin,zmax,nx,ny,nz,g);
1150  result->setName(input);
1151  result->setBoundaryTolerance(input);
1152  if (!err5) {
1153  result->setMedium(medium);
1154  }
1155  result->setXYZLabels(input);
1156  result->setLabels(input);
1157  return result;
1158  }
1159  else if (!is_xyz && input->compare("EGS_XYZGeometry",type)) {
1160  string dens_file, ramp_file, egsphant_file, interfile_file;
1161  int ierr1 = input->getInput("density matrix",dens_file);
1162  int ierr2 = input->getInput("ct ramp",ramp_file);
1163  int ierr3 = input->getInput("egsphant file",egsphant_file);
1164  int ierr4 = input->getInput("interfile header",interfile_file);
1165  int dens_or_egsphant_or_interfile = -1;
1166  if (!ierr1) {
1167  dens_or_egsphant_or_interfile = 0;
1168  }
1169  else if (!ierr3) {
1170  dens_or_egsphant_or_interfile = 1;
1171  dens_file = egsphant_file;
1172  }
1173  else if (!ierr4) {
1174  dens_or_egsphant_or_interfile = 2;
1175  dens_file = interfile_file;
1176  }
1177  if (dens_or_egsphant_or_interfile >= 0 || !ierr2) {
1178  if (dens_or_egsphant_or_interfile < 0) {
1179  egsWarning("%s: no 'density matrix', 'egsphant file' or 'interfile header' input\n",func);
1180  return 0;
1181  }
1182  if (ierr2) {
1183  egsWarning("%s: no 'ct ramp' input\n",func);
1184  return 0;
1185  }
1186  EGS_XYZGeometry *result =
1187  EGS_XYZGeometry::constructGeometry(dens_file.c_str(),ramp_file.c_str(),dens_or_egsphant_or_interfile);
1188 
1189  if (result) {
1190  result->setName(input);
1191  result->setBoundaryTolerance(input);
1192  result->setBScaling(input);
1193  }
1194 
1195  return result;
1196  }
1197  vector<EGS_Float> xpos, ypos, zpos, xslab, yslab, zslab;
1198  int ix = input->getInput("x-planes",xpos);
1199  int iy = input->getInput("y-planes",ypos);
1200  int iz = input->getInput("z-planes",zpos);
1201  int ix1 = input->getInput("x-slabs",xslab);
1202  int iy1 = input->getInput("y-slabs",yslab);
1203  int iz1 = input->getInput("z-slabs",zslab);
1204  int nx=0, ny=0, nz=0;
1205  if (!ix1) {
1206  if (xslab.size() != 3) {
1207  egsWarning("createGeometry(XYZ): exactly 3 inputs are required"
1208  " when using 'x-slabs' input method\n");
1209  ix1 = 1;
1210  }
1211  else {
1212  nx = (int)(xslab[2]+0.1);
1213  if (nx < 1) {
1214  egsWarning("createGeometry(XYZ): number of slabs must be"
1215  " positive!\n");
1216  ix1 = 1;
1217  }
1218  if (xslab[1] <= 0) {
1219  egsWarning("createGeometry(XYZ): slab thickness must be"
1220  " positive!\n");
1221  ix1 = 1;
1222  }
1223  }
1224  }
1225  if (ix && ix1) {
1226  egsWarning("createGeometry(XYZ): wrong/missing 'x-planes' "
1227  "and 'x-slabs' input\n");
1228  return 0;
1229  }
1230  if (!iy1) {
1231  if (yslab.size() != 3) {
1232  egsWarning("createGeometry(XYZ): exactly 3 inputs are required"
1233  " when using 'y-slabs' input method\n");
1234  iy1 = 1;
1235  }
1236  else {
1237  ny = (int)(yslab[2]+0.1);
1238  if (ny < 1) {
1239  egsWarning("createGeometry(XYZ): number of slabs must be"
1240  " positive!\n");
1241  iy1 = 1;
1242  }
1243  if (yslab[1] <= 0) {
1244  egsWarning("createGeometry(XYZ): slab thickness must be"
1245  " positive!\n");
1246  iy1 = 1;
1247  }
1248  }
1249  }
1250  if (iy && iy1) {
1251  egsWarning("createGeometry(XYZ): wrong/missing 'y-planes' "
1252  "and 'y-slabs' input\n");
1253  return 0;
1254  }
1255  if (!iz1) {
1256  if (zslab.size() != 3) {
1257  egsWarning("createGeometry(XYZ): exactly 3 inputs are required"
1258  " when using 'z-slabs' input method\n");
1259  iz1 = 1;
1260  }
1261  else {
1262  nz = (int)(zslab[2]+0.1);
1263  if (nz < 1) {
1264  egsWarning("createGeometry(XYZ): number of slabs must be"
1265  " positive!\n");
1266  iz1 = 1;
1267  }
1268  if (zslab[1] <= 0) {
1269  egsWarning("createGeometry(XYZ): slab thickness must be"
1270  " positive!\n");
1271  iz1 = 1;
1272  }
1273  }
1274  }
1275  if (iz && iz1) {
1276  egsWarning("createGeometry(XYZ): wrong/missing 'z-planes' "
1277  "and 'z-slabs' input\n");
1278  return 0;
1279  }
1280  EGS_PlanesX *xp = !ix1 ?
1281  new EGS_PlanesX(xslab[0],xslab[1],nx,"",EGS_XProjector("x-planes")) :
1282  new EGS_PlanesX(xpos,"",EGS_XProjector("x-planes"));
1283  EGS_PlanesY *yp = !iy1 ?
1284  new EGS_PlanesY(yslab[0],yslab[1],ny,"",EGS_YProjector("y-planes")) :
1285  new EGS_PlanesY(ypos,"",EGS_YProjector("y-planes"));
1286  EGS_PlanesZ *zp = !iz1 ?
1287  new EGS_PlanesZ(zslab[0],zslab[1],nz,"",EGS_ZProjector("z-planes")) :
1288  new EGS_PlanesZ(zpos,"",EGS_ZProjector("z-planes"));
1289  EGS_XYZGeometry *result = new EGS_XYZGeometry(xp,yp,zp);
1290 
1291  if (result) {
1292  egsWarning("**********************************************\n");
1293  EGS_BaseGeometry *g = result;
1294  result->setName(input);
1295  result->setBoundaryTolerance(input);
1296  g->setMedia(input);
1297  result->voxelizeGeometry(input);
1298 
1299  // labels
1300  result->setXYZLabels(input);
1301  g->setLabels(input);
1302  }
1303  return result;
1304  }
1305 #endif
1306  vector<EGS_BaseGeometry *> dims;
1307  EGS_Input *ij;
1308  int error = 0;
1309  while ((ij = input->takeInputItem("geometry",false)) != 0) {
1310  EGS_BaseGeometry *g = EGS_BaseGeometry::createSingleGeometry(ij);
1311  if (g) {
1312  dims.push_back(g);
1313  g->ref();
1314  }
1315  else {
1316  error++;
1317  }
1318  delete ij;
1319  }
1320  vector<string> gnames;
1321  int err1 = input->getInput("dimensions",gnames);
1322  if (!err1) {
1323  for (unsigned int j=0; j<gnames.size(); j++) {
1324  EGS_BaseGeometry *g = EGS_BaseGeometry::getGeometry(gnames[j]);
1325  if (g) {
1326  dims.push_back(g);
1327  g->ref();
1328  }
1329  else {
1330  egsWarning("Geometry %s does not exist\n",gnames[j].c_str());
1331  error++;
1332  }
1333  }
1334  }
1335  if (error) {
1336  egsWarning("createGeometry(ND_Geometry): %d errors while "
1337  "creating/geting geometries defining individual dimensions\n",error);
1338  return 0;
1339  }
1340  if (dims.size() < 2) {
1341  egsWarning("createGeometry(ND_Geometry): why do you want to "
1342  "construct a ND geometry with a single dimension?\n");
1343  input->print(0,cerr);
1344  for (int j=0; j<gnames.size(); j++) egsWarning("dimension %d: %s\n",
1345  j+1,gnames[j].c_str());
1346  }
1347  int n_concav = 0;
1348  for (int j=0; j<dims.size(); j++) if (!dims[j]->isConvex()) {
1349  n_concav++;
1350  }
1351  bool is_ok = true;
1352  if (n_concav > 1) {
1353  egsWarning("createGeometry(ND_Geometry): a ND geometry can not have "
1354  " more than one non-convex dimension, yours has %d\n",n_concav);
1355  is_ok = false;
1356  }
1357  else if (n_concav == 1) {
1358  if (dims[dims.size()-1]->isConvex()) {
1359  egsWarning("createGeometry(ND_Geometry): the non-convex "
1360  "dimension must be the last dimension\n");
1361  is_ok = false;
1362  }
1363  }
1364  if (!is_ok) {
1365  for (int j=0; j<dims.size(); j++)
1366  if (dims[j]->deref()) {
1367  delete dims[j];
1368  }
1369  return 0;
1370  }
1371  int hn_method=0;
1372  err1 = input->getInput("hownear method",hn_method);
1373  EGS_BaseGeometry *result;
1374  if (!err1 && hn_method == 1) {
1375  result = new EGS_NDGeometry(dims,"",false);
1376  }
1377  else {
1378  result = new EGS_NDGeometry(dims);
1379  }
1380  result->setName(input);
1381  result->setBoundaryTolerance(input);
1382  result->setMedia(input);
1383  result->setLabels(input);
1384  return result;
1385  }
1386 
1387 
1388  void EGS_XYZGeometry::setXYZLabels(EGS_Input *input) {
1389 
1390  // x,y,z labels
1391  string inp;
1392  int err;
1393 
1394  err = input->getInput("set x label", inp);
1395  if (!err) {
1396  xp->setLabels(inp);
1397  }
1398 
1399  err = input->getInput("set y label", inp);
1400  if (!err) {
1401  yp->setLabels(inp);
1402  }
1403 
1404  err = input->getInput("set z label", inp);
1405  if (!err) {
1406  zp->setLabels(inp);
1407  }
1408  }
1409 
1410 
1411  void EGS_NDGeometry::ndRegions(int r, int dim, int dimk, int k, vector<int> &regs) {
1412 
1413  // skip looping over selected dimension
1414  if (dim == dimk) {
1415  r += k*n[dimk];
1416  if (dim == N-1) {
1417  regs.push_back(r);
1418  }
1419  else {
1420  ndRegions(r, dim+1, dimk, k, regs);
1421  }
1422  }
1423 
1424  // last dimension: end recursion and record global region number
1425  else if (dim == N-1) {
1426  for (int j=0; j<g[dim]->regions(); j++) {
1427  regs.push_back(r+j*n[dim]);
1428  }
1429  }
1430 
1431  // keep collecting by recursion
1432  else {
1433  for (int j=0; j<g[dim]->regions(); j++) {
1434  ndRegions(r + j*n[dim], dim+1, dimk, k, regs);
1435  }
1436  }
1437  }
1438 
1439 
1440  void EGS_NDGeometry::getLabelRegions(const string &str, vector<int> &regs) {
1441 
1442  // label defined in the sub-geometries
1443  vector<int> local_regs;
1444  for (int i=0; i<N; i++) {
1445  local_regs.clear();
1446  if (g[i]) {
1447  g[i]->getLabelRegions(str, local_regs);
1448  }
1449  if (local_regs.size() == 0) {
1450  continue;
1451  }
1452 
1453  // recurse to collect all global regions comprised in each local region
1454  for (int j=0; j<local_regs.size(); j++) {
1455  ndRegions(0, 0, i, local_regs[j], regs);
1456  }
1457  }
1458 
1459  // label defined in self (nd input block)
1460  EGS_BaseGeometry::getLabelRegions(str, regs);
1461 
1462  }
1463 
1464 
1465  void EGS_XYZGeometry::getLabelRegions(const string &str, vector<int> &regs) {
1466 
1467  vector<int> local_regs;
1468 
1469  // x plane labels
1470  local_regs.clear();
1471  xp->getLabelRegions(str, local_regs);
1472  for (int i=0; i<local_regs.size(); i++) {
1473  for (int j=0; j<ny; j++) {
1474  for (int k=0; k<nz; k++) {
1475  regs.push_back(local_regs[i] + nx*j + nxy*k);
1476  }
1477  }
1478  }
1479 
1480  // y plane labels
1481  local_regs.clear();
1482  yp->getLabelRegions(str, local_regs);
1483  for (int j=0; j<local_regs.size(); j++) {
1484  for (int i=0; i<nx; i++) {
1485  for (int k=0; k<nz; k++) {
1486  regs.push_back(i + nx*local_regs[j] + nxy*k);
1487  }
1488  }
1489  }
1490 
1491  // z plane labels
1492  local_regs.clear();
1493  zp->getLabelRegions(str, local_regs);
1494  for (int k=0; k<local_regs.size(); k++) {
1495  for (int i=0; i<nx; i++) {
1496  for (int j=0; j<ny; j++) {
1497  regs.push_back(i + nx*j + nxy*local_regs[k]);
1498  }
1499  }
1500  }
1501 
1502  // label defined in self (xyz input block)
1503  EGS_BaseGeometry::getLabelRegions(str, regs);
1504 
1505  }
1506 
1507 
1508  void EGS_XYZRepeater::getLabelRegions(const string &str, vector<int> &regs) {
1509 
1510  vector<int> local_regs;
1511 
1512  // label in repeated geometry
1513  local_regs.clear();
1514  g->getLabelRegions(str, local_regs);
1515  for (int i=0; i<nx; i++) {
1516  for (int j=0; j<ny; j++) {
1517  for (int k=0; k<nz; k++) {
1518  for (int r=0; r<local_regs.size(); r++) {
1519  regs.push_back(ng*(i + nx*j + nxy*k) + local_regs[r]);
1520  }
1521  }
1522  }
1523  }
1524 
1525  // x plane labels
1526  local_regs.clear();
1527  xyz->getXLabelRegions(str, local_regs);
1528  for (int i=0; i<local_regs.size(); i++) {
1529  for (int j=0; j<ny; j++) {
1530  for (int k=0; k<nz; k++) {
1531  for (int r=0; r<ng; r++) {
1532  regs.push_back(ng*(local_regs[i] + nx*j + nxy*k) + r);
1533  }
1534  }
1535  }
1536  }
1537 
1538  // y plane labels
1539  local_regs.clear();
1540  xyz->getYLabelRegions(str, local_regs);
1541  for (int j=0; j<local_regs.size(); j++) {
1542  for (int i=0; i<nx; i++) {
1543  for (int k=0; k<nz; k++) {
1544  for (int r=0; r<ng; r++) {
1545  regs.push_back(ng*(i + nx*local_regs[j] + nxy*k) + r);
1546  }
1547  }
1548  }
1549  }
1550 
1551  // z plane labels
1552  local_regs.clear();
1553  xyz->getZLabelRegions(str, local_regs);
1554  for (int k=0; k<local_regs.size(); k++) {
1555  for (int i=0; i<nx; i++) {
1556  for (int j=0; j<ny; j++) {
1557  for (int r=0; r<ng; r++) {
1558  regs.push_back(ng*(i + nx*j + nxy*local_regs[k]) + r);
1559  }
1560  }
1561  }
1562  }
1563 
1564  // label defined in self (repeater input block)
1565  EGS_BaseGeometry::getLabelRegions(str, regs);
1566 
1567  }
1568 
1569 }
loopMax
const EGS_I64 loopMax
The maximum number of iterations for near-infinite loops.
Definition: egs_functions.h:95
EGS_BaseGeometry::regions
int regions() const
Returns the number of local regions in this geometry.
Definition: egs_base_geometry.h:267
EGS_BaseGeometry::deref
int deref()
Decrease the reference count to this geometry.
Definition: egs_base_geometry.h:690
EGS_BaseGeometry::printInfo
virtual void printInfo() const
Print information about this geometry.
Definition: egs_base_geometry.cpp:726
EGS_BaseGeometry::bfactor
EGS_Float * bfactor
Array with B field scaling factors.
Definition: egs_base_geometry.h:785
EGS_NDGeometry::N
int N
Number of dimensions.
Definition: egs_nd_geometry.h:568
EGS_AffineTransform
A class providing affine transformations.
Definition: egs_transformations.h:417
EGS_BaseGeometry::name
string name
Name of this geometry.
Definition: egs_base_geometry.h:748
egs_input.h
EGS_Input class header file.
EGS_DeformedXYZ::setDeformations
int setDeformations(const char *defFile)
Definition: egs_nd_geometry.cpp:828
EGS_Vector
A class representing 3D vectors.
Definition: egs_vector.h:56
EGS_BaseGeometry::setLabels
int setLabels(EGS_Input *input)
Set the labels from an input block.
Definition: egs_base_geometry.cpp:1098
egs_functions.h
Global egspp functions header file.
EGS_NDGeometry::n
int * n
Used for calculating region indeces.
Definition: egs_nd_geometry.h:570
EGS_Input::print
void print(int nind, ostream &)
Used for debugging purposes.
Definition: egs_input.cpp:1174
EGS_BaseGeometry::med
int med
Medium index.
Definition: egs_base_geometry.h:764
createGeometry
EGS_GLIB_EXPORT EGS_BaseGeometry * createGeometry(EGS_Input *input)
Definition: egs_glib.cpp:57
egsSwapBytes
void egsSwapBytes(int *n)
Swap the bytes of 32 bit integers.
Definition: egs_functions.cpp:180
EGS_BaseGeometry
Base geometry class. Every geometry class must be derived from EGS_BaseGeometry.
Definition: egs_base_geometry.h:93
EGS_BaseGeometry::setRelativeRho
virtual void setRelativeRho(int start, int end, EGS_Float rho)
Set the relative mass density in regions.
Definition: egs_base_geometry.cpp:813
EGS_BaseGeometry::setBoundaryTolerance
void setBoundaryTolerance(EGS_Input *inp)
Set the value of the boundary tolerance from the input inp.
Definition: egs_base_geometry.cpp:717
EGS_BaseGeometry::isConvex
bool isConvex() const
Is the geometry convex?
Definition: egs_base_geometry.h:124
egs_nd_geometry.h
N-dimensional geometries: header.
EGS_BaseGeometry::setBScaling
virtual void setBScaling(int start, int end, EGS_Float bf)
Set the B field scaling factor in regions.
Definition: egs_base_geometry.cpp:861
EGS_BaseGeometry::getGeometry
static EGS_BaseGeometry * getGeometry(const string &Name)
Get a pointer to the geometry named Name.
Definition: egs_base_geometry.cpp:468
egsFatal
EGS_InfoFunction EGS_EXPORT egsFatal
Always use this function for reporting fatal errors.
Definition: egs_functions.cpp:136
EGS_BaseGeometry::getLabelRegions
virtual void getLabelRegions(const string &str, vector< int > &regs)
Get the list of all regions labeled with str.
Definition: egs_base_geometry.cpp:1069
EGS_BaseGeometry::setMedia
void setMedia(EGS_Input *inp)
Set the media in the geometry from the input pointed to by inp.
Definition: egs_base_geometry.cpp:743
EGS_BaseGeometry::addMedium
static int addMedium(const string &medname)
Add a medium or get the index of an existing medium.
Definition: egs_base_geometry.cpp:488
egsInformation
EGS_InfoFunction EGS_EXPORT egsInformation
Always use this function for reporting the progress of a simulation and any other type of information...
Definition: egs_functions.cpp:134
EGS_PlanesT
A set of parallel planes.
Definition: egs_planes.h:166
EGS_BaseGeometry::is_convex
bool is_convex
Is this geometry convex?
Definition: egs_base_geometry.h:822
egsGetEndian
int egsGetEndian()
Get the endianess of the machine.
Definition: egs_functions.cpp:368
EGS_BaseGeometry::setMedium
void setMedium(const string &Name)
Set all regions to a medium with name Name.
Definition: egs_base_geometry.cpp:480
EGS_BaseGeometry::setName
void setName(EGS_Input *inp)
Set the name of the geometry from the input inp.
Definition: egs_base_geometry.cpp:609
egs_transformations.h
EGS_AffineTransform and EGS_RotationMatrix class header file.
EGS_NDGeometry::setMedia
void setMedia(EGS_Input *inp, int nmed, const int *med_ind)
Define media.
Definition: egs_nd_geometry.cpp:111
EGS_XYZGeometry
An XYZ-geometry.
Definition: egs_nd_geometry.h:740
EGS_XYZRepeater
A geometry repeated on a regular XYZ grid.
Definition: egs_nd_geometry.h:1698
EGS_ZProjector
A projector into the z-plane.
Definition: egs_projectors.h:203
EGS_AffineTransform::transform
void transform(EGS_Vector &v) const
Transforms the vector v.
Definition: egs_transformations.h:546
EGS_BaseGeometry::region_media
short * region_media
Array of media indeces.
Definition: egs_base_geometry.h:756
EGS_YProjector
A projector into the y-plane.
Definition: egs_projectors.h:161
EGS_AffineTransform::getTransformation
static EGS_AffineTransform * getTransformation(EGS_Input *inp)
Constructs an affine transformation object from the input pointed to by inp and returns a pointer to ...
Definition: egs_transformations.cpp:43
EGS_BaseGeometry::has_rho_scaling
bool has_rho_scaling
Does this geometry have relative mass density scvaling?
Definition: egs_base_geometry.h:770
EGS_BaseGeometry::ref
int ref()
Increase the reference count to this geometry.
Definition: egs_base_geometry.h:678
EGS_XProjector
A projector into the x-plane.
Definition: egs_projectors.h:119
EGS_BaseGeometry::createSingleGeometry
static EGS_BaseGeometry * createSingleGeometry(EGS_Input *inp)
Create a single geometry from the input inp.
Definition: egs_base_geometry.cpp:543
epsilon
const EGS_Float epsilon
The epsilon constant for floating point comparisons.
Definition: egs_functions.h:61
EGS_NDGeometry
A class modeling a N-dimensional geometry.
Definition: egs_nd_geometry.h:312
EGS_Input
A class for storing information in a tree-like structure of key-value pairs. This class is used throu...
Definition: egs_input.h:182
EGS_BaseGeometry::getMediumName
static const char * getMediumName(int ind)
Get the name of medium with index ind.
Definition: egs_base_geometry.cpp:531
EGS_BaseGeometry::getName
const string & getName() const
Get the name of this geometry.
Definition: egs_base_geometry.h:496
EGS_Input::compare
static bool compare(const string &s1, const string &s2)
Definition: egs_input.cpp:1170
EGS_NDGeometry::g
EGS_BaseGeometry ** g
The dimensions.
Definition: egs_nd_geometry.h:569
EGS_BaseGeometry::nreg
int nreg
Number of local regions in this geometry.
Definition: egs_base_geometry.h:741
EGS_Input::takeInputItem
EGS_Input * takeInputItem(const string &key, bool self=true)
Get the property named key.
Definition: egs_input.cpp:226
EGS_NDGeometry::EGS_NDGeometry
EGS_NDGeometry(int ng, EGS_BaseGeometry **G, const string &Name="", bool O=true)
Definition: egs_nd_geometry.cpp:55
EGS_BaseGeometry::rhor
EGS_Float * rhor
Array with relative mass densities.
Definition: egs_base_geometry.h:775
EGS_Input::getInput
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_NDGeometry::type
static string type
The geometry type.
Definition: egs_nd_geometry.h:571
EGS_BaseGeometry::has_B_scaling
bool has_B_scaling
Does this geometry has B field scaling factor?
Definition: egs_base_geometry.h:780
egsWarning
EGS_InfoFunction EGS_EXPORT egsWarning
Always use this function for reporting warnings.
Definition: egs_functions.cpp:135
Generated by   doxygen 1.8.5