src/Surfaces/TrapezoidBounds.cpp

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 #include "Acts/Surfaces/TrapezoidBounds.hpp"
0010
0011 #include "Acts/Surfaces/ConvexPolygonBounds.hpp"
0012
0013 #include <iomanip>
0014 #include <iostream>
0015
0016 namespace Acts {
0017
0018 TrapezoidBounds::TrapezoidBounds(double halfXnegY, double halfXposY,
0019                                  double halfY, double rotAngle) noexcept(false)
0020     : m_values({halfXnegY, halfXposY, halfY, rotAngle}),
0021       m_boundingBox(std::max(halfXnegY, halfXposY), halfY) {
0022   rotateBoundingBox();
0023   checkConsistency();
0024 }
0025
0026 TrapezoidBounds::TrapezoidBounds(
0027     const std::array<double, eSize>& values) noexcept(false)
0028     : m_values(values),
0029       m_boundingBox(
0030           std::max(values[eHalfLengthXnegY], values[eHalfLengthXposY]),
0031           values[eHalfLengthY]) {
0032   rotateBoundingBox();
0033   checkConsistency();
0034 }
0035
0036 std::vector<double> TrapezoidBounds::values() const {
0037   return {m_values.begin(), m_values.end()};
0038 }
0039
0040 bool TrapezoidBounds::inside(const Vector2& lposition) const {
0041   const double hlXnY = get(TrapezoidBounds::eHalfLengthXnegY);
0042   const double hlXpY = get(TrapezoidBounds::eHalfLengthXposY);
0043   const double hlY = get(TrapezoidBounds::eHalfLengthY);
0044   const double rotAngle = get(TrapezoidBounds::eRotationAngle);
0045
0046   const Vector2 extPosition = Eigen::Rotation2Dd(rotAngle) * lposition;
0047   const double x = extPosition[0];
0048   const double y = extPosition[1];
0049
0050   if (std::abs(y) - hlY > 0) {
0051     // outside y range
0052     return false;
0053   }
0054
0055   if (std::abs(x) - std::max(hlXnY, hlXpY) > 0) {
0056     // outside x range
0057     return false;
0058   }
0059
0060   if (std::abs(x) - std::min(hlXnY, hlXpY) <= 0) {
0061     // inside x range
0062     return true;
0063   }
0064
0065   // at this stage, the point can only be in the triangles
0066   // run slow-ish polygon check
0067   std::array<Vector2, 4> vertices{
0068       {{-hlXnY, -hlY}, {hlXnY, -hlY}, {hlXpY, hlY}, {-hlXpY, hlY}}};
0069   return detail::VerticesHelper::isInsidePolygon(extPosition, vertices);
0070 }
0071
0072 Vector2 TrapezoidBounds::closestPoint(const Vector2& lposition,
0073                                       const SquareMatrix2& metric) const {
0074   const double hlXnY = get(TrapezoidBounds::eHalfLengthXnegY);
0075   const double hlXpY = get(TrapezoidBounds::eHalfLengthXposY);
0076   const double hlY = get(TrapezoidBounds::eHalfLengthY);
0077   const double rotAngle = get(TrapezoidBounds::eRotationAngle);
0078
0079   const Vector2 extPosition = Eigen::Rotation2Dd(rotAngle) * lposition;
0080
0081   std::array<Vector2, 4> vertices{
0082       {{-hlXnY, -hlY}, {hlXnY, -hlY}, {hlXpY, hlY}, {-hlXpY, hlY}}};
0083
0084   Vector2 extClosest = detail::VerticesHelper::computeClosestPointOnPolygon(
0085       extPosition, vertices, metric);
0086
0087   return Eigen::Rotation2Dd(-rotAngle) * extClosest;
0088 }
0089
0090 std::vector<Vector2> TrapezoidBounds::vertices(
0091     unsigned int /*ignoredSegments*/) const {
0092   const double hlXnY = get(TrapezoidBounds::eHalfLengthXnegY);
0093   const double hlXpY = get(TrapezoidBounds::eHalfLengthXposY);
0094   const double hlY = get(TrapezoidBounds::eHalfLengthY);
0095   const double rotAngle = get(TrapezoidBounds::eRotationAngle);
0096
0097   std::vector<Vector2> vertices = {
0098       {-hlXnY, -hlY}, {hlXnY, -hlY}, {hlXpY, hlY}, {-hlXpY, hlY}};
0099   for (auto& v : vertices) {
0100     v = Eigen::Rotation2Dd(-rotAngle) * v;
0101   }
0102   return vertices;
0103 }
0104
0105 const RectangleBounds& TrapezoidBounds::boundingBox() const {
0106   return m_boundingBox;
0107 }
0108
0109 Vector2 TrapezoidBounds::center() const {
0110   return Vector2::Zero();
0111 }
0112
0113 std::ostream& TrapezoidBounds::toStream(std::ostream& sl) const {
0114   sl << std::setiosflags(std::ios::fixed);
0115   sl << std::setprecision(7);
0116   sl << "Acts::TrapezoidBounds:  (halfXnegY, halfXposY, halfY, rotAngle) = "
0117      << "(" << get(eHalfLengthXnegY) << ", " << get(eHalfLengthXposY) << ", "
0118      << get(eHalfLengthY) << ", " << get(eRotationAngle) << ")";
0119   sl << std::setprecision(-1);
0120   return sl;
0121 }
0122
0123 void TrapezoidBounds::rotateBoundingBox() noexcept(false) {
0124   const double rotAngle = get(eRotationAngle);
0125
0126   if (rotAngle != 0.) {
0127     m_boundingBox = ConvexPolygonBounds<4>(vertices()).boundingBox();
0128   }
0129 }
0130
0131 void TrapezoidBounds::checkConsistency() noexcept(false) {
0132   if (get(eHalfLengthXnegY) <= 0. || get(eHalfLengthXposY) <= 0.) {
0133     throw std::invalid_argument("TrapezoidBounds: invalid local x setup");
0134   }
0135   if (get(eHalfLengthY) <= 0.) {
0136     throw std::invalid_argument("TrapezoidBounds: invalid local y setup");
0137   }
0138 }
0139
0140 }  // namespace Acts