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 // Simple wrapper around boost::disjoint_sets. In theory, could use
0023 // boost::vector_property_map and use boost::disjoint_sets without
0024 // wrapping, but it's way slower
0025 class DisjointSets {
0026  public:
0027   explicit DisjointSets(std::size_t initial_size = 128)
0028       : m_defaultSize(initial_size),
0029         m_size(initial_size),
0030         m_rank(m_size),
0031         m_parent(m_size),
0032         m_ds(&m_rank[0], &m_parent[0]) {}
0033
0034   Acts::Ccl::Label makeSet() {
0035     // Empirically, m_size = 128 seems to be good default. If we
0036     // exceed this, take a performance hit and do the right thing.
0037     while (m_globalId >= m_size) {
0038       m_size *= 2;
0039       m_rank.resize(m_size);
0040       m_parent.resize(m_size);
0041       m_ds = boost::disjoint_sets<std::size_t*, std::size_t*>(&m_rank[0],
0042                                                               &m_parent[0]);
0043     }
0044     m_ds.make_set(m_globalId);
0045     return static_cast<Acts::Ccl::Label>(m_globalId++);
0046   }
0047
0048   void unionSet(std::size_t x, std::size_t y) { m_ds.union_set(x, y); }
0049   Acts::Ccl::Label findSet(std::size_t x) {
0050     return static_cast<Acts::Ccl::Label>(m_ds.find_set(x));
0051   }
0052
0053   void clear() {
0054     m_size = m_defaultSize;
0055     m_rank.clear();
0056     m_parent.clear();
0057     m_rank.resize(m_size);
0058     m_parent.resize(m_size);
0059     m_globalId = 1;
0060     m_ds = boost::disjoint_sets<std::size_t*, std::size_t*>(&m_rank[0],
0061                                                             &m_parent[0]);
0062   }
0063
0064  private:
0065   std::size_t m_defaultSize{128};
0066   std::size_t m_globalId = 1;
0067   std::size_t m_size{m_defaultSize};
0068   std::vector<std::size_t> m_rank;
0069   std::vector<std::size_t> m_parent;
0070   boost::disjoint_sets<std::size_t*, std::size_t*> m_ds;
0071 };
0072
0073 struct ClusteringData {
0074   void clear() {
0075     labels.clear();
0076     nClusters.clear();
0077     ds.clear();
0078   }
0079
0080   std::vector<Acts::Ccl::Label> labels{};
0081   std::vector<std::size_t> nClusters{};
0082   Acts::Ccl::DisjointSets ds{};
0083 };
0084
0085 template <typename Cell>
0086 concept HasRetrievableColumnInfo = requires(Cell cell) {
0087   { getCellColumn(cell) } -> std::same_as<int>;
0088 };
0089
0090 template <typename Cell>
0091 concept HasRetrievableRowInfo = requires(Cell cell) {
0092   { getCellRow(cell) } -> std::same_as<int>;
0093 };
0094
0095 template <typename Cell, typename Cluster>
0096 concept CanAcceptCell = requires(Cell cell, Cluster cluster) {
0097   { clusterAddCell(cluster, cell) } -> std::same_as<void>;
0098 };
0099
0100 template <typename Cluster>
0101 concept CanReserve = requires(Cluster cluster, std::size_t n) {
0102   { clusterReserve(cluster, n) } -> std::same_as<void>;
0103 };
0104
0105 // When looking for a cell connected to a reference cluster, the code
0106 // always loops backward, starting from the reference cell. Since
0107 // the cells are globally sorted column-wise, the connection function
0108 // can therefore tell when the search should be stopped.
0109 enum class ConnectResult {
0110   eNoConn,      // No connections, keep looking
0111   eNoConnStop,  // No connections, stop looking
0112   eConn,        // Found connection
0113   eDuplicate    // Found duplicate cell, throw an exception
0114 };
0115
0116 // Default connection type for 2-D grids: 4- or 8-cell connectivity
0117 template <typename Cell>
0118   requires(Acts::Ccl::HasRetrievableColumnInfo<Cell> &&
0119            Acts::Ccl::HasRetrievableRowInfo<Cell>)
0120 struct Connect2D {
0121   bool conn8{true};
0122   Connect2D() = default;
0123   explicit Connect2D(bool commonCorner) : conn8{commonCorner} {}
0124   virtual ConnectResult operator()(const Cell& ref, const Cell& iter) const;
0125   virtual ~Connect2D() = default;
0126 };
0127
0128 // Default connection type for 1-D grids: 2-cell connectivity
0129 template <Acts::Ccl::HasRetrievableColumnInfo Cell>
0130 struct Connect1D {
0131   virtual ConnectResult operator()(const Cell& ref, const Cell& iter) const;
0132   virtual ~Connect1D() = default;
0133 };
0134
0135 // Default connection type based on GridDim
0136 template <typename Cell, std::size_t GridDim = 2>
0137 struct DefaultConnect {
0138   static_assert(GridDim != 1 && GridDim != 2,
0139                 "Only grid dimensions of 1 or 2 are supported");
0140 };
0141
0142 template <typename Cell>
0143 struct DefaultConnect<Cell, 1> : public Connect1D<Cell> {
0144   ~DefaultConnect() override = default;
0145 };
0146
0147 template <typename Cell>
0148 struct DefaultConnect<Cell, 2> : public Connect2D<Cell> {
0149   explicit DefaultConnect(bool commonCorner) : Connect2D<Cell>(commonCorner) {}
0150   DefaultConnect() = default;
0151   ~DefaultConnect() override = default;
0152 };
0153
0154 /// @brief labelClusters
0155 ///
0156 /// In-place connected component labelling using the Hoshen-Kopelman algorithm.
0157 /// The `Cell` type must have the following functions defined:
0158 ///   int  getCellRow(const Cell&),
0159 ///   int  getCellColumn(const Cell&)
0160 ///
0161 /// @param [in] data collection of quantities for clusterization
0162 /// @param [in] cells the cell collection to be labeled
0163 /// @param [in] connect the connection type (see DefaultConnect)
0164 /// @throws std::invalid_argument if the input contains duplicate cells.
0165 template <typename CellCollection, std::size_t GridDim = 2,
0166           typename Connect =
0167               DefaultConnect<typename CellCollection::value_type, GridDim>>
0168 void labelClusters(Acts::Ccl::ClusteringData& data, CellCollection& cells,
0169                    Connect&& connect = Connect());
0170
0171 /// @brief mergeClusters
0172 ///
0173 /// Merge a set of cells previously labeled (for instance with `labelClusters`)
0174 /// into actual clusters. The Cluster type must have the following function
0175 /// defined:
0176 ///   void clusterAddCell(Cluster&, const Cell&)
0177 template <typename CellCollection, typename ClusterCollection>
0178   requires(Acts::Ccl::CanAcceptCell<typename CellCollection::value_type,
0179                                     typename ClusterCollection::value_type>)
0180 void mergeClusters(Acts::Ccl::ClusteringData& data, const CellCollection& cells,
0181                    ClusterCollection& outv);
0182
0183 /// @brief createClusters
0184 /// Convenience function which runs both labelClusters and createClusters.
0185 template <typename CellCollection, typename ClusterCollection,
0186           std::size_t GridDim = 2,
0187           typename Connect =
0188               DefaultConnect<typename CellCollection::value_type, GridDim>>
0189 [[deprecated]]
0190 ClusterCollection createClusters(CellCollection& cells,
0191                                  Connect&& connect = Connect());
0192
0193 /// @brief createClusters
0194 /// Alternative convenience function which runs both labelClusters and
0195 /// createClusters.
0196 ///
0197 /// @throws std::invalid_argument if the input contains duplicate cells.
0198 template <typename CellCollection, typename ClusterCollection,
0199           std::size_t GridDim = 2,
0200           typename Connect =
0201               DefaultConnect<typename CellCollection::value_type, GridDim>>
0202   requires(GridDim == 1 || GridDim == 2)
0203 void createClusters(Acts::Ccl::ClusteringData& data, CellCollection& cells,
0204                     ClusterCollection& clusters, Connect&& connect = Connect());
0205
0206 }  // namespace Acts::Ccl
0207
0208 #include "Acts/Clusterization/Clusterization.ipp"