Math/GenVector/Plane3D.h

0001 // @(#)root/mathcore:$Id$
0002 // Authors: L. Moneta    12/2005
0003
0004 /**********************************************************************
0005  *                                                                    *
0006  * Copyright (c) 2005 , LCG ROOT MathLib Team                         *
0007  *                                                                    *
0008  *                                                                    *
0009  **********************************************************************/
0010
0011 // Header file for class LorentzVector
0012 //
0013 // Created by:    moneta   at Fri Dec 02   2005
0014 //
0015 // Last update: $Id$
0016 //
0017 #ifndef ROOT_Math_GenVector_Plane3D
0018 #define ROOT_Math_GenVector_Plane3D  1
0019
0020 #include <type_traits>
0021
0022 #include "Math/GenVector/DisplacementVector3D.h"
0023 #include "Math/GenVector/PositionVector3D.h"
0024
0025
0026
0027 namespace ROOT {
0028
0029 namespace Math {
0030
0031 namespace Impl {
0032
0033 //_______________________________________________________________________________
0034 /**
0035    Class describing a geometrical plane in 3 dimensions.
0036    A Plane3D is a 2 dimensional surface spanned by two linearly independent vectors.
0037    The plane is described by the equation
0038    \f$ a*x + b*y + c*z + d = 0 \f$ where (a,b,c) are the components of the
0039    normal vector to the plane \f$ n = (a,b,c)  \f$ and \f$ d = - n \dot x \f$, where x is any point
0040    belonging to plane.
0041    More information on the mathematics describing a plane in 3D is available on
0042    <A HREF=http://mathworld.wolfram.com/Plane.html>MathWord</A>.
0043    The Plane3D class contains the 4 scalar values in T which represent the
0044    four coefficients, fA, fB, fC, fD. fA, fB, fC are the normal components normalized to 1,
0045    i.e. fA**2 + fB**2 + fC**2 = 1
0046
0047    @ingroup GenVector
0048
0049    @sa Overview of the @ref GenVector "physics vector library"
0050 */
0051
0052 template <typename T = double>
0053 class Plane3D {
0054
0055 public:
0056    // ------ ctors ------
0057
0058    typedef T Scalar;
0059
0060    typedef DisplacementVector3D<Cartesian3D<T>, DefaultCoordinateSystemTag> Vector;
0061    typedef PositionVector3D<Cartesian3D<T>, DefaultCoordinateSystemTag>     Point;
0062
0063    /**
0064       default constructor create plane z = 0
0065    */
0066    Plane3D() : fA(0), fB(0), fC(1), fD(0) {}
0067
0068    /**
0069     generic constructors from the four scalar values describing the plane
0070     according to the equation ax + by + cz + d = 0
0071       \param a scalar value
0072       \param b scalar value
0073       \param c scalar value
0074       \param d sxcalar value
0075    */
0076    Plane3D(const Scalar &a, const Scalar &b, const Scalar &c, const Scalar &d) : fA(a), fB(b), fC(c), fD(d)
0077    {
0078       // renormalize a,b,c to unit
0079       Normalize();
0080    }
0081
0082    /**
0083     constructor a Plane3D from a normal vector and a point coplanar to the plane
0084     \param n normal expressed as a ROOT::Math::DisplacementVector3D<Cartesian3D<T> >
0085     \param p point  expressed as a  ROOT::Math::PositionVector3D<Cartesian3D<T> >
0086    */
0087    Plane3D(const Vector &n, const Point &p) { BuildFromVecAndPoint(n, p); }
0088
0089    /**
0090     Construct from a generic DisplacementVector3D (normal vector) and PositionVector3D (point coplanar to
0091     the plane)
0092     \param n normal expressed as a generic ROOT::Math::DisplacementVector3D
0093     \param p point  expressed as a generic ROOT::Math::PositionVector3D
0094    */
0095    template <class T1, class T2, class U>
0096    Plane3D(const DisplacementVector3D<T1, U> &n, const PositionVector3D<T2, U> &p)
0097    {
0098       BuildFromVecAndPoint(Vector(n), Point(p));
0099    }
0100
0101    /**
0102     constructor from three Cartesian point belonging to the plane
0103     \param p1 point1  expressed as a generic ROOT::Math::PositionVector3D
0104     \param p2 point2  expressed as a generic ROOT::Math::PositionVector3D
0105     \param p3 point3  expressed as a generic ROOT::Math::PositionVector3D
0106    */
0107    Plane3D(const Point &p1, const Point &p2, const Point &p3) { BuildFrom3Points(p1, p2, p3); }
0108
0109    /**
0110     constructor from three generic point belonging to the plane
0111     \param p1 point1 expressed as  ROOT::Math::DisplacementVector3D<Cartesian3D<T> >
0112     \param p2 point2 expressed as  ROOT::Math::DisplacementVector3D<Cartesian3D<T> >
0113     \param p3 point3 expressed as  ROOT::Math::DisplacementVector3D<Cartesian3D<T> >
0114    */
0115    template <class T1, class T2, class T3, class U>
0116    Plane3D(const PositionVector3D<T1, U> &p1, const PositionVector3D<T2, U> &p2, const PositionVector3D<T3, U> &p3)
0117    {
0118       BuildFrom3Points(Point(p1.X(), p1.Y(), p1.Z()), Point(p2.X(), p2.Y(), p2.Z()), Point(p3.X(), p3.Y(), p3.Z()));
0119    }
0120
0121    // compiler-generated copy ctor and dtor are fine.
0122    Plane3D(const Plane3D &) = default;
0123
0124    // ------ assignment ------
0125
0126    /**
0127       Assignment operator from other Plane3D class
0128    */
0129    Plane3D &operator=(const Plane3D &) = default;
0130
0131    /**
0132       Return the a coefficient of the plane equation \f$ a*x + b*y + c*z + d = 0 \f$. It is also the
0133       x-component of the vector perpendicular to the plane.
0134    */
0135    Scalar A() const { return fA; }
0136
0137    /**
0138       Return the b coefficient of the plane equation \f$ a*x + b*y + c*z + d = 0 \f$. It is also the
0139       y-component of the vector perpendicular to the plane
0140    */
0141    Scalar B() const { return fB; }
0142
0143    /**
0144       Return the c coefficient of the plane equation \f$ a*x + b*y + c*z + d = 0 \f$. It is also the
0145       z-component of the vector perpendicular to the plane
0146    */
0147    Scalar C() const { return fC; }
0148
0149    /**
0150       Return the d coefficient of the plane equation \f$ a*x + b*y + c*z + d = 0 \f$. It is also
0151       the distance from the origin (HesseDistance)
0152    */
0153    Scalar D() const { return fD; }
0154
0155    /**
0156       Return normal vector to the plane as Cartesian DisplacementVector
0157    */
0158    Vector Normal() const { return Vector(fA, fB, fC); }
0159
0160    /**
0161     Return the Hesse Distance (distance from the origin) of the plane or
0162     the d coefficient expressed in normalize form
0163    */
0164    Scalar HesseDistance() const { return fD; }
0165
0166    /**
0167     Return the signed distance to a Point.
0168     The distance is signed positive if the Point is in the same side of the
0169     normal vector to the plane.
0170     \param p Point expressed in Cartesian Coordinates
0171     */
0172    Scalar Distance(const Point &p) const { return fA * p.X() + fB * p.Y() + fC * p.Z() + fD; }
0173
0174    /**
0175     Return the distance to a Point described with generic coordinates
0176     \param p Point expressed as generic ROOT::Math::PositionVector3D
0177     */
0178    template <class T1, class U>
0179    Scalar Distance(const PositionVector3D<T1, U> &p) const
0180    {
0181       return Distance(Point(p.X(), p.Y(), p.Z()));
0182    }
0183
0184    /**
0185     Return the projection of a Cartesian point to a plane
0186     \param p Point expressed as PositionVector3D<Cartesian3D<T> >
0187     */
0188    Point ProjectOntoPlane(const Point &p) const
0189    {
0190       const Scalar d = Distance(p);
0191       return XYZPoint(p.X() - fA * d, p.Y() - fB * d, p.Z() - fC * d);
0192    }
0193
0194    /**
0195     Return the projection of a point to a plane
0196     \param p Point expressed as generic ROOT::Math::PositionVector3D
0197     */
0198    template <class T1, class U>
0199    PositionVector3D<T1, U> ProjectOntoPlane(const PositionVector3D<T1, U> &p) const
0200    {
0201       const Point pxyz = ProjectOntoPlane(Point(p.X(), p.Y(), p.Z()));
0202       return PositionVector3D<T, U>(pxyz.X(), pxyz.Y(), pxyz.Z());
0203    }
0204
0205    // ------------------- Equality -----------------
0206
0207    /**
0208       Exact equality
0209    */
0210    bool operator==(const Plane3D &rhs) const { return (fA == rhs.fA && fB == rhs.fB && fC == rhs.fC && fD == rhs.fD); }
0211    bool operator!=(const Plane3D &rhs) const { return !(operator==(rhs)); }
0212
0213 protected:
0214    /**
0215       Normalize the normal (a,b,c) plane components
0216    */
0217    template <typename SCALAR = T, typename std::enable_if<std::is_arithmetic<SCALAR>::value>::type * = nullptr>
0218    void Normalize()
0219    {
0220       // normalize the plane
0221       using std::sqrt;
0222       const SCALAR s = sqrt(fA * fA + fB * fB + fC * fC);
0223       // what to do if s = 0 ?
0224       if (s == SCALAR(0)) {
0225          fD = SCALAR(0);
0226       } else {
0227          const SCALAR w = Scalar(1) / s;
0228          fA *= w;
0229          fB *= w;
0230          fC *= w;
0231          fD *= w;
0232       }
0233    }
0234
0235    /**
0236      Normalize the normal (a,b,c) plane components
0237    */
0238    template <typename SCALAR = T, typename std::enable_if<!std::is_arithmetic<SCALAR>::value>::type * = nullptr>
0239    void Normalize()
0240    {
0241       // normalize the plane
0242       using std::sqrt;
0243       SCALAR s = sqrt(fA * fA + fB * fB + fC * fC);
0244       // what to do if s = 0 ?
0245       const auto m = (s == SCALAR(0));
0246       // set zero entries to 1 in the vector to avoid /0 later on
0247       s(m)           = SCALAR(1);
0248       fD(m)          = SCALAR(0);
0249       const SCALAR w = SCALAR(1) / s;
0250       fA *= w;
0251       fB *= w;
0252       fC *= w;
0253       fD *= w;
0254    }
0255
0256 private:
0257    // internal method to construct class from a vector and a point
0258    void BuildFromVecAndPoint(const Vector &n, const Point &p)
0259    {
0260       // build from a normal vector and a point
0261       fA = n.X();
0262       fB = n.Y();
0263       fC = n.Z();
0264       fD = -n.Dot(p);
0265       Normalize();
0266    }
0267
0268    // internal method to construct class from 3 points
0269    void BuildFrom3Points(const Point &p1, const Point &p2, const Point &p3)
0270    {
0271       // plane from three points
0272       // normal is (x3-x1) cross (x2 -x1)
0273       const Vector n = (p2 - p1).Cross(p3 - p1);
0274       fA             = n.X();
0275       fB             = n.Y();
0276       fC             = n.Z();
0277       fD             = -n.Dot(p1);
0278       Normalize();
0279    }
0280
0281    // plane data members the four scalar which  satisfies fA*x + fB*y + fC*z + fD = 0
0282    // for every point (x,y,z) belonging to the plane.
0283    // fA**2 + fB**2 + fC** =1 plane is stored in normalized form
0284    Scalar fA;
0285    Scalar fB;
0286    Scalar fC;
0287    Scalar fD;
0288
0289    };  // Plane3D<>
0290
0291    /**
0292       Stream Output and Input
0293    */
0294    // TODO - I/O should be put in the manipulator form
0295    template <typename T>
0296    std::ostream &operator<<(std::ostream &os, const Plane3D<T> &p)
0297    {
0298       os << "\n"
0299          << p.Normal().X() << "  " << p.Normal().Y() << "  " << p.Normal().Z() << "  " << p.HesseDistance() << "\n";
0300       return os;
0301    }
0302
0303    } // end namespace Impl
0304
0305    // typedefs for double and float versions
0306    typedef Impl::Plane3D<double> Plane3D;
0307    typedef Impl::Plane3D<float>  Plane3DF;
0308
0309 } // end namespace Math
0310
0311 } // end namespace ROOT
0312
0313
0314 #endif