Utilities/detail/interpolation_impl.hpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2017-2018 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 #pragma once
0010
0011 #include <array>
0012
0013 namespace Acts::detail {
0014
0015 /// @brief check types for requirements needed by interpolation
0016 ///
0017 /// @tparam Point1 type for specifying geometric positions
0018 /// @tparam Point2 type for specifying geometric positions
0019 /// @tparam Point3 type for specifying geometric positions
0020 /// @tparam Value  type of values to be interpolated
0021 ///
0022 /// This helper struct provides compile-time information whether the provided
0023 /// @c Point and @c Value types can be used in the Acts::interpolate function.
0024 ///
0025 /// The following boolean variable
0026 /// @code{.cpp}
0027 /// Acts::detail::can_interpolate<Point1,Point2,Point3,Value>::value
0028 /// @endcode
0029 ///
0030 /// is @c true if all @c Point types and @c Value fulfill the type
0031 /// requirements for being used in the interpolation function, otherwise it is
0032 /// @c false. This expression can be employed in @c std::enable_if_t to use
0033 /// SFINAE patterns to enable/disable (member) functions.
0034 template <typename Point1, typename Point2, typename Point3, typename Value>
0035 struct can_interpolate {
0036   template <typename C>
0037   static auto value_type_test(C* c)
0038       -> decltype(static_cast<C>(std::declval<double>() * std::declval<C>() +
0039                                  std::declval<double>() * std::declval<C>()),
0040                   std::true_type());
0041   template <typename C>
0042   static std::false_type value_type_test(...);
0043
0044   template <typename C>
0045   static auto point_type_test(C* c)
0046       -> decltype(static_cast<double>(std::declval<C>()[0]), std::true_type());
0047   template <typename C>
0048   static std::false_type point_type_test(...);
0049
0050   static const bool value =
0051       std::is_same<std::true_type,
0052                    decltype(value_type_test<Value>(nullptr))>::value &&
0053       std::is_same<std::true_type,
0054                    decltype(point_type_test<Point1>(nullptr))>::value &&
0055       std::is_same<std::true_type,
0056                    decltype(point_type_test<Point2>(nullptr))>::value &&
0057       std::is_same<std::true_type,
0058                    decltype(point_type_test<Point3>(nullptr))>::value;
0059 };
0060
0061 /// @brief determine number of dimension from power of 2
0062 ///
0063 /// @tparam N power of 2
0064 template <std::size_t N>
0065 struct get_dimension {
0066   /// exponent @c d such that \f$2^d = N \f$
0067   static constexpr std::size_t value = get_dimension<(N >> 1)>::value + 1u;
0068 };
0069
0070 /// @cond
0071 template <>
0072 struct get_dimension<2u> {
0073   static constexpr std::size_t value = 1u;
0074 };
0075 /// @endcond
0076
0077 /// @brief helper struct for performing multi-dimensional linear interpolation
0078 ///
0079 /// @tparam T      type of values to be interpolated
0080 /// @tparam Point1 type specifying geometric positions
0081 /// @tparam Point2 type specifying geometric positions
0082 /// @tparam Point3 type specifying geometric positions
0083 /// @tparam D      current dimension on which to perform reduction
0084 /// @tparam N      number of hyper box corners
0085 ///
0086 /// @note
0087 /// - Given @c U and @c V of value type @c T as well as two @c double @c a and
0088 /// @c b, then the following must be a valid expression <tt>a * U + b * V</tt>
0089 /// yielding an object which is (implicitly) convertible to @c T.
0090 /// - The @c Point types must represent d-dimensional positions and support
0091 /// coordinate access using @c operator[]. Coordinate indices must start at 0.
0092 /// - @c N is the number of hyper box corners which is \f$2^d\f$ where \f$d\f$
0093 /// is the dimensionality of the hyper box. The dimensionality must be
0094 /// consistent with the provided @c Point types.
0095 template <typename T, class Point1, class Point2, class Point3, std::size_t D,
0096           std::size_t N>
0097 struct interpolate_impl;
0098
0099 /// @cond
0100 // recursive implementation of linear interpolation in multiple dimensions
0101 template <typename T, class Point1, class Point2, class Point3, std::size_t D,
0102           std::size_t N>
0103 struct interpolate_impl {
0104   static T run(const Point1& pos, const Point2& lowerLeft,
0105                const Point3& upperRight, const std::array<T, N>& fields) {
0106     // get distance to lower boundary relative to total bin width
0107     const double f = (pos[D] - lowerLeft[D]) / (upperRight[D] - lowerLeft[D]);
0108
0109     std::array<T, (N >> 1)> newFields{};
0110     for (std::size_t i = 0; i < N / 2; ++i) {
0111       newFields.at(i) = (1 - f) * fields.at(2 * i) + f * fields.at(2 * i + 1);
0112     }
0113
0114     return interpolate_impl<T, Point1, Point2, Point3, D - 1, (N >> 1)>::run(
0115         pos, lowerLeft, upperRight, newFields);
0116   }
0117 };
0118
0119 // simple linear interpolation in 1D
0120 template <typename T, class Point1, class Point2, class Point3, std::size_t D>
0121 struct interpolate_impl<T, Point1, Point2, Point3, D, 2u> {
0122   static T run(const Point1& pos, const Point2& lowerLeft,
0123                const Point3& upperRight, const std::array<T, 2u>& fields) {
0124     // get distance to lower boundary relative to total bin width
0125     const double f = (pos[D] - lowerLeft[D]) / (upperRight[D] - lowerLeft[D]);
0126
0127     return (1 - f) * fields.at(0) + f * fields.at(1);
0128   }
0129 };
0130 /// @endcond
0131 }  // namespace Acts::detail