|
|
|||
File indexing completed on 2026-01-10 09:16:04
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 "Acts/Definitions/Algebra.hpp" 0012 #include "Acts/Material/IVolumeMaterial.hpp" 0013 #include "Acts/Material/Material.hpp" 0014 #include "Acts/Utilities/BinUtility.hpp" 0015 #include "Acts/Utilities/Interpolation.hpp" 0016 0017 #include <functional> 0018 #include <optional> 0019 0020 namespace Acts { 0021 0022 /// @addtogroup material 0023 /// @{ 0024 0025 /// @brief Struct for mapping global 3D positions to material values 0026 /// 0027 /// Global 3D positions are transformed into a @c DIM_POS Dimensional vector 0028 /// which is used to look up the material classification value in the 0029 /// underlying material map. 0030 template <typename G> 0031 struct MaterialMapLookup { 0032 public: 0033 /// Type alias for material grid 0034 using Grid_t = G; 0035 /// Dimensionality of the position space for material interpolation 0036 static constexpr std::size_t DIM_POS = Grid_t::DIM; 0037 0038 /// @brief Struct representing smallest grid unit in material grid 0039 /// 0040 /// This type encapsulate all required information to perform linear 0041 /// interpolation of material classification values within a 3D volume. 0042 struct MaterialCell { 0043 /// Number of corner points defining the confining hyper-box 0044 static constexpr unsigned int N = 1 << DIM_POS; 0045 0046 /// @brief Default constructor 0047 /// 0048 /// @param [in] transformPos Mapping of global 3D coordinates onto grid 0049 /// space 0050 /// @param [in] lowerLeft Generalized lower-left corner of hyper box 0051 /// (containing the minima of the hyper box along 0052 /// each Dimension) 0053 /// @param [in] upperRight Generalized upper-right corner of hyper box 0054 /// (containing the maxima of the hyper box along 0055 /// each Dimension) 0056 /// @param [in] materialValues Material classification values at the hyper 0057 /// box corners sorted in the canonical order defined in Acts::interpolate 0058 MaterialCell( 0059 std::function<ActsVector<DIM_POS>(const Vector3&)> transformPos, 0060 std::array<double, DIM_POS> lowerLeft, 0061 std::array<double, DIM_POS> upperRight, 0062 std::array<Material::ParametersVector, N> materialValues) 0063 : m_transformPos(std::move(transformPos)), 0064 m_lowerLeft(std::move(lowerLeft)), 0065 m_upperRight(std::move(upperRight)), 0066 m_materialValues(std::move(materialValues)) {} 0067 0068 /// @brief Retrieve material at given position 0069 /// 0070 /// @param [in] position Global 3D position 0071 /// @return Material at the given position 0072 /// 0073 /// @pre The given @c position must lie within the current cell. 0074 Material getMaterial(const Vector3& position) const { 0075 // defined in Interpolation.hpp 0076 return Material(interpolate(m_transformPos(position), m_lowerLeft, 0077 m_upperRight, m_materialValues)); 0078 } 0079 0080 /// @brief Check whether given 3D position is inside this cell 0081 /// 0082 /// @param [in] position Global 3D position 0083 /// @return @c true if position is inside the current cell, 0084 /// otherwise @c false 0085 bool isInside(const Vector3& position) const { 0086 const auto& gridCoordinates = m_transformPos(position); 0087 for (unsigned int i = 0; i < DIM_POS; ++i) { 0088 if (gridCoordinates[i] < m_lowerLeft.at(i) || 0089 gridCoordinates[i] >= m_upperRight.at(i)) { 0090 return false; 0091 } 0092 } 0093 return true; 0094 } 0095 0096 private: 0097 /// Geometric transformation applied to global 3D positions 0098 std::function<ActsVector<DIM_POS>(const Vector3&)> m_transformPos; 0099 0100 /// Generalized lower-left corner of the confining hyper-box 0101 std::array<double, DIM_POS> m_lowerLeft; 0102 0103 /// Generalized upper-right corner of the confining hyper-box 0104 std::array<double, DIM_POS> m_upperRight; 0105 0106 /// @brief Material component vectors at the hyper-box corners 0107 /// 0108 /// @note These values must be order according to the prescription detailed 0109 /// in Acts::interpolate. 0110 std::array<Material::ParametersVector, N> m_materialValues; 0111 }; 0112 0113 /// @brief Default constructor 0114 /// 0115 /// @param [in] transformPos Mapping of global 3D coordinates (cartesian) 0116 /// onto grid space 0117 /// @param [in] grid Grid storing material classification values 0118 MaterialMapLookup( 0119 std::function<ActsVector<DIM_POS>(const Vector3&)> transformPos, 0120 Grid_t grid) 0121 : m_transformPos(std::move(transformPos)), m_grid(std::move(grid)) {} 0122 0123 /// @brief Retrieve binned material at given position 0124 /// 0125 /// @param [in] position Global 3D position 0126 /// @return Material at the given position 0127 /// 0128 /// @pre The given @c position must lie within the range of the underlying 0129 /// map. 0130 Material material(const Vector3& position) const { 0131 return Material(m_grid.atLocalBins( 0132 m_grid.localBinsFromLowerLeftEdge(m_transformPos(position)))); 0133 } 0134 0135 /// @brief Retrieve interpolated material at given position 0136 /// 0137 /// @param [in] position Global 3D position 0138 /// @return Material at the given position 0139 /// 0140 /// @pre The given @c position must lie within the range of the underlying 0141 /// map. 0142 Material getMaterial(const Vector3& position) const { 0143 return Material(m_grid.interpolate(m_transformPos(position))); 0144 } 0145 0146 /// @brief Retrieve material cell for given position 0147 /// 0148 /// @param [in] position Global 3D position 0149 /// @return material cell containing the given global position 0150 /// 0151 /// @pre The given @c position must lie within the range of the underlying 0152 /// map. 0153 MaterialCell getMaterialCell(const Vector3& position) const { 0154 const auto& gridPosition = m_transformPos(position); 0155 std::size_t bin = m_grid.globalBinFromPosition(gridPosition); 0156 const auto& indices = m_grid.localBinsFromPosition(bin); 0157 const auto& lowerLeft = m_grid.lowerLeftBinEdge(indices); 0158 const auto& upperRight = m_grid.upperRightBinEdge(indices); 0159 0160 // Loop through all corner points 0161 constexpr std::size_t nCorners = 1 << DIM_POS; 0162 std::array<Material::ParametersVector, nCorners> neighbors{}; 0163 const auto& cornerIndices = m_grid.closestPointsIndices(gridPosition); 0164 0165 std::size_t i = 0; 0166 for (std::size_t index : cornerIndices) { 0167 neighbors.at(i++) = m_grid.at(index); 0168 } 0169 0170 return MaterialCell(m_transformPos, lowerLeft, upperRight, 0171 std::move(neighbors)); 0172 } 0173 0174 /// @brief Get the number of bins for all axes of the map 0175 /// 0176 /// @return Vector returning number of bins for all map axes 0177 std::vector<std::size_t> getNBins() const { 0178 auto nBinsArray = m_grid.numLocalBins(); 0179 return std::vector<std::size_t>(nBinsArray.begin(), nBinsArray.end()); 0180 } 0181 0182 /// @brief Get the minimum value of all axes of the map 0183 /// 0184 /// @return Vector returning the minima of all map axes 0185 std::vector<double> getMin() const { 0186 auto minArray = m_grid.minPosition(); 0187 return std::vector<double>(minArray.begin(), minArray.end()); 0188 } 0189 0190 /// @brief Get the maximum value of all axes of the map 0191 /// 0192 /// @return Vector returning the maxima of all map axes 0193 std::vector<double> getMax() const { 0194 auto maxArray = m_grid.maxPosition(); 0195 return std::vector<double>(maxArray.begin(), maxArray.end()); 0196 } 0197 0198 /// @brief Check whether given 3D position is inside look-up domain 0199 /// 0200 /// @param [in] position Global 3D position 0201 /// @return @c true if position is inside the defined look-up grid, 0202 /// otherwise @c false 0203 bool isInside(const Vector3& position) const { 0204 return m_grid.isInside(m_transformPos(position)); 0205 } 0206 0207 /// @brief Get a const reference on the underlying grid structure 0208 /// 0209 /// @return Grid reference 0210 const Grid_t& getGrid() const { return m_grid; } 0211 0212 private: 0213 /// Geometric transformation applied to global 3D positions 0214 std::function<ActsVector<DIM_POS>(const Vector3&)> m_transformPos; 0215 /// Grid storing material values 0216 Grid_t m_grid; 0217 }; 0218 0219 /// @brief Interpolate material classification values from material values on a 0220 /// given grid 0221 /// 0222 /// This class implements a material service which is initialized by a 0223 /// material map defined by: 0224 /// - a list of material values on a regular grid in some n-Dimensional space, 0225 /// - a transformation of global 3D coordinates onto this n-Dimensional 0226 /// space. 0227 /// - a transformation of local n-Dimensional material coordinates into 0228 /// global (cartesian) 3D coordinates 0229 /// 0230 /// The material value for a given global position is then determined by: 0231 /// - mapping the position onto the grid, 0232 /// - looking up the material classification values on the closest grid points, 0233 /// - doing a linear interpolation of these values. 0234 /// @warning Each classification number of the material is interpolated 0235 /// independently and thus does not consider any correlations that exists 0236 /// between these values. This might work out since the used material is already 0237 /// a mean of the materials in a certain bin and can therewith be treated as a 0238 /// collection of numbers. 0239 /// @tparam G Type of the grid 0240 template <typename Mapper_t> 0241 class InterpolatedMaterialMap : public IVolumeMaterial { 0242 public: 0243 /// @brief Temporary storage of a certain cell to improve material access 0244 struct Cache { 0245 /// Stored material cell 0246 std::optional<typename Mapper_t::MaterialCell> matCell; 0247 /// Boolean statement if the cell is initialized 0248 bool initialized = false; 0249 }; 0250 0251 /// @brief Create interpolated map 0252 /// 0253 /// @param [in] mapper Material map 0254 explicit InterpolatedMaterialMap(Mapper_t&& mapper) 0255 : m_mapper(std::move(mapper)) {} 0256 0257 /// @brief Create interpolated map 0258 /// 0259 /// @param [in] mapper Material map 0260 /// @param [in] bu @c BinUtility for build from 0261 InterpolatedMaterialMap(Mapper_t&& mapper, BinUtility bu) 0262 : m_mapper(std::move(mapper)), m_binUtility(std::move(bu)) {} 0263 0264 /// @brief Retrieve the binned material 0265 /// 0266 /// @param [in] position Global 3D position 0267 /// 0268 /// @return Material at given position 0269 const Material material(const Vector3& position) const override { 0270 return m_mapper.material(position); 0271 } 0272 0273 /// @brief Retrieve the interpolated material 0274 /// 0275 /// @param [in] position Global 3D position 0276 /// 0277 /// @return material at given position 0278 Material getMaterial(const Vector3& position) const { 0279 return m_mapper.getMaterial(position); 0280 } 0281 0282 /// @brief Retrieve material 0283 /// 0284 /// @param [in] position Global 3D position 0285 /// @param [in,out] cache Cache object. Contains material cell used for 0286 /// interpolation 0287 /// 0288 /// @return material at given position 0289 Material getMaterial(const Vector3& position, Cache& cache) const { 0290 if (!cache.initialized || !(*cache.matCell).isInside(position)) { 0291 cache.matCell = getMaterialCell(position); 0292 cache.initialized = true; 0293 } 0294 return (*cache.matCell).getMaterial(position); 0295 } 0296 0297 /// @brief Retrieve material value & its "gradient" 0298 /// 0299 /// @param [in] position Global 3D position 0300 /// @return Material 0301 /// 0302 /// @note Currently the derivative is not calculated 0303 /// @todo return derivative 0304 Material getMaterialGradient(const Vector3& position, 0305 ActsMatrix<5, 5>& /*derivative*/) const { 0306 return m_mapper.getMaterial(position); 0307 } 0308 0309 /// @brief Retrieve material value & its "gradient" 0310 /// 0311 /// @param [in] position Global 3D position 0312 /// @return Material 0313 /// 0314 /// @note Currently the derivative is not calculated 0315 /// @note Cache is not used currently 0316 /// @todo return derivative 0317 Material getMaterialGradient(const Vector3& position, 0318 ActsMatrix<5, 5>& /*derivative*/, 0319 Cache& /*cache*/) const { 0320 return m_mapper.getMaterial(position); 0321 } 0322 0323 /// @brief Convenience method to access underlying material mapper 0324 /// 0325 /// @return The material mapper 0326 const Mapper_t& getMapper() const { return m_mapper; } 0327 0328 /// @brief Check whether given 3D position is inside look-up domain 0329 /// 0330 /// @param [in] position Global 3D position 0331 /// @return @c true if position is inside the defined map, otherwise @c false 0332 bool isInside(const Vector3& position) const { 0333 return m_mapper.isInside(position); 0334 } 0335 0336 /// Return the BinUtility 0337 /// @return Const reference to the bin utility for the material map 0338 const BinUtility& binUtility() const { return m_binUtility; } 0339 0340 /// Output Method for std::ostream 0341 /// 0342 /// @param sl The outoput stream 0343 /// @return Reference to the output stream for method chaining 0344 std::ostream& toStream(std::ostream& sl) const override { 0345 sl << "Acts::InterpolatedMaterialMap : " << std::endl; 0346 sl << " - Number of Material bins [0,1] : " << m_binUtility.max(0) + 1 0347 << " / " << m_binUtility.max(1) + 1 << std::endl; 0348 sl << " - Parse full update material : " << std::endl; 0349 return sl; 0350 } 0351 0352 private: 0353 /// @brief Retrieve cell for given position 0354 /// 0355 /// @param [in] position Global 3D position 0356 /// @return Material cell containing the given global position 0357 /// 0358 /// @pre The given @c position must lie within the range of the underlying 0359 /// map. 0360 typename Mapper_t::MaterialCell getMaterialCell( 0361 const Vector3& position) const { 0362 return m_mapper.getMaterialCell(position); 0363 } 0364 0365 /// @brief object for global coordinate transformation and interpolation 0366 /// 0367 /// This object performs the mapping of the global 3D coordinates onto the 0368 /// material grid and the interpolation of the material component values on 0369 /// close-by grid points. 0370 Mapper_t m_mapper; 0371 0372 BinUtility m_binUtility{}; 0373 }; 0374 0375 /// @} 0376 } // namespace Acts
| [ Source navigation ] | [ Diff markup ] | [ Identifier search ] | [ general search ] |
|
This page was automatically generated by the 2.3.7 LXR engine. The LXR team |
|