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