egs_shapes.h

Go to the documentation of this file.

/*
###############################################################################
#
#  EGSnrc egs++ shapes 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:    Frederic Tessier
#                   Marc Chamberland
#                   Reid Townson
#
###############################################################################
*/
 
 
#ifndef EGS_SHAPES_
#define EGS_SHAPES_
 
#include "egs_vector.h"
#include "egs_transformations.h"
#include "egs_rndm.h"
#include "egs_object_factory.h"
#include "egs_input_struct.h"
 
#include <string>
using std::string;
 
class EGS_Input;
 
inline void setShapeInputs(shared_ptr<EGS_BlockInput> shapePtr) {
    auto libPtr = shapePtr->addSingleInput("library", false, "The type of shape, loaded by shared library in egs++/dso.");
    auto typePtr = shapePtr->addSingleInput("type", false, "The type of shape - this input includes only a small set of simple shapes. For more options, use the 'library' input instead.", {"point", "box", "sphere", "cylinder"});
 
    // Only one of "library" or "type" are allowed
    libPtr->addDependency(typePtr, "", true);
    typePtr->addDependency(libPtr, "", true);
 
    // Point
    shapePtr->addSingleInput("position", true, "The x, y, z position that the source will emit particles from.")->addDependency(typePtr, "point");
 
    // Box
    shapePtr->addSingleInput("box size", true, "The side lengths of the box, in cm. Enter 1 number for a cube, or 3 numbers to denote the x, y, and z side lengths.")->addDependency(typePtr, "box");
 
    // Sphere
    auto radiusPtr = shapePtr->addSingleInput("radius", true, "The radius of the sphere or cylinder, in cm.");
    radiusPtr->addDependency(typePtr, "sphere");
    auto midPtr = shapePtr->addSingleInput("midpoint", false, "The x, y and z coordinates of the midpoint of the sphere or cylinder, in cm. Defaults to 0, 0, 0.");
    midPtr->addDependency(typePtr, "sphere");
 
    // Cylinder
    radiusPtr->addDependency(typePtr, "cylinder");
    midPtr->addDependency(typePtr, "cylinder");
    shapePtr->addSingleInput("height", true, "The height of the cylinder, in cm.")->addDependency(typePtr, "cylinder");
    shapePtr->addSingleInput("phi range", false, "The minimum and maximum phi values, in degrees. This allows you restrict the cylinder to a shape like a slice of pie!")->addDependency(typePtr, "cylinder");
    shapePtr->addSingleInput("axis", true, "A unit vector that defines the axis of the cylinder.")->addDependency(typePtr, "cylinder");
 
    addTransformationBlock(shapePtr);
}
 
class EGS_EXPORT EGS_BaseShape : public EGS_Object {
 
public:
 
    EGS_BaseShape(const string &Name="",EGS_ObjectFactory *f=0) :
        EGS_Object(Name,f), T(0) {
        otype = "base_shape";
    };
    virtual ~EGS_BaseShape() {
        if (T) {
            delete T;
        }
    };
 
    virtual EGS_Vector getRandomPoint(EGS_RandomGenerator *rndm) {
        if (T) {
            return (*T)*getPoint(rndm);
        }
        else {
            return getPoint(rndm);
        }
    };
 
    virtual EGS_Vector getPoint(EGS_RandomGenerator *rndm) {
        (void)rndm;
        egsFatal("You need to implement the getPoint function in your "
                 "derived class\n");
        return EGS_Vector();
    };
 
    void setTransformation(EGS_Input *inp);
 
    void setTransformation(EGS_AffineTransform *t) {
        if (T) {
            delete T;
        }
        T = new EGS_AffineTransform(*t);
    };
 
    const EGS_AffineTransform *getTransform() const {
        return T;
    };
 
    static EGS_BaseShape *createShape(EGS_Input *inp);
 
    static EGS_BaseShape *getShape(const string &Name);
 
    virtual bool supportsDirectionMethod() const {
        return false;
    };
 
    /* getNextShapePosition is the equivalent of getNextParticle but for a shape object. Its goal is to determine the next state of the geometry, either by synchronizing itself to the
     * source time parameter, or by sampling it's own time parameter and updating itself accordingly if the source has provided no time index.
     *
     *This function has a non-empty implementation in 2 cases.
     *1) it is re implemented in any composite shape, where it will call getNextShapePosition on all of its components
     *2) it is re implemented in the dynamic shape class. This is where the code will find the current (non static) state of the shape. */
    virtual void getNextShapePosition(EGS_RandomGenerator *rndm) {
        (void)rndm;
    };
 
    virtual void getPointSourceDirection(const EGS_Vector &xo,
                                         EGS_RandomGenerator *rndm, EGS_Vector &u, EGS_Float &wt) {
        (void)xo;
        (void)rndm;
        (void)u;
        (void)wt;
        egsFatal("getPointSourceDirection: you have to implement this "
                 "method for the %s shape if you want to use it\n",otype.c_str());
    };
 
    virtual EGS_Float area() const {
        return 1;
    };
 
    virtual void updatePosition(EGS_Float time) {
        (void)time;
    };
 
protected:
 
    EGS_AffineTransform *T; 
 
};
 
class  EGS_EXPORT EGS_SurfaceShape : public EGS_BaseShape {
 
public:
 
    EGS_SurfaceShape(const string &Name="",EGS_ObjectFactory *f=0) :
        EGS_BaseShape(Name,f), A(1) {};
    ~EGS_SurfaceShape() {};
    bool supportsDirectionMethod() const {
        return true;
    };
    EGS_Float area() const {
        return A;
    };
    void getPointSourceDirection(const EGS_Vector &Xo,
                                 EGS_RandomGenerator *rndm, EGS_Vector &u, EGS_Float &wt) {
        EGS_Vector xo = T ? Xo*(*T) : Xo;
        EGS_Vector x = getPoint(rndm);
        u = x - xo;
        EGS_Float d2i = 1/u.length2(), di = sqrt(d2i);
        u *= di;
        wt = A*fabs(u.z)*d2i;
        if (T) {
            T->rotate(u);
        }
    };
 
protected:
 
    EGS_Float A;
 
};
 
class EGS_EXPORT EGS_PointShape : public EGS_BaseShape {
 
public:
 
    EGS_PointShape(const EGS_Vector &Xo = EGS_Vector(),
                   const string &Name="",EGS_ObjectFactory *f=0) :
        EGS_BaseShape(Name,f), xo(Xo) {
        otype = "point";
    };
    ~EGS_PointShape() { };
    EGS_Vector getPoint(EGS_RandomGenerator *) {
        return xo;
    };
    EGS_Object *createObject(EGS_Input *inp);
 
protected:
 
    EGS_Vector  xo; 
 
};
 
class EGS_EXPORT EGS_BoxShape : public EGS_BaseShape {
 
protected:
 
    EGS_Float  ax, ay, az; 
 
public:
 
    EGS_BoxShape(const string &Name="",EGS_ObjectFactory *f=0) :
        EGS_BaseShape(Name,f), ax(1), ay(1), az(1) {
        otype="box";
    };
    EGS_BoxShape(EGS_Float A, const EGS_AffineTransform *t = 0,
                 const string &Name="",EGS_ObjectFactory *f=0) :
        EGS_BaseShape(Name,f), ax(A), ay(A), az(A) {
        if (t) {
            T = new EGS_AffineTransform(*t);
        }
        otype="box";
    };
    EGS_BoxShape(EGS_Float Ax, EGS_Float Ay, EGS_Float Az,
                 const EGS_AffineTransform *t = 0,
                 const string &Name="",EGS_ObjectFactory *f=0) :
        EGS_BaseShape(Name,f), ax(Ax), ay(Ay), az(Az) {
        if (t) {
            T = new EGS_AffineTransform(*t);
        }
        otype="box";
    };
    ~EGS_BoxShape() { };
 
    EGS_Vector getPoint(EGS_RandomGenerator *rndm) {
        EGS_Vector v(ax*(rndm->getUniform()-0.5),
                     ay*(rndm->getUniform()-0.5),
                     az*(rndm->getUniform()-0.5));
        return v;
    };
 
    EGS_Object *createObject(EGS_Input *);
 
    bool supportsDirectionMethod() const {
        return true;
    };
 
    void getPointSourceDirection(const EGS_Vector &Xo,
                                 EGS_RandomGenerator *rndm, EGS_Vector &u, EGS_Float &wt) {
        EGS_Vector xo = T ? Xo*(*T) : Xo;
        EGS_Float eta = rndm->getUniform()*area();
        if (eta < 2*ax*ay) {
            u.x = ax*(rndm->getUniform()-0.5);
            u.y = ay*(rndm->getUniform()-0.5);
            if (eta < ax*ay) {
                u.z = az/2;
                wt = u.z - xo.z;
            }
            else {
                u.z = -az/2;
                wt = xo.z - u.z;
            }
        }
        else if (eta < 2*(ax*ay + ax*az)) {
            u.x = ax*(rndm->getUniform()-0.5);
            u.z = az*(rndm->getUniform()-0.5);
            if (eta < 2*ax*ay + ax*az) {
                u.y = ay/2;
                wt = u.y - xo.y;
            }
            else {
                u.y = -ay/2;
                wt = xo.y - u.y;
            }
        }
        else {
            eta -= 2*(ax*ay + ax*az);
            u.y = ay*(rndm->getUniform()-0.5);
            u.z = az*(rndm->getUniform()-0.5);
            if (eta < ay*az) {
                u.x = ax/2;
                wt = u.x - xo.x;
            }
            else {
                u.x = -ax/2;
                wt = xo.x - u.x;
            }
        }
        u -= xo;
        EGS_Float d2 = u.length2(), d = sqrt(d2);
        u *= (1/d);
        wt *= (area()/(d2*d));
        if (T) {
            T->rotate(u);
        }
    };
    EGS_Float area() const {
        return 2*(ax*ay + ax*az + ay*az);
    };
 
};
 
class EGS_EXPORT EGS_SphereShape : public EGS_BaseShape {
 
protected:
 
    EGS_Float  R;     
    EGS_Vector xo;    
 
public:
 
    EGS_SphereShape(const string &Name="",EGS_ObjectFactory *f=0) :
        EGS_BaseShape(Name,f), R(1), xo() {
        otype="sphere";
    };
    EGS_SphereShape(EGS_Float r, const EGS_Vector &Xo = EGS_Vector(0,0,0),
                    const string &Name="",EGS_ObjectFactory *f=0) :
        EGS_BaseShape(Name,f), R(r), xo(Xo) {
        otype = "sphere";
    };
    ~EGS_SphereShape() {};
 
    EGS_Vector getPoint(EGS_RandomGenerator *rndm) {
        EGS_Float r = rndm->getUniform(), r1 = rndm->getUniform(),
                  r2 = rndm->getUniform();
        if (r1 > r) {
            r = r1;
        }
        if (r2 > r) {
            r = r2;
        }
        EGS_Float cost = 2*rndm->getUniform()-1;
        EGS_Float sint = sqrt(1-cost*cost);
        r1 = R*r*sint;
        EGS_Float cphi, sphi;
        rndm->getAzimuth(cphi,sphi);
        return xo + EGS_Vector(r1*cphi,r1*sphi,R*r*cost);
    };
 
    EGS_Object *createObject(EGS_Input *inp);
 
    bool supportsDirectionMethod() const {
        return true;
    };
 
    void getPointSourceDirection(const EGS_Vector &Xo,
                                 EGS_RandomGenerator *rndm, EGS_Vector &u, EGS_Float &wt) {
        EGS_Vector xo = T ? Xo*(*T) : Xo;
        EGS_Float cost = 2*rndm->getUniform()-1;
        EGS_Float sint = 1-cost*cost;
        EGS_Vector x;
        if (sint > epsilon) {
            EGS_Float cphi, sphi;
            rndm->getAzimuth(cphi,sphi);
            sint = R*sqrt(sint);
            x.x = sint*cphi;
            x.y = sint*sphi;
            x.z = R*cost;
        }
        else {
            x.z = R*cost;
        }
        u = (x + this->xo) - xo;
        EGS_Float di = 1/u.length();
        u *= di;
        wt = u*x*4*M_PI*R*di*di;
    };
    EGS_Float area() const {
        return 4*M_PI*R*R;
    };
};
 
class EGS_EXPORT EGS_CylinderShape : public EGS_BaseShape {
 
protected:
 
    EGS_Float  R;       
    EGS_Float  h;       
    EGS_Vector xo;      
    EGS_Vector a;       
    EGS_Float  phi_min;
    EGS_Float  phi_max;
    bool       has_phi; 
 
    inline void getPointInCircle(EGS_RandomGenerator *rndm, EGS_Float &x,
                                 EGS_Float &y) {
        if (!has_phi) {
            do {
                x = 2*rndm->getUniform()-1;
                y = 2*rndm->getUniform()-1;
            }
            while (x*x + y*y > 1);
            x *= R;
            y *= R;
        }
        else {
            EGS_Float r = R*sqrt(rndm->getUniform());
            EGS_Float eta = rndm->getUniform();
            EGS_Float phi = phi_min*(1-eta) + phi_max*eta;
            x = r*cos(phi);
            y = r*sin(phi);
        }
    };
 
 
public:
 
    EGS_CylinderShape(const string &Name="",EGS_ObjectFactory *f=0) :
        EGS_BaseShape(), R(1), h(1), xo(), a(0,0,1),
        phi_min(0), phi_max(2*M_PI), has_phi(false) {
        (void)Name;
        (void)f;
        otype="cylinder";
    };
    EGS_CylinderShape(EGS_Float r, EGS_Float H,
                      const EGS_Vector &Xo = EGS_Vector(0,0,0),
                      const EGS_Vector &A = EGS_Vector(0,0,1),
                      const string &Name="",EGS_ObjectFactory *f=0) :
        EGS_BaseShape(Name,f), R(r), h(H), xo(Xo), a(A),
        phi_min(0), phi_max(2*M_PI), has_phi(false) {
        EGS_RotationMatrix rmat(a);
        if (xo.length2() > epsilon || !rmat.isI()) {
            T = new EGS_AffineTransform(rmat.inverse(),xo);
        }
        otype="cylinder";
    };
    EGS_CylinderShape(EGS_Float r, EGS_Float H, const EGS_AffineTransform *t,
                      const string &Name="",EGS_ObjectFactory *f=0) :
        EGS_BaseShape(Name,f), R(r), h(H),
        phi_min(0), phi_max(2*M_PI), has_phi(false) {
        if (t) {
            T = new EGS_AffineTransform(*t);
        }
        otype="cylinder";
    };
    ~EGS_CylinderShape() { };
 
    void setPhiRange(EGS_Float Phi_min, EGS_Float Phi_max) {
        if (Phi_min < Phi_max) {
            phi_min = Phi_min;
            phi_max = Phi_max;
        }
        else                    {
            phi_min = Phi_max;
            phi_max = Phi_min;
        }
        if (phi_max - phi_min < 1.99999*M_PI) {
            has_phi = true;
        }
        else {
            has_phi = false;
        }
    };
 
    EGS_Vector getPoint(EGS_RandomGenerator *rndm) {
        EGS_Float x,y;
        getPointInCircle(rndm,x,y);
        EGS_Float z = h*(rndm->getUniform()-0.5);
        return EGS_Vector(x,y,z);
    };
 
    EGS_Object *createObject(EGS_Input *);
 
    EGS_Float getRadius() const {
        return R;
    };
    EGS_Float getHeight() const {
        return h;
    };
 
    bool supportsDirectionMethod() const {
        return true;
    };
 
    void getPointSourceDirection(const EGS_Vector &Xo,
                                 EGS_RandomGenerator *rndm, EGS_Vector &u, EGS_Float &wt) {
        EGS_Vector xo = T ? Xo*(*T) : Xo;
        EGS_Float eta = rndm->getUniform()*(R+h);
        EGS_Vector x;                               // point on cylinder with respect to midpoint
        EGS_Vector n;                               // normal to cylinder at point x
        if (eta < R) {
            getPointInCircle(rndm,x.x,x.y);
            if (2*eta < R) {
                x.z = h/2;    // top face: normal is up
                n.z = 1;
            }
            else {
                x.z = -h/2;    // bottom face: normal is down
                n.z = -1;
            }
        }
        else {
            EGS_Float cphi,sphi;
            rndm->getAzimuth(cphi,sphi);
            x.x = R*cphi;
            x.y = R*sphi;
            x.z = h*(rndm->getUniform()-0.5);
            n.x = x.x;
            n.y = x.y;                   // side face: normal is (x,y)
        }
        u = (x+this->xo) - xo;                      // direction vector from origin to cylinder point
        EGS_Float d2 = u.length2(), d = sqrt(d2);
        u *= (1/d);                                 // normalize direction vectors
        n.normalize();                              // normalize normal
        wt = u*n*area()/d2;
        if (T) {
            T->rotate(u);
        }
    };
    EGS_Float area() const {
        return 2*M_PI*R*(R+h);
    };
};
 
#endif