Acts/Utilities/Interpolation.hpp

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/Utilities/detail/interpolation_impl.hpp"
0013
0014 #include <array>
0015 #include <type_traits>
0016
0017 namespace Acts {
0018
0019 /// @brief performs linear interpolation inside a hyper box
0020 ///
0021 /// @tparam T      type of values to be interpolated
0022 /// @tparam N      number of hyper box corners
0023 /// @tparam Point1 type specifying geometric positions
0024 /// @tparam Point2 type specifying geometric positions
0025 /// @tparam Point3 type specifying geometric positions
0026 ///
0027 /// @param [in] position    position to which to interpolate
0028 /// @param [in] lowerCorner generalized lower-left corner of hyper box
0029 ///                         (containing the minima of the hyper box along each
0030 ///                         dimension)
0031 /// @param [in] upperCorner generalized upper-right corner of hyper box
0032 ///                         (containing the maxima of the hyper box along each
0033 ///                         dimension)
0034 /// @param [in] values      field values at the hyper box corners sorted in the
0035 ///                         canonical order defined below.
0036 ///
0037 /// @return interpolated value at given position
0038 ///
0039 /// @pre @c position must describe a position inside the given hyper box, that
0040 ///      is \f$\text{lowerCorner}[i] \le \text{position}[i] \le
0041 ///      \text{upperCorner}[i] \quad \forall i=0, \dots, d-1\f$.
0042 ///
0043 /// @note
0044 /// - Given @c U and @c V of value type @c T as well as two @c double @c a and
0045 /// @c b, then the following must be a valid expression <tt>a * U + b * V</tt>
0046 /// yielding an object which is (implicitly) convertible to @c T.
0047 /// - All @c Point types must represent d-dimensional positions and support
0048 /// coordinate access using @c operator[] which should return a @c double (or a
0049 /// value which is implicitly convertible). Coordinate indices must start at 0.
0050 /// - @c N is the number of hyper box corners which is \f$2^d\f$ where \f$d\f$
0051 /// is the dimensionality of the hyper box. The dimensionality must be
0052 /// consistent with the provided @c Point types.
0053 /// - Definition of the canonical order for sorting the field values: The hyper
0054 /// box corners are numbered according to the following scheme. Each corner is
0055 /// defined by the set of lower/upper boundary limits in each dimension @c i.
0056 /// This can be represented by a binary code (from left to right) where a @c 0
0057 /// stands for a lower bound along this axis and a @c 1 stand for the upper
0058 /// bound along this axis. The left most bit corresponds to the first dimension
0059 /// and the bits to the left correspond to the 2nd, 3rd... dimension. The binary
0060 /// code can be interpreted as integer which gives the number of the
0061 /// corresponding hyper box corner. The field values are ordered according to
0062 /// ascending hyper box corner numbers.<br />
0063 /// As an example assume we have a 3D box with @c lowerCorner = (1,2,3) and @c
0064 /// upperCorner = (4,5,6). The eight corners with their bit patterns and corner
0065 /// numbers are:
0066 ///    - (1,2,3): 000 = 0
0067 ///    - (1,2,6): 001 = 1
0068 ///    - (1,5,3): 010 = 2
0069 ///    - (1,5,6): 011 = 3
0070 ///    - (4,2,3): 100 = 4
0071 ///    - (4,2,6): 101 = 5
0072 ///    - (4,5,3): 110 = 6
0073 ///    - (4,5,6): 111 = 7
0074 template <typename T, std::size_t N, class Point1, class Point2 = Point1,
0075           class Point3 = Point2>
0076 inline T interpolate(const Point1& position, const Point2& lowerCorner,
0077                      const Point3& upperCorner, const std::array<T, N>& values)
0078   requires(Concepts::interpolatable<T, Point1, Point2, Point3>)
0079 {
0080   return detail::interpolate_impl<T, Point1, Point2, Point3,
0081                                   detail::get_dimension<N>::value - 1,
0082                                   N>::run(position, lowerCorner, upperCorner,
0083                                           values);
0084 }
0085
0086 }  // namespace Acts