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 <type_traits> // std::declval()
#include <utility>     // std::move()

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 vect::rounded01(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(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 vect::dot(v, MainDir()); }

    auto VectorSecondaryComponent(Vector_t const& v) const { return 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 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