egs_polygon.h

Go to the documentation of this file.

/*
###############################################################################
#
#  EGSnrc egs++ polygon 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
#                   Reid Townson
#
###############################################################################
*/


#ifndef EGS_POLYGON_
#define EGS_POLYGON_

#include "egs_vector.h"
#include "egs_projectors.h"
#include "egs_math.h"
#include "egs_functions.h"

#include <vector>

using namespace std;

class EGS_EXPORT EGS_2DPolygon {

public:

    EGS_2DPolygon(vector<EGS_2DVector> &points, bool Open=false);

    ~EGS_2DPolygon();

    int getN() const {
        return np-1;
    };

    bool isConvex() const {
        return is_convex;
    };

    EGS_2DVector getNormal(int j) const {
        if (j >= 0 && j < np-1) {
            return a[j];
        }
        else {
            return EGS_2DVector();
        }
    };

    EGS_Float hownear(bool in, const EGS_2DVector &x) const {
        if (!open) {
            EGS_Float tperp = veryFar;
            bool do_it = true;
            for (int j=0; j<np-1; j++) {
                EGS_2DVector v(x - p[j]);
                EGS_Float lam = uj[j]*v;
                if (lam >= 0 && lam <= uj[j].length2()) {
                    do_it = false;
                    EGS_Float t = fabs(d[j] - x*a[j]);
                    if (t < tperp) {
                        tperp = t;
                    }
                }
                else if (lam < 0 && do_it) {
                    EGS_Float t = v.length();
                    if (t < tperp) {
                        tperp = t;
                    }
                }
                else {
                    do_it = true;
                }
            }
            return tperp;
        }
        EGS_2DVector v(x - p[0]);
        EGS_Float lam = uj[0]*v;
        EGS_Float tperp;
        bool do_it;
        if (lam <= uj[0].length2()) {
            do_it = false;
            tperp = fabs(d[0] - x*a[0]);
        }
        else {
            do_it = true;
            tperp = veryFar;
        }
        v = x - p[1];
        lam = uj[1]*v;
        if (lam >= 0) {
            EGS_Float t = fabs(d[1] - x*a[1]);
            if (t < tperp) {
                tperp = t;
            }
        }
        else if (do_it) {
            EGS_Float t = v.length();
            if (t < tperp) {
                tperp = t;
            }
        }
        return tperp;
    };

    bool isInside(const EGS_2DVector &x) const {
        if (!open &&
                (x.x<xmin+epsilon || x.x+epsilon>xmax || x.y<ymin+epsilon || x.y+epsilon>ymax)) {
            return false;
        }
        if (is_convex) {
            int nn = open ? np-2 : np-1;
            for (int j=0; j<nn; j++)
                if (!inside(j,x)) {
                    return false;
                }
            return true;
        }
        if (!cpol->isInside(x)) {
            return false;
        }
        for (int j=0; j<ncut; j++) if (cut[j]->isInside(x)) {
                return false;
            }
        return true;
    };

    bool howfar(bool in, const EGS_2DVector &x, const EGS_2DVector &u,
                EGS_Float &t, EGS_2DVector *normal = 0) {
        EGS_Float xp, up;
        bool res = false;
        int nn = open ? np-2 : np-1;
        if (in) {
            int jhit=0;
            for (int j=0; j<nn; j++)  {
                if ((up = u*a[j]) < 0 && (xp = x*a[j])+epsilon > d[j]) {
                    EGS_Float tt = d[j] - xp;
                    if (tt+epsilon >= t*up) {
                        tt /= up;
                        bool ok = is_convex || pc[j];
                        if (!ok) {
                            EGS_Float lam = uj[j]*(x-p[j]+u*tt);
                            if (lam >= 0 && lam < uj[j].length2()) {
                                ok = true;
                            }
                        }
                        if (ok) {
                            t = tt;
                            res = true;
                            jhit = j;
                            //if( normal ) *normal = a[j];
                        }
                    }
                }
            }
            if (res && normal) {
                *normal = a[jhit];
            }
        }
        else {
            int jhit=0;
            if (open) {
                if ((up = u*a[0]) > 0 && (xp = x*a[0]) < d[0]+epsilon) {
                    EGS_Float tt = (d[0] - xp)/up;
                    if (tt <= t+epsilon) {
                        EGS_Float lam = uj[0]*(x-p[0]+u*tt);
                        if (lam < uj[0].length2()) {
                            t = tt;
                            res = true;
                            jhit = 0;
                        }
                    }
                }
                if ((up = u*a[1]) > 0 && (xp = x*a[1]) < d[1]+epsilon) {
                    EGS_Float tt = (d[1] - xp)/up;
                    if (tt <= t+epsilon) {
                        EGS_Float lam = uj[1]*(x-p[1]+u*tt);
                        if (lam > 0) {
                            t = tt;
                            res = true;
                            jhit = 1;
                        }
                    }
                }
            }
            else {
                for (int j=0; j<nn; j++)  {
                    if ((up = u*a[j]) > 0 && (xp = x*a[j]) < d[j]+epsilon) {
                        EGS_Float tt = (d[j] - xp)/up;
                        if (tt <= t+epsilon) {
                            EGS_Float lam = uj[j]*(x-p[j]+u*tt);
                            if (lam >= 0 && lam < uj[j].length2()) {
                                if (tt < 0) {
                                    t = 0;
                                }
                                else {
                                    t = tt;
                                }
                                res = true;
                                jhit = j;
                                //if( normal ) *normal = a[j]*(-1);
                            }
                        }
                    }
                }
            }
            if (res && normal) {
                *normal = a[jhit]*(-1);
            }
        }
        return res;
    };

    EGS_2DVector getPoint(int j) const {
        if (j >= 0 && j < np) {
            return p[j];
        }
        else {
            return EGS_2DVector();
        }
    };


private:

    EGS_2DVector *p;   
    EGS_2DVector *a;   
    EGS_2DVector *uj;  
    EGS_Float    *d;   
    bool         *pc;  

    EGS_Float    xmin, xmax, ymin, ymax;
    int          np;   
    bool         is_convex; 
    bool         open; 
    int          ncut;  
    EGS_2DPolygon  **cut; 
    EGS_2DPolygon  *cpol; 

    static bool checkCCW(const vector<EGS_2DVector> &points);
    bool inside(int j, const EGS_2DVector &x) const {
        if (x*a[j]+epsilon >= d[j]) {
            return true;
        }
        else {
            return false;
        }
    };

};

template <class T>
class EGS_EXPORT EGS_PolygonT {

public:

    EGS_PolygonT(vector<EGS_2DVector> &points, const T &projector,
                 bool Open=false) : p(new EGS_2DPolygon(points,Open)), a(projector) {};

    ~EGS_PolygonT() {
        delete p;
    };

    inline bool isConvex() const {
        return p->isConvex();
    };

    inline int getN() const {
        return p->getN();
    };

    inline EGS_Vector getPoint(int j) const {
        return a.getPoint(p->getPoint(j));
    };

    inline EGS_Vector getNormal() const {
        return a.normal();
    };

    inline EGS_Vector getNormal(const EGS_2DVector &x) const {
        return a.normal(x);
    };

    inline EGS_Vector getNormal(int j) const {
        return a.normal(p->getNormal(j));
    };

    inline EGS_2DVector getProjection(const EGS_Vector &x) const {
        return a.getProjection(x);
    };

    inline EGS_Float distance(const EGS_Vector &x) const {
        return a.distance(x);
    };

    inline bool isInside2D(const EGS_2DVector &xp) const {
        return p->isInside(xp);
    };

    inline bool isInside2D(const EGS_Vector &x) const {
        return p->isInside(a.getProjection(x));
    };

    inline bool isInside(const EGS_Vector &x) const {
        return (a.distance(x) >= 0);
    };

    inline EGS_Float hownear2D(bool in, const EGS_Vector &x) const {
        return p->hownear(in,a.getProjection(x));
    };

    inline EGS_Float hownear(bool in, const EGS_Vector &x) const {
        EGS_2DVector pos(a.getProjection(x));
        EGS_Float t1 = fabs(a.distance(x));
        if (p->isInside(pos)) {
            return t1;
        }
        EGS_Float t2 = p->hownear(true,pos);
        return sqrt(t1*t1+t2*t2);
    };

    inline bool howfar2D(bool in, const EGS_Vector &x, const EGS_Vector &u,
                         EGS_Float &t, EGS_2DVector *normal = 0) const {
        EGS_2DVector dir(a.getProjection(u));
        if (u.length2() < epsilon) {
            return false;
        }
        return p->howfar(in,a.getProjection(x),dir,t,normal);
    };

    inline bool howfar(bool in, const EGS_Vector &x, const EGS_Vector &u,
                       EGS_Float &t) const {
        EGS_Float up = a*u;
        if ((in && up >= 0)  || (!in && up <= 0)) {
            return false;
        }
        EGS_Float tt = -a.distance(x)/up;
        if (tt <= t+epsilon) {
            EGS_Vector xp(x + u*tt);
            if (p->isInside(a.getProjection(xp))) {
                t = tt;
                return true;
            }
        }
        return false;
    };

    const string &getType() const {
        return a.getType();
    };

private:

    EGS_2DPolygon *p;
    T             a;

};

typedef EGS_PolygonT<EGS_XProjector> EGS_PolygonYZ;
typedef EGS_PolygonT<EGS_YProjector> EGS_PolygonXZ;
typedef EGS_PolygonT<EGS_ZProjector> EGS_PolygonXY;
typedef EGS_PolygonT<EGS_Projector>  EGS_Polygon;

EGS_EXPORT EGS_Polygon *makePolygon(const vector<EGS_Vector> &points);

#endif