egs_dynamic_shape.cpp

Go to the documentation of this file.

/*
###############################################################################
#
#  EGSnrc egs++ dynamic shape
#  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:          Reid Townson
#
#  Contributors:
#
###############################################################################
*/
 
 
#include "egs_dynamic_shape.h"
#include "egs_input.h"
#include "egs_functions.h"
#include <sstream>
 
static bool EGS_DYNAMIC_SHAPE_LOCAL inputSet = false;
static shared_ptr<EGS_BlockInput> EGS_DYNAMIC_SHAPE_LOCAL shapeBlockInput = make_shared<EGS_BlockInput>("shape");
 
extern "C" {
 
    static void setInputs() {
        inputSet = true;
 
        setShapeInputs(shapeBlockInput);
        shapeBlockInput->getSingleInput("library")->setValues({"egs_dynamic_shape"});
 
        auto shapePtr = shapeBlockInput->addBlockInput("shape");
        setShapeInputs(shapePtr);
 
        shared_ptr<EGS_BlockInput> motionBlock = shapeBlockInput->addBlockInput("motion", true);
        motionBlock->addSingleInput("control point", true, "Parameters to define motion: timeIndex xtrans ytrans ztrans xrot yrot zrot");
    }
 
    EGS_DYNAMIC_SHAPE_EXPORT string getExample() {
        string example;
        example = {
            R"(
    # Example of egs_dynamic_shape
    #:start shape:
        library = egs_dynamic_shape
        :start shape:
            definition of the shape to be 'dynamic'
        :stop shape:
        :start motion: # units of cm and degrees
            control point = timeIndex(1) xtrans(1) ytrans(1) ztrans(1) xrot(1) yrot(1) zrot(1)
            control point = timeIndex(2) xtrans(2) ytrans(2) ztrans(2) xrot(2) yrot(2) zrot(2)
            .
            .
            .
            control point = timeIndex(N) xtrans(N) ytrans(N) ztrans(N) xrot(N) yrot(N) zrot(N)
        :stop motion:
    :stop shape:
)"};
        return example;
    }
 
    EGS_DYNAMIC_SHAPE_EXPORT shared_ptr<EGS_BlockInput> getInputs() {
        if(!inputSet) {
            setInputs();
        }
        return shapeBlockInput;
    }
 
    EGS_DYNAMIC_SHAPE_EXPORT EGS_BaseShape *createShape(EGS_Input *input,
            EGS_ObjectFactory *f) {
        if (!input) {
            egsWarning("createShape(dynamic shape): null input?\n");
            return 0;
        }
        EGS_Input *ishape = input->takeInputItem("shape",false);
        EGS_BaseShape *shape=0;;
        if (ishape) {
            shape = EGS_BaseShape::createShape(ishape);
            delete ishape;
        }
        if (!shape) {
            string shape_name;
            int err = input->getInput("shape name",shape_name);
            if (err) {
                egsWarning("createShape(dynamic shape): no inline shape definition"
                           " and no 'shape name' keyword\n");
                return 0;
            }
            shape = EGS_BaseShape::getShape(shape_name);
            if (!shape) {
                egsWarning("createShape(dynamic shape): no shape named %s "
                           "exists\n",shape_name.c_str());
                return 0;
            }
        }
 
        // Now getting motion information for dynamic component
        EGS_Input *dyninp = input->takeInputItem("motion");
 
        if (dyninp) {
            EGS_DynamicShape *s = new EGS_DynamicShape(shape, dyninp, "", f);
            s->setName(input);
            return s;
        }
        else {
            egsWarning("EGS_DynamicShape: no control points input.\n");
            return 0;
        }
    }
 
    void EGS_DynamicShape::buildDynamicShape(EGS_Input *dyninp) {
        stringstream itos;
        ncpts=0;
        vector<EGS_Float> point;
        EGS_ControlPoint cpt;
        int err;
        int icpts=1;
        itos << icpts;
        string inputTag = "control point";
        string inputTag_backCompat = "control point " + itos.str();
        EGS_Input *currentInput;
        int rotsize=0;
 
        // Control points read one by one from motion block in dynamic geometry definition, and saved as a vector to points
        while (true) {
            currentInput = dyninp->takeInputItem(inputTag);
            if (!currentInput || currentInput->getInput(inputTag, point)) {
                currentInput = dyninp->takeInputItem(inputTag_backCompat);
                if (!currentInput || currentInput->getInput(inputTag_backCompat, point)) {
                    delete currentInput;
                    break;
                }
            }
            delete currentInput;
 
            // Checking the size to make sure it is a valid control point input
            if (point.size()!=6 && point.size()!=7) {
                egsFatal("EGS_DynamicShape: Control point %i must specify either 6 or 7 values.\n",icpts);
            }
            else {
                if (ncpts == 0) {
                    rotsize=point.size();// Variable to make sure all control point definitions have consistent formats
                }
                if (ncpts>0 && point[0] < cpts[ncpts-1].time) {// Make sure each time index is larger than the last
                    egsFatal("EGS_DynamicShape: Time index of control point %i < time index of control point %i\n",icpts,ncpts);
                }
                else if (point[0] < 0.) {// Checks that time index is valid (larger than zero)
                    egsFatal("EGS_DynamicShape: Time index of control point %i < 0.0\n",icpts);
                }
                else if (ncpts> 0 && point.size() != rotsize) { // Checks that control point formats follow the first
                    egsFatal("EGS_DynamicShape: Rotation definition inconsistent \n");
                }
                else {
                    ncpts++;
 
                    if (ncpts ==1 && point[0] > 0.0) {
                        egsWarning("EGS_DynamicShape: Time index of control point 1 > 0.0.  This will generate many warning messages.\n");
                    }
 
                    vector<EGS_Float> T_vect;// Vector storing translation information
                    vector<EGS_Float> R_vect;// Vector storing rotation information
                    // Add translation coordinates to the translation vector
                    T_vect.push_back(point[1]);
                    T_vect.push_back(point[2]);
                    T_vect.push_back(point[3]);
                    // Add rotation coordinates to rotation vector (cases differentiate cpt formats, 6 is 2 rotation parameters, 7 is 3)
                    // In each case vector order is determined by format of EGS_rotationMatrix constructor (EGS_Float arguments)
                    if (point.size()==6) {
                        R_vect.push_back(point[6]);//rotation about z
                        R_vect.push_back(point[4]);//rotation about x
                    }//2 number case
                    if (point.size()==7) {
                        R_vect.push_back(point[4]);//rotation about x
                        R_vect.push_back(point[5]);//rotation about y
                        R_vect.push_back(point[6]);//rotation about z
                    }//3 number case
 
                    // Add it to the vector of control points
                    cpt.time = point[0];
                    cpt.trnsl = T_vect;
                    cpt.rot = R_vect;
                    cpts.push_back(cpt);// Add created control point to list of control points
                    icpts++;
                    itos.str("");
                    itos << icpts;
                    inputTag_backCompat = "control point " + itos.str();// Define next control point i string for getInput in while condition
                }
            }
        }
        if (ncpts<=1) {
            egsFatal("EGS_DynamicShape: not enough or missing control points.\n");
        }
        if (cpts[ncpts-1].time == 0.0) {
            egsFatal("EGS_DynamicShape: time index of last control point = 0.  Something's wrong.\n");
        }
        else {
            // Normalize time index to max. value
            for (int i=0; i<=ncpts-1; i++) {// I changed the normalization here (and in dynamic source) to have <= instead of <. Otherwise last cpt never gets normalized, which is not an issue for final cpt time>1 but is an issue if final cpt time<1.
                cpts[i].time /= cpts[ncpts-1].time;
            }
        }
        updatePosition(0); // Sets position to initial time in egs_view upon opening
    };
 
    void EGS_DynamicShape::getNextShapePosition(EGS_RandomGenerator *rndm) {
        int errg = 1;
        EGS_ControlPoint gipt;
 
        // Here get source from activeapplication in order to extract time
        EGS_Application *app = EGS_Application::activeApplication();
        while (errg) {
            // Gets time if it's already set (otherwise gives -1).
            ptime = app->getTimeIndex();
 
            if (ptime<0) {
                // If no time is given by the source the shape will randomly sample from 0 to 1.
                ptime = rndm->getUniform();
 
                // Set randomly sampled time index for all objects in the simulation
                app->setTimeIndex(ptime);
            }
 
            // Now run the get coord method that will sample the cpt given to find the transformation that will be applied for the current history
            errg = getCoord(ptime,gipt);
        }
 
        // Create and set the current shape transformation using the sampled coordinates from getCoord. This is where overloaded EGS_AffineTransform is used
        EGS_AffineTransform *tDG = EGS_AffineTransform::getTransformation(gipt.trnsl, gipt.rot);
        shape->setTransformation(tDG);
        delete tDG;
 
        // Call getNextShapePosition on base shape in case there are lower level dynamic shapes
        shape->getNextShapePosition(rndm);
    };
 
    int EGS_DynamicShape::getCoord(EGS_Float rand, EGS_ControlPoint &gipt) {
        int iindex=0;
        int i;
 
        // The following loop determines which 2 control points the current time index falls between
        for (i=0; i<ncpts; i++) {
            if (rand < cpts[i].time-epsilon) { // If the control point i time index is larger than the current time index, we know this is the upper bound control point (represented by iindex)
                iindex = i;
                break;
            }
        }
 
        if (i==ncpts) {
            egsWarning("EGS_DynamicShape: could not locate control point.\n");
            return 1;
        }
 
        // Below 3 vectors are defined, a vector containing the lower bound translation coordinates, a vector containing the upper bound translation coordinates, and a vector for the sampled translation coordinates
        vector<EGS_Float> translation_LB=cpts[iindex-1].trnsl;
        vector<EGS_Float> translation_UB=cpts[iindex].trnsl;
        vector<EGS_Float> translation_samp;
 
        // The following is a factor (between 0 and 1) used for sampling. Essentially, it represents the fractional position within the interval
        EGS_Float factor = (rand-cpts[iindex-1].time)/(cpts[iindex].time-cpts[iindex-1].time);
 
        // Translations are given as the lower bound plus the length of the interval multiplied by the fractional position factor. So its essentially lower bound + n% of the interval length
        translation_samp.push_back(translation_LB[0]+(translation_UB[0]-translation_LB[0])*factor);
        translation_samp.push_back(translation_LB[1]+(translation_UB[1]-translation_LB[1])*factor);
        translation_samp.push_back(translation_LB[2]+(translation_UB[2]-translation_LB[2])*factor);
        // Update the translation coordinates in the current state control point object
        gipt.trnsl=translation_samp;
 
        // Again 3 vectors defined, lower bound, upper bound, and sampled, this time for rotation coordinates
        vector<EGS_Float> rotation_LB=cpts[iindex-1].rot;
        vector<EGS_Float> rotation_UB=cpts[iindex].rot;
        vector<EGS_Float> rotation_samp;
 
        // Now set rotations. Current coordinates computed as before (lowerbound+n% of interval length), but now we also convert degrees to radians
        // These must be done case by case (since order of arguments matters)
        if (cpts[iindex].rot.size()==2) {
            rotation_samp.push_back((rotation_LB[0]+(rotation_UB[0]-rotation_LB[0])*factor)*(M_PI/180));
            rotation_samp.push_back((rotation_LB[1]+(rotation_UB[1]-rotation_LB[1])*factor)*(M_PI/180));
        }
        else if (cpts[iindex].rot.size()==3) {
            rotation_samp.push_back((rotation_LB[0]+(rotation_UB[0]-rotation_LB[0])*factor)*(M_PI/180));
            rotation_samp.push_back((rotation_LB[1]+(rotation_UB[1]-rotation_LB[1])*factor)*(M_PI/180));
            rotation_samp.push_back((rotation_LB[2]+(rotation_UB[2]-rotation_LB[2])*factor)*(M_PI/180));
        }
        else {
            egsWarning("EGS_DynamicShape: Invalid number of rotation parameters\n");
        }
        gipt.rot=rotation_samp;
 
        // Update the rotation coordinates in the current state control point object
        return 0;
    };
 
}  // extern "C"