Decomposer.h

Go to the documentation of this file.

#ifndef LARCOREALG_GEOMETRY_DECOMPOSER_H
#define LARCOREALG_GEOMETRY_DECOMPOSER_H

// LArSoft libraries
#include "larcorealg/Geometry/geo_vectors_utils.h" // geo::vect

// C/C++ standard libraries
#include <cmath> // std::abs()
#include <utility> // std::move()
#include <type_traits> // std::declval()


namespace geo {
  
  // --- BEGIN Decomposition objects -------------------------------------------
  
  template <typename Vector>
  class PlaneBase {
      public:
    using Vector_t = Vector; 
    
    static constexpr double RoundingTol = 1e-4;
    
    PlaneBase(Vector_t const& main, Vector_t const& secondary)
      : fMain(PastorizeUnitVector(main))
      , fSecondary(PastorizeUnitVector(secondary))
      , fNormal(ComputeNormal())
      {}
    
    Vector_t const& MainDir() const { return fMain; }
    
    Vector_t const& SecondaryDir() const { return fSecondary; }
    
    Vector_t const& NormalDir() const { return fNormal; }
    
    void SetMainDir(Vector_t const& dir) { fMain = dir; ResetNormal(); }
    
    void SetSecondaryDir(Vector_t const& dir)
      { fSecondary = dir; ResetNormal(); }
    
    
    static Vector_t PastorizeUnitVector(Vector_t dir)
      { return geo::vect::rounded01(geo::vect::normalize(dir), RoundingTol); }
    
      private:
    Vector_t fMain;      
    Vector_t fSecondary; 
    Vector_t fNormal;    
    
    Vector_t ComputeNormal() const
      { return PastorizeUnitVector(MainDir().Cross(SecondaryDir())); }
    
    void ResetNormal() { fNormal = ComputeNormal(); }
    
  }; // class PlaneBase<>
  
  
  template <typename Vector, typename Point = Vector>
  class AffinePlaneBase {
    using PlaneBase_t = PlaneBase<Vector>;
    
      public:
    using Vector_t = typename PlaneBase_t::Vector_t; 
    using Point_t  = Point; 
    
    AffinePlaneBase
      (Point_t const& origin, Vector_t const& main, Vector_t const& secondary)
      : fOrigin(origin)
      , fBase(main, secondary)
      {}
    
    Vector_t const& MainDir() const { return fBase.MainDir(); }
    
    Vector_t const& SecondaryDir() const { return fBase.SecondaryDir(); }
    
    Vector_t const& NormalDir() const { return fBase.NormalDir(); }
    
    Point_t Origin() const { return fOrigin; }
    
    Vector_t ToVector(Point_t const& point) const { return point - Origin(); }
    
    
    void SetOrigin(Point_t const& point) { fOrigin = point; }
    
    void SetMainDir(Vector_t const& dir) { fBase.SetMainDir(dir); }
    
    void SetSecondaryDir(Vector_t const& dir) { fBase.SetSecondaryDir(dir); }
    
      private:
    
    Point_t     fOrigin; 
    PlaneBase_t fBase;   
    
  }; // class AffinePlaneBase<>
  
  
  template <typename ProjVector>
  struct DecomposedVector {
    
    using Projection_t = ProjVector; 
    
    using Distance_t = decltype(geo::vect::mag2(std::declval<Projection_t>()));
    
    
    Projection_t projection; 
    Distance_t   distance;   
    
    
    DecomposedVector() = default;
    
    DecomposedVector(Projection_t const& projection, Distance_t distance)
      : projection(projection), distance(distance) {}
    
    DecomposedVector(Distance_t distance, Projection_t const& projection)
      : projection(projection), distance(distance) {}
    
  }; // struct DecomposedVector
  
  
  
  template <typename Vector, typename Point, typename ProjVector>
  class PlaneDecomposer {
    
      public:
    
    using AffinePlaneBase_t = AffinePlaneBase<Vector, Point>;
    
    using DecomposedVector_t = DecomposedVector<ProjVector>;
    
    using Point_t = typename AffinePlaneBase_t::Point_t;
    
    using Vector_t = typename AffinePlaneBase_t::Vector_t;
    
    using Projection_t = typename DecomposedVector_t::Projection_t;
    
    using Distance_t = typename DecomposedVector_t::Distance_t;
    
    
    PlaneDecomposer()
      : fPlaneBase(
        { 0.0, 0.0, 0.0 }, // origin
        { 1.0, 0.0, 0.0 }, // x axis
        { 0.0, 1.0, 0.0 }  // y axis
        )
      {}
    
    PlaneDecomposer(AffinePlaneBase_t&& base): fPlaneBase(std::move(base)) {}
    
    PlaneDecomposer(AffinePlaneBase_t const& base): fPlaneBase(base) {}
    
    
    void SetBase(AffinePlaneBase_t&& base) { fPlaneBase = std::move(base); }
    
    void SetBase(AffinePlaneBase_t const& base) { fPlaneBase = base; }
    
    void SetOrigin(Point_t const& point) { fPlaneBase.SetOrigin(point); }
    
    void SetMainDir(Vector_t const& dir) { fPlaneBase.SetMainDir(dir); }
    
    void SetSecondaryDir(Vector_t const& dir)
      { fPlaneBase.SetSecondaryDir(dir); }
    
    
    
    Point_t ReferencePoint() const { return Base().Origin(); }
    
    Vector_t const& MainDir() const { return Base().MainDir(); }
    
    Vector_t const& SecondaryDir() const { return Base().SecondaryDir(); }
    
    AffinePlaneBase_t const& Base() const { return fPlaneBase; }
    
    
    
    
    auto MainComponent(Projection_t const& v) const { return v.X(); }
    
    auto SecondaryComponent(Projection_t const& v) const { return v.Y(); }
    
    
    
    
    auto PointMainComponent(Point_t const& point) const
      { return VectorMainComponent(Base().ToVector(point)); }
    
    auto PointSecondaryComponent(Point_t const& point) const
      { return VectorSecondaryComponent(Base().ToVector(point)); }
    
    Projection_t PointProjection(Point_t const& point) const
      { return VectorProjection(Base().ToVector(point)); }
    
    
    auto VectorMainComponent(Vector_t const& v) const
      { return geo::vect::dot(v, MainDir()); }
    
    auto VectorSecondaryComponent(Vector_t const& v) const
      { return geo::vect::dot(v, SecondaryDir()); }
    
    Projection_t VectorProjection(Vector_t const& v) const
      { return { VectorMainComponent(v), VectorSecondaryComponent(v) }; }
    
    double Angle(Vector_t const& v) const
      {
        double const a
          = std::atan2(VectorSecondaryComponent(v), VectorMainComponent(v));
        return (a >= M_PI)? -M_PI: a;
      }
    
    
    
    
    Vector_t ComposeVector(Projection_t const& projection) const
      {
        return MainComponent(projection) * MainDir()
          + SecondaryComponent(projection) * SecondaryDir()
          ;
      }
    
    Point_t ComposePoint(Projection_t const& projection) const
      { return ReferencePoint() + ComposeVector(projection); }
    
    
      private:
    AffinePlaneBase_t fPlaneBase; 
    
  }; // class PlaneDecomposer<>
  
  
  
  template <typename Vector, typename Point, typename ProjVector>
  class Decomposer {
    
    using PlaneDecomposer_t = PlaneDecomposer<Vector, Point, ProjVector>;
    
    PlaneDecomposer_t fPlaneDecomp; 
    
    PlaneDecomposer_t const& Plane() const { return fPlaneDecomp; }
    
      public:
    using Point_t = typename PlaneDecomposer_t::Point_t;
    
    using Vector_t = typename PlaneDecomposer_t::Vector_t;
    
    using Projection_t = typename PlaneDecomposer_t::Projection_t;
    
    using Distance_t = typename PlaneDecomposer_t::Distance_t;
    
    using DecomposedVector_t  = typename PlaneDecomposer_t::DecomposedVector_t;
    
    using AffinePlaneBase_t = typename PlaneDecomposer_t::AffinePlaneBase_t;
    
    
    Decomposer() = default;
    
    Decomposer(AffinePlaneBase_t&& base): fPlaneDecomp(std::move(base)) {}
    
    Decomposer(AffinePlaneBase_t const& base): fPlaneDecomp(base) {}
    
    
    void SetBase(AffinePlaneBase_t&& base)
      { fPlaneDecomp.SetBase(std::move(base)); }
    
    void SetBase(AffinePlaneBase_t const& base)
      { fPlaneDecomp.SetBase(base); }
    
    void SetOrigin(Point_t const& point) { fPlaneDecomp.SetOrigin(point); }
    
    void SetMainDir(Vector_t const& dir) { fPlaneDecomp.SetMainDir(dir); }
    
    void SetSecondaryDir(Vector_t const& dir)
      { fPlaneDecomp.SetSecondaryDir(dir); }
    
    
    
    
    Point_t ReferencePoint() const { return Plane().ReferencePoint(); }
    
    AffinePlaneBase_t const& Base() const { return Plane().Base(); }
    
    Vector_t const& MainDir() const { return Plane().MainDir(); }
    
    Vector_t const& SecondaryDir() const { return Plane().SecondaryDir(); }
    
    Vector_t const& NormalDir() const { return Base().NormalDir(); }
    
    
    
    
    auto PointMainComponent(Point_t const& point) const
      { return VectorMainComponent(Base().ToVector(point)); }
    
    auto PointSecondaryComponent(Point_t const& point) const
      { return VectorSecondaryComponent(Base().ToVector(point)); }
    
    auto PointNormalComponent(Point_t const& point) const
      { return VectorNormalComponent(Base().ToVector(point)); }
    
    Projection_t ProjectPointOnPlane(Point_t const& point) const
      { return Plane().PointProjection(point); }
      
    DecomposedVector_t DecomposePoint(Point_t const& point) const
      { return DecomposeVector(Base().ToVector(point)); }
    
    
    
    
    auto VectorMainComponent(Vector_t const& v) const
      { return Plane().VectorMainComponent(v); }
    
    auto VectorSecondaryComponent(Vector_t const& v) const
      { return Plane().VectorSecondaryComponent(v); }
    
    auto VectorNormalComponent(Vector_t const& v) const
      { return geo::vect::dot(v, NormalDir()); }
    
    Projection_t ProjectVectorOnPlane(Vector_t const& v) const
      { return Plane().VectorProjection(v); }
      
    DecomposedVector_t DecomposeVector(Vector_t const& v) const
      { return { VectorNormalComponent(v), ProjectVectorOnPlane(v) }; }
    
    double Angle(Vector_t const& v) const
      { return Plane().Angle(v); }
    
    
    
    
    auto MainComponent(Projection_t const& v) const
      { return Plane().MainComponent(v); }
    
    auto SecondaryComponent(Projection_t const& v) const
      { return Plane().SecondaryComponent(v); }
    
    
    
    
    Point_t ComposePoint(DecomposedVector_t const& decomp) const
      { return ComposePoint(decomp.distance, decomp.projection); }
    
    Point_t ComposePoint(double distance, Projection_t const& proj) const
      { return ReferencePoint() + ComposeVector(distance, proj); }
    
    
    
    Vector_t ComposeVector(DecomposedVector_t const& decomp) const
      { return ComposeVector(decomp.distance, decomp.projection); }
    
    Vector_t ComposeVector(double distance, Projection_t const& proj) const
      { return Plane().ComposeVector(proj) + distance * NormalDir(); }
    
    
    
  }; // class Decomposer<>
  
  // --- END Decomposition objects ---------------------------------------------
  
} // namespace geo


#endif // LARCOREALG_GEOMETRY_DECOMPOSER_H