detray/utils/axis_rotation.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 // Project include(s).
0012 #include "detray/definitions/algebra.hpp"
0013 #include "detray/definitions/detail/qualifiers.hpp"
0014 #include "detray/definitions/math.hpp"
0015
0016 namespace detray {
0017
0018 /// @brief Helper struct to rotate a vector around a given axis and angle
0019 /// counterclockwisely
0020 template <concepts::algebra algebra_t>
0021 struct axis_rotation {
0022  public:
0023   using algebra_type = algebra_t;
0024   using scalar_type = dscalar<algebra_t>;
0025   using vector3_type = dvector3D<algebra_t>;
0026   template <std::size_t ROWS, std::size_t COLS>
0027   using matrix_type = dmatrix<algebra_t, ROWS, COLS>;
0028
0029   /// @brief Constructor for axis rotation
0030   ///
0031   /// @param axis rotation axis
0032   /// @param theta rotation angle
0033   DETRAY_HOST_DEVICE
0034   axis_rotation(const vector3_type& axis, const scalar_type theta) {
0035     // normalize the axis
0036     const vector3_type U = vector::normalize(axis);
0037
0038     scalar_type cos_theta{math::cos(theta)};
0039
0040     auto I = matrix::identity<matrix_type<3, 3>>();
0041     matrix_type<3, 3> axis_cross = matrix::cross_matrix(U);
0042     matrix_type<3, 3> axis_outer = matrix::outer_product(U, U);
0043
0044     R = cos_theta * I + math::sin(theta) * axis_cross +
0045         (1.f - cos_theta) * axis_outer;
0046   }
0047
0048   /// @param v vector to be rotated
0049   /// @returns Get the counterclockwisely-rotated vector
0050   template <typename vector3_t>
0051     requires(concepts::vector3D<vector3_t> ||
0052              concepts::column_matrix3D<vector3_t>)
0053   DETRAY_HOST_DEVICE vector3_t operator()(const vector3_t& v) const {
0054     return R * v;
0055   }
0056
0057  private:
0058   /// Rotation matrix
0059   matrix_type<3, 3> R{matrix::identity<matrix_type<3, 3>>()};
0060 };
0061
0062 /// @brief Helper struct to perform an euler rotation for a given vector
0063 /// All rotation operations are counterclockwise
0064 template <concepts::algebra algebra_t>
0065 struct euler_rotation {
0066  public:
0067   using algebra_type = algebra_t;
0068   using scalar_type = dscalar<algebra_t>;
0069   using vector3_type = dvector3D<algebra_t>;
0070
0071   // Following the z-x-z convention
0072   vector3_type x{1.0f, 0.f, 0.f};
0073   vector3_type z{0.f, 0.f, 1.f};
0074   scalar_type alpha{0.f};
0075   scalar_type beta{0.f};
0076   scalar_type gamma{0.f};
0077
0078   /// @returns Get the new x and z axis
0079   DETRAY_HOST_DEVICE std::pair<vector3_type, vector3_type> operator()() const {
0080     // alpha around z axis
0081     axis_rotation<algebra_t> axis_rot_alpha(z, alpha);
0082     auto new_x = axis_rot_alpha(x);
0083
0084     // beta around x' axis
0085     axis_rotation<algebra_t> axis_rot_beta(new_x, beta);
0086     const auto new_z = axis_rot_beta(z);
0087
0088     axis_rotation<algebra_t> axis_rot_gamma(new_z, gamma);
0089     new_x = axis_rot_gamma(new_x);
0090
0091     return std::make_pair(new_x, new_z);
0092   }
0093
0094   /// @param v vector to be rotated
0095   /// @returns Get the rotated vector
0096   DETRAY_HOST_DEVICE vector3_type operator()(const vector3_type& v0) const {
0097     // alpha around z axis
0098     axis_rotation<algebra_t> axis_rot_alpha(z, alpha);
0099     const auto new_x = axis_rot_alpha(x);
0100     const auto v1 = axis_rot_alpha(v0);
0101
0102     // beta around x' axis
0103     axis_rotation<algebra_t> axis_rot_beta(new_x, beta);
0104     const auto new_z = axis_rot_beta(z);
0105     const auto v2 = axis_rot_beta(v1);
0106
0107     // gamma around z' axis
0108     axis_rotation<algebra_t> axis_rot_gamma(new_z, gamma);
0109     const auto v3 = axis_rot_gamma(v2);
0110
0111     return v3;
0112   }
0113 };
0114
0115 }  // namespace detray