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