interpolation_impl.hpp

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file interpolation_impl.hpp

// This file is part of the Acts project.
//
// Copyright (C) 2017-2018 CERN for the benefit of the Acts project
//
// This Source Code Form is subject to the terms of the Mozilla Public
// License, v. 2.0. If a copy of the MPL was not distributed with this
// file, You can obtain one at http://mozilla.org/MPL/2.0/.

#pragma once

#include <array>

namespace Acts {

namespace detail {

template <typename Point1, typename Point2, typename Point3, typename Value>
struct can_interpolate {
  template <typename C>
  static auto value_type_test(C* c)
      -> decltype(C(std::declval<double>() * std::declval<C>() +
                    std::declval<double>() * std::declval<C>()),
                  std::true_type());
  template <typename C>
  static std::false_type value_type_test(...);

  template <typename C>
  static auto point_type_test(C* c)
      -> decltype(double(std::declval<C>()[0]), std::true_type());
  template <typename C>
  static std::false_type point_type_test(...);

  static const bool value =
      std::is_same<std::true_type,
                   decltype(value_type_test<Value>(nullptr))>::value and
      std::is_same<std::true_type,
                   decltype(point_type_test<Point1>(nullptr))>::value and
      std::is_same<std::true_type,
                   decltype(point_type_test<Point2>(nullptr))>::value and
      std::is_same<std::true_type,
                   decltype(point_type_test<Point3>(nullptr))>::value;
};

template <size_t N>
struct get_dimension {
  static constexpr size_t value = get_dimension<(N >> 1)>::value + 1u;
};

template <>
struct get_dimension<2u> {
  static constexpr size_t value = 1u;
};

template <typename T, class Point1, class Point2, class Point3, size_t D,
          size_t N>
struct interpolate_impl;

// recursive implementation of linear interpolation in multiple dimensions
template <typename T, class Point1, class Point2, class Point3, size_t D,
          size_t N>
struct interpolate_impl {
  static T run(const Point1& pos, const Point2& lowerLeft,
               const Point3& upperRight, const std::array<T, N>& fields) {
    // get distance to lower boundary relative to total bin width
    const double f = (pos[D] - lowerLeft[D]) / (upperRight[D] - lowerLeft[D]);

    std::array<T, (N >> 1)> newFields;
    for (size_t i = 0; i < N / 2; ++i) {
      newFields.at(i) = (1 - f) * fields.at(2 * i) + f * fields.at(2 * i + 1);
    }

    return interpolate_impl<T, Point1, Point2, Point3, D - 1, (N >> 1)>::run(
        pos, lowerLeft, upperRight, newFields);
  }
};

// simple linear interpolation in 1D
template <typename T, class Point1, class Point2, class Point3, size_t D>
struct interpolate_impl<T, Point1, Point2, Point3, D, 2u> {
  static T run(const Point1& pos, const Point2& lowerLeft,
               const Point3& upperRight, const std::array<T, 2u>& fields) {
    // get distance to lower boundary relative to total bin width
    const double f = (pos[D] - lowerLeft[D]) / (upperRight[D] - lowerLeft[D]);

    return (1 - f) * fields.at(0) + f * fields.at(1);
  }
};
}  // namespace detail

}  // namespace Acts