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 };
0114
0115 // Default connection type for 2-D grids: 4- or 8-cell connectivity
0116 template <typename Cell>
0117   requires(Acts::Ccl::HasRetrievableColumnInfo<Cell> &&
0118            Acts::Ccl::HasRetrievableRowInfo<Cell>)
0119 struct Connect2D {
0120   bool conn8{true};
0121   Connect2D() = default;
0122   explicit Connect2D(bool commonCorner) : conn8{commonCorner} {}
0123   virtual ConnectResult operator()(const Cell& ref, const Cell& iter) const;
0124   virtual ~Connect2D() = default;
0125 };
0126
0127 // Default connection type for 1-D grids: 2-cell connectivity
0128 template <Acts::Ccl::HasRetrievableColumnInfo Cell>
0129 struct Connect1D {
0130   virtual ConnectResult operator()(const Cell& ref, const Cell& iter) const;
0131   virtual ~Connect1D() = default;
0132 };
0133
0134 // Default connection type based on GridDim
0135 template <typename Cell, std::size_t GridDim = 2>
0136 struct DefaultConnect {
0137   static_assert(GridDim != 1 && GridDim != 2,
0138                 "Only grid dimensions of 1 or 2 are supported");
0139 };
0140
0141 template <typename Cell>
0142 struct DefaultConnect<Cell, 1> : public Connect1D<Cell> {
0143   ~DefaultConnect() override = default;
0144 };
0145
0146 template <typename Cell>
0147 struct DefaultConnect<Cell, 2> : public Connect2D<Cell> {
0148   explicit DefaultConnect(bool commonCorner) : Connect2D<Cell>(commonCorner) {}
0149   DefaultConnect() = default;
0150   ~DefaultConnect() override = default;
0151 };
0152
0153 /// @brief labelClusters
0154 ///
0155 /// In-place connected component labelling using the Hoshen-Kopelman algorithm.
0156 /// The `Cell` type must have the following functions defined:
0157 ///   int  getCellRow(const Cell&),
0158 ///   int  getCellColumn(const Cell&)
0159 ///
0160 /// @param [in] data collection of quantities for clusterization
0161 /// @param [in] cells the cell collection to be labeled
0162 /// @param [in] connect the connection type (see DefaultConnect)
0163 template <typename CellCollection, std::size_t GridDim = 2,
0164           typename Connect =
0165               DefaultConnect<typename CellCollection::value_type, GridDim>>
0166 void labelClusters(Acts::Ccl::ClusteringData& data, CellCollection& cells,
0167                    Connect&& connect = Connect());
0168
0169 /// @brief mergeClusters
0170 ///
0171 /// Merge a set of cells previously labeled (for instance with `labelClusters`)
0172 /// into actual clusters. The Cluster type must have the following function
0173 /// defined:
0174 ///   void clusterAddCell(Cluster&, const Cell&)
0175 ///
0176 /// @return nothing
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 template <typename CellCollection, typename ClusterCollection,
0197           std::size_t GridDim = 2,
0198           typename Connect =
0199               DefaultConnect<typename CellCollection::value_type, GridDim>>
0200   requires(GridDim == 1 || GridDim == 2)
0201 void createClusters(Acts::Ccl::ClusteringData& data, CellCollection& cells,
0202                     ClusterCollection& clusters, Connect&& connect = Connect());
0203
0204 }  // namespace Acts::Ccl
0205
0206 #include "Acts/Clusterization/Clusterization.ipp"