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