60 #define S_STREAM std::istringstream
63 #define S_STREAM std::istrstream
70 EGS_VoxelInfo *EGS_VoxelGeometry::v = 0;
71 int EGS_VoxelGeometry::nv = 0;
73 VHP_OrganData::VHP_OrganData(
const char *fname,
int slice_min,
int slice_max)
74 : nslice(0), ok(false) {
75 ifstream data(fname,ios::binary);
77 egsWarning(
"VHP_OrganData: failed to open file %s for reading\n",
82 data.read(&endian,
sizeof(endian));
84 if (endian != my_endian) {
85 egsWarning(
"VHP_OrganData: wrong endianess. \n"
86 "This machine has endianess %d, data was written on %d\n"
87 "This requires byte swaping, which has not been implemented yet\n",
88 (
int)my_endian,(
int)endian);
91 data.read((
char *)&dx,
sizeof(dx));
92 data.read((
char *)&dy,
sizeof(dy));
93 data.read((
char *)&dz,
sizeof(dz));
94 data.read((
char *)&nslice,
sizeof(nslice));
98 if (slice_max > nslice) {
101 int nwant = slice_max - slice_min;
102 slices =
new VHP_SliceInfo [nwant];
105 for (
int islice=0; islice<nslice; ++islice) {
106 if (islice >= slice_min) {
107 slices[is].readData(data);
109 slices[is++].getSliceDimensions(ni,nj);
121 if (islice == slice_max - 1) {
126 nslice = slice_max - slice_min;
128 egsWarning(
"VHP_OrganData: I/O error while reading data\n");
135 VHP_OrganData::~VHP_OrganData() {
141 void VHP_RowInfo::readData(istream &in) {
143 in.read((
char *)&n,
sizeof(n));
156 pix =
new unsigned short [nbin+1];
157 org =
new unsigned char [nbin];
159 in.read((
char *)pix,(nbin+1)*
sizeof(
unsigned short));
160 in.read((
char *)org,nbin*
sizeof(
unsigned char));
163 void VHP_SliceInfo::readData(istream &in) {
168 in.read((
char *)&firstr,
sizeof(firstr));
169 in.read((
char *)&lastr,
sizeof(lastr));
170 int nrow = lastr-firstr+1;
172 rinfo =
new VHP_RowInfo [nrow];
173 for (
int j=0; j<nrow; ++j) {
174 rinfo[j].readData(in);
180 EGS_VHPGeometry::EGS_VHPGeometry(
const char *phantom_file,
181 const char *media_file,
int slice_min,
int slice_max,
182 const string &Name) :
184 organs =
new VHP_OrganData(phantom_file,slice_min,slice_max);
185 if (!organs->isOK()) {
186 egsWarning(
"EGS_VHPGeometry: failed to construct organ data\n");
191 ifstream med_data(media_file);
193 egsWarning(
"EGS_VHPGeometry: failed to open media data file %s\n",
202 egsWarning(
"EGS_VHPGeometry: %d organs in media data file %s?\n",
209 for (j=0; j<256; ++j) {
211 organ_names[j] =
"Undefined";
215 med_data.getline(buf,1023);
216 for (j=0; j<norg; ++j) {
219 med_data >> iorg >> imed >> c;
220 med_data.getline(buf,1023);
221 if (iorg < 0 || iorg > 255)
222 egsWarning(
"EGS_VHPGeometry: bad organ number %d on line %d\n",
225 organ_media[iorg] = imed;
227 organ_names[iorg] = c;
228 organ_names[iorg] += buf;
232 med_data.getline(buf,1023);
233 vector<string> media_names;
234 vector<int> media_indeces;
235 for (j=0; j<nmed; ++j) {
236 med_data.getline(buf,1023);
238 media_names.push_back(
med);
242 for (j=0; j<nmed; ++j) {
244 media_names[j].c_str(),media_indeces[j]);
246 for (j=0; j<256; ++j) {
247 if (organ_media[j] >= 0) {
248 if (organ_media[j] >= nmed) {
249 egsWarning(
"EGS_VHPGeometry: invalid medium index %d\n",
253 organ_media[j] = media_indeces[organ_media[j]];
257 EGS_Float dx, dy, dz;
259 organs->getVoxelSizes(dx,dy,dz);
260 organs->getPhantomDimensions(nx,ny,nz);
261 vg =
new EGS_VoxelGeometry(nx,ny,nz,dx,dy,dz);
264 egsInformation(
"Phantom size in voxels: x=%d y=%d z=%d\n",nx,ny,nz);
268 void EGS_VHPGeometry::setMicros(
EGS_Input *input) {
269 if (!vg || !organs) {
270 egsWarning(
"EGS_VHPGeometry::setMicros: called for an invalid "
275 egsWarning(
"EGS_VHPGeometry::setMicros: no micro matrix definitions\n");
278 vector<EGS_MicroMatrixCluster *> mv;
280 input->
getInput(
"micro matrix folder",dir);
283 const static char *err_msg =
284 "EGS_VHPGeometry::setMicros: wrong/missing '%s' input\n";
285 int err = input->
getInput(
"BSC thickness",t_bsc);
286 if (err || t_bsc < 0) {
290 string tb_medium, bm_medium;
291 err = input->
getInput(
"TB medium",tb_medium);
296 err = input->
getInput(
"BM medium",bm_medium);
302 egsWarning(
"EGS_VHPGeometry::setMicros: wrong/missing inputs found\n"
303 " => no micro matrices set for use\n");
307 vector<string> minput;
308 while (!input->
getInput(
"micro matrix",minput)) {
310 if (minput.size() < 2) {
311 egsWarning(
"EGS_VHPGeometry::setMicros: 'micro matrix' input"
312 " requires at least two inputs => ignored\n");
315 string fname = dir.size() ?
egsJoinPath(dir,minput[0]) : minput[0];
316 EGS_MicroMatrixCluster *mcluster =
new EGS_MicroMatrixCluster(
317 vg->dx,vg->dy,vg->dz,t_bsc,tb_med,bm_med,fname.c_str());
318 if (!mcluster->isValid()) {
319 egsWarning(
"failed to construct micro cluster from %s\n",
326 for (
int i=1; i<minput.size(); ++i) {
327 S_STREAM s(minput[i].c_str());
331 egsWarning(
"EGS_VHPGeometry::setMicros: invalid input %s "
332 "for micro matrix %s\n",minput[i].c_str(),
337 if (org < 0 || org > 255) {
338 egsWarning(
"EGS_VHPGeometry::setMicros: invalid input"
339 " %d for micro matrix %s\n",org,minput[0].c_str());
348 mv.push_back(mcluster);
350 " organs:",fname.c_str());
351 for (
int i=0; i<orgs.size(); ++i) {
353 organ_micro[orgs[i]] = mv.size()-1;
364 for (
int i=0; i<nmicro; ++i) {
370 micros =
new EGS_MicroMatrixCluster* [nmicro];
372 for (
int i=0; i<nmicro; ++i) {
374 int nr = micros[i]->nxy*micros[i]->nz;
375 int nm = micros[i]->mxy*micros[i]->mz;
380 nmacro = vg->nxy*vg->nz;
384 if (ntot > 2147483647)
egsFatal(
"EGS_VHPGeometry::setMicros: "
385 "too many micro regions\n");
386 nreg = nmacro + nmax*nmicro;
389 egsInformation(
"Using regions %d...%d to identify micro voxels\n",
394 EGS_VHPGeometry::~EGS_VHPGeometry() {
402 for (
int j=0; j<nmicro; ++j) {
410 void EGS_VHPGeometry::setMedia(
EGS_Input *,
int,
const int *) {
411 egsWarning(
"EGS_VHPGeometry::setMedia: don't use this method.\n"
412 " Media are fixed by the organ data defining an VHP geometry\n");
415 void EGS_VHPGeometry::printInfo()
const {
426 string EGS_VHPGeometry::type =
"EGS_VHPGeometry";
429 VHPBox EGS_MicroMatrixCluster::bm_boxes[64];
430 bool EGS_MicroMatrixCluster::bm_boxes_initialized =
false;
432 EGS_MicroMatrixCluster::EGS_MicroMatrixCluster(
433 EGS_Float Dx, EGS_Float Dy, EGS_Float Dz, EGS_Float bsc_thickness,
434 int tb_med,
int bm_med,
const char *micro_file) : med_tb(tb_med), med_bm(bm_med),
439 ifstream in(micro_file,ios::binary);
441 egsWarning(
"EGS_MicroMatrixCluster: failed to open file %s\n",
449 unsigned char Mx, My, Mz;
450 in.read((
char *)&Mx,
sizeof(Mx));
451 in.read((
char *)&My,
sizeof(My));
452 in.read((
char *)&Mz,
sizeof(Mz));
458 micros =
new unsigned char *[nmic];
459 in.read((
char *)&Mx,
sizeof(Mx));
460 in.read((
char *)&My,
sizeof(My));
461 in.read((
char *)&Mz,
sizeof(Mz));
466 int nxyz = nxy*nz, j;
469 for (
int j=0; j<nmic; ++j) {
470 micros[j] =
new unsigned char [nxyz];
471 in.read((
char *)micros[j],nxyz*
sizeof(
unsigned char));
473 egsWarning(
"EGS_MicroMatrixCluster: I/O error while reading "
474 "%d'th micto-matrix from file %s\n",j+1,micro_file);
475 for (
int i=0; i<=j; ++i) {
494 EGS_Float ddx = Dx/nx, ddy = Dy/ny, ddz = Dz/nz;
503 vg =
new EGS_VoxelGeometry(nx,ny,nz,ddx,ddy,ddz);
509 if (!bm_boxes_initialized) {
510 bm_boxes_initialized =
true;
511 for (
int ixm=0; ixm<2; ++ixm) {
512 EGS_Float xmin = bsc_thickness*ixm;
513 for (
int ixp=0; ixp<2; ++ixp) {
514 EGS_Float xmax = ddx - bsc_thickness*ixp;
515 for (
int iym=0; iym<2; ++iym) {
516 EGS_Float ymin = bsc_thickness*iym;
517 for (
int iyp=0; iyp<2; ++iyp) {
518 EGS_Float ymax = ddy - bsc_thickness*iyp;
519 for (
int izm=0; izm<2; ++izm) {
520 EGS_Float zmin = bsc_thickness*izm;
521 for (
int izp=0; izp<2; ++izp) {
522 EGS_Float zmax = ddz - bsc_thickness*izp;
523 int ibox = ixm+2*ixp+4*iym+8*iyp+16*izm+32*izp;
524 bm_boxes[ibox].xmin = xmin;
525 bm_boxes[ibox].xmax = xmax;
526 bm_boxes[ibox].ymin = ymin;
527 bm_boxes[ibox].ymax = ymax;
528 bm_boxes[ibox].zmin = zmin;
529 bm_boxes[ibox].zmax = zmax;
547 v_tb =
new EGS_Float [nmic];
548 v_bm =
new EGS_Float [nmic];
549 v_bsc =
new EGS_Float [nmic];
550 double tot_tb = 0, tot_bm = 0, tot_bsc = 0;
551 for (j=0; j<nmic; ++j) {
553 int iy = (j-iz*mxy)/mx;
554 int ix = j-iz*mxy-iy*mx;
579 EGS_Float bsc_vol = 0, bm_vol = 0, tb_vol = 0;
580 for (
int jx=0; jx<nx; ++jx) {
581 for (
int jy=0; jy<ny; ++jy) {
582 for (
int jz=0; jz<nz; ++jz) {
583 int ireg = jx + jy*nx + jz*nxy;
584 if (micros[j][ireg]) {
585 int micxm = !(jx > 0 ? micros[j][ireg-1] :
586 micros[ixm+iy*mx+iz*mxy][nx-1+jy*nx+jz*nxy]);
587 int micxp = !(jx < nx-1 ? micros[j][ireg+1] :
588 micros[ixp+iy*mx+iz*mxy][jy*nx+jz*nxy]);
589 int micym = !(jy > 0 ? micros[j][ireg-nx] :
590 micros[ix+iym*mx+iz*mxy][jx+(ny-1)*nx+jz*nxy]);
591 int micyp = !(jy < ny-1 ? micros[j][ireg+nx] :
592 micros[ix+iyp*mx+iz*mxy][jx+jz*nxy]);
593 int miczm = !(jz > 0 ? micros[j][ireg-nxy] :
594 micros[ix+iy*mx+izm*mxy][jx+jy*nx+(nz-1)*nxy]);
595 int miczp = !(jz < nz-1 ? micros[j][ireg+nxy] :
596 micros[ix+iy*mx+izp*mxy][jx+jy*nx]);
597 if ((micxm && micxp && 2*bsc_thickness>ddx) ||
598 (micym && micyp && 2*bsc_thickness>ddy) ||
599 (miczm && miczp && 2*bsc_thickness>ddz)) {
604 int ibox=micxm+2*micxp+4*micym+8*micyp+
606 micros[j][ireg] = ibox + 2;
612 (bm_boxes[ibox].xmax-bm_boxes[ibox].xmin)*
613 (bm_boxes[ibox].ymax-bm_boxes[ibox].ymin)*
614 (bm_boxes[ibox].zmax-bm_boxes[ibox].zmin);
632 v_tb[j] = tb_vol*ddx*ddy*ddz;
633 v_bm[j] = bm_vol*ddx*ddy*ddz;
634 v_bsc[j] = bsc_vol*ddx*ddy*ddz;
645 EGS_MicroMatrixCluster::~EGS_MicroMatrixCluster() {
649 EGS_TestMicro::EGS_TestMicro(EGS_Float Dx, EGS_Float Dy, EGS_Float Dz,
650 EGS_Float bsc_t,
int tb_med,
int bm_med,
const char *micro_file,
652 micro =
new EGS_MicroMatrixCluster(Dx,Dy,Dz,bsc_t,tb_med,bm_med,micro_file);
653 vg =
new EGS_VoxelGeometry(micro->mx,micro->my,micro->mz,Dx,Dy,Dz);
654 nr = micro->nx*micro->ny*micro->nz;
655 nreg = nr*micro->mx*micro->my*micro->mz;
658 void EGS_TestMicro::setMedia(
EGS_Input *,
int,
const int *) {
659 egsFatal(
"EGS_TestMicro::setMedia(EGS_Input*,int,const int *):\n"
660 " Don't use this method. Media are only set via the micro matrix"
664 EGS_TestMicro::~EGS_TestMicro() {
667 string EGS_TestMicro::type =
"EGS_TestMicro";
669 const static EGS_VHP_LOCAL
char *vhp_error_msg1 =
670 "createGeometry(VHP): wrong/missing %s input for micro_test geometry\n";
676 int err = input->
getInput(
"type",type);
677 if (!err && type ==
"micro_test") {
679 vector<EGS_Float> vs;
680 err = input->
getInput(
"voxel sizes",vs);
681 if (err || vs.size() != 3) {
686 err = input->
getInput(
"BSC thickness",bsc_t);
691 string tb_medium, bm_medium;
692 int err1 = input->
getInput(
"TB medium",tb_medium);
693 int err2 = input->
getInput(
"BM medium",bm_medium);
704 err = input->
getInput(
"micro data file",micro_data);
706 egsWarning(vhp_error_msg1,
"'micro data file'");
715 tb_med,bm_med,micro_data.c_str());
721 string phantom, media;
722 int err1 = input->
getInput(
"phantom data",phantom);
723 int err2 = input->
getInput(
"media data",media);
724 if (err1)
egsWarning(
"createGeometry(EGS_VHPGeometry): "
725 "missing 'phantom data' input\n");
726 if (err2)
egsWarning(
"createGeometry(EGS_VHPGeometry): "
727 "missing 'media data' input\n");
732 err1 = input->
getInput(
"slice range",srange);
734 if (!err1 && srange.size() == 2 && srange[1] > srange[0]) result =
735 new EGS_VHPGeometry(phantom.c_str(),media.c_str(),srange[0],srange[1]);
739 if (!result->isOK()) {
740 egsWarning(
"createGeometry(EGS_VHPGeometry): failed to construct "
747 result->setMicros(input);
virtual void printInfo() const
Print information about this geometry.
string egsJoinPath(const string &first, const string &second)
Join two path variables (or a path and a file name) using the platform specific directory separator a...
int setLabels(EGS_Input *input)
Set the labels from an input block.
Global egspp functions header file.
Voxelized Human Phantom (VHP) geometry: header.
EGS_GLIB_EXPORT EGS_BaseGeometry * createGeometry(EGS_Input *input)
Base geometry class. Every geometry class must be derived from EGS_BaseGeometry.
void setBoundaryTolerance(EGS_Input *inp)
Set the value of the boundary tolerance from the input inp.
EGS_InfoFunction EGS_EXPORT egsFatal
Always use this function for reporting fatal errors.
static int addMedium(const string &medname)
Add a medium or get the index of an existing medium.
A Voxelized Human Phantom (VHP) geometry.
EGS_InfoFunction EGS_EXPORT egsInformation
Always use this function for reporting the progress of a simulation and any other type of information...
int egsGetEndian()
Get the endianess of the machine.
void setName(EGS_Input *inp)
Set the name of the geometry from the input inp.
int nreg
Number of local regions in this geometry.
EGS_InfoFunction EGS_EXPORT egsWarning
Always use this function for reporting warnings.