Acts/Clusterization/Clusterization.hpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #pragma once
0010
0011 #include <memory>
0012 #include <vector>
0013
0014 #include <boost/pending/disjoint_sets.hpp>
0015
0016 namespace Acts::Ccl {
0017 using Label = int;
0018 constexpr Label NO_LABEL = 0;
0019 }  // namespace Acts::Ccl
0020
0021 namespace Acts::Ccl {
0022 /// Wrapper around `boost::disjoint_sets` used by the clustering helpers.
0023 ///
0024 /// Avoids the overhead of `boost::vector_property_map` while keeping the
0025 /// interface small for the clustering workflow.
0026 class DisjointSets {
0027  public:
0028   /// @brief Constructor
0029   /// @param initial_size Initial size of the disjoint sets
0030   explicit DisjointSets(std::size_t initial_size = 128)
0031       : m_defaultSize(initial_size),
0032         m_size(initial_size),
0033         m_rank(m_size),
0034         m_parent(m_size),
0035         m_ds(&m_rank[0], &m_parent[0]) {}
0036
0037   /// @brief Create a new set
0038   /// @return The label of the new set
0039   Acts::Ccl::Label makeSet() {
0040     // Empirically, m_size = 128 seems to be good default. If we
0041     // exceed this, take a performance hit and do the right thing.
0042     while (m_globalId >= m_size) {
0043       m_size *= 2;
0044       m_rank.resize(m_size);
0045       m_parent.resize(m_size);
0046       m_ds = boost::disjoint_sets<std::size_t*, std::size_t*>(&m_rank[0],
0047                                                               &m_parent[0]);
0048     }
0049     m_ds.make_set(m_globalId);
0050     return static_cast<Acts::Ccl::Label>(m_globalId++);
0051   }
0052
0053   /// @brief Union two sets
0054   /// @param x First set
0055   /// @param y Second set
0056   void unionSet(std::size_t x, std::size_t y) { m_ds.union_set(x, y); }
0057   /// @brief Find the representative of a set
0058   /// @param x Element to find
0059   /// @return The label of the set containing x
0060   Acts::Ccl::Label findSet(std::size_t x) {
0061     return static_cast<Acts::Ccl::Label>(m_ds.find_set(x));
0062   }
0063
0064   /// @brief Clear all sets and reset to initial state
0065   void clear() {
0066     m_size = m_defaultSize;
0067     m_rank.clear();
0068     m_parent.clear();
0069     m_rank.resize(m_size);
0070     m_parent.resize(m_size);
0071     m_globalId = 1;
0072     m_ds = boost::disjoint_sets<std::size_t*, std::size_t*>(&m_rank[0],
0073                                                             &m_parent[0]);
0074   }
0075
0076  private:
0077   std::size_t m_defaultSize{128};
0078   std::size_t m_globalId = 1;
0079   std::size_t m_size{m_defaultSize};
0080   std::vector<std::size_t> m_rank;
0081   std::vector<std::size_t> m_parent;
0082   boost::disjoint_sets<std::size_t*, std::size_t*> m_ds;
0083 };
0084
0085 /// Collection of mutable data used by clustering utilities.
0086 struct ClusteringData {
0087   /// Clear all clustering data
0088   void clear() {
0089     labels.clear();
0090     nClusters.clear();
0091     ds.clear();
0092   }
0093
0094   /// Cluster labels for each cell
0095   std::vector<Acts::Ccl::Label> labels{};
0096   /// Number of cells per cluster
0097   std::vector<std::size_t> nClusters{};
0098   /// Disjoint sets data structure for clustering
0099   Acts::Ccl::DisjointSets ds{};
0100 };
0101
0102 template <typename Cell>
0103 concept HasRetrievableColumnInfo = requires(Cell cell) {
0104   { getCellColumn(cell) } -> std::same_as<int>;
0105 };
0106
0107 template <typename Cell>
0108 concept HasRetrievableRowInfo = requires(Cell cell) {
0109   { getCellRow(cell) } -> std::same_as<int>;
0110 };
0111
0112 template <typename Cell, typename Cluster>
0113 concept CanAcceptCell = requires(Cell cell, Cluster cluster) {
0114   { clusterAddCell(cluster, cell) } -> std::same_as<void>;
0115 };
0116
0117 template <typename Cluster>
0118 concept CanReserve = requires(Cluster cluster, std::size_t n) {
0119   { clusterReserve(cluster, n) } -> std::same_as<void>;
0120 };
0121
0122 // When looking for a cell connected to a reference cluster, the code
0123 // always loops backward, starting from the reference cell. Since
0124 // the cells are globally sorted column-wise, the connection function
0125 // can therefore tell when the search should be stopped.
0126 enum class ConnectResult {
0127   eNoConn,      // No connections, keep looking
0128   eNoConnStop,  // No connections, stop looking
0129   eConn,        // Found connection
0130   eDuplicate    // Found duplicate cell, throw an exception
0131 };
0132
0133 /// Default connection type for 2-D grids with 4- or 8-cell connectivity.
0134 template <typename Cell>
0135   requires(Acts::Ccl::HasRetrievableColumnInfo<Cell> &&
0136            Acts::Ccl::HasRetrievableRowInfo<Cell>)
0137 struct Connect2D {
0138   /// Whether to use 8-cell connectivity (true) or 4-cell connectivity (false)
0139   bool conn8{true};
0140   Connect2D() = default;
0141   /// @brief Constructor
0142   /// @param commonCorner Whether cells sharing only a corner are considered connected
0143   explicit Connect2D(bool commonCorner) : conn8{commonCorner} {}
0144   /// @brief Check if two cells are connected
0145   /// @param ref Reference cell
0146   /// @param iter Candidate cell
0147   /// @return Connection result
0148   virtual ConnectResult operator()(const Cell& ref, const Cell& iter) const;
0149   virtual ~Connect2D() = default;
0150 };
0151
0152 /// Default connection type for 1-D grids with 2-cell connectivity.
0153 template <Acts::Ccl::HasRetrievableColumnInfo Cell>
0154 struct Connect1D {
0155   /// @brief Check if two cells are connected
0156   /// @param ref Reference cell
0157   /// @param iter Candidate cell
0158   /// @return Connection result
0159   virtual ConnectResult operator()(const Cell& ref, const Cell& iter) const;
0160   virtual ~Connect1D() = default;
0161 };
0162
0163 // Default connection type based on GridDim
0164 template <typename Cell, std::size_t GridDim = 2>
0165 struct DefaultConnect {
0166   static_assert(GridDim != 1 && GridDim != 2,
0167                 "Only grid dimensions of 1 or 2 are supported");
0168 };
0169
0170 /// Default connection type for 1-D grids.
0171 template <typename Cell>
0172 struct DefaultConnect<Cell, 1> : public Connect1D<Cell> {
0173   ~DefaultConnect() override = default;
0174 };
0175
0176 /// Default connection type for 2-D grids.
0177 template <typename Cell>
0178 struct DefaultConnect<Cell, 2> : public Connect2D<Cell> {
0179   /// @brief Constructor
0180   /// @param commonCorner Whether cells sharing only a corner are considered connected
0181   explicit DefaultConnect(bool commonCorner) : Connect2D<Cell>(commonCorner) {}
0182   DefaultConnect() = default;
0183   ~DefaultConnect() override = default;
0184 };
0185
0186 /// @brief labelClusters
0187 ///
0188 /// In-place connected component labelling using the Hoshen-Kopelman algorithm.
0189 /// The `Cell` type must have the following functions defined:
0190 ///   int  getCellRow(const Cell&),
0191 ///   int  getCellColumn(const Cell&)
0192 ///
0193 /// @param [in] data collection of quantities for clusterization
0194 /// @param [in] cells the cell collection to be labeled
0195 /// @param [in] connect the connection type (see DefaultConnect)
0196 /// @throws std::invalid_argument if the input contains duplicate cells.
0197 template <typename CellCollection, std::size_t GridDim = 2,
0198           typename Connect =
0199               DefaultConnect<typename CellCollection::value_type, GridDim>>
0200 void labelClusters(Acts::Ccl::ClusteringData& data, CellCollection& cells,
0201                    Connect&& connect = Connect());
0202
0203 /// @brief mergeClusters
0204 ///
0205 /// Merge a set of cells previously labeled (for instance with `labelClusters`)
0206 /// into actual clusters. The Cluster type must have the following function
0207 /// defined:
0208 ///   void clusterAddCell(Cluster&, const Cell&)
0209 template <typename CellCollection, typename ClusterCollection>
0210   requires(Acts::Ccl::CanAcceptCell<typename CellCollection::value_type,
0211                                     typename ClusterCollection::value_type>)
0212 void mergeClusters(Acts::Ccl::ClusteringData& data, const CellCollection& cells,
0213                    ClusterCollection& outv);
0214
0215 /// @brief createClusters
0216 /// Convenience function which runs both labelClusters and mergeClusters.
0217 ///
0218 /// @throws std::invalid_argument if the input contains duplicate cells.
0219 template <typename CellCollection, typename ClusterCollection,
0220           std::size_t GridDim = 2,
0221           typename Connect =
0222               DefaultConnect<typename CellCollection::value_type, GridDim>>
0223   requires(GridDim == 1 || GridDim == 2)
0224 void createClusters(Acts::Ccl::ClusteringData& data, CellCollection& cells,
0225                     ClusterCollection& clusters, Connect&& connect = Connect());
0226
0227 }  // namespace Acts::Ccl
0228
0229 #include "Acts/Clusterization/Clusterization.ipp"