EGSnrc C++ class library  Report PIRS-898 (2021)
Iwan Kawrakow, Ernesto Mainegra-Hing, Frederic Tessier, Reid Townson and Blake Walters
egs_voxelized_shape.cpp
Go to the documentation of this file.
1 /*
2 ###############################################################################
3 #
4 # EGSnrc egs++ voxelized shape
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, 2009
25 #
26 # Contributors: Frederic Tessier
27 # Reid Townson
28 #
29 ###############################################################################
30 */
31 
32 
38 #include "egs_voxelized_shape.h"
39 #include "egs_input.h"
40 #include "egs_functions.h"
41 
42 #include <fstream>
43 using namespace std;
44 
45 void EGS_VoxelizedShape::EGS_VoxelizedShapeFormat0(const char *fname,
46  const string &Name,EGS_ObjectFactory *f) {
47  prob=0;
48  xpos=0;
49  ypos=0;
50  zpos=0;
51  map=0;
52  nx=0;
53  ny=0;
54  nz=0;
55  nxy=0;
56  nreg=0;
57  type=-1;
58  const static char *func = "EGS_VoxelizedShape::EGS_VoxelizedShape";
59  otype = "voxelized_shape";
60  if (!fname) {
61  return;
62  }
63  ifstream data(fname,ios::binary);
64  if (!data) {
65  egsWarning("%s: failed to open file %s\n",func,fname);
66  return;
67  }
68  char endian, form;
69  data.read(&endian,1);
70  data.read(&form,1);
71  if (data.fail()) {
72  egsWarning("%s: failed to read endianess and format from %s\n",func,fname);
73  return;
74  }
75  if (form < 0 || form > 1) {
76  egsWarning("%s: unknwon format %d found in file %s\n",func,(int)form,fname);
77  return;
78  }
79  if (endian != egsGetEndian()) {
80  egsWarning("%s: data in %s from a machine with different endianess.\n"
81  " Byte swaping not implemented yet.\n",func,fname);
82  return;
83  }
84  short snx, sny, snz;
85  data.read((char *)&snx,sizeof(short));
86  data.read((char *)&sny,sizeof(short));
87  data.read((char *)&snz,sizeof(short));
88  if (data.fail()) {
89  egsWarning("%s: failed to read number of voxels from %s\n",func,fname);
90  return;
91  }
92  if (snx < 1 || snx > 10000 || sny < 1 || sny > 10000 || snz < 1 || snz > 10000) {
93  egsWarning("%s: number of voxels seems strange: nx=%d, ny=%d, nz=%d\n",
94  func,snx,sny,snz);
95  return;
96  }
97  nx = snx;
98  ny = sny;
99  nz = snz;
100  nxy = nx*ny;
101  nreg = nxy*nz;
102  egsInformation("Distribution has %d x %d x %d voxels\n",nx,ny,nz);
103  float *x = new float [nx+1], *y = new float [ny+1], *z = new float [nz+1];
104  data.read((char *)x, (nx+1)*sizeof(float));
105  data.read((char *)y, (ny+1)*sizeof(float));
106  data.read((char *)z, (nz+1)*sizeof(float));
107  if (data.fail()) {
108  egsWarning("%s: failed to read voxel boundaries\n",func);
109  delete [] x;
110  delete [] y;
111  delete [] z;
112  return;
113  }
114  int nmap;
115  float *p;
116  if (form == 0) {
117  p = new float [nreg];
118  egsInformation("Format 0, reading %d values\n",nreg);
119  data.read((char *)p, nreg*sizeof(float));
120  if (data.fail()) {
121  egsWarning("%s: failed to read probabilities\n",func);
122  delete [] p;
123  delete [] x;
124  delete [] y;
125  delete [] z;
126  return;
127  }
128  nmap = nreg;
129  }
130  else {
131  egsInformation("Format 1, reading data\n");
132  data.read((char *)&nmap, sizeof(int));
133  if (data.fail() || nmap < 1) {
134  egsWarning("%s: failed to read nmap\n",func);
135  delete [] x;
136  delete [] y;
137  delete [] z;
138  return;
139  }
140  map = new int [nmap];
141  p = new float [nmap];
142  data.read((char *)map, nmap*sizeof(int));
143  data.read((char *)p, nmap*sizeof(float));
144  if (data.fail()) {
145  egsWarning("%s: failed to read probabilities and bin numbers\n",func);
146  delete [] p;
147  delete [] map;
148  delete [] x;
149  delete [] y;
150  delete [] z;
151  return;
152  }
153  }
154  EGS_Float *p1 = new EGS_Float [nmap];
155  int j;
156  for (j=0; j<nmap; ++j) {
157  p1[j] = p[j];
158  }
159  delete [] p;
160  egsInformation("Making alias table\n");
161  prob = new EGS_SimpleAliasTable(nmap,p1);
162  xpos = new EGS_Float [nx+1];
163  ypos = new EGS_Float [ny+1];
164  zpos = new EGS_Float [nz+1];
165  for (j=0; j<=nx; ++j) {
166  xpos[j] = x[j];
167  }
168  for (j=0; j<=ny; ++j) {
169  ypos[j] = y[j];
170  }
171  for (j=0; j<=nz; ++j) {
172  zpos[j] = z[j];
173  }
174  delete [] x;
175  delete [] y;
176  delete [] z;
177  type = form;
178  egsInformation("Done\n");
179 }
180 
181 EGS_VoxelizedShape::EGS_VoxelizedShape(int file_format, const char *fname,
182  const string &Name,EGS_ObjectFactory *f) : EGS_BaseShape(Name,f),
183  prob(0), xpos(0), ypos(0), zpos(0), map(0), nx(0), ny(0), nz(0), nxy(0),
184  nreg(0), type(-1) {
185  const static char *func = "EGS_VoxelizedShape::EGS_VoxelizedShape";
186  if (file_format == 0) { // binary file format -> call original constructor
187  string s(fname);
188  EGS_VoxelizedShapeFormat0(s.c_str());
189  return;
190  }
191  if (file_format != 1) { // unknown file format
192  egsWarning("%s: unknown file format = %d\n",func,file_format);
193  return;
194  }
195  // interfile format -> continue
196  otype = "voxelized_shape";
197  if (!fname) {
198  return;
199  }
200  ifstream h_file(fname);
201  if (!h_file) {
202  egsWarning("%s: failed to open file %s\n",func,fname);
203  return;
204  }
205  string data_file("");
206  int data_type = -1;
207  int Nx, Ny, Nz;
208  float scale_x=0.0, scale_y=0.0, scale_z=0.0;
209  for (EGS_I64 loopCount=0; loopCount<=loopMax; ++loopCount) {
210  if (loopCount == loopMax) {
211  egsFatal("EGS_VoxelizedShape::EGS_VoxelizedShape: Too many iterations were required! Input may be invalid, or consider increasing loopMax.");
212  return;
213  }
214  string line, key, value;
215  size_t pos;
216  getline(h_file, line);
217  if (h_file.eof() || h_file.fail() || !h_file.good()) {
218  break;
219  }
220  pos = line.find(":=");
221  if (pos == string::npos) {
222  continue;
223  }
224  key = line.substr(0, int(pos));
225  value = line.substr(int(pos)+2);
226  while ((key[0] == '!') || (key[0] == ' ')) {
227  key.erase(0,1);
228  }
229  while (key[key.length()-1] == ' ') {
230  key.erase(key.length()-1, 1);
231  }
232  while (value[0] == ' ') {
233  value.erase(0,1);
234  }
235  while (value[value.length()-1] == ' ') {
236  value.erase(value.length()-1, 1);
237  }
238  if (key == "matrix size [1]") {
239  sscanf(value.c_str(), "%u", &Nx);
240  }
241  else if (key == "matrix size [2]") {
242  sscanf(value.c_str(), "%u", &Ny);
243  }
244  else if ((key == "number of slices") || (key == "number of images")) {
245  sscanf(value.c_str(), "%u", &Nz);
246  }
247  else if (key == "scaling factor (mm/pixel) [1]") {
248  sscanf(value.c_str(), "%f", &scale_x);
249  }
250  else if (key == "scaling factor (mm/pixel) [2]") {
251  sscanf(value.c_str(), "%f", &scale_y);
252  }
253  else if (key == "slice thickness (pixels)") {
254  sscanf(value.c_str(), "%f", &scale_z); // ????????????????
255  }
256  else if (key == "name of data file") {
257  data_file = value;
258  }
259  else if (key == "number format") {
260  if ((value == "float") || (value == "FLOAT")) {
261  data_type = 0;
262  }
263  else if ((value == "unsigned integer") || (value == "UNSIGNED INTEGER")) {
264  data_type = 1;
265  }
266  else if ((value == "signed integer") || (value == "SIGNED INTEGER")) {
267  data_type = 2;
268  }
269  else {
270  egsWarning("%s: unrecognised 'number format' type: %s \n",func,value.c_str());
271  }
272  }
273  }
274  if (Nx < 1 || Ny < 1 || Nz < 1 || scale_x <= 0.0 || scale_y <= 0.0 || scale_z <= 0.0 || data_file == "" || data_type == -1) {
275  egsWarning("%s: invalid interfile header information: "
276  "Nx=%d Ny=%d Nz=%d scale_x=%f scale_y=%f scale_z=%f "
277  "data_file='%s' number_format=%d\n",func,Nx,Ny,Nz,scale_x,scale_y,scale_z,data_file.c_str(),data_type);
278  return;
279  }
280  nx = Nx;
281  ny = Ny;
282  nz = Nz;
283  nxy = nx*ny;
284  nreg = nxy*nz;
285  egsInformation("Distribution has %d x %d x %d voxels\n",nx,ny,nz);
286  float *x = new float [nx+1], *y = new float [ny+1], *z = new float [nz+1];
287  scale_x /= 10.0;
288  scale_y /= 10.0;
289  scale_z /= 10.0; // convert to cm
290  {
291  int j;
292  for (j=0; j<=nx; j++) {
293  x[j] = -(float)nx*scale_x/2.0 + j*scale_x;
294  }
295  for (j=0; j<=ny; j++) {
296  y[j] = -(float)ny*scale_y/2.0 + j*scale_y;
297  }
298  for (j=0; j<=nz; j++) {
299  z[j] = -(float)nz*scale_z/2.0 + j*scale_z;
300  }
301  }
302  ifstream i_file(data_file.c_str(),ios::binary);
303  if (!i_file) {
304  egsWarning("%s: failed to open interfile data "
305  "%s\n",func,data_file.c_str());
306  return;
307  }
308  int nmap = nreg;
309  float *p = new float [nreg];
310  if (data_type == 0) {
311  i_file.read((char *)p, nreg*sizeof(float));
312  }
313  else if (data_type == 1) {
314  unsigned short int *p_tmp = new unsigned short int [nreg];
315  i_file.read((char *)p_tmp, nreg*sizeof(unsigned short int));
316  for (int cc = 0; cc<nreg; cc++) {
317  p[cc] = (float)(p_tmp[cc]);
318  }
319  delete [] p_tmp;
320  }
321  else {
322  short int *p_tmp = new short int [nreg];
323  i_file.read((char *)p_tmp, nreg*sizeof(short int));
324  for (int cc = 0; cc<nreg; cc++) {
325  p[cc] = (float)(p_tmp[cc]);
326  }
327  delete [] p_tmp;
328  }
329  EGS_Float *p1 = new EGS_Float [nmap];
330  int j;
331  for (j=0; j<nmap; ++j) {
332  p1[j] = p[j];
333  }
334  delete [] p;
335  egsInformation("Making alias table\n");
336  prob = new EGS_SimpleAliasTable(nmap,p1);
337  xpos = new EGS_Float [nx+1];
338  ypos = new EGS_Float [ny+1];
339  zpos = new EGS_Float [nz+1];
340  for (j=0; j<=nx; ++j) {
341  xpos[j] = x[j];
342  }
343  for (j=0; j<=ny; ++j) {
344  ypos[j] = y[j];
345  }
346  for (j=0; j<=nz; ++j) {
347  zpos[j] = z[j];
348  egsInformation("%f ", zpos[j]);
349  }
350  delete [] x;
351  delete [] y;
352  delete [] z;
353  type = 0;
354  egsInformation("Done\n");
355 }
356 
357 EGS_VoxelizedShape::~EGS_VoxelizedShape() {
358  if (isValid()) {
359  delete prob;
360  delete [] xpos;
361  delete [] ypos;
362  delete [] zpos;
363  if (type == 1) {
364  delete [] map;
365  }
366  }
367 }
368 
369 
370 extern "C" {
371 
372  EGS_VOXELIZED_SHAPE_EXPORT EGS_BaseShape *createShape(EGS_Input *input,
373  EGS_ObjectFactory *f) {
374  static const char *func = "createShape(voxelized shape)";
375  if (!input) {
376  egsWarning("%s: null input?\n",func);
377  return 0;
378  }
379  string fname;
380  int err = input->getInput("file name",fname);
381  int file_format;
382  int err2 = input->getInput("file format",file_format);
383  if (err) {
384  egsWarning("%s: missing 'file name' input\n",func);
385  return 0;
386  }
387  if (err2) {
388  egsInformation("%s: 'file format' input missing. Using default 'binary'"
389  "file format \n",func);
390  file_format = 0;
391  }
392  EGS_VoxelizedShape *shape = new EGS_VoxelizedShape(file_format, fname.c_str());
393  if (!shape->isValid()) {
394  delete shape;
395  return 0;
396  }
397  shape->setName(input);
398  shape->setTransformation(input);
399  return shape;
400  }
401 
402 }
Base shape class. All shapes in the EGSnrc C++ class library are derived from EGS_BaseShape.
Definition: egs_shapes.h:114
void setTransformation(EGS_Input *inp)
Set the transformation attached to this shape.
Definition: egs_shapes.cpp:69
A class for storing information in a tree-like structure of key-value pairs. This class is used throu...
Definition: egs_input.h:182
int 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
An object factory.
void setName(EGS_Input *inp)
Set the name of the object from the information provided by inp.
string otype
The object type.
A class for sampling random bins from a given probability distribution using the alias table techniqu...
A "voxelized shape".
EGS_Float * zpos
! The y-positions of the grid
EGS_Float * xpos
! The alias table for randomly picking voxels
int * map
! The z-positions of the grid
EGS_Float * ypos
! The x-positions of the grid
int nx
! Voxel map (for type=1)
EGS_VoxelizedShape(int file_format, const char *fname, const string &Name="", EGS_ObjectFactory *f=0)
Conctructor.
Global egspp functions header file.
EGS_Input class header file.
A "voxelized shape": header.
int egsGetEndian()
Get the endianess of the machine.
EGS_InfoFunction EGS_EXPORT egsInformation
Always use this function for reporting the progress of a simulation and any other type of information...
EGS_InfoFunction EGS_EXPORT egsFatal
Always use this function for reporting fatal errors.
const EGS_I64 loopMax
The maximum number of iterations for near-infinite loops.
Definition: egs_functions.h:96
EGS_InfoFunction EGS_EXPORT egsWarning
Always use this function for reporting warnings.