egs_transformations.h

Go to the documentation of this file.

/*
###############################################################################
#
#  EGSnrc egs++ transformations 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:    Marc Chamberland
#                   Frederic Tessier
#                   Reid Townson
#                   Jan Weidner
#                   Alexandre Demelo
#
###############################################################################
*/
 
 
#ifndef EGS_TRANSFORMATIONS_
 
#define EGS_TRANSFORMATIONS_
 
#include "egs_vector.h"
#include "egs_libconfig.h"
#include "egs_math.h"
#include "egs_functions.h"
#include "egs_input_struct.h"
 
#include <iostream>
#include <vector>
 
using namespace std;
 
class EGS_Input;
 
inline void addTransformationBlock(shared_ptr<EGS_BlockInput> blockPtr) {
    shared_ptr<EGS_BlockInput> transBlock = blockPtr->addBlockInput("transformation");
    transBlock->addSingleInput("translation", false, "The x, y, z translation offsets in cm.");
    auto vecPtr = transBlock->addSingleInput("rotation vector", false, "Defines a rotation which, when applied to the 3D vector defined by this input, transforms it into a vector along the positive z-axis.");
    auto rotPtr = transBlock->addSingleInput("rotation", false, "2, 3 or 9 floating point numbers define a rotation. See the documentation for details.");
    vecPtr->addDependency(rotPtr, "", true);
    rotPtr->addDependency(vecPtr, "", true);
}
 
class EGS_EXPORT EGS_RotationMatrix {
 
protected:
 
    EGS_Float   rxx, rxy, rxz,
                ryx, ryy, ryz,
                rzx, rzy, rzz;
 
public:
 
    EGS_RotationMatrix() {
        rxx=1;
        rxy=0;
        rxz=0;
        ryx=0;
        ryy=1;
        ryz=0;
        rzx=0;
        rzy=0;
        rzz=1;
    };
 
    //
    // hopefully no one will use this constructor
    // with a non-rotation matrix elements
    //
    EGS_RotationMatrix(EGS_Float xx, EGS_Float xy, EGS_Float xz,
                       EGS_Float yx, EGS_Float yy, EGS_Float yz,
                       EGS_Float zx, EGS_Float zy, EGS_Float zz) {
        rxx=xx;
        rxy=xy;
        rxz=xz;
        ryx=yx;
        ryy=yy;
        ryz=yz;
        rzx=zx;
        rzy=zy;
        rzz=zz;
    };
 
    EGS_RotationMatrix(const EGS_RotationMatrix &m) {
        rxx=m.rxx;
        rxy=m.rxy;
        rxz=m.rxz;
        ryx=m.ryx;
        ryy=m.ryy;
        ryz=m.ryz;
        rzx=m.rzx;
        rzy=m.rzy;
        rzz=m.rzz;
    };
 
    bool isRotation() const {
        EGS_Float d = det() - 1;
        EGS_RotationMatrix t((*this)*inverse());
        if (fabs(d) > epsilon || !t.isI()) {
            return false;
        }
        return true;
    };
 
    // creates a matrix which, when applied to the vector v,
    // transforms it into a vector along the z-axis.
    // why z-axis?
    // well, it has to be one of the axis and in physics
    // things usually happen along the z-axis!
    EGS_RotationMatrix(const EGS_Vector &v) {
        EGS_Float sinz = v.x*v.x + v.y*v.y;
        EGS_Float norm = sinz + v.z*v.z;
        if (norm < epsilon) egsFatal("EGS_RotationMatrix::EGS_RotationMatrix: \n"
                                         "  no construction from a zero vector possible!\n");
        norm = sqrt(norm);
        if (sinz > epsilon) {
            sinz = sqrt(sinz);
            EGS_Float cphi = v.x/sinz;
            EGS_Float sphi = v.y/sinz;
            EGS_Float cost = v.z/norm;
            EGS_Float sint = -sinz/norm;
            *this = rotY(cost,sint)*rotZ(cphi,sphi);
        }
        else if (v.z < 0.) {
            // v along negative z is degenerate: a rotation of pi around any axis in the
            // xy plane satisfies the condition; we pick the x axis
            *this = rotX(M_PI);
        }
        else { // v is along the z-axis => matrix is the unit transformation
            rxx = 1;
            rxy = 0;
            rxz = 0;
            ryx = 0;
            ryy = 1;
            ryz = 0;
            rzx = 0;
            rzy = 0;
            rzz = 1;
        }
    };
 
    // R_x(alpha) * R_y(beta) * R_z(gamma)
    EGS_RotationMatrix(EGS_Float alpha, EGS_Float beta, EGS_Float gamma) {
        *this = rotX(alpha)*rotY(beta)*rotZ(gamma);
    };
 
    // R_z(phi) * R_x(theta)
    EGS_RotationMatrix(EGS_Float phi, EGS_Float theta) {
        *this = rotZ(phi)*rotX(theta);
    };
 
    EGS_RotationMatrix &operator=(const EGS_RotationMatrix &m) {
        rxx=m.rxx;
        rxy=m.rxy;
        rxz=m.rxz;
        ryx=m.ryx;
        ryy=m.ryy;
        ryz=m.ryz;
        rzx=m.rzx;
        rzy=m.rzy;
        rzz=m.rzz;
        return *this;
    };
 
    bool operator==(const EGS_RotationMatrix &m) {
        return
            ((rxx == m.rxx) && (rxy == m.rxy) && (rxz == m.rxz) &&
             (ryx == m.ryx) && (ryy == m.ryy) && (ryz == m.ryz) &&
             (rzx == m.rxz) && (rzy == m.rzy) && (rzz == m.rzz)) ? true : false;
    };
 
    bool isI() const {
 
        // compare with unit matrix components
        if (fabs(rxy)   < epsilon &&
                fabs(rxz)   < epsilon &&
                fabs(ryx)   < epsilon &&
                fabs(ryz)   < epsilon &&
                fabs(rzx)   < epsilon &&
                fabs(rzy)   < epsilon &&
                fabs(rxx-1) < epsilon &&
                fabs(ryy-1) < epsilon &&
                fabs(rzz-1) < epsilon) {
            return true;
        }
        return false;
    };
 
    EGS_Vector operator*(const EGS_Vector &v) const {
        return EGS_Vector(rxx*v.x + rxy*v.y + rxz*v.z,
                          ryx*v.x + ryy*v.y + ryz*v.z,
                          rzx*v.x + rzy*v.y + rzz*v.z);
    };
 
    EGS_RotationMatrix operator*(const EGS_RotationMatrix &m) const {
        return EGS_RotationMatrix(
                   rxx*m.rxx+rxy*m.ryx+rxz*m.rzx, rxx*m.rxy+rxy*m.ryy+rxz*m.rzy,
                   rxx*m.rxz+rxy*m.ryz+rxz*m.rzz,
                   ryx*m.rxx+ryy*m.ryx+ryz*m.rzx, ryx*m.rxy+ryy*m.ryy+ryz*m.rzy,
                   ryx*m.rxz+ryy*m.ryz+ryz*m.rzz,
                   rzx*m.rxx+rzy*m.ryx+rzz*m.rzx, rzx*m.rxy+rzy*m.ryy+rzz*m.rzy,
                   rzx*m.rxz+rzy*m.ryz+rzz*m.rzz);
    };
 
    // r *= m means r = r*m
    EGS_RotationMatrix &operator*=(const EGS_RotationMatrix &m) {
        return *this = operator * (m);
    };
 
    // r.multiply(m) means r = m*r;
    void multiply(const EGS_RotationMatrix &m) {
        *this = m.operator * (*this);
    };
 
    EGS_RotationMatrix inverse() const {
        return EGS_RotationMatrix(rxx,ryx,rzx,rxy,ryy,rzy,rxz,ryz,rzz);
    };
 
    EGS_RotationMatrix T() {
        return EGS_RotationMatrix(rxx,ryx,rzx,rxy,ryy,rzy,rxz,ryz,rzz);
    };
 
    EGS_RotationMatrix &invert() {
        return *this = inverse();
    };
 
    static EGS_RotationMatrix rotX(EGS_Float cphi,EGS_Float sphi) {
        return EGS_RotationMatrix(
                   (EGS_Float)1,(EGS_Float)0,(EGS_Float)0,
                   (EGS_Float)0,        cphi,        sphi,
                   (EGS_Float)0,       -sphi,        cphi);
    };
 
    static EGS_RotationMatrix rotY(EGS_Float cphi,EGS_Float sphi) {
        return EGS_RotationMatrix(cphi, (EGS_Float)0,         -sphi,
                                  (EGS_Float)0, (EGS_Float)1, (EGS_Float)0,
                                  sphi, (EGS_Float)0,         cphi);
    };
 
    static EGS_RotationMatrix rotZ(EGS_Float cphi,EGS_Float sphi) {
        return EGS_RotationMatrix(cphi,        sphi, (EGS_Float)0,
                                  -sphi,        cphi, (EGS_Float)0,
                                  (EGS_Float)0,(EGS_Float)0, (EGS_Float)1);
    };
 
    static EGS_RotationMatrix rotX(EGS_Float phi) {
        return rotX(cos(phi),sin(phi));
    };
 
    static EGS_RotationMatrix rotY(EGS_Float phi) {
        return rotY(cos(phi),sin(phi));
    };
 
    static EGS_RotationMatrix rotZ(EGS_Float phi) {
        return rotZ(cos(phi),sin(phi));
    };
 
    static EGS_RotationMatrix rotV(EGS_Float phi, const EGS_Vector &v) {
        return rotV(cos(phi),sin(phi),v);
    };
 
    static EGS_RotationMatrix rotV(EGS_Float cphi, EGS_Float sphi,
                                   const EGS_Vector &v) {
        EGS_RotationMatrix m(v);
        return (m.inverse())*(rotZ(cphi,sphi))*(m);
    };
 
    EGS_Float det() const {
        return rxx*ryy*rzz + rxy*ryz*rzx + ryx*rzy*rxz -
               rxz*ryy*rzx - rxy*ryx*rzz - rzy*ryz*rxx;
    };
 
    friend EGS_Vector operator*(const EGS_Vector &v,
                                const EGS_RotationMatrix &m) {
        return EGS_Vector(v.x*m.rxx+v.y*m.ryx+v.z*m.rzx,
                          v.x*m.rxy+v.y*m.ryy+v.z*m.rzy,
                          v.x*m.rxz+v.y*m.ryz+v.z*m.rzz);
    };
 
    friend EGS_Vector &operator*=(EGS_Vector &v, const EGS_RotationMatrix &m) {
        v = v*m;
        return v;
    };
 
    inline EGS_Float xx() const {
        return rxx;
    };
    inline EGS_Float xy() const {
        return rxy;
    };
    inline EGS_Float xz() const {
        return rxz;
    };
    inline EGS_Float yx() const {
        return ryx;
    };
    inline EGS_Float yy() const {
        return ryy;
    };
    inline EGS_Float yz() const {
        return ryz;
    };
    inline EGS_Float zx() const {
        return rzx;
    };
    inline EGS_Float zy() const {
        return rzy;
    };
    inline EGS_Float zz() const {
        return rzz;
    };
 
};
 
class EGS_EXPORT EGS_AffineTransform {
 
protected:
 
    EGS_RotationMatrix R;
    EGS_Vector         t;
    bool               has_t, has_R;
 
public:
 
    EGS_AffineTransform() : R(),t(),has_t(false),has_R(false) {};
 
    EGS_AffineTransform(const EGS_AffineTransform &tr) :
        R(tr.R),t(tr.t),has_t(tr.has_t),has_R(tr.has_R) {};
 
    EGS_AffineTransform(const EGS_RotationMatrix &m, const EGS_Vector &v) :
        R(m),t(v) {
        if (t.length2() > 0) {
            has_t = true;
        }
        else {
            has_t = false;
        }
        if (R.isI()) {
            has_R = false;
        }
        else {
            has_R = true;
        }
    };
 
    EGS_AffineTransform(const EGS_RotationMatrix &m) : R(m),t() {
        has_t = false;
        if (R.isI()) {
            has_R = false;
        }
        else {
            has_R = true;
        }
    };
 
    EGS_AffineTransform(const EGS_Vector &v) : R(),t(v) {
        has_R = false;
        if (t.length2() > 0) {
            has_t = true;
        }
        else {
            has_t = false;
        }
    };
 
    // Explicitly declare the copy assignment operator
    EGS_AffineTransform &operator=(const EGS_AffineTransform &other) {
        if (this != &other) { // Protect against self-assignment
            // Copy all member variables manually:
            this->R = other.R;
            this->t = other.t;
            this->has_t = other.has_t;
            this->has_R = other.has_R;
        }
        return *this;
    };
 
    EGS_AffineTransform operator*(const EGS_AffineTransform &tr) const {
        return EGS_AffineTransform(R*tr.R,R*tr.t+t);
    };
 
    EGS_AffineTransform operator*(const EGS_RotationMatrix &m) const {
        return EGS_AffineTransform(R*m,t);
    };
 
    EGS_AffineTransform &operator*=(const EGS_AffineTransform &tr) {
        return *this = operator * (tr);
    };
 
    EGS_AffineTransform &operator*=(const EGS_RotationMatrix &m) {
        return *this = operator * (m);
    };
 
    EGS_AffineTransform operator+(const EGS_Vector &v) const {
        return EGS_AffineTransform(R,t+v);
    };
 
    EGS_AffineTransform &operator+=(const EGS_Vector &v) {
        t += v;
        return *this;
    };
 
    EGS_Vector operator*(const EGS_Vector &v) const {
        return (R*v + t);
    };
 
    friend EGS_Vector operator*(const EGS_Vector &v,
                                const EGS_AffineTransform &tr) {
        return ((v-tr.t)*tr.R);
    };
 
    friend EGS_Vector &operator*=(EGS_Vector &v,
                                  const EGS_AffineTransform &tr) {
        v = v*tr;
        return v;
    };
 
    void transform(EGS_Vector &v) const {
        if (has_R) {
            v = R*v;
        }
        if (has_t) {
            v += t;
        }
    };
 
    void inverseTransform(EGS_Vector &v) const {
        if (has_t) {
            v -= t;
        }
        if (has_R) {
            v *= R;
        }
        //v -= t; v *= R;
    };
 
    EGS_AffineTransform inverse() const {
        EGS_Vector tmp;
        tmp -= t*R;
        return EGS_AffineTransform(R.inverse(),tmp);
    };
 
    void rotate(EGS_Vector &v) const {
        if (has_R) {
            v = R*v;
        }
    };
    void rotateInverse(EGS_Vector &v) const {
        if (has_R) {
            v *= R;
        }
    };
    void translate(EGS_Vector &v) const {
        v += t;
    };
 
    const EGS_Vector &getTranslation() const {
        return t;
    };
    const EGS_RotationMatrix &getRotation() const {
        return R;
    };
 
    bool isI() const {
        return (!has_R && !has_t);
    };
 
    bool hasTranslation() const {
        return has_t;
    };
 
    bool hasRotation() const {
        return has_R;
    };
 
    static EGS_AffineTransform *getTransformation(EGS_Input *inp);
 
    static EGS_AffineTransform *getTransformation(vector<EGS_Float> trnsl, vector<EGS_Float> rot);
};
 
#endif