SpacePartition.h

Go to the documentation of this file.

#ifndef LAREXAMPLES_ALGORITHMS_REMOVEISOLATEDSPACEPOINTS_SPACEPARTITION_H
#define LAREXAMPLES_ALGORITHMS_REMOVEISOLATEDSPACEPOINTS_SPACEPARTITION_H

// LArSoft libraries
#include "lardata/Utilities/GridContainers.h"

// C/C++ standard libraries
#include <cstddef> // std::ptrdiff_t
#include <cmath> // std::ceil()
#include <vector>
#include <array>
#include <string>
#include <stdexcept> // std::runtime_error


namespace lar {
  namespace example {
    
    // BEGIN RemoveIsolatedSpacePoints group -----------------------------------

    template <typename Point>
    struct PositionExtractor;
    
    namespace details {
      template <typename Point>
      struct PointTraits_t;
      
      template <typename Point>
      using ExtractCoordType_t
        = typename details::PointTraits_t<Point>::Coord_t;
      
    } // namespace details
    
    
    template <typename Coord>
    struct CoordRange {
      using Range_t = CoordRange<Coord>; 
      using Coord_t = Coord; 
      
      Coord_t lower; 
      Coord_t upper; 
      
      bool contains(Coord_t c) const;
      
      bool empty() const;
      
      bool valid() const;
      
      Coord_t size() const;
      
      Coord_t offset(Coord_t c) const { return c - lower; }
      
      bool operator== (const Range_t& as) const;
      
      bool operator!= (const Range_t& than) const;
      
    }; // CoordRange<>
    
    
    template <typename Coord>
    struct CoordRangeCells: public CoordRange<Coord> {
      using Base_t = CoordRange<Coord>;
      
      using Coord_t = typename Base_t::Coord_t;
      
      Coord_t cellSize; 
      
      CoordRangeCells(Coord_t low, Coord_t high, Coord_t cs);
      
      CoordRangeCells(Base_t const& range, Coord_t cs);
      
      
      std::ptrdiff_t findCell(Coord_t c) const;
      
    }; // CoordRangeCells<>
    
    
    template <typename PointIter>
    class SpacePartition {
      using Point_t = decltype(*(PointIter())); 
      using Grid_t = util::GridContainer3D<PointIter>; 
      
        public:
      using Coord_t = details::ExtractCoordType_t<Point_t>;
      using Range_t = CoordRangeCells<Coord_t>; 
      
      using Indexer_t = typename Grid_t::Indexer_t;
      
      using CellIndexOffset_t = typename Indexer_t::CellIndexOffset_t;
      
      using CellIndex_t = typename Indexer_t::CellIndex_t; 
      
      using CellID_t = typename Indexer_t::CellID_t;
      
      using Cell_t = typename Grid_t::Cell_t;
      
      SpacePartition
        (Range_t rangeX, Range_t rangeY, Range_t rangeZ);
      
      void fill(PointIter begin, PointIter end);
      
      CellIndexOffset_t pointIndex(Point_t const& point) const;
      
      Indexer_t const& indexManager() const
        { return data.indexManager(); }
      
      bool has(CellIndexOffset_t ofs) const
        { return data.has(ofs); }
      
      Cell_t const& operator[] (CellIndex_t index) const
        { return data[index]; }
      
      typename Grid_t::const_iterator begin() const;
      
      typename Grid_t::const_iterator end() const;
      
        protected:
      using CellDimIndex_t = typename Grid_t::CellDimIndex_t;
      
      Coord_t cellSize; 
      
      Range_t xRange; 
      Range_t yRange; 
      Range_t zRange; 
      
      Grid_t data; 
      
    }; // SpacePartition<>
    
    
    
    //--------------------------------------------------------------------------
    
    namespace details {
      template <typename Cont, typename Data>
      struct PositionExtractorFromArray {
        static Data x(Cont const& p) { return p[0]; }
        static Data y(Cont const& p) { return p[1]; }
        static Data z(Cont const& p) { return p[2]; }
      }; // PositionExtractorFromArray<T*>
      
      template <typename Point>
      auto extractPositionX(Point const& point)
        { return PositionExtractor<Point>::x(point); }
      template <typename Point>
      auto extractPositionY(Point const& point)
        { return PositionExtractor<Point>::y(point); }
      template <typename Point>
      auto extractPositionZ(Point const& point)
        { return PositionExtractor<Point>::z(point); }
      
      template <typename Point>
      struct PointTraits_t {
        using Coord_t = decltype(extractPositionX(std::decay_t<Point>()));
      }; // ExtractCoordTypeHelper_t
    
    } // namespace details
    
    template <typename T>
    struct PositionExtractor<T*>:
      public details::PositionExtractorFromArray<T*, T>
    {};
    
    template <typename T>
    struct PositionExtractor<std::array<T, 3U>>:
      public details::PositionExtractorFromArray<std::array<T, 3U>, T>
    {};
    
    template <typename T>
    struct PositionExtractor<std::vector<T>>: 
      public details::PositionExtractorFromArray<std::vector<T>, T>
    {};
    
    
    // END RemoveIsolatedSpacePoints group -------------------------------------
    
  } // namespace example
} // namespace lar


//------------------------------------------------------------------------------
//--- lar::example::details
//---
namespace lar {
  namespace example {
    namespace details {
      
      template <typename Coord>
      std::array<size_t, 3> diceVolume(
        CoordRangeCells<Coord> const& rangeX,
        CoordRangeCells<Coord> const& rangeY,
        CoordRangeCells<Coord> const& rangeZ
        )
        {
          return {{
            size_t(std::ceil(rangeX.size() / rangeX.cellSize)),
            size_t(std::ceil(rangeY.size() / rangeY.cellSize)),
            size_t(std::ceil(rangeZ.size() / rangeZ.cellSize))
            }};
        } // diceVolume()
      
    } // namespace details
  } // namespace example
} // namespace lar


//------------------------------------------------------------------------------
//--- template implementation
//------------------------------------------------------------------------------
//--- lar::example::CoordRange
//---
template <typename Coord>
bool lar::example::CoordRange<Coord>::contains(Coord_t c) const
  { return (lower <= c) && (upper >= c); }

template <typename Coord>
bool lar::example::CoordRange<Coord>::Range_t::empty() const
  { return lower == upper; }


template <typename Coord>
bool lar::example::CoordRange<Coord>::valid() const
  { return lower <= upper; }

template <typename Coord>
Coord lar::example::CoordRange<Coord>::size() const
  { return upper - lower; }

template <typename Coord>
bool lar::example::CoordRange<Coord>::operator== (const Range_t& as) const
  { return (upper == as.upper) && (lower == as.lower); }

template <typename Coord>
bool lar::example::CoordRange<Coord>::operator!= (const Range_t& than) const
  { return (upper != than.upper) || (lower != than.lower); }


//------------------------------------------------------------------------------
//--- lar::example::CoordRange
//---
template <typename Coord>
lar::example::CoordRangeCells<Coord>::CoordRangeCells
  (Coord_t low, Coord_t high, Coord_t cs)
  : Base_t(low, high), cellSize(cs)
  {}

template <typename Coord>
lar::example::CoordRangeCells<Coord>::CoordRangeCells
  (Base_t const& range, Coord_t cs)
  : Base_t(range), cellSize(cs)
  {}

//------------------------------------------------------------------------------
template <typename Coord>
std::ptrdiff_t lar::example::CoordRangeCells<Coord>::findCell(Coord_t c) const
  { return std::ptrdiff_t(Base_t::offset(c) / cellSize); }
      

//------------------------------------------------------------------------------
//--- lar::example::SpacePartition
//---
template <typename PointIter>
lar::example::SpacePartition<PointIter>::SpacePartition
  (Range_t rangeX, Range_t rangeY, Range_t rangeZ)
  : xRange(rangeX)
  , yRange(rangeY)
  , zRange(rangeZ)
  , data(details::diceVolume(xRange, yRange, zRange))
{
  /*
    std::cout << "Grid: "
      << indexManager().sizeX() << " x "
      << indexManager().sizeY() << " x " << indexManager().sizeZ()
      << " (" << indexManager().size() << " cells)"
      << "\n  range X: " << xRange.lower << " -- " << xRange.upper << " [/" << xRange.cellSize << "]"
      << "\n  range Y: " << yRange.lower << " -- " << yRange.upper << " [/" << yRange.cellSize << "]"
      << "\n  range Z: " << zRange.lower << " -- " << zRange.upper << " [/" << zRange.cellSize << "]"
      << std::endl;
  */
} // lar::example::SpacePartition<>::SpacePartition


//--------------------------------------------------------------------------
template <typename PointIter>
void lar::example::SpacePartition<PointIter>::fill
  (PointIter begin, PointIter end)
{
  
  PointIter it = begin;
  while (it != end) {
    // if the point is outside the volume, pointIndex will throw an exception
    data.insert(pointIndex(*it), it);
    ++it;
  } // while
  
} // lar::example::SpacePartition<>::fill()


//--------------------------------------------------------------------------
template <typename PointIter>
typename lar::example::SpacePartition<PointIter>::CellIndexOffset_t
lar::example::SpacePartition<PointIter>::pointIndex(Point_t const& point) const
{
  // compute the cell ID coordinates
  Coord_t const x = details::extractPositionX(point);
  CellDimIndex_t const xc = xRange.findCell(x);
  if (!data.hasX(xc)) {
    throw std::runtime_error
      ("Point out of the volume (x = " + std::to_string(x) + ")");
  }
  
  Coord_t const y = details::extractPositionY(point);
  CellDimIndex_t const yc = yRange.findCell(y);
  if (!data.hasY(yc)) {
    throw std::runtime_error
      ("Point out of the volume (y = " + std::to_string(y) + ")");
  }
  
  Coord_t const z = details::extractPositionZ(point);
  CellDimIndex_t const zc = zRange.findCell(z);
  if (!data.hasZ(zc)) {
    throw std::runtime_error
      ("Point out of the volume (z = " + std::to_string(z) + ")");
  }
  
  // return its index
  return data.index(CellID_t{{ xc, yc, zc }});
  
} // lar::example::SpacePartition<>::pointIndex()


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

#endif // LAREXAMPLES_ALGORITHMS_REMOVEISOLATEDSPACEPOINTS_SPACEPARTITION_H