Back to home page

EIC code displayed by LXR

 
 

    


Warning, file /include/Acts/Seeding/HoughTransformUtils.hpp was not indexed or was modified since last indexation (in which case cross-reference links may be missing, inaccurate or erroneous).

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2022 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 http://mozilla.org/MPL/2.0/.
0008 
0009 /// This file implements the tools for a hough transform.
0010 
0011 #pragma once
0012 #include "Acts/Utilities/Delegate.hpp"
0013 #include "Acts/Utilities/Logger.hpp"
0014 #include "Acts/Utilities/Result.hpp"
0015 
0016 #include <array>
0017 #include <map>
0018 #include <optional>
0019 #include <set>
0020 #include <unordered_set>
0021 
0022 #include "HoughVectors.hpp"
0023 
0024 namespace Acts::HoughTransformUtils {
0025 
0026 /// this type is responsible for encoding the parameters of our hough space
0027 using CoordType = double;
0028 
0029 // this type is used to encode hit counts.
0030 // Floating point to allow hit weights to be applied
0031 using YieldType = float;
0032 
0033 /// @brief this function represents a mapping of a coordinate point in detector space to a line in
0034 /// hough space. Given the value of the first hough coordinate, it shall return
0035 /// the corresponding second coordinate according to the line parametrisation.
0036 /// Should be implemented by the user.
0037 /// @tparam PointType: The class representing a point in detector space (can differ between implementations)
0038 template <class PointType>
0039 using LineParametrisation =
0040     std::function<CoordType(CoordType, const PointType&)>;
0041 
0042 /// @brief struct to define the ranges of the hough histogram.
0043 /// Used to move between parameter and bin index coordinates.
0044 /// Disconnected from the hough plane binning to be able to re-use
0045 /// a plane with a given binning for several parameter ranges
0046 struct HoughAxisRanges {
0047   CoordType xMin = 0.0f;  // minimum value of the first coordinate
0048   CoordType xMax = 0.0f;  // maximum value of the first coordinate
0049   CoordType yMin = 0.0f;  // minimum value of the second coordinate
0050   CoordType yMax = 0.0f;  // maximum value of the second coordinate
0051 };
0052 
0053 /// convenience functions to link bin indices to axis coordinates
0054 
0055 /// @brief find the bin index corresponding to a certain abscissa
0056 /// of the coordinate axis, based on the axis limits and binning.
0057 /// @param min: Start of axis range
0058 /// @param max: End of axis range
0059 /// @param nSteps: Number of bins in axis
0060 /// @param val: value to find the corresponding bin for
0061 /// @return the bin number.
0062 /// No special logic to prevent over-/underflow, checking these is
0063 /// left to the caller
0064 inline int binIndex(double min, double max, unsigned nSteps, double val) {
0065   return static_cast<int>((val - min) / (max - min) * nSteps);
0066 }
0067 // Returns the lower bound of the bin specified by step
0068 /// @param min: Start of axis range
0069 /// @param max: End of axis range
0070 /// @param nSteps: Number of bins in axis
0071 /// @param binIndex: The index of the bin
0072 /// @return the parameter value at the lower bin edge.
0073 /// No special logic to prevent over-/underflow, checking these is
0074 /// left to the caller
0075 inline double lowerBinEdge(double min, double max, unsigned nSteps,
0076                            std::size_t binIndex) {
0077   return min + (max - min) * binIndex / nSteps;
0078 }
0079 // Returns the lower bound of the bin specified by step
0080 /// @param min: Start of axis range
0081 /// @param max: End of axis range
0082 /// @param nSteps: Number of bins in axis
0083 /// @param binIndex: The index of the bin
0084 /// @return the parameter value at the bin center.
0085 /// No special logic to prevent over-/underflow, checking these is
0086 /// left to the caller
0087 inline double binCenter(double min, double max, unsigned nSteps,
0088                         std::size_t binIndex) {
0089   return min + (max - min) * 0.5 * (2 * binIndex + 1) / nSteps;
0090 }
0091 
0092 /// @brief data class for the information to store for each
0093 /// cell of the hough histogram.
0094 /// @tparam identifier_t: Type of the identifier to associate to the hits
0095 ///                       Should be sortable. Used to uniquely identify each
0096 ///                       hit and to eventually return the list of hits per cell
0097 template <class identifier_t>
0098 class HoughCell {
0099  public:
0100   /// @brief construct the cell as empty
0101   HoughCell() = default;
0102   /// @brief add an entry to this cell
0103   /// @param identifier: Identifier of the hit (used to distinguish hits from another)
0104   /// @param layer: Layer of the hit (used when counting layers)
0105   /// @param weight: Optional weight to assign to the hit
0106   void fill(const identifier_t& identifier, unsigned int layer,
0107             YieldType weight = 1.);
0108   /// @brief access the number of layers with hits compatible with this cell
0109   YieldType nLayers() const { return m_nLayers; }
0110   /// @brief access the number of unique hits compatible with this cell
0111   YieldType nHits() const { return m_nHits; }
0112   /// @brief access the set of layers compatible with this cell
0113   const std::unordered_set<unsigned>& layers() const { return m_layers; }
0114   /// @brief access the set of unique hits compatible with this cell
0115   const std::unordered_set<identifier_t>& hits() const { return m_hits; }
0116   /// @brief reset this cell, removing any existing content.
0117   void reset();
0118 
0119  private:
0120   /// data members
0121 
0122   YieldType m_nLayers =
0123       0;                  // (weighted) number of layers with hits on this cell
0124   YieldType m_nHits = 0;  // (weighted) number of unique hits on this cell
0125   std::unordered_set<unsigned> m_layers =
0126       {};  // set of layers with hits on this cell
0127   std::unordered_set<identifier_t> m_hits =
0128       {};  // set of unique hits on this cell
0129 };
0130 
0131 /// @brief Configuration - number of bins in each axis.
0132 /// The Hough plane is agnostic of how the bins map to
0133 /// coordinates, allowing to re-use a plane for several
0134 /// (sub) detectors of different dimensions if the bin number
0135 /// remains applicable
0136 struct HoughPlaneConfig {
0137   std::size_t nBinsX = 0;  // number of bins in the first coordinate
0138   std::size_t nBinsY = 0;  // number of bins in the second coordinate
0139 };
0140 
0141 /// @brief Representation of the hough plane - the histogram used
0142 /// for the hough transform with methods to fill and evaluate
0143 /// the histogram. Templated to a class used as identifier for the hits
0144 template <class identifier_t>
0145 class HoughPlane {
0146  public:
0147   /// @brief hough histogram representation as a 2D-indexable vector of hough cells
0148   using HoughHist = MultiIndexedVector2D<HoughCell<identifier_t>>;
0149 
0150   /// @brief instantiate the (empty) hough plane
0151   /// @param cfg: configuration
0152   HoughPlane(const HoughPlaneConfig& cfg);
0153 
0154   /// fill and reset methods to modify the grid content
0155 
0156   /// @brief add one measurement to the hough plane
0157   /// @tparam PointType: Type of the objects to use when adding measurements (e.g. experiment EDM object)
0158   /// @param measurement: The measurement to add
0159   /// @param axisRanges: Ranges of the hough axes, used to map the bin numbers to parameter values
0160   /// @param linePar: The function y(x) parametrising the hough space line for a given measurement
0161   /// @param widthPar: The function dy(x) parametrising the width of the y(x) curve
0162   ///                   for a given measurement
0163   /// @param identifier: The unique identifier for the given hit
0164   /// @param layer: A layer index for this hit
0165   /// @param weight: An optional weight to assign to this hit
0166   template <class PointType>
0167   void fill(const PointType& measurement, const HoughAxisRanges& axisRanges,
0168             LineParametrisation<PointType> linePar,
0169             LineParametrisation<PointType> widthPar,
0170             const identifier_t& identifier, unsigned layer = 0,
0171             YieldType weight = 1.0f);
0172   /// @brief resets the contents of the grid. Can be used to avoid reallocating the histogram
0173   /// when switching regions / (sub)detectors
0174   void reset();
0175 
0176   //// user-facing accessors
0177 
0178   /// @brief get the layers with hits in one cell of the histogram
0179   /// @param xBin: bin index in the first coordinate
0180   /// @param yBin: bin index in the second coordinate
0181   /// @return the set of layer indices that have hits for this cell
0182   const std::unordered_set<unsigned>& layers(std::size_t xBin,
0183                                              std::size_t yBin) const {
0184     return m_houghHist(xBin, yBin).layers();
0185   }
0186 
0187   /// @brief get the (weighted) number of layers  with hits in one cell of the histogram
0188   /// @param xBin: bin index in the first coordinate
0189   /// @param yBin: bin index in the second coordinate
0190   /// @return the (weighed) number of layers that have hits for this cell
0191   YieldType nLayers(std::size_t xBin, std::size_t yBin) const {
0192     return m_houghHist(xBin, yBin).nLayers();
0193   }
0194 
0195   /// @brief get the identifiers of all hits in one cell of the histogram
0196   /// @param xBin: bin index in the first coordinate
0197   /// @param yBin: bin index in the second coordinate
0198   /// @return the set of identifiers of the hits for this cell
0199   const std::unordered_set<identifier_t>& hitIds(std::size_t xBin,
0200                                                  std::size_t yBin) const {
0201     return m_houghHist(xBin, yBin).hits();
0202   }
0203   /// @brief get the (weighted) number of hits in one cell of the histogram
0204   /// @param xBin: bin index in the first coordinate
0205   /// @param yBin: bin index in the second coordinate
0206   /// @return the (weighted) number of hits for this cell
0207   YieldType nHits(std::size_t xBin, std::size_t yBin) const {
0208     return m_houghHist(xBin, yBin).nHits();
0209   }
0210 
0211   /// @brief get the number of bins on the first coordinate
0212   std::size_t nBinsX() const { return m_cfg.nBinsX; }
0213   /// @brief get the number of bins on the second coordinate
0214   std::size_t nBinsY() const { return m_cfg.nBinsY; }
0215 
0216   /// @brief get the maximum number of (weighted) hits seen in a single
0217   /// cell across the entire histrogram.
0218   YieldType maxHits() const { return m_maxHits; }
0219 
0220   /// @brief get the list of cells with non-zero content.
0221   /// Useful for peak-finders in sparse data
0222   /// to avoid looping over all cells
0223   const std::set<std::pair<std::size_t, std::size_t>>& getNonEmptyBins() const {
0224     return m_touchedBins;
0225   }
0226 
0227   /// @brief get the bin indices of the cell containing the largest number
0228   /// of (weighted) hits across the entire histogram
0229   std::pair<std::size_t, std::size_t> locMaxHits() const {
0230     return m_maxLocHits;
0231   }
0232 
0233   /// @brief get the maximum number of (weighted) layers with hits  seen
0234   /// in a single cell across the entire histrogram.
0235   YieldType maxLayers() const { return m_maxLayers; }
0236 
0237   /// @brief get the bin indices of the cell containing the largest number
0238   /// of (weighted) layers with hits across the entire histogram
0239   std::pair<std::size_t, std::size_t> locMaxLayers() const {
0240     return m_maxLocLayers;
0241   }
0242 
0243  private:
0244   /// @brief Helper method to fill a bin of the hough histogram.
0245   /// Updates the internal helper data structures (maximum tracker etc).
0246   /// @param binX: bin number along x
0247   /// @param binY: bin number along y
0248   /// @param identifier: hit identifier
0249   /// @param layer: layer index
0250   /// @param w: optional hit weight
0251   void fillBin(std::size_t binX, std::size_t binY,
0252                const identifier_t& identifier, unsigned layer, double w = 1.0f);
0253 
0254   YieldType m_maxHits = 0.0f;    // track the maximum number of hits seen
0255   YieldType m_maxLayers = 0.0f;  // track the maximum number of layers seen
0256   std::pair<std::size_t, std::size_t> m_maxLocHits = {
0257       0, 0};  // track the location of the maximum in hits
0258   std::pair<std::size_t, std::size_t> m_maxLocLayers = {
0259       0, 0};  // track the location of the maximum in layers
0260   std::set<std::pair<std::size_t, std::size_t>> m_touchedBins =
0261       {};                  // track the bins with non-trivial content
0262   HoughPlaneConfig m_cfg;  // the configuration object
0263   HoughHist m_houghHist;   // the histogram data object
0264 };
0265 
0266 /// example peak finders.
0267 namespace PeakFinders {
0268 /// configuration for the LayerGuidedCombinatoric peak finder
0269 struct LayerGuidedCombinatoricConfig {
0270   YieldType threshold = 3.0f;  // min number of layers to obtain a maximum
0271   int localMaxWindowSize = 0;  // Only create candidates from a local maximum
0272 };
0273 
0274 /// @brief Peak finder inspired by ATLAS ITk event filter developments.
0275 /// Builds peaks based on layer counts and allows for subsequent resolution
0276 /// of the combinatorics by building multiple candidates per peak if needed.
0277 /// In flat regions, peak positions are moved towards smaller values of the
0278 /// second and first coordinate.
0279 /// @tparam identifier_t: The identifier type to use. Should match the one used in the hough plane.
0280 template <class identifier_t>
0281 class LayerGuidedCombinatoric {
0282  public:
0283   /// @brief data class representing the found maxima.
0284   /// Here, we just have a list of cluster identifiers
0285   struct Maximum {
0286     std::unordered_set<identifier_t> hitIdentifiers =
0287         {};  // identifiers of contributing hits
0288   };
0289   /// @brief constructor
0290   /// @param cfg: Configuration object
0291   LayerGuidedCombinatoric(const LayerGuidedCombinatoricConfig& cfg);
0292 
0293   /// @brief main peak finder method.
0294   /// @param plane: Filled hough plane to search
0295   /// @return vector of found maxima
0296   std::vector<Maximum> findPeaks(const HoughPlane<identifier_t>& plane) const;
0297 
0298  private:
0299   /// @brief check if a given bin is a local maximum.
0300   /// @param plane: The filled hough plane
0301   /// @param xBin: x bin index
0302   /// @param yBin: y bin index
0303   /// @return true if a maximum, false otherwise
0304   bool passThreshold(const HoughPlane<identifier_t>& plane, std::size_t xBin,
0305                      std::size_t yBin) const;  // did we pass extensions?
0306 
0307   LayerGuidedCombinatoricConfig m_cfg;  // configuration data object
0308 };
0309 /// @brief Configuration for the IslandsAroundMax peak finder
0310 struct IslandsAroundMaxConfig {
0311   YieldType threshold =
0312       3.0f;  // min number of weigted hits required in a maximum
0313   YieldType fractionCutoff =
0314       0;  // Fraction of the global maximum at which to cut off maxima
0315   std::pair<CoordType, CoordType> minSpacingBetweenPeaks = {
0316       0.0f, 0.0f};  // minimum distance of a new peak from existing peaks in
0317                     // parameter space
0318 };
0319 /// @brief Peak finder inspired by ATLAS muon reconstruction.
0320 /// Looks for regions above a given fraction of the global maximum
0321 /// hit count and connects them into islands comprising adjacent bins
0322 /// above the threshold. Peak positions are averaged across cells in the island,
0323 /// weighted by hit counts
0324 /// @tparam identifier_t: The identifier type to use. Should match the one used in the hough plane.
0325 template <class identifier_t>
0326 class IslandsAroundMax {
0327  public:
0328   /// @brief data struct representing a local maximum.
0329   /// Comes with a position estimate and a list of hits within the island
0330   struct Maximum {
0331     CoordType x = 0;   // x value of the maximum
0332     CoordType y = 0;   // y value of the maximum
0333     CoordType wx = 0;  // x width of the maximum
0334     CoordType wy = 0;  // y width of the maximum
0335     std::unordered_set<identifier_t> hitIdentifiers =
0336         {};  // identifiers of contributing hits
0337   };
0338   /// @brief constructor.
0339   /// @param cfg: configuration object
0340   IslandsAroundMax(const IslandsAroundMaxConfig& cfg);
0341 
0342   /// @brief main peak finder method.
0343   /// @param plane: The filled hough plane to search
0344   /// @param ranges: The axis ranges used for mapping between parameter space and bins.
0345   /// @return List of the found maxima
0346   std::vector<Maximum> findPeaks(const HoughPlane<identifier_t>& plane,
0347                                  const HoughAxisRanges& ranges);
0348 
0349  private:
0350   /// @brief method to incrementally grow an island by adding adjacent cells
0351   /// Performs a breadth-first search for neighbours above threshold and adds
0352   /// them to candidate. Stops when no suitable neighbours are left.
0353   /// @param inMaximum: List of cells found in the island. Incrementally populated by calls to the method
0354   /// @param toExplore: List of neighbour cell candidates left to explore. Method will not do anything once this is empty
0355   /// @param threshold: the threshold to apply to check if a cell should be added to an island
0356   /// @param yieldMap: A map of the hit content of above-threshold cells. Used cells will be set to empty content to avoid re-use by subsequent calls
0357   void extendMaximum(
0358       std::vector<std::pair<std::size_t, std::size_t>>& inMaximum,
0359       std::vector<std::pair<std::size_t, std::size_t>>& toExplore,
0360       YieldType threshold,
0361       std::map<std::pair<std::size_t, std::size_t>, YieldType>& yieldMap);
0362 
0363   IslandsAroundMaxConfig m_cfg;  // configuration data object
0364 
0365   /// @brief array of steps to consider when exploring neighbouring cells.
0366   const std::array<std::pair<int, int>, 8> m_stepDirections{
0367       std::make_pair(-1, -1), std::make_pair(0, -1), std::make_pair(1, -1),
0368       std::make_pair(-1, 0),  std::make_pair(1, 0),  std::make_pair(-1, 1),
0369       std::make_pair(0, 1),   std::make_pair(1, 1)};
0370 };
0371 }  // namespace PeakFinders
0372 }  // namespace Acts::HoughTransformUtils
0373 
0374 #include "HoughTransformUtils.ipp"