Geom.h

Go to the documentation of this file.

#ifndef GEOM_H_INCLUDED
#define GEOM_H_INCLUDED
/*
 * Copyright 2011-2016 Universiteit Antwerpen
 *
 * Licensed under the EUPL, Version 1.1 or  as soon they will be approved by
 * the European Commission - subsequent versions of the EUPL (the "Licence");
 * You may not use this work except in compliance with the Licence.
 * You may obtain a copy of the Licence at: http://ec.europa.eu/idabc/eupl5
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the Licence is distributed on an "AS IS" basis,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the Licence for the specific language governing
 * permissions and limitations under the Licence.
 */
 /*
 *  Some of the procedures are adapted from Dan Sunday's  Geometry Algorithms
 *  code. See http://geomalgorithms.com/index.html. When this is the case,
 *  the appropriate url is mentioned in comment.
 *
 *  Copyright 2000 softSurfer, 2012 Dan Sunday
 *  This code may be freely used, distributed and modified for any purpose
 *  providing that this copyright notice is included with it.
 *  SoftSurfer makes no warranty for this code, and cannot be held
 *  liable for any real or imagined damage resulting from its use.
 *  Users of this code must verify correctness for their application.
 */
#include "array3.h"
#include "bio/Node.h"
#include "util/container/CircularIterator.h"
#include "util/container/ConstCircularIterator.h"

#include <array>
#include <tuple>

namespace SimPT_Sim {

using namespace SimPT_Sim::Util;

class Geom
{
public:
        static std::tuple<double, double, double>
        CircumCircle(double x1, double y1, double x2, double y2, double x3, double y3);

        static inline std::array<double, 3> Intersection(
                        const std::array<double, 3>& v1, const std::array<double, 3>& v2,
                        const std::array<double, 3>& v3, const std::array<double, 3>& v4)
        {
                const double denom  = (v4[1] - v3[1]) * (v2[0] - v1[0]) - (v4[0] - v3[0]) * (v2[1] - v1[1]);
                const double ua     = ((v4[0] - v3[0]) * (v1[1] - v3[1]) - (v4[1] - v3[1]) * (v1[0] - v3[0])) / denom;
                return v1 + ua * (v2 - v1);
        }

        static inline int
        isLeft(const std::array<double, 3>& P0, const std::array<double, 3>& P1, const std::array<double, 3>& P2 )
        {
            return (P1[0] - P0[0])*(P2[1] - P0[1]) - (P2[0] -  P0[0])*(P1[1] - P0[1]);
        }

        static inline int
        wn_PnPoly(const std::array<double, 3>& p, const std::vector<Node*>& nodes)
        {
                using std::array;
                using namespace SimPT_Sim::Container;

                int  wn = 0;                                              // winding number counter
                auto it = make_const_circular(nodes);

                // loop through all edges of the polygon
                do {   // edge from e_a to e_b
                        array<double, 3> e_a  = *(*it);
                        array<double, 3> e_b  = *(*next(it));

                        if (e_a[1] <= p[1]) {                             // start y <= p.y
                                if (e_b[1]  > p[1]) {                     // an upward crossing
                                        if (isLeft( e_a, e_b, p) > 0) {   // p left of  edge
                                                ++wn;                     // have valid up intersect
                                        }
                                }
                        } else {                                          // start y > p.y (no test needed)
                                if (e_b[1]  <= p[1]) {                    // a downward crossing
                                        if (isLeft( e_a, e_b, p) < 0) {   // p right of  edge
                                                --wn;                     // have valid down intersect
                                        }
                                }
                        }

                } while(++it != nodes.begin());
                return wn;
        }

        };
}

#endif // end-of-include-guard