egs_nd_geometry.h

Go to the documentation of this file.

/*
###############################################################################
#
#  EGSnrc egs++ nd geometry headers
#  Copyright (C) 2015 National Research Council Canada
#
#  This file is part of EGSnrc.
#
#  EGSnrc is free software: you can redistribute it and/or modify it under
#  the terms of the GNU Affero General Public License as published by the
#  Free Software Foundation, either version 3 of the License, or (at your
#  option) any later version.
#
#  EGSnrc is distributed in the hope that it will be useful, but WITHOUT ANY
#  WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
#  FOR A PARTICULAR PURPOSE.  See the GNU Affero General Public License for
#  more details.
#
#  You should have received a copy of the GNU Affero General Public License
#  along with EGSnrc. If not, see <http://www.gnu.org/licenses/>.
#
###############################################################################
#
#  Author:          Iwan Kawrakow, 2005
#
#  Contributors:    Blake Walters
#                   Ernesto Mainegra-Hing
#                   Frederic Tessier
#                   Reid Townson
#                   Randle Taylor
#
###############################################################################
*/


#ifndef EGS_ND_GEOMETRY_
#define EGS_ND_GEOMETRY_

#ifdef WIN32

    #ifdef BUILD_NDG_DLL
        #define EGS_NDG_EXPORT __declspec(dllexport)
    #else
        #define EGS_NDG_EXPORT __declspec(dllimport)
    #endif
    #define EGS_NDG_LOCAL

#else

    #ifdef HAVE_VISIBILITY
        #define EGS_NDG_EXPORT __attribute__ ((visibility ("default")))
        #define EGS_NDG_LOCAL  __attribute__ ((visibility ("hidden")))
    #else
        #define EGS_NDG_EXPORT
        #define EGS_NDG_LOCAL
    #endif

#endif


#include "egs_base_geometry.h"

#include<vector>
#include <iomanip>
using std::vector;

class EGS_NDG_EXPORT EGS_NDGeometry : public EGS_BaseGeometry {

public:

    EGS_NDGeometry(int ng, EGS_BaseGeometry **G, const string &Name = "",
                   bool O=true);

    EGS_NDGeometry(vector<EGS_BaseGeometry *> &G, const string &Name = "",
                   bool O=true);

    ~EGS_NDGeometry();

    bool isInside(const EGS_Vector &x) {
        for (int j=0; j<N; j++) if (!g[j]->isInside(x)) {
                return false;
            }
        return true;
    };

    int isWhere(const EGS_Vector &x) {
        int ireg = 0;
        for (int j=0; j<N; j++) {
            int ij = g[j]->isWhere(x);
            if (ij < 0) {
                return -1;
            }
            ireg += ij*n[j];
        }
        return ireg;
    };

    int inside(const EGS_Vector &x) {
        return isWhere(x);
    };

    EGS_Float howfarToOutside(int ireg, const EGS_Vector &x,
                              const EGS_Vector &u) {
        if (ireg < 0) {
            return 0;
        }
        int itmp = ireg;
        EGS_Float d = veryFar;
        for (int j=N-1; j>=0; j--) {
            int l = itmp/n[j];
            EGS_Float t = g[j]->howfarToOutside(l,x,u);
            if (t <= 0) {
                return t;
            }
            if (t < d) {
                d = t;
            }
        }
        return d;
    };

    int howfar(int ireg, const EGS_Vector &x, const EGS_Vector &u,
               EGS_Float &t, int *newmed=0, EGS_Vector *normal=0) {
        if (ireg >= 0) {
            int itmp = ireg;
            int inext = -1, idelta, lnew_j=0;
            for (int j=N-1; j>=0; j--) {
                int l = itmp/n[j];
                int lnew = g[j]->howfar(l,x,u,t,0,normal);
                if (lnew != l) {
                    inext = j;
                    idelta = lnew - l;
                    lnew_j = lnew;
                }
                itmp -= l*n[j];
            }
            if (inext < 0) {
                return ireg;
            }
            int inew = lnew_j >= 0 ? ireg + idelta*n[inext] : lnew_j;
            //if( lnew_j < 0 ) return lnew_j;
            //int inew = ireg + idelta*n[inext];
            if (newmed) {
                *newmed = inew >= 0 ? medium(inew) : -1;
            }
            return inew;
        }
        for (int j=0; j<N; j++) {
            EGS_Float tj = t;
            bool is_in = g[j]->isInside(x);
            // Note: this logic fails if the geometry is not convex!
            if (!is_in) {
                EGS_Vector nn, *np=0;
                if (normal) {
                    np = &nn;
                }
                int i = g[j]->howfar(ireg,x,u,tj,0,np);
                if (i >= 0) {
                    EGS_Vector tmp(x + u*tj);
                    bool is_ok = true;
                    int res = 0;
                    for (int k=0; k<N; k++) {
                        if (k != j) {
                            int check = g[k]->isWhere(tmp);
                            if (check < 0) {
                                is_ok = false;
                                break;
                            }
                            res += check*n[k];
                        }
                    }
                    if (is_ok) {
                        t = tj;
                        res += i*n[j];
                        if (newmed) {
                            *newmed = medium(res);
                        }
                        if (normal) {
                            *normal = nn;
                        }
                        return res;
                    }
                }
            }
            // and so, to fix it, we do the following for non-convex
            // dimensions. But this only works if and only if there
            // is not more than 1 non-convex dimension and this
            // dimension is last in the list.
            else if (!g[j]->isConvex()) {
                EGS_Vector tmp(x);
                EGS_Float tleft = t;
                int ii = g[j]->isWhere(x);
                EGS_Float ttot = 0;
                for (EGS_I64 loopCount=0; loopCount<=loopMax; ++loopCount) {
                    if (loopCount == loopMax) {
                        egsFatal("EGS_NDGeometry::howfar: Too many iterations were required! Input may be invalid, or consider increasing loopMax.");
                        return -1;
                    }
                    EGS_Float tt = tleft;
                    int inew = g[j]->howfar(ii,tmp,u,tt);
                    if (inew == ii) {
                        break;
                    }
                    tleft -= tt;
                    tmp += u*tt;
                    ii = inew;
                    ttot += tt;
                    //egsWarning(" inew = %d tt = %g ttot = %g tmp = (%g,%g,%g)\n",inew,tt,ttot,tmp.x,tmp.y,tmp.z);
                    if (ii < 0) {
                        break;
                    }
                }
                //egsWarning("doing non-convex part: x = (%g,%g,%g) ii = %d t = %g ttot = %g\n",x.x,x.y,x.z,ii,t,ttot);
                if (ii < 0) {
                    EGS_Vector nn, *np = normal ? &nn : 0;
                    EGS_Float tt = tleft;
                    int inew = g[j]->howfar(ii,tmp,u,tt,0,np);
                    //egsWarning(" exited and tried to reenter: %d %g\n",inew,ttot+tt);
                    if (inew >= 0) {
                        tmp += u*tt;
                        bool is_ok = true;
                        int res = 0;
                        for (int k=0; k<N; k++) {
                            if (k != j) {
                                int check = g[k]->isWhere(tmp);
                                if (check < 0) {
                                    is_ok = false;
                                    break;
                                }
                                res += check*n[k];
                            }
                        }
                        if (is_ok) {
                            t = ttot + tt;
                            res += inew*n[j];
                            if (newmed) {
                                *newmed = medium(res);
                            }
                            if (normal) {
                                *normal = nn;
                            }
                            return res;
                        }
                    }
                }
            }
        }
        return -1;
    };

    EGS_Float hownear(int ireg, const EGS_Vector &x) {
        EGS_Float tmin = veryFar;
        if (ireg >= 0) {
            int itmp = ireg;
            for (int j=N-1; j>=0; j--) {
                int l = itmp/n[j];
                EGS_Float t = g[j]->hownear(l,x);
                if (t < tmin) {
                    tmin = t;
                    if (tmin <= 0) {
                        return 0;
                    }
                }
                itmp -= l*n[j];
            }
            return tmin;
        }
        if (ortho) {
            int nc = 0;
            EGS_Float s1=0, s2=0;
            for (int j=0; j<N; j++) {
                if (!g[j]->isInside(x)) {
                    EGS_Float t = g[j]->hownear(-1,x);
                    nc++;
                    s1 += t;
                    s2 += t*t;
                }
            }
            if (nc == 1) {
                return s1;
            }
            return sqrt(s2);
        }
        else {
            EGS_Float tmin = veryFar;
            for (int j=0; j<N; j++) {
                EGS_Float t;
                int i = g[j]->isWhere(x);
                t = g[j]->hownear(i,x);
                if (t < tmin) {
                    tmin = t;
                    if (tmin <= 0) {
                        return 0;
                    }
                }
            }
            return tmin;
        }
    };

    int getMaxStep() const {
        int nstep = 1;
        for (int j=0; j<N; ++j) {
            nstep += g[j]->getMaxStep();
        }
        return nstep;
    };

    const string &getType() const {
        return type;
    };

    void printInfo() const;

    virtual void getLabelRegions(const string &str, vector<int> &regs);
    void ndRegions(int r, int dim, int dimk, int k, vector<int> &regs);

protected:

    int              N;       
    EGS_BaseGeometry **g;     
    int              *n;      
    static string    type;    
    bool ortho;               

    void setup();

    void setMedia(EGS_Input *inp, int nmed, const int *med_ind);

private:

    void setM(int ib, int idim, const vector<int> &ranges, int medium);
};

#ifdef EXPLICIT_XYZ
#include "../egs_planes/egs_planes.h"

class EGS_NDG_EXPORT EGS_XYZGeometry : public EGS_BaseGeometry {

public:

    EGS_XYZGeometry(EGS_PlanesX *Xp, EGS_PlanesY *Yp, EGS_PlanesZ *Zp,
                    const string &Name = "");

    ~EGS_XYZGeometry();

    bool isInside(const EGS_Vector &x) {
        if (!xp->isInside(x)) {
            return false;
        }
        if (!yp->isInside(x)) {
            return false;
        }
        if (!zp->isInside(x)) {
            return false;
        }
        return true;
    };

    int isWhere(const EGS_Vector &x) {
        int ix = xp->isWhere(x);
        if (ix < 0) {
            return ix;
        }
        int iy = yp->isWhere(x);
        if (iy < 0) {
            return iy;
        }
        int iz = zp->isWhere(x);
        if (iz < 0) {
            return iz;
        }
        return ix + iy*nx + iz*nxy;
    };

    int inside(const EGS_Vector &x) {
        return isWhere(x);
    };

    int computeIntersections(int ireg, int n, const EGS_Vector &X,
                             const EGS_Vector &u, EGS_GeometryIntersections *isections) {
        if (n < 1) {
            return -1;
        }
        int ifirst = 0;
        EGS_Float t, t_ini = 0;
        EGS_Vector x(X);
        int imed;
        if (ireg < 0) {
            t = veryFar;
            ireg = howfar(ireg,x,u,t,&imed);
            if (ireg < 0) {
                return 0;
            }
            isections[0].t = t;
            isections[0].rhof = 1;
            isections[0].ireg = -1;
            isections[0].imed = -1;
            t_ini = t;
            ++ifirst;
            x += u*t;
        }
        else {
            imed = medium(ireg);
        }
        EGS_Float *px = xp->getPositions(),
                   *py = yp->getPositions(),
                    *pz = zp->getPositions();
        int iz = ireg/nxy;
        int ir = ireg - iz*nxy;
        int iy = ir/nx;
        int ix = ir - iy*nx;
        EGS_Float uxi, uyi, uzi, nextx, nexty, nextz;
        int dirx, icx, diry, icy, dirz, icz;
        if (u.x > 0)      {
            uxi = 1/u.x;
            dirx =  1;
            icx = 1;
        }
        else if (u.x < 0) {
            uxi = 1/u.x;
            dirx = -1;
            icx = 0;
        }
        else               {
            uxi = veryFar*1e3;
            dirx =  1;
            icx = 1;
        }
        if (u.y > 0)      {
            uyi = 1/u.y;
            diry =  1;
            icy = 1;
        }
        else if (u.y < 0) {
            uyi = 1/u.y;
            diry = -1;
            icy = 0;
        }
        else               {
            uyi = veryFar*1e3;
            diry =  1;
            icy = 1;
        }
        if (u.z > 0)      {
            uzi = 1/u.z;
            dirz =  1;
            icz = 1;
        }
        else if (u.z < 0) {
            uzi = 1/u.z;
            dirz = -1;
            icz = 0;
        }
        else               {
            uzi = veryFar*1e3;
            dirz =  1;
            icz = 1;
        }
        nextx = (px[ix+icx] - x.x)*uxi + t_ini;
        nexty = (py[iy+icy] - x.y)*uyi + t_ini;
        nextz = (pz[iz+icz] - x.z)*uzi + t_ini;
        for (int j=ifirst; j<n; j++) {
            isections[j].imed = imed;
            isections[j].ireg = ireg;
            isections[j].rhof = getRelativeRho(ireg);
            int inew;
            if (nextx < nexty && nextx < nextz) {
                t = nextx;
                ix += dirx;
                if (ix < 0 || ix >= nx) {
                    inew = -1;
                }
                else {
                    inew = ireg + dirx;
                    nextx = (px[ix+icx] - x.x)*uxi + t_ini;
                }
            }
            else if (nexty < nextz) {
                t = nexty;
                iy += diry;
                if (iy < 0 || iy >= ny) {
                    inew = -1;
                }
                else {
                    inew = ireg + nx*diry;
                    nexty = (py[iy+icy] - x.y)*uyi + t_ini;
                }
            }
            else {
                t = nextz;
                iz += dirz;
                if (iz < 0 || iz >= nz) {
                    inew = -1;
                }
                else {
                    inew = ireg + nxy*dirz;
                    nextz = (pz[iz+icz] - x.z)*uzi + t_ini;
                }
            }
            isections[j].t = t;
            if (inew < 0) {
                return j+1;
            }
            ireg = inew;
            imed = medium(ireg);
        }
        return ireg >= 0 ? -1 : n;
    };

    EGS_Float howfarToOutside(int ireg, const EGS_Vector &x,
                              const EGS_Vector &u) {
        if (ireg < 0) {
            return 0;
        }
        int iz = ireg/nxy;
        int ir = ireg - iz*nxy;
        int iy = ir/nx;
        int ix = ir - iy*nx;
        EGS_Float  d = xp->howfarToOutside(ix,x,u);
        EGS_Float ty = yp->howfarToOutside(iy,x,u);
        if (ty < d) {
            d = ty;
        }
        EGS_Float tz = zp->howfarToOutside(iz,x,u);
        if (tz < d) {
            d = tz;
        }
        return d;
    };

    int howfar(int ireg, const EGS_Vector &x, const EGS_Vector &u,
               EGS_Float &t, int *newmed=0, EGS_Vector *normal=0) {
        if (ireg >= 0) {
            int iz = ireg/nxy;
            int ir = ireg - iz*nxy;
            int iy = ir/nx;
            int ix = ir - iy*nx;
            int inew = ireg;
            if (u.x > 0) {
                EGS_Float d = (xpos[ix+1]-x.x)/u.x;
                if (d <= t) {
                    if (d <= boundaryTolerance) {
                        // t=0 works on most cases but can result in
                        // getting stuck on an edge, so use boundaryTolerance
                        t = halfBoundaryTolerance;
                    }
                    else {
                        t = d;
                    }
                    if ((++ix) < nx) {
                        inew = ireg + 1;
                    }
                    else {
                        inew = -1;
                    }
                    if (normal) {
                        *normal = EGS_Vector(-1,0,0);
                    }
                }
            }
            else if (u.x < 0) {
                EGS_Float d = (xpos[ix]-x.x)/u.x;
                if (d <= t) {
                    if (d <= boundaryTolerance) {
                        t = halfBoundaryTolerance;
                    }
                    else {
                        t = d;
                    }
                    if ((--ix) >= 0) {
                        inew = ireg - 1;
                    }
                    else {
                        inew = -1;
                    }
                    if (normal) {
                        *normal = EGS_Vector(1,0,0);
                    }
                }
            }
            if (u.y > 0) {
                EGS_Float d = (ypos[iy+1]-x.y)/u.y;
                if (d <= t) {
                    if (d <= boundaryTolerance) {
                        t = halfBoundaryTolerance;
                    }
                    else {
                        t = d;
                    }
                    if ((++iy) < ny) {
                        inew = ireg + nx;
                    }
                    else {
                        inew = -1;
                    }
                    if (normal) {
                        *normal = EGS_Vector(0,-1,0);
                    }
                }
            }
            else if (u.y < 0) {
                EGS_Float d = (ypos[iy]-x.y)/u.y;
                if (d <= t) {
                    if (d <= boundaryTolerance) {
                        t = halfBoundaryTolerance;
                    }
                    else {
                        t = d;
                    }
                    if ((--iy) >= 0) {
                        inew = ireg - nx;
                    }
                    else {
                        inew = -1;
                    }
                    if (normal) {
                        *normal = EGS_Vector(0,1,0);
                    }
                }
            }
            if (u.z > 0) {
                EGS_Float d = (zpos[iz+1]-x.z)/u.z;
                if (d <= t) {
                    if (d <= boundaryTolerance) {
                        t = halfBoundaryTolerance;
                    }
                    else {
                        t = d;
                    }
                    if ((++iz) < nz) {
                        inew = ireg+nxy;
                    }
                    else {
                        inew = -1;
                    }
                    if (normal) {
                        *normal = EGS_Vector(0,0,-1);
                    }
                }
            }
            else if (u.z < 0) {
                EGS_Float d = (zpos[iz]-x.z)/u.z;
                if (d <= t) {
                    if (d <= boundaryTolerance) {
                        t = halfBoundaryTolerance;
                    }
                    else {
                        t = d;
                    }
                    if ((--iz) >= 0) {
                        inew = ireg-nxy;
                    }
                    else {
                        inew = -1;
                    }
                    if (normal) {
                        *normal = EGS_Vector(0,0,1);
                    }
                }
            }
            if (newmed && inew >= 0) {
                *newmed = medium(inew);
            }
            return inew;
        }
        int face,ix,iy,iz;
        int inew = howfarFromOut(x,u,t,ix,iy,iz,face,normal);
        if (inew >= 0 && newmed) {
            *newmed = medium(inew);
        }
        return inew;
    };

    EGS_Float hownear(int ireg, const EGS_Vector &x) {
        if (ireg >= 0) {
            int iz = ireg/nxy;
            int ir = ireg - iz*nxy;
            int iy = ir/nx;
            int ix = ir - iy*nx;
            EGS_Float t = xp->hownear(ix,x);
            EGS_Float t1 = yp->hownear(iy,x);
            if (t1 < t) {
                t = t1;
            }
            t1 = zp->hownear(iz,x);
            if (t1 < t) {
                t = t1;
            }
            return t;
        }
        int nc = 0;
        EGS_Float s1=0, s2=0;
        if (!xp->isInside(x)) {
            EGS_Float t = xp->hownear(-1,x);
            nc++;
            s1 += t;
            s2 += t*t;
        }
        if (!yp->isInside(x)) {
            EGS_Float t = yp->hownear(-1,x);
            nc++;
            s1 += t;
            s2 += t*t;
        }
        if (!zp->isInside(x)) {
            EGS_Float t = zp->hownear(-1,x);
            nc++;
            s1 += t;
            s2 += t*t;
        }
        if (nc == 1) {
            return s1;
        }
        return sqrt(s2);
    };


    EGS_Float getVolume(int ireg) {
        int iz = ireg/nxy;
        int ir = ireg - iz*nxy;
        int iy = ir/nx;
        int ix = ir - iy*nx;
        return (xpos[ix+1]-xpos[ix])*(ypos[iy+1]-ypos[iy])*(zpos[iz+1]-zpos[iz]);
    }

    EGS_Float getBound(int idir, int ind) {
        EGS_Float bound;
        if (idir == 0) {
            if (ind>=0 && ind<=nx) {
                bound = xpos[ind];
            }
            else {
                egsWarning("Error in getBound: Looking for X bound out of range %d\n"
                           "Will set this to zero.\n",ind);
                bound = 0.0;
            }
        }
        else if (idir == 1) {
            if (ind>=0 && ind<=ny) {
                bound = ypos[ind];
            }
            else {
                egsWarning("Error in getBound: Looking for Y bound out of range %d\n"
                           "Will set this to zero.\n",ind);
                bound = 0.0;
            }
        }
        else if (idir == 2) {
            if (ind>=0 && ind<=nz) {
                bound = zpos[ind];
            }
            else {
                egsWarning("Error in getBound: Looking for Z bound out of range %d\n"
                           "Will set this to zero.\n",ind);
                bound = 0.0;
            }
        }
        else {
            egsWarning("Error in getBound: Dimension %d undefined in EGS_XYZGeometry.\n",idir);
            bound = 0.0;
        }
        return bound;
    }

    int getNRegDir(int idir) {
        if (idir==0) {
            return nx;
        }
        else if (idir==1) {
            return ny;
        }
        else if (idir==2) {
            return nz;
        }
        else {
            egsWarning("Error in getNregDir: Dimension %d undefined in EGS_XYZGeometry.\n",
                       idir);
            return 0;
        }
    }


    int getMaxStep() const {
        return nx+ny+nz+1;
    };

    int getNx() const {
        return nx;
    };
    int getNy() const {
        return ny;
    };
    int getNz() const {
        return nz;
    };

    EGS_Float *getXPositions() {
        return xp->getPositions();
    };
    EGS_Float *getYPositions() {
        return yp->getPositions();
    };
    EGS_Float *getZPositions() {
        return zp->getPositions();
    };

    static int getDigits(int i);

    const string &getType() const {
        return type;
    };

    void printInfo() const;

    static EGS_XYZGeometry *constructGeometry(const char *dens_or_egphant_file,
            const char *ramp_file, int dens_or_egphant=0);
    static EGS_XYZGeometry *constructCTGeometry(const char *dens_or_egphant_file);

    void voxelizeGeometry(EGS_Input *input);

    void setXYZLabels(EGS_Input *input);

    virtual void getLabelRegions(const string &str, vector<int> &regs);
    void getXLabelRegions(const string &str, vector<int> &regs) {
        xp->getLabelRegions(str, regs);
    }
    void getYLabelRegions(const string &str, vector<int> &regs) {
        yp->getLabelRegions(str, regs);
    }
    void getZLabelRegions(const string &str, vector<int> &regs) {
        zp->getLabelRegions(str, regs);
    }

protected:

    EGS_PlanesX      *xp;
    EGS_PlanesY      *yp;
    EGS_PlanesZ      *zp;
    EGS_Float        *xpos;
    EGS_Float        *ypos;
    EGS_Float        *zpos;
    int              nx, ny, nz, nxy;
    static string    type;

    void setup();

    void setMedia(EGS_Input *inp, int nmed, const int *med_ind);

    int howfarFromOut(const EGS_Vector &x, const EGS_Vector &u,
                      EGS_Float &t, int &ix, int &iy, int &iz,
                      int &face, EGS_Vector *normal=0) {
        ix = -1;
        EGS_Float t1;
        if (x.x <= xpos[0] && u.x > 0) {
            t1 = (xpos[0]-x.x)/u.x;
            ix = 0;
            face = 0;
        }
        else if (x.x >= xpos[nx] && u.x < 0) {
            t1 = (xpos[nx]-x.x)/u.x;
            ix = nx-1;
            face = 1;
        }
        if (ix >= 0 && t1 <= t) {
            EGS_Vector tmp(x + u*t1);
            iy = yp->isWhere(tmp);
            if (iy >= 0) {
                iz = zp->isWhere(tmp);
                if (iz >= 0) {
                    int inew =  ix + iy*nx + iz*nxy;
                    if (t1 < boundaryTolerance) {
                        int itest = howfar(inew, tmp, u, t1);
                        if (itest == -1 && t1 < boundaryTolerance) {
                            return -1;
                        }
                    }
                    if (normal) *normal = face == 0 ? EGS_Vector(-1,0,0) :
                                              EGS_Vector(1,0,0);
                    t = t1;
                    return inew;
                }
            }
        }
        iy = -1;
        if (x.y <= ypos[0] && u.y > 0) {
            t1 = (ypos[0]-x.y)/u.y;
            iy = 0;
            face = 2;
        }
        else if (x.y >= ypos[ny] && u.y < 0) {
            t1 = (ypos[ny]-x.y)/u.y;
            iy = ny-1;
            face = 3;
        }
        if (iy >= 0 && t1 <= t) {
            EGS_Vector tmp(x + u*t1);
            ix = xp->isWhere(tmp);
            if (ix >= 0) {
                iz = zp->isWhere(tmp);
                if (iz >= 0) {
                    int inew =  ix + iy*nx + iz*nxy;
                    if (t1 < boundaryTolerance) {
                        int itest = howfar(inew, tmp, u, t1);
                        if (itest == -1 && t1 < boundaryTolerance) {
                            return -1;
                        }
                    }
                    if (normal) *normal = face == 2 ? EGS_Vector(0,-1,0) :
                                              EGS_Vector(0, 1,0);
                    t = t1;
                    return inew;
                }
            }
        }
        iz = -1;
        if (x.z <= zpos[0] && u.z > 0) {
            t1 = (zpos[0]-x.z)/u.z;
            iz = 0;
            face = 4;
        }
        else if (x.z >= zpos[nz] && u.z < 0) {
            t1 = (zpos[nz]-x.z)/u.z;
            iz = nz-1;
            face = 5;
        }
        if (iz >= 0 && t1 <= t) {
            EGS_Vector tmp(x + u*t1);
            ix = xp->isWhere(tmp);
            if (ix >= 0) {
                iy = yp->isWhere(tmp);
                if (iy >= 0) {
                    int inew =  ix + iy*nx + iz*nxy;
                    if (t1 < boundaryTolerance) {
                        int itest = howfar(inew, tmp, u, t1);
                        if (itest == -1 && t1 < boundaryTolerance) {
                            return -1;
                        }
                    }
                    if (normal) *normal = face == 4 ? EGS_Vector(0,0,-1) :
                                              EGS_Vector(0,0, 1);
                    t = t1;
                    return inew;
                }
            }
        }
        return -1;
    };
};

class EGS_NDG_EXPORT EGS_DeformedXYZ : public EGS_XYZGeometry {

public:

    EGS_DeformedXYZ(EGS_PlanesX *Xp, EGS_PlanesY *Yp, EGS_PlanesZ *Zp,
                    const char *defFile, const string &Name = "");

    ~EGS_DeformedXYZ();

    int isWhere(const EGS_Vector &x) {
        int ireg = EGS_XYZGeometry::isWhere(x);
        if (ireg >= 0) egsFatal("EGS_DeformedXYZ::isWhere(): this geometry "
                                    "does not allow the use of this method with position inside\n");
        return ireg;
    };

    int inside(const EGS_Vector &x) {
        return isWhere(x);
    };

    int medium(int ireg) const {
        return EGS_XYZGeometry::medium(ireg/6);
    };

    int howfar(int ireg, const EGS_Vector &x, const EGS_Vector &u,
               EGS_Float &t, int *newmed=0, EGS_Vector *normal=0) {
        if (ireg >= 0) {
            int ir = ireg/6, tetra = ireg - 6*ir, tetra4 = 4*tetra;
            int ind[3];
            ind[2] = ir/nxy;
            ind[1] = (ir - ind[2]*nxy)/nx;
            ind[0] = ir - ind[2]*nxy - ind[1]*nx;
            int edge = ir + ind[1] + ind[2]*nxyp1;
            EGS_Vector n0(vectors[edge+tnodes[tetra4]]),
                       n1(vectors[edge+tnodes[tetra4+1]]-n0),
                       n2(vectors[edge+tnodes[tetra4+2]]-n0),
                       n3(vectors[edge+tnodes[tetra4+3]]-n0);
            EGS_Vector n0x(n0-x);
            EGS_Float disti[4],dpi[4];
            EGS_Vector normvec;
            normvec = n1%n2;
            dpi[0]=normvec*u;
            disti[0]=normvec*n0x;
            normvec = n3%n1;
            dpi[2]=normvec*u;
            disti[2]=normvec*n0x;
            normvec = n2%n3;
            dpi[1]=normvec*u;
            disti[1]=normvec*n0x;
            normvec = (n3-n1)%(n2-n1);
            dpi[3] = normvec*u;
            disti[3] = normvec*(n0x+n1);
            int nextplane = 4;
            EGS_Float mindist = t, mindp=1;
            for (int i=0; i<4; ++i) {
                if (dpi[i] > 0 && disti[i]*mindp < mindist*dpi[i]) {
                    mindist = disti[i];
                    mindp=dpi[i];
                    nextplane=i;
                }
            }
            if (nextplane == 4) {
                return ireg;
            }
            t = mindist/mindp;
            int j = 12*tetra + 3*nextplane;
            int next = tetra_data[j];
            int irnew = ir;
            if (next >= 0) {
                ind[next] += tetra_data[j+1];
                if (ind[next] < 0 || ind[next] >= np[next]) {
                    return -1;
                }
                irnew = ind[0] + ind[1]*nx + ind[2]*nxy;
            }
            if (newmed) {
                *newmed = EGS_XYZGeometry::medium(irnew);
            }
            if (normal) {
                switch (nextplane) {
                case 0:
                    *normal = n2%n1;
                    break;
                case 1:
                    *normal = n1%n3;
                    break;
                case 2:
                    *normal = n3%n2;
                    break;
                case 3:
                    *normal = (n2-n1)%(n3-n1);
                }
            }
            return 6*irnew + tetra_data[j+2];
        }
        int ix,iy,iz,face;
        int inew = howfarFromOut(x,u,t,ix,iy,iz,face,normal);
        if (inew < 0) {
            return inew;
        }
        if (newmed) {
            *newmed = EGS_XYZGeometry::medium(inew);
        }
        // need to determine the tetrahedron we are entering.
        int edge = inew + iy + iz*nxyp1;
        EGS_Vector tmp(x + u*t);
        int tetra1 = enter_tetra1[face], plane1 = enter_plane1[face];
        int node1 = tnodes[4*tetra1 + plane_order[3*plane1]],
            node2 = tnodes[4*tetra1 + plane_order[3*plane1+1]],
            node3 = tnodes[4*tetra1 + plane_order[3*plane1+2]];
        if (inTriangle(vectors[edge+node1],vectors[edge+node2],
                       vectors[edge+node3],tmp)) {
            return 6*inew + tetra1;
        }
        int tetra2 = enter_tetra2[face], plane2 = enter_plane2[face];
        int node1a = tnodes[4*tetra2 + plane_order[3*plane2]],
            node2a = tnodes[4*tetra2 + plane_order[3*plane2+1]],
            node3a = tnodes[4*tetra2 + plane_order[3*plane2+2]];
        if (inTriangle(vectors[edge+node1a],vectors[edge+node2a],
                       vectors[edge+node3a],tmp)) {
            return 6*inew + tetra2;
        }
        // roundoff problem
        egsWarning("deformed geometry: entering from face %d but inTriangle()"
                   " fails for both triangles\n",face);
        egsWarning("x_enter=(%16.10f,%16.10f,%16.10f)\n",tmp.x,tmp.y,tmp.z);
        return -1;
    };

    EGS_Float hownear(int ireg, const EGS_Vector &x) {
        if (ireg < 0) {
            return EGS_XYZGeometry::hownear(ireg,x);
        }
        int ir = ireg/6, tetra = ireg - 6*ir, tetra4 = 4*tetra;
        int ind[3];
        ind[2] = ir/nxy;
        ind[1] = (ir - ind[2]*nxy)/nx;
        ind[0] = ir - ind[2]*nxy - ind[1]*nx;
        int edge = ir + ind[1] + ind[2]*nxyp1;
        EGS_Vector n0(vectors[edge+tnodes[tetra4]]),
                   n1(vectors[edge+tnodes[tetra4+1]]-n0),
                   n2(vectors[edge+tnodes[tetra4+2]]-n0),
                   n3(vectors[edge+tnodes[tetra4+3]]-n0);
        EGS_Vector n0x(n0-x);
        EGS_Float disti[4],dpi[4];
        EGS_Vector normvec;
        normvec = n1%n2;
        dpi[0]=normvec.length2();
        disti[0]=normvec*n0x;
        normvec = n3%n1;
        dpi[2]=normvec.length2();
        disti[2]=normvec*n0x;
        normvec = n2%n3;
        dpi[1]=normvec.length2();
        disti[1]=normvec*n0x;
        normvec = (n3-n1)%(n2-n1);
        dpi[3] = normvec.length2();
        disti[3] = normvec*(n0x+n1);
        EGS_Float mindist = veryFar, mindp=1;
        for (int i=0; i<4; ++i) {
            if (disti[i]*mindp < mindist*dpi[i]) {
                mindist = disti[i];
                mindp=dpi[i];
            }
        }
        return mindist/mindp;
    };

    int computeIntersections(int ireg, int n, const EGS_Vector &X,
                             const EGS_Vector &u, EGS_GeometryIntersections *isections) {
        return EGS_BaseGeometry::computeIntersections(ireg,n,X,u,isections);
    };

    EGS_Float howfarToOutside(int ireg, const EGS_Vector &x,
                              const EGS_Vector &u) {
        if (ireg < 0) {
            return 0;
        }
        int ir = ireg/6;
        int iz = ir/nxy;
        int ir1 = ir - iz*nxy;
        int iy = ir1/nx;
        int ix = ir1 - iy*nx;
        EGS_Float  d = xp->howfarToOutside(ix,x,u);
        EGS_Float ty = yp->howfarToOutside(iy,x,u);
        if (ty < d) {
            d = ty;
        }
        EGS_Float tz = zp->howfarToOutside(iz,x,u);
        if (tz < d) {
            d = tz;
        }
        return d;
    };

    int getMaxStep() const {
        return 6*(nx+ny+nz) + 1;
    };

    const string &getType() const {
        return def_type;
    };

    void printInfo() const {};

    int  setDeformations(const char *defFile);

    bool inTriangle(const EGS_Vector &n0, const EGS_Vector &n1,
                    const EGS_Vector &n2, const EGS_Vector &x) {
        EGS_Vector u1(n1-n0), u2(n2-n0), xp(x-n0);
        EGS_Float a=u1.length2(), b=u1*u2, c=u2.length2(), b1=xp*u1, b2=xp*u2;
        EGS_Float D=a*c-b*b;
        EGS_Float u=b1*c-b2*b;
        if (u<0||u>D) {
            return false;
        }
        EGS_Float v=b2*a-b1*b;
        if (v<0||u+v>D) {
            return false;
        }
        return true;
    };

protected:

    EGS_Vector   *vectors;

    int          tnodes[24];

    int          np[3];

    int          nxyp1;

    static string def_type;

    static char tetrahedra[];

    static char plane_order[];

    static char tetra_data[];

    static char enter_tetra1[];
    static char enter_plane1[];
    static char enter_tetra2[];
    static char enter_plane2[];


};

class EGS_NDG_EXPORT EGS_XYZRepeater : public EGS_BaseGeometry {

public:

    EGS_XYZRepeater(EGS_Float Xmin, EGS_Float Xmax,
                    EGS_Float Ymin, EGS_Float Ymax,
                    EGS_Float Zmin, EGS_Float Zmax,
                    int Nx, int Ny, int Nz, EGS_BaseGeometry *G,
                    const string &Name = "");

    ~EGS_XYZRepeater();

    int isWhere(const EGS_Vector &x) {
        int cell = xyz->isWhere(x);
        //egsWarning("isInside(%g,%g,%g): cell=%d\n",x.x,x.y,x.z,cell);
        if (cell < 0) {
            return cell;
        }
        EGS_Vector xp(x - translation[cell]);
        int ir = g->isWhere(xp);
        //egsWarning("xp=(%g,%g,%g) ir=%d\n",xp.x,xp.y,xp.z,ir);
        return ir >= 0 ? cell*ng + ir : nreg-1;
    };

    bool isInside(const EGS_Vector &x) {
        return xyz->isInside(x);
    };

    int inside(const EGS_Vector &x) {
        return isWhere(x);
    };

    int medium(int ireg) const {
        if (ireg == nreg-1) {
            return med;
        }
        int cell = ireg/ng;
        int ilocal = ireg - ng*cell;
        return g->medium(ilocal);
    };

    EGS_Float howfarToOutside(int ireg, const EGS_Vector &x,
                              const EGS_Vector &u) {
        if (ireg < 0) {
            return 0;
        }
        return xyz->howfarToOutside(0,x,u);
    };

    int howfar(int ireg, const EGS_Vector &x, const EGS_Vector &u,
               EGS_Float &t, int *newmed=0, EGS_Vector *normal=0) {
        //egsWarning("\n***howfar(reg=%d,x=(%g,%g,%g),u=(%g,%g,%g),t=%g)\n",
        //        ireg,x.x,x.y,x.z,u.x,u.y,u.z,t);
        if (ireg < 0) {
            // outside of the box.
            int cell = xyz->howfar(ireg,x,u,t,0,normal);
            //egsWarning("entering: cell=%d t=%g\n",cell,t);
            if (cell >= 0) {
                if (newmed) {
                    *newmed = med;
                }
                return nreg-1;
            }
            return -1;
        }
        if (ireg >= 0) {
            if (ireg < nreg-1) {  // in one of the repetitions
                int cell = ireg/ng;
                int ilocal = ireg - cell*ng;
                EGS_Vector xp(x - translation[cell]);
                int inew = g->howfar(ilocal,xp,u,t,newmed,normal);
                //egsWarning("cell=%d il=%d inew=%d t=%g xp=(%g,%g,%g)\n",
                //        cell,ilocal,inew,t,xp.x,xp.y,xp.z);
                if (inew >= 0) {
                    return cell*ng + inew;
                }
                if (newmed) {
                    *newmed = med;
                }
                return nreg-1;
            }
        }
        // if here, we are in the box outside of all repetitions
        // find the cell we are in
        int ix = (int)(x.x*dxi + ax);
        if (ix >= nx) {
            ix = nx-1;
        }
        int iy = (int)(x.y*dyi + ay);
        if (iy >= ny) {
            iy = ny-1;
        }
        int iz = (int)(x.z*dzi + az);
        if (iz >= nz) {
            iz = nz-1;
        }
        int cell = ix + iy*nx + iz*nxy;
        //egsWarning("in box: ix=%d iy=%d iz=%d\n",ix,iy,iz);
        EGS_Float t_left = t;
        EGS_Vector xtmp(x);
        EGS_Float ttot = 0;
        for (EGS_I64 loopCount=0; loopCount<=loopMax; ++loopCount) {
            if (loopCount == loopMax) {
                egsFatal("EGS_XYZRepeater::howfar: Too many iterations were required! Input may be invalid, or consider increasing loopMax.");
                return -1;
            }
            EGS_Float this_t = t_left;
            EGS_Vector xp(xtmp - translation[cell]);
            int inew = g->howfar(-1,xp,u,this_t,newmed,normal);
            //egsWarning("cell=%d tleft=%g xp=(%g,%g,%g) inew=%d this_t=%g\n",
            //        cell,t_left,xp.x,xp.y,xp.z,inew,this_t);
            if (inew >= 0) {
                t = ttot + this_t;
                return cell*ng + inew;
            }
            int next_cell = xyz->howfar(cell,xtmp,u,this_t,0,normal);
            //egsWarning("next: %d %g\n",next_cell,this_t);
            if (next_cell == cell) {
                return ireg;
            }
            if (next_cell < 0) {
                t = ttot + this_t;
                return -1;
            }
            ttot += this_t;
            t_left -= this_t;
            xtmp += u*this_t;
            cell = next_cell;
        }

        return ireg;
    };

    EGS_Float hownear(int ireg, const EGS_Vector &x) {
        if (ireg < 0) {
            return xyz->hownear(ireg,x);
        }
        if (ireg < nreg-1) {
            int cell = ireg/ng;
            int ilocal = ireg - cell*ng;
            EGS_Vector xp(x - translation[cell]);
            return g->hownear(ilocal,xp);
        }
        // if here, we are in the box outside of all repetitions
        // find the cell we are in
        int ix = (int)(x.x*dxi + ax);
        if (ix >= nx) {
            ix = nx-1;
        }
        int iy = (int)(x.y*dyi + ay);
        if (iy >= ny) {
            iy = ny-1;
        }
        int iz = (int)(x.z*dzi + az);
        if (iz >= nz) {
            iz = nz-1;
        }
        int cell = ix + iy*nx + iz*nxy;
        EGS_Vector xp(x - translation[cell]);
        EGS_Float t = g->hownear(-1,xp);
        EGS_Float tcell = xyz->hownear(cell,x);
        if (t < tcell) {
            return t;
        }
        for (int iix=ix-1; iix<=ix+1; ++iix) {
            if (iix >= 0 && iix < nx) {
                for (int iiy=iy-1; iiy<=iy+1; ++iiy) {
                    if (iiy >= 0 && iiy < ny) {
                        for (int iiz=iz-1; iiz<=iz+1; ++iiz) {
                            if (iiz >= 0 && iiz < nz) {
                                if (iix != ix || iiy != iy || iiz != iz) {
                                    int cell1 = iix+iiy*nx+iiz*nz;
                                    EGS_Vector tmp(x-translation[cell1]);
                                    EGS_Float t1 = g->hownear(-1,tmp);
                                    if (t1 < t) {
                                        t = t1;
                                    }
                                }
                            }
                        }
                    }
                }
            }
        }
        return t;
    };

    int getMaxStep() const {
        return nxyz*(g->getMaxStep() + 1) + 1;
    };

    const string &getType() const {
        return type;
    };

    void printInfo() const;

    void setXYZLabels(EGS_Input *input) {
        xyz->setXYZLabels(input);
    }

    virtual void getLabelRegions(const string &str, vector<int> &regs);

protected:

    EGS_XYZGeometry  *xyz;
    EGS_BaseGeometry *g;
    EGS_Vector       *translation;
    EGS_Float        dx, dxi, ax;
    EGS_Float        dy, dyi, ay;
    EGS_Float        dz, dzi, az;
    int              nx, ny, nz, nxy, nxyz, ng;

    static string    type;

};

#endif

#endif