SimpleGeo.h

Go to the documentation of this file.

#ifndef LARCOREALG_GEOMETRY_SIMPLEGEO_H
#define LARCOREALG_GEOMETRY_SIMPLEGEO_H


// C/C++ standard library
#include <array>
#include <algorithm> // std::min(), std::max()
#include <stdexcept> // std::runtime_error


namespace lar {
  namespace util {
    
    // --- BEGIN Geometry group ------------------------------------------------

    namespace simple_geo {
      
      
      template <typename Data = double>
      struct Point2D {
        using Data_t = Data;
        Data_t x = 0.;
        Data_t y = 0.;
        
        Point2D() = default;
        Point2D(Data_t x, Data_t y): x(x), y(y) {}
        Point2D(Data_t const* p): x(p[0]), y(p[1]) {}
        
      }; // struct Point2D<>
      
      template <typename T>
      bool operator== (Point2D<T> const& a, Point2D<T> const& b)
        { return (a.x == b.x) && (a.y == b.y); }
      template <typename T>
      bool operator!= (Point2D<T> const& a, Point2D<T> const& b)
        { return (a.x != b.x) || (a.y != b.y); }
      template <typename T>
      Point2D<T> operator+ (Point2D<T> const& a, Point2D<T> const& b)
        { return { a.x + b.x, a.y + b.y }; }
      template <typename T>
      Point2D<T> operator* (Point2D<T> const& p, typename Point2D<T>::Data_t f)
        { return { p.x * f, p.y * f }; }
      template <typename T>
      Point2D<T> operator/ (Point2D<T> const& p, typename Point2D<T>::Data_t f)
        { return { p.x / f, p.y / f }; }
      template <typename Stream, typename T>
      Stream& operator<< (Stream&& out, Point2D<T> const& p)
        { out << "( " << p.x << " ; " << p.y << " )"; return out; }
      
      
      template <typename Data = double>
      struct Point3D {
        using Data_t = Data;
        Data_t x = 0.;
        Data_t y = 0.;
        Data_t z = 0.;
        
        Point3D() = default;
        Point3D(Data_t x, Data_t y, Data_t z): x(x), y(y), z(z) {}
        Point3D(Data_t const* p): x(p[0]), y(p[1]), z(p[2]) {}
      }; // struct Point3D<>
      
      template <typename T>
      bool operator== (Point3D<T> const& a, Point3D<T> const& b)
        { return (a.x == b.x) && (a.y == b.y) && (a.z == b.z); }
      template <typename T>
      bool operator!= (Point3D<T> const& a, Point3D<T> const& b)
        { return (a.x != b.x) || (a.y != b.y) || (a.z != b.z); }
      template <typename T>
      Point3D<T> operator+ (Point3D<T> const& a, Point3D<T> const& b)
        { return { a.x + b.x, a.y + b.y, a.z + b.z }; }
      template <typename T>
      Point3D<T> operator* (Point3D<T> const& p, typename Point2D<T>::Data_t f)
        { return { p.x * f, p.y * f, p.z * f }; }
      template <typename T>
      Point3D<T> operator/ (Point3D<T> const& p, typename Point2D<T>::Data_t f)
        { return { p.x / f, p.y / f, p.z / f }; }
      template <typename Stream, typename T>
      Stream& operator<< (Stream&& out, Point3D<T> const& p)
        {
          out << "( " << p.x << " ; " << p.y << " ; " << p.z << " )";
          return out;
        }
      
      
      
      template <typename Point>
      class AreaBase {
          public:
        using Point_t = Point;
        using Area_t = AreaBase<Point>;
        
        class NullIntersection: public std::runtime_error
          { using std::runtime_error::runtime_error; };
        
        
        AreaBase() = default;
        
        AreaBase(Point_t const& a, Point_t const& b)
          { set_sorted(min.x, max.x, a.x, b.x); }
        
        Point_t const& Min() const { return min; }
        Point_t const& Max() const { return max; }
        Point_t Center() const { return (min + max) * 0.5; }
        
        auto DeltaX() const { return Max().x - Min().x; }
        bool isEmptyX() const { return (DeltaX() == 0); }
        bool isNullX() const { return Max().x < Min().x; }
        bool isNull() const { return isNullX(); }
        unsigned int nonEmptyDims() const { return isEmptyX()? 0: 1; }
        bool isEmpty() const { return nonEmptyDims() == 0; }
        bool isLine() const { return nonEmptyDims() == 1; }
        
        auto thinnestSize() const { return DeltaX(); }
        
        unsigned int thinnestSide() const { return 0; }
        
        void IncludePoint(Point_t const& point)
          { set_min_max(min.x, max.x, point.x); }
        
        void Include(Area_t const& area)
          { IncludePoint(area.min); IncludePoint(area.max); }
        
        void Intersect(Area_t const& area)
          {
            set_max(min.x, area.min.x);
            set_min(max.x, area.max.x);
            if (isNullX()) throw NullIntersection("null x dimension");
          }
        
        bool operator== (Area_t const& as) const
          { return (min == as.min) && (max == as.max); }
        bool operator!= (Area_t const& as) const
          { return (min != as.min) || (max != as.max); }
        
          protected:
        using Data_t = typename Point_t::Data_t;
        
        Point_t min, max;
        
        static void set_min(Data_t& var, Data_t val)
          { if (val < var) var = val; }
        
        static void set_max(Data_t& var, Data_t val)
          { if (val > var) var = val; }
          
        static void set_min_max(Data_t& min_var, Data_t& max_var, Data_t val)
          { set_min(min_var, val); set_max(max_var, val); }
          
        static void set_sorted
          (Data_t& min_var, Data_t& max_var, Data_t a, Data_t b)
          {
            if (a > b) { min_var = b; max_var = a; }
            else       { min_var = a; max_var = b; }
          }
        
      }; // class AreaBase<>
      
      
      template <typename Point = Point2D<double>>
      class Area: public AreaBase<Point> {
        using Base_t = AreaBase<Point>;
        
          public:
        using Point_t = typename Base_t::Point_t;
        using Area_t = Area<Point_t>;
        
        Area() = default;
        
        Area(Point_t const& a, Point_t const& b)
          : Base_t(a, b)
          { Base_t::set_sorted(Base_t::min.y, Base_t::max.y, a.y, b.y); }
        
        auto DeltaY() const { return Base_t::Max().y - Base_t::Min().y; }
        bool isEmptyY() const { return DeltaY() == 0; }
        unsigned int nonEmptyDims() const
          { return Base_t::nonEmptyDims() + (isEmptyY()? 0: 1); }
        bool isNullY() const { return Base_t::Max().y < Base_t::Min().y; }
        bool isNull() const { return Base_t::isNull() || isNullY(); }
        bool isEmpty() const { return nonEmptyDims() == 0; }
        bool isLine() const { return nonEmptyDims() == 1; }
        bool isPlane() const { return nonEmptyDims() == 2; }
        
        auto thinnestSize() const
          { return std::min(DeltaY(), Base_t::thinnestSize()); }
        
        unsigned int thinnestSide() const
          {
            return
              (DeltaY() < Base_t::thinnestSize())? 1: Base_t::thinnestSide();
          }
        
        void IncludePoint(Point_t const& point)
          {
            Base_t::IncludePoint(point);
            Base_t::set_min_max(Base_t::min.y, Base_t::max.y, point.y);
          }
        
        void Include(Area_t const& area)
          { IncludePoint(area.min); IncludePoint(area.max); }
        
        void Intersect(Area_t const& area)
          {
            Base_t::Intersect(area);
            Base_t::set_max(Base_t::min.y, area.min.y);
            Base_t::set_min(Base_t::max.y, area.max.y);
            if (isNullY())
              throw typename Base_t::NullIntersection("null y dimension");
          }
        
      }; // class Area<>
      
      template <typename Stream, typename Point>
      Stream& operator<< (Stream&& out, AreaBase<Point> const& area)
        { out << area.Min() << " - " << area.Max(); return out; }
      
      
      template <typename Point = Point3D<double>>
      class Volume: public Area<Point> {
        using Base_t = Area<Point>;
        
          public:
        using Point_t = typename Base_t::Point_t;
        using Volume_t = Volume<Point_t>;
        
        Volume() = default;
        
        Volume(Point_t const& a, Point_t const& b)
          : Base_t(a, b)
          { Base_t::set_sorted(Base_t::min.z, Base_t::max.z, a.z, b.z); }
        
        auto DeltaZ() const { return Base_t::Max().z - Base_t::Min().z; }
        bool isEmptyZ() const { return DeltaZ() == 0; }
        unsigned int nonEmptyDims() const
          { return Base_t::nonEmptyDims() + (isEmptyZ()? 0: 1); }
        bool isNullZ() const { return Base_t::Max().z < Base_t::Min().z; }
        bool isNull() const { return Base_t::isNull() || isNullZ(); }
        bool isEmpty() const { return nonEmptyDims() == 0; }
        bool isLine() const { return nonEmptyDims() == 1; }
        bool isPlane() const { return nonEmptyDims() == 2; }
        bool isVolume() const { return nonEmptyDims() == 3; }
        
        auto thinnestSize() const
          { return std::min(DeltaZ(), Base_t::thinnestSize()); }
        
        unsigned int thinnestSide() const
          {
            return
              (DeltaZ() < Base_t::thinnestSize())? 2: Base_t::thinnestSide();
          }
        
        void IncludePoint(Point_t const& point)
          {
            Base_t::IncludePoint(point);
            Base_t::set_min_max(Base_t::min.z, Base_t::max.z, point.z);
          }
        
        void Include(Volume_t const& volume)
          { IncludePoint(volume.min); IncludePoint(volume.max); }
        
        void Intersect(Volume_t const& volume)
          {
            Base_t::Intersect(volume);
            Base_t::set_max(Base_t::min.z, volume.min.z);
            Base_t::set_min(Base_t::max.z, volume.max.z);
            if (isNullZ())
              throw typename Base_t::NullIntersection("null z dimension");
          }
        
      }; // class Volume<>
      
      
      
      template <typename Data = double>
      struct Range {
        using Data_t = Data; 
        using Range_t = Range<Data>; 
        
        Data_t lower = 1.0;  
        Data_t upper = 0.0; 
        
        Range() = default;
        
        Range(Data_t lower, Data_t upper, bool doSort = false)
          : lower(lower), upper(upper)
          { if (doSort) sort(); }
        
        bool isNull() const { return lower >= upper; }
        
        Data_t length() const { return std::max(upper - lower, Data_t(0.0)); }
        
        bool contains(Data_t v) const { return (v >= lower) && (v <= upper); }
        
        bool overlaps(Range_t const& r) const;
        
        Data_t delta(Data_t v, Data_t margin = 0.0) const;
        
        void extendToInclude(Data_t);
        
        void extendToInclude(Range_t const& r);
        
        void intersect(Range_t const& r);
        
          private:
        void sort() { if (lower > upper) std::swap(lower, upper); }
        
        void makeNull() { *this = Range_t{}; }
        
      }; // struct Range<>
      
      template <typename Stream, typename Data>
      Stream& operator<< (Stream&& out, Range<Data> const& range);
      
      template <typename Data = double>
      struct Rectangle {
        using Data_t = Data; 
        using Rectangle_t = Rectangle<Data>; 
        using Range_t = Range<Data_t>; 
        
        Range_t width; 
        Range_t depth; 
        
        Rectangle() = default;
        
        Rectangle(Range_t const& width, Range_t const& depth)
          : width(width), depth(depth)
          {}
        
        bool isNull() const { return width.isNull() || depth.isNull(); }
        
        bool contains(Data_t w, Data_t d) const
          { return width.contains(w) && depth.contains(d); }
        
        bool overlaps(Rectangle_t const& r) const;
        
        void extendToInclude(Rectangle_t const& r);
        
      }; // Rectangle<>
      
      template <typename Stream, typename Data>
      Stream& operator<< (Stream&& out, Rectangle<Data> const& rect);
      
      
      
    } // namespace simple_geo
    
    // --- END Geometry group --------------------------------------------------

  } // namespace util
} // namespace lar


//==============================================================================
//--- Template implementation
//--- 
//------------------------------------------------------------------------------
//---  geo::PlaneGeo::Range<>
//---  
template <typename Data>
auto lar::util::simple_geo::Range<Data>::delta
  (Data_t v, Data_t margin /* = 0 */) const
  -> Data_t
{
  
  if (v < (lower + margin)) return lower + margin - v; // always positive
  if (v > (upper - margin)) return upper - margin - v; // always negative
  return 0.0;                                          // always zero
  
} // lar::util::simple_geo::Range<Data>::delta()


//------------------------------------------------------------------------------
template <typename Data>
void lar::util::simple_geo::Range<Data>::extendToInclude(Data_t v) {
  
  if (lower > upper) lower = upper = v;
  else if (lower > v) lower = v;
  else if (upper < v) upper = v;
  
} // lar::util::simple_geo::Range<Data>::extendToInclude()


//------------------------------------------------------------------------------
template <typename Data>
void lar::util::simple_geo::Range<Data>::extendToInclude(Range_t const& r) {
  
  if (r.isNull()) return;
  if (isNull()) {
    *this = r;
    return;
  }
  extendToInclude(r.lower);
  extendToInclude(r.upper);
  
} // lar::util::simple_geo::Range<Data>::extendToInclude()


//------------------------------------------------------------------------------
template <typename Data>
void lar::util::simple_geo::Range<Data>::intersect(Range_t const& r) {
  // this implementation explicitly makes the range default-constructed-null
  // if the intersection results empty
  if (isNull()) return;
  if (r.isNull()) {
    makeNull();
    return;
  }
  if (lower < r.lower) lower = r.lower;
  if (upper > r.upper) upper = r.upper;
  if (lower > upper) makeNull();
} // lar::util::simple_geo::Range<Data>::intersect()


//------------------------------------------------------------------------------
template <typename Data>
bool lar::util::simple_geo::Range<Data>::overlaps(Range_t const& r) const {
  if (isNull() || r.isNull()) return false;
  return (r.lower < upper) && (lower < r.upper);
} // lar::util::simple_geo::Range<Data>::overlaps()


//------------------------------------------------------------------------------
template <typename Stream, typename Data>
Stream& lar::util::simple_geo::operator<<
  (Stream&& out, Range<Data> const& range)
{
  out << "( " << range.lower << " -- " << range.upper << " )";
  return out;
} // operator<< (Range)


//------------------------------------------------------------------------------
template <typename Data>
void lar::util::simple_geo::Rectangle<Data>::extendToInclude
  (Rectangle_t const& r)
{
  width.extendToInclude(r.width);
  depth.extendToInclude(r.depth);
} // lar::util::simple_geo::Rectangle<Data>::extendToInclude()


//------------------------------------------------------------------------------
template <typename Data>
bool lar::util::simple_geo::Rectangle<Data>::overlaps
  (Rectangle_t const& r) const
{
  if (isNull() || r.isNull()) return false;
  return width.overlap(r.width) && depth.overlap(r.depth);
} // lar::util::simple_geo::Rectangle<Data>::overlaps()


//------------------------------------------------------------------------------
template <typename Stream, typename Data>
Stream& lar::util::simple_geo::operator<<
  (Stream&& out, Rectangle<Data> const& rect)
{
  out << "w=" << rect.width << " d=" << rect.depth;
  return out;
} // operator<< (Rectangle)


//------------------------------------------------------------------------------


#endif // LARCOREALG_GEOMETRY_SIMPLEGEO_H