PointIsolationAlg.h

Go to the documentation of this file.

#ifndef LAREXAMPLES_ALGORITHMS_REMOVEISOLATEDSPACEPOINTS_POINTISOLATIONALG_H
#define LAREXAMPLES_ALGORITHMS_REMOVEISOLATEDSPACEPOINTS_POINTISOLATIONALG_H

// LArSoft libraries
#include "larexamples/Algorithms/RemoveIsolatedSpacePoints/SpacePartition.h"

// infrastructure and utilities
#include "cetlib/pow.h" // cet::sum_squares()

// C/C++ standard libraries
#include <cassert> // assert()
#include <cmath> // std::sqrt()
#include <vector>
#include <array>
#include <string>
#include <type_traits> // std::add_const_t<>
#include <iterator> // std::cbegin(), std::cend(), std::distance()
#include <stdexcept> // std::runtime_error


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

    template <typename Coord = double>
    class PointIsolationAlg {
      
        public:
      using Coord_t = Coord;
      using Range_t = CoordRange<Coord_t>;
      
      struct Configuration_t {
        Range_t rangeX;   
        Range_t rangeY;   
        Range_t rangeZ;   
        Coord_t radius2;  
        size_t maxMemory = 100 * 1048576;
      }; // Configuration_t
      
      
      
      PointIsolationAlg(Configuration_t const& first_config)
        : config(first_config)
        {}
      
      void reconfigure(Configuration_t const& new_config)
        { config = new_config; }
      
      
      Configuration_t& configuration() const { return config; }
      
      
      
      template <typename PointIter>
      std::vector<size_t> removeIsolatedPoints
        (PointIter begin, PointIter end) const;
      
      
      template <typename Cont>
      std::vector<size_t> removeIsolatedPoints (Cont const& points) const
        { return removeIsolatedPoints(std::cbegin(points), std::cend(points)); }
      
      
      template <typename PointIter>
      std::vector<size_t> bruteRemoveIsolatedPoints
        (PointIter begin, PointIter end) const;
      
      
      
      static void validateConfiguration(Configuration_t const& config);
      
      
      
      static Coord_t maximumOptimalCellSize(Coord_t radius)
        { return radius / std::sqrt(3.); }
      
      
        private:
      using Indexer_t = util::GridContainer3DIndices; // same in GridContainer
      
      using NeighAddresses_t = std::vector<Indexer_t::CellIndexOffset_t>;
      
      template <typename PointIter>
      using Partition_t = SpacePartition<PointIter>;
      
      template <typename PointIter>
      using Point_t = decltype(*PointIter());
      
      
      Configuration_t config; 
      
      
      template <typename PointIter = std::array<double, 3> const*>
      Coord_t computeCellSize() const;
      
      
      NeighAddresses_t buildNeighborhood
        (Indexer_t const& indexer, unsigned int neighExtent) const;
      
      template <typename PointIter>
      bool isPointIsolatedFrom(
        Point_t<PointIter> const& point,
        typename Partition_t<PointIter>::Cell_t const& otherPoints
        ) const;
        
      template <typename PointIter>
      bool isPointIsolatedWithinNeighborhood(
        Partition_t<PointIter> const& partition,
        Indexer_t::CellIndex_t cellIndex,
        Point_t<PointIter> const& point,
        NeighAddresses_t const& neighList
        ) const;

      template <typename Point>
      bool closeEnough(Point const& A, Point const& B) const;
      
      
      static std::string rangeString(Coord_t from, Coord_t to);
      
      static std::string rangeString(Range_t range)
        { return rangeString(range.lower, range.upper); }
      
    }; // class PointIsolationAlg
    
    
    //--------------------------------------------------------------------------
    // END RemoveIsolatedSpacePoints group -------------------------------------
    
  } // namespace example
} // namespace lar



//------------------------------------------------------------------------------
//--- template implementation
//------------------------------------------------------------------------------
template <typename Coord>
template <typename PointIter>
std::vector<size_t> lar::example::PointIsolationAlg<Coord>::removeIsolatedPoints
  (PointIter begin, PointIter end) const
{
  
  std::vector<size_t> nonIsolated;
  
  Coord_t const R = std::sqrt(config.radius2);
  
  //
  // determine space partition settings: cell size
  // 
  // maximum: the volume of a single cell must be contained in a sphere with
  // radius equal to the isolation radius R
  // 
  // minimum: needs tuning
  // 
  
  Coord_t cellSize = computeCellSize<PointIter>();
  assert(cellSize > 0);
  Partition_t<PointIter> partition(
    { config.rangeX, cellSize },
    { config.rangeY, cellSize },
    { config.rangeZ, cellSize }
    );
  
  // if a cell is contained in a sphere with 
  bool const cellContainedInIsolationSphere
    = (cellSize <= maximumOptimalCellSize(R));
  
  // 
  // determine neighbourhood
  // the neighbourhood is the number of cells that might contain points closer
  // than R to a cell; it is equal to R in cell size units, rounded up;
  // it's expressed as a list of coordinate shifts from a base cell to all the
  // others in the neighbourhood; it is contained in a cube
  //
  unsigned int const neighExtent = (int) std::ceil(R / cellSize);
  NeighAddresses_t neighList
    = buildNeighborhood(partition.indexManager(), neighExtent);
  
  // if a cell is not fully contained in a isolation radius, we need to check
  // the points of the cell with each other: their cell becomes part of the
  // neighbourhood
  if (!cellContainedInIsolationSphere)
    neighList.insert(neighList.begin(), { 0, 0, 0 });
  
  //
  // populate the partition
  //
  partition.fill(begin, end);
  
  //
  // for each cell in the partition:
  //
  size_t const nCells = partition.indexManager().size();
  for (Indexer_t::CellIndex_t cellIndex = 0; cellIndex < nCells; ++cellIndex) {
    auto const& cellPoints = partition[cellIndex];
    
    //
    // if the cell has more than one element, mark all points as non-isolated;
    // true only if the cell is completely contained within a R rßadius
    //
    if (cellContainedInIsolationSphere && (cellPoints.size() > 1)) {
      for (auto const& pointPtr: cellPoints)
        nonIsolated.push_back(std::distance(begin, pointPtr));
      continue;
    } // if all non-isolated
    
    //
    // brute force approach: try all the points in this cell against all the
    // points in the neighbourhood
    //
    for (auto const pointPtr: cellPoints) {
      //
      // optimisation (speed): mark the points from other cells as non-isolated
      // when they trigger non-isolation in points of the current one
      //
      
      // TODO
      
      if (!isPointIsolatedWithinNeighborhood
        (partition, cellIndex, *pointPtr, neighList)
        )
      {
        nonIsolated.push_back(std::distance(begin, pointPtr));
      }
    } // for points in cell
    
  } // for cell
  
  return nonIsolated;
} // lar::example::PointIsolationAlg::removeIsolatedPoints()


//--------------------------------------------------------------------------
template <typename Coord>
void lar::example::PointIsolationAlg<Coord>::validateConfiguration
  (Configuration_t const& config)
{
  std::vector<std::string> errors;
  if (config.radius2 < Coord_t(0)) {
    errors.push_back
      ("invalid radius squared (" + std::to_string(config.radius2) + ")");
  }
  if (!config.rangeX.valid()) {
    errors.push_back("invalid x range " + rangeString(config.rangeX));
  }
  if (!config.rangeY.valid()) {
    errors.push_back("invalid y range " + rangeString(config.rangeY));
  }
  if (!config.rangeZ.valid()) {
    errors.push_back("invalid z range " + rangeString(config.rangeZ));
  }
  
  if (errors.empty()) return;
  
  // compose the full error message as concatenation:
  std::string message
    (std::to_string(errors.size()) + " configuration errors found:");
  
  for (auto const& error: errors) message += "\n * " + error;
  throw std::runtime_error(message);
  
} // lar::example::PointIsolationAlg::validateConfiguration()


//--------------------------------------------------------------------------
template <typename Coord>
template <typename PointIter /* = std::array<double, 3> const* */>
Coord lar::example::PointIsolationAlg<Coord>::computeCellSize() const {
  
  Coord_t const R = std::sqrt(config.radius2);
  
  // maximum: the maximum distance between two points in the cell (that is,
  //   the diagonal of the cell) must be no larger than the isolation radius R;
  // minimum: needs tuning
  Coord_t cellSize = maximumOptimalCellSize(R); // try the minimum for now

  //
  // optimisation (memory): determine minimum size of box
  //
  
  // TODO
  
  if (config.maxMemory == 0) return cellSize;
  
  do {
    std::array<size_t, 3> partition = details::diceVolume(
      CoordRangeCells<Coord_t>{ config.rangeX, cellSize },
      CoordRangeCells<Coord_t>{ config.rangeY, cellSize },
      CoordRangeCells<Coord_t>{ config.rangeZ, cellSize }
      );
    
    size_t const nCells = partition[0] * partition[1] * partition[2];
    if (nCells <= 1) break; // we can't reduce it any further
    
    // is memory low enough?
    size_t const memory
      = nCells * sizeof(typename SpacePartition<PointIter>::Cell_t);
    if (memory < config.maxMemory) break;
    
    cellSize *= 2;
  } while (true);
  
  return cellSize;
} // lar::example::PointIsolationAlg<Coord>::computeCellSize()


//------------------------------------------------------------------------------
template <typename Coord>
typename lar::example::PointIsolationAlg<Coord>::NeighAddresses_t
lar::example::PointIsolationAlg<Coord>::buildNeighborhood
  (Indexer_t const& indexer, unsigned int neighExtent) const
{
  unsigned int const neighSize = 1 + 2 * neighExtent;
  NeighAddresses_t neighList;
  neighList.reserve(neighSize * neighSize * neighSize - 1);
  
  using CellID_t = Indexer_t::CellID_t;
  using CellDimIndex_t = Indexer_t::CellDimIndex_t;
  
  //
  // optimisation (speed): reshape the neighbourhood
  // neighbourhood might cut out cells close to the vertices
  //
  
  // TODO
  
  CellDimIndex_t const ext = neighExtent; // convert into the right signedness
  
  CellID_t center{{ 0, 0, 0 }}, cellID;
  for (CellDimIndex_t ixOfs = -ext; ixOfs <= ext; ++ixOfs) {
    cellID[0] = ixOfs;
    for (CellDimIndex_t iyOfs = -ext; iyOfs <= ext; ++iyOfs) {
      cellID[1] = iyOfs;
      for (CellDimIndex_t izOfs = -ext; izOfs <= ext; ++izOfs) {
        if ((ixOfs == 0) && (iyOfs == 0) && (izOfs == 0)) continue;
        cellID[2] = izOfs;
        
        neighList.push_back(indexer.offset(center, cellID));
        
      } // for ixOfs
    } // for iyOfs
  } // for izOfs
  
  return neighList;
} // lar::example::PointIsolationAlg<Coord>::buildNeighborhood()

  
//--------------------------------------------------------------------------
template <typename Coord>
template <typename PointIter>
bool lar::example::PointIsolationAlg<Coord>::isPointIsolatedFrom(
  Point_t<PointIter> const& point,
  typename Partition_t<PointIter>::Cell_t const& otherPoints
) const
{
  
  for (auto const& otherPointPtr: otherPoints) {
    // make sure that we did not compare the point with itself
    if (closeEnough(point, *otherPointPtr) && (&point != &*otherPointPtr))
      return false;
  }
  
  return true;
  
} // lar::example::PointIsolationAlg<Coord>::isPointIsolatedFrom()


//--------------------------------------------------------------------------
template <typename Coord>
template <typename PointIter>
bool lar::example::PointIsolationAlg<Coord>::isPointIsolatedWithinNeighborhood(
  Partition_t<PointIter> const& partition,
  Indexer_t::CellIndex_t cellIndex,
  Point_t<PointIter> const& point,
  NeighAddresses_t const& neighList
) const
{
  
  // check in all cells of the neighbourhood
  for (Indexer_t::CellIndexOffset_t neighOfs: neighList) {
    
    //
    // optimisation (speed): have neighbour offsets so that the invalid ones
    // are all at the beginning and at the end, so that skipping is faster
    //
    
    if (!partition.has(cellIndex + neighOfs)) continue;
    auto const& neighCellPoints = partition[cellIndex + neighOfs];

    if (!isPointIsolatedFrom<PointIter>(point, neighCellPoints)) return false;
    
  } // for neigh cell
  
  return true;
  
} // lar::example::PointIsolationAlg<Coord>::isPointIsolatedWithinNeighborhood()


//--------------------------------------------------------------------------
template <typename Coord>
template <typename PointIter>
std::vector<size_t>
lar::example::PointIsolationAlg<Coord>::bruteRemoveIsolatedPoints
  (PointIter begin, PointIter end) const
{
  //
  // reference implementation: brute force dumb one
  //
  
  std::vector<size_t> nonIsolated;
  
  size_t i = 0;
  for (auto it = begin; it != end; ++it, ++i) {
    
    for (auto ioth = begin; ioth != end; ++ioth) {
      if (it == ioth) continue;
      
      if (closeEnough(*it, *ioth)) {
        nonIsolated.push_back(i);
        break;
      }
      
    } // for oth
    
  } // for (it)
  
  return nonIsolated;
} // lar::example::PointIsolationAlg::bruteRemoveIsolatedPoints()


//--------------------------------------------------------------------------
template <typename Coord>
template <typename Point>
bool lar::example::PointIsolationAlg<Coord>::closeEnough
  (Point const& A, Point const& B) const
{
  return cet::sum_of_squares(
    details::extractPositionX(A) - details::extractPositionX(B),
    details::extractPositionY(A) - details::extractPositionY(B),
    details::extractPositionZ(A) - details::extractPositionZ(B)
    ) <= config.radius2;
} // lar::example::PointIsolationAlg<Point>::closeEnough()


//--------------------------------------------------------------------------
template <typename Coord>
std::string lar::example::PointIsolationAlg<Coord>::rangeString
  (Coord_t from, Coord_t to)
  { return "(" + std::to_string(from) + " to " + std::to_string(to) + ")"; }


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

#endif // LAREXAMPLES_ALGORITHMS_REMOVEISOLATEDSPACEPOINTS_POINTISOLATIONALG_H