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