lar::example::PointIsolationAlg< Coord > Class Template Reference

Algorithm to detect isolated space points. More...

#include <PointIsolationAlg.h>

Classes
struct	Configuration_t
	Type containing all configuration parameters of the algorithm. More...

Public Types
using	Coord_t = Coord
	Type of coordinate. More...
using	Range_t = CoordRange< Coord_t >

Public Member Functions
template<typename PointIter >
std::vector< size_t >	removeIsolatedPoints (PointIter begin, PointIter end) const
	Returns the set of points that are not isolated. More...
template<typename Cont >
std::vector< size_t >	removeIsolatedPoints (Cont const &points) const
	Returns the set of points that are not isolated. More...
template<typename PointIter >
std::vector< size_t >	bruteRemoveIsolatedPoints (PointIter begin, PointIter end) const
	Brute-force reference algorithm. More...
template<typename PointIter >
Coord	computeCellSize () const

Static Public Member Functions
static Coord_t	maximumOptimalCellSize (Coord_t radius)
	Returns the maximum optimal cell size when using a isolation radius. More...

Private Types
using	Indexer_t = ::util::GridContainer3DIndices
	type managing cell indices More...
using	NeighAddresses_t = std::vector< Indexer_t::CellIndexOffset_t >
	type of neighbourhood cell offsets More...
template<typename PointIter >
using	Partition_t = SpacePartition< PointIter >
template<typename PointIter >
using	Point_t = decltype(*PointIter())

Private Member Functions
template<typename PointIter = std::array<double, 3> const*>
Coord_t	computeCellSize () const
	Computes the cell size to be used. More...
NeighAddresses_t	buildNeighborhood (Indexer_t const &indexer, unsigned int neighExtent) const
	Returns a list of cell offsets for the neighbourhood of given radius. More...
template<typename PointIter >
bool	isPointIsolatedFrom (Point_t< PointIter > const &point, typename Partition_t< PointIter >::Cell_t const &otherPoints) const
	Returns whether a point is isolated with respect to all the others. More...
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
	Returns whether a point is isolated in the specified neighbourhood. More...
template<typename Point >
bool	closeEnough (Point const &A, Point const &B) const
	Returns whether A and B are close enough to be considered non-isolated. More...

Static Private Member Functions
static std::string	rangeString (Coord_t from, Coord_t to)
	Helper function. Returns a string "(\<from\> to \<to\>)" More...
static std::string	rangeString (Range_t range)
	Helper function. Returns a string "(\<from\> to \<to\>)" More...

Private Attributes
Configuration_t	config
	all configuration data More...

Configuration
	PointIsolationAlg (Configuration_t const &first_config)
	Constructor with configuration validation. More...
void	reconfigure (Configuration_t const &new_config)
Configuration_t &	configuration () const
static void	validateConfiguration (Configuration_t const &config)

Detailed Description

template<typename Coord = double>
class lar::example::PointIsolationAlg< Coord >

Algorithm to detect isolated space points.

Template Parameters

Coord

type of the coordinate

See also

RemoveIsolatedSpacePoints example overview

This algorithm returns a selection of the input points which are not isolated. Point is defined as isolated if: where describes the position of the point in space and is the isolation radius.

This class must be configured by providing a complete Configuration_t object. Configuration can be changed at any time after that.

The configuration information (Configuration_t) defines the volume the points span and the square of the isolation radius. The information on the volume may be used to optimise the algorithm, and it is not checked. If that information is wrong (that means input points lie outside that volume), the result is undefined. No check is automatically performed to assess if the configuration is valid.

The algorithm can be run on any collection of points, as long as the point class supports the PositionExtractor class. A typical cycle of use is:

// creation and configuration
lar::examples::PointIsolationAlg<float>::Configuration_t config;
config.rangeX = { -1., 1. };
config.rangeY = { -1., 1. };
config.rangeZ = { -5., 5. };
config.radius2 = 0.25;

lar::examples::PointIsolationAlg<float> algo(config);

// preparation/retrieval of input
std::vector<std::array<float, 3>> points;
// points are filled here

// execution
std::vector<size_t> indices
  = algo.removeIsolatedPoints(points.begin(), points.end());

// utilization of the result;
// - e.g., create a collection of non-isolated points...
std::vector<std::array<float, 3>> nonIsolatedPoints;
nonIsolatedPoints.reserve(indices.size());
for (size_t index: indices)
  nonIsolatedPoints.push_back(points[index]);

// - ... or their pointers
std::vector<std::array<float, 3> const*> nonIsolatedPointPtrs;
nonIsolatedPointPtrs.reserve(indices.size());
for (size_t index: indices)
  nonIsolatedPointPtrs.push_back(&(points[index]));

The point type here is std::array<float, 3>, for which a lar::examples::PositionExtractor<std::array<float, 3>> is defined in this same library. The algorithm can be executed multiple times, and the configuration can be changed at any time (reconfigure()).

Validation of the configuration is optional, and needs to be explicitly called if desired (validateConfiguration()).

Description of the algorithm

The basic method to determine the isolation of a point is by brute force, by computing the distance with all others and, as soon as one of them is found too close, declare the point non-isolated.

A refinement is implemented: the points are grouped in cubic "cells" and points in cells that are farther than isolation radius are not checked against each other. This requires some memory to allocate the structure, that can become huge. The maximum memory parameter keeps this sane.

Other refinements are not implemented. When a point is found non-isolated also the point that makes it non-isolated should also be marked so. Cell radius might be tuned to be smaller. Some of the neighbour cells may be too far and should not be checked. The grid allocates a vector for each cell, whether it's empty or not; using a sparse structure might reduce the memory; also if the grid contains pointers to vectors instead of vectors, and the grid is very sparse, there should still be some memory saving.

Definition at line 127 of file PointIsolationAlg.h.

Member Typedef Documentation

template<typename Coord = double>

using lar::example::PointIsolationAlg< Coord >::Coord_t = Coord

Type of coordinate.

Definition at line 131 of file PointIsolationAlg.h.

template<typename Coord = double>

using lar::example::PointIsolationAlg< Coord >::Indexer_t = ::util::GridContainer3DIndices

private

type managing cell indices

Definition at line 241 of file PointIsolationAlg.h.

template<typename Coord = double>

using lar::example::PointIsolationAlg< Coord >::NeighAddresses_t = std::vector<Indexer_t::CellIndexOffset_t>

private

type of neighbourhood cell offsets

Definition at line 244 of file PointIsolationAlg.h.

template<typename Coord = double>

template<typename PointIter >

using lar::example::PointIsolationAlg< Coord >::Partition_t = SpacePartition<PointIter>

private

Definition at line 247 of file PointIsolationAlg.h.

template<typename Coord = double>

template<typename PointIter >

using lar::example::PointIsolationAlg< Coord >::Point_t = decltype(*PointIter())

private

Definition at line 250 of file PointIsolationAlg.h.

template<typename Coord = double>

using lar::example::PointIsolationAlg< Coord >::Range_t = CoordRange<Coord_t>

Definition at line 132 of file PointIsolationAlg.h.

Constructor & Destructor Documentation

template<typename Coord = double>

lar::example::PointIsolationAlg< Coord >::PointIsolationAlg ( Configuration_t const &  first_config )

inline

Constructor with configuration validation.

Parameters

first_config

configuration parameter structure

For the configuration, see SpacePointIsolationAlg documentation. No validation is performed on the configuration.

Definition at line 155 of file PointIsolationAlg.h.

        : config(first_config)
        {}

Member Function Documentation

template<typename Coord >

template<typename PointIter >

std::vector< size_t > lar::example::PointIsolationAlg< Coord >::bruteRemoveIsolatedPoints	(	PointIter	begin,
		PointIter	end
	)		const

Brute-force reference algorithm.

Template Parameters

PointIter

random access iterator to a point type

Parameters

begin	iterator to the first point to be considered
end	iterator after the last point to be considered

Returns

a list of indices of non-isolated points in the input range

See also

removeIsolatedPoints

This algorithm executes the task in a way, slow and supposedly reliable. The interface is the same as removeIsolatedPoints. Use this only for tests.

Definition at line 575 of file PointIsolationAlg.h.

{
  //
  // 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 >

lar::example::PointIsolationAlg< Coord >::NeighAddresses_t lar::example::PointIsolationAlg< Coord >::buildNeighborhood ( Indexer_t const &  indexer, unsigned int  neighExtent  ) const

private

Returns a list of cell offsets for the neighbourhood of given radius.

Definition at line 482 of file PointIsolationAlg.h.

{
  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 Point >

bool lar::example::PointIsolationAlg< Coord >::closeEnough ( Point const &  A, Point const &  B  ) const

private

Returns whether A and B are close enough to be considered non-isolated.

Definition at line 606 of file PointIsolationAlg.h.

{
  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 = double>

template<typename PointIter = std::array<double, 3> const*>

Coord_t lar::example::PointIsolationAlg< Coord >::computeCellSize ( ) const

private

Computes the cell size to be used.

template<typename Coord = double>

template<typename PointIter >

Coord lar::example::PointIsolationAlg< Coord >::computeCellSize

(

)

const

Definition at line 439 of file PointIsolationAlg.h.

                                                                {

  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 = double>

Configuration_t& lar::example::PointIsolationAlg< Coord >::configuration ( ) const

inline

Returns a constant reference to the current configuration

See also

reconfigure()

Definition at line 168 of file PointIsolationAlg.h.

{ return config; }

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

private

Returns whether a point is isolated with respect to all the others.

Definition at line 522 of file PointIsolationAlg.h.

{

  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

private

Returns whether a point is isolated in the specified neighbourhood.

Definition at line 542 of file PointIsolationAlg.h.

{

  // 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 = double>

static Coord_t lar::example::PointIsolationAlg< Coord >::maximumOptimalCellSize ( Coord_t  radius )

inlinestatic

Returns the maximum optimal cell size when using a isolation radius.

Definition at line 235 of file PointIsolationAlg.h.

        { return radius / std::sqrt(3.); }

template<typename Coord >

std::string lar::example::PointIsolationAlg< Coord >::rangeString ( Coord_t  from, Coord_t  to  )

staticprivate

Helper function. Returns a string "(\<from\> to \<to\>)"

Definition at line 619 of file PointIsolationAlg.h.

  { return "(" + std::to_string(from) + " to " + std::to_string(to) + ")"; }

template<typename Coord = double>

static std::string lar::example::PointIsolationAlg< Coord >::rangeString ( Range_t  range )

inlinestaticprivate

Helper function. Returns a string "(\<from\> to \<to\>)"

Definition at line 290 of file PointIsolationAlg.h.

        { return rangeString(range.lower, range.upper); }

template<typename Coord = double>

void lar::example::PointIsolationAlg< Coord >::reconfigure ( Configuration_t const &  new_config )

inline

Reconfigures the algorithm with the specified configuration (no validation is performed)

See also

configuration()

Definition at line 162 of file PointIsolationAlg.h.

        { config = new_config; }

template<typename Coord >

template<typename PointIter >

std::vector< size_t > lar::example::PointIsolationAlg< Coord >::removeIsolatedPoints	(	PointIter	begin,
		PointIter	end
	)		const

Returns the set of points that are not isolated.

Template Parameters

PointIter

random access iterator to a point type

Parameters

begin	iterator to the first point to be considered
end	iterator after the last point to be considered

Returns

a list of indices of non-isolated points in the input range

This method is the operating core of the algorithm.

The input is two iterators. The output is a collection of the indices of the elements that are not isolated. The index is equivalent to std::distance(begin, point). The order of the elements in the collection is not specified.

This method can use any collection of input data, as long as a PositionExtractor object is available for it.

Definition at line 311 of file PointIsolationAlg.h.

{

  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 = double>

template<typename Cont >

std::vector<size_t> lar::example::PointIsolationAlg< Coord >::removeIsolatedPoints ( Cont const &  points ) const

inline

Returns the set of points that are not isolated.

Parameters

points

list of the reconstructed space points

Returns

a list of indices of non-isolated points in the vector

See also

removeIsolatedPoints(PointIter begin, PointIter end) const

Definition at line 202 of file PointIsolationAlg.h.

        { return removeIsolatedPoints(std::cbegin(points), std::cend(points)); }

template<typename Coord >

void lar::example::PointIsolationAlg< Coord >::validateConfiguration ( Configuration_t const &  config )

static

Throws an exception if the configuration is invalid

Exceptions

std::runtime_error

if configuration is invalid

Definition at line 407 of file PointIsolationAlg.h.

{
  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()

Member Data Documentation

template<typename Coord = double>

Configuration_t lar::example::PointIsolationAlg< Coord >::config

private

all configuration data

Definition at line 253 of file PointIsolationAlg.h.

---

The documentation for this class was generated from the following file:

• larexamples/larexamples/Algorithms/RemoveIsolatedSpacePoints/PointIsolationAlg.h