src/Geometry/Polyhedron.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/Geometry/Polyhedron.hpp"
0010
0011 #include "Acts/Surfaces/detail/VerticesHelper.hpp"
0012 #include "Acts/Utilities/BinningType.hpp"
0013 #include "Acts/Visualization/IVisualization3D.hpp"
0014
0015 #include <algorithm>
0016 #include <cmath>
0017 #include <limits>
0018 #include <numbers>
0019 #include <utility>
0020
0021 void Acts::Polyhedron::merge(const Acts::Polyhedron& other) {
0022   std::size_t cvert = vertices.size();
0023   vertices.insert(vertices.end(), other.vertices.begin(), other.vertices.end());
0024   /// Add the new faces with offsets
0025   auto join = [&](std::vector<FaceType>& existing,
0026                   const std::vector<FaceType>& additional) -> void {
0027     for (const auto& aface : additional) {
0028       FaceType nface = aface;
0029       std::transform(nface.begin(), nface.end(), nface.begin(),
0030                      [&](std::size_t x) { return (x + cvert); });
0031       existing.push_back(nface);
0032     }
0033   };
0034   // For faces and triangular mesh
0035   join(faces, other.faces);
0036   join(triangularMesh, other.triangularMesh);
0037 }
0038
0039 void Acts::Polyhedron::move(const Transform3& transform) {
0040   for_each(vertices.begin(), vertices.end(),
0041            [&](auto& v) { v = transform * v; });
0042 }
0043
0044 Acts::Extent Acts::Polyhedron::extent(const Transform3& transform) const {
0045   Extent extent;
0046   auto vtxs = vertices;
0047   std::transform(vtxs.begin(), vtxs.end(), vtxs.begin(), [&](auto& v) {
0048     auto vt = (transform * v);
0049     extent.extend(vt);
0050     return (vt);
0051   });
0052
0053   // Special checks of AxisDirection::AxisR for hyper plane surfaces
0054   if (detail::VerticesHelper::onHyperPlane(vtxs)) {
0055     // Check inclusion of origin (i.e. convex around origin)
0056     Vector3 origin =
0057         transform * Vector3(0., 0., extent.medium(AxisDirection::AxisZ));
0058     for (const auto& face : faces) {
0059       std::vector<Vector3> tface;
0060       tface.reserve(face.size());
0061       for (auto f : face) {
0062         tface.push_back(vtxs[f]);
0063       }
0064       if (detail::VerticesHelper::isInsidePolygon(origin, tface)) {
0065         extent.range(AxisDirection::AxisR).setMin(0.);
0066         extent.range(AxisDirection::AxisPhi)
0067             .set(-std::numbers::pi, std::numbers::pi);
0068         break;
0069       }
0070     }
0071     if (exact) {
0072       // Check for radial extend in 2D
0073       auto radialDistance = [&](const Vector3& pos1,
0074                                 const Vector3& pos2) -> double {
0075         Vector2 O(0, 0);
0076         Vector2 p1p2 = (pos2.block<2, 1>(0, 0) - pos1.block<2, 1>(0, 0));
0077         double L = p1p2.norm();
0078         Vector2 p1O = (O - pos1.block<2, 1>(0, 0));
0079
0080         // Don't try parallel lines
0081         if (L < 1e-7) {
0082           return std::numeric_limits<double>::max();
0083         }
0084         double f = p1p2.dot(p1O) / L;
0085
0086         // Clamp to [0, |p1p2|]
0087         f = std::min(L, std::max(0., f));
0088         Vector2 closest = f * p1p2.normalized() + pos1.block<2, 1>(0, 0);
0089         double dist = (closest - O).norm();
0090         return dist;
0091       };
0092
0093       for (std::size_t iv = 1; iv < vtxs.size() + 1; ++iv) {
0094         std::size_t fpoint = iv < vtxs.size() ? iv : 0;
0095         double testR = radialDistance(vtxs[fpoint], vtxs[iv - 1]);
0096         extent.range(AxisDirection::AxisR).expandMin(testR);
0097       }
0098     }
0099   }
0100   return extent;
0101 }
0102
0103 void Acts::Polyhedron::visualize(IVisualization3D& helper,
0104                                  const ViewConfig& viewConfig) const {
0105   if (viewConfig.visible) {
0106     if (!viewConfig.triangulate) {
0107       helper.faces(vertices, faces, viewConfig.color);
0108     } else {
0109       helper.faces(vertices, triangularMesh, viewConfig.color);
0110     }
0111   }
0112 }