|
||||
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"
[ Source navigation ] | [ Diff markup ] | [ Identifier search ] | [ general search ] |
This page was automatically generated by the 2.3.7 LXR engine. The LXR team |