src/Geometry/GenericCuboidVolumeBounds.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/GenericCuboidVolumeBounds.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/Direction.hpp"
0013 #include "Acts/Geometry/Volume.hpp"
0014 #include "Acts/Surfaces/ConvexPolygonBounds.hpp"
0015 #include "Acts/Surfaces/PlaneSurface.hpp"
0016 #include "Acts/Surfaces/Surface.hpp"
0017 #include "Acts/Utilities/BoundingBox.hpp"
0018 #include "Acts/Visualization/IVisualization3D.hpp"
0019
0020 #include <algorithm>
0021 #include <array>
0022 #include <cmath>
0023 #include <cstddef>
0024 #include <memory>
0025 #include <ostream>
0026 #include <stdexcept>
0027 #include <type_traits>
0028 #include <utility>
0029
0030 Acts::GenericCuboidVolumeBounds::GenericCuboidVolumeBounds(
0031     const std::array<Acts::Vector3, 8>& vertices) noexcept(false)
0032     : m_vertices(vertices) {
0033   construct();
0034 }
0035
0036 Acts::GenericCuboidVolumeBounds::GenericCuboidVolumeBounds(
0037     const std::array<double, GenericCuboidVolumeBounds::BoundValues::eSize>&
0038         values) noexcept(false)
0039     : m_vertices() {
0040   for (std::size_t iv = 0; iv < 8; ++iv) {
0041     m_vertices[iv] =
0042         Vector3(values[iv * 3], values[iv * 3 + 1], values[iv * 3 + 2]);
0043   }
0044   construct();
0045 }
0046
0047 bool Acts::GenericCuboidVolumeBounds::inside(const Acts::Vector3& gpos,
0048                                              double tol) const {
0049   constexpr std::array<std::size_t, 6> vtxs = {0, 4, 0, 1, 2, 1};
0050   // needs to be on same side, get ref
0051   bool ref = std::signbit((gpos - m_vertices[vtxs[0]]).dot(m_normals[0]));
0052   for (std::size_t i = 1; i < 6; i++) {
0053     double dot = (gpos - m_vertices[vtxs[i]]).dot(m_normals[i]);
0054     if (std::signbit(dot) != ref) {
0055       // technically outside, but how far?
0056       if (std::abs(dot) > tol) {
0057         // distance greater than tol
0058         return false;
0059       }
0060       // distance smaller than tol, ignore
0061     }
0062   }
0063   return true;
0064 }
0065
0066 std::vector<Acts::OrientedSurface>
0067 Acts::GenericCuboidVolumeBounds::orientedSurfaces(
0068     const Transform3& transform) const {
0069   std::vector<OrientedSurface> oSurfaces;
0070
0071   // approximate cog of the volume
0072   Vector3 cog(0, 0, 0);
0073
0074   for (std::size_t i = 0; i < 8; i++) {
0075     cog += m_vertices[i];
0076   }
0077
0078   cog *= 0.125;  // 1/8.
0079
0080   auto make_surface = [&](const auto& a, const auto& b, const auto& c,
0081                           const auto& d) -> void {
0082     // calculate centroid of these points
0083     Vector3 ctrd = (a + b + c + d) / 4.;
0084     // create normal
0085     const Vector3 ab = b - a, ac = c - a;
0086     Vector3 normal = ab.cross(ac).normalized();
0087
0088     Direction dir = Direction::fromScalar((cog - d).dot(normal));
0089
0090     // build transform from z unit to normal
0091     // z is normal in local coordinates
0092     // Volume local to surface local
0093     Transform3 vol2srf;
0094
0095     // GCC13+ Complains about maybe uninitialized memory inside Eigen's SVD code
0096     // This warning is ignored in this compilation unit by using the pragma at
0097     // the top of this file.
0098     vol2srf = (Eigen::Quaternion<Transform3::Scalar>().setFromTwoVectors(
0099         normal, Vector3::UnitZ()));
0100
0101     vol2srf = vol2srf * Translation3(-ctrd);
0102
0103     // now calculate position of vertices in surface local frame
0104     Vector3 a_l, b_l, c_l, d_l;
0105     a_l = vol2srf * a;
0106     b_l = vol2srf * b;
0107     c_l = vol2srf * c;
0108     d_l = vol2srf * d;
0109
0110     std::vector<Vector2> vertices({{a_l.x(), a_l.y()},
0111                                    {b_l.x(), b_l.y()},
0112                                    {c_l.x(), c_l.y()},
0113                                    {d_l.x(), d_l.y()}});
0114
0115     auto polyBounds = std::make_shared<const ConvexPolygonBounds<4>>(vertices);
0116     auto srfTrf = transform * vol2srf.inverse();
0117     auto srf = Surface::makeShared<PlaneSurface>(srfTrf, polyBounds);
0118
0119     oSurfaces.push_back(OrientedSurface{std::move(srf), dir});
0120   };
0121
0122   make_surface(m_vertices[0], m_vertices[1], m_vertices[2], m_vertices[3]);
0123   make_surface(m_vertices[4], m_vertices[5], m_vertices[6], m_vertices[7]);
0124   make_surface(m_vertices[0], m_vertices[3], m_vertices[7], m_vertices[4]);
0125   make_surface(m_vertices[1], m_vertices[2], m_vertices[6], m_vertices[5]);
0126   make_surface(m_vertices[2], m_vertices[3], m_vertices[7], m_vertices[6]);
0127   make_surface(m_vertices[1], m_vertices[0], m_vertices[4], m_vertices[5]);
0128
0129   return oSurfaces;
0130 }
0131
0132 std::ostream& Acts::GenericCuboidVolumeBounds::toStream(
0133     std::ostream& sl) const {
0134   sl << "Acts::GenericCuboidVolumeBounds: vertices (x, y, z) =\n";
0135   for (std::size_t i = 0; i < 8; i++) {
0136     if (i > 0) {
0137       sl << ",\n";
0138     }
0139     sl << "[" << m_vertices[i].transpose() << "]";
0140   }
0141   return sl;
0142 }
0143
0144 void Acts::GenericCuboidVolumeBounds::construct() noexcept(false) {
0145   // calculate approximate center of gravity first, so we can make sure
0146   // the normals point inwards
0147   Vector3 cog(0, 0, 0);
0148
0149   for (std::size_t i = 0; i < 8; i++) {
0150     cog += m_vertices[i];
0151   }
0152
0153   cog *= 0.125;  // 1/8.
0154
0155   std::size_t idx = 0;
0156
0157   auto handle_face = [&](const auto& a, const auto& b, const auto& c,
0158                          const auto& d) {
0159     // we assume a b c d to be counter clockwise
0160     const Vector3 ab = b - a, ac = c - a;
0161     Vector3 normal = ab.cross(ac).normalized();
0162
0163     if ((cog - a).dot(normal) < 0) {
0164       // normal points outwards, flip normal
0165       normal *= -1.;
0166     }
0167
0168     // get rid of -0 values if present
0169     normal += Vector3::Zero();
0170
0171     // check if d is on the surface
0172     if (std::abs((a - d).dot(normal)) > 1e-6) {
0173       throw(std::invalid_argument(
0174           "GenericCuboidBounds: Four points do not lie on the same plane!"));
0175     }
0176
0177     m_normals[idx] = normal;
0178     idx++;
0179   };
0180
0181   // handle faces
0182   handle_face(m_vertices[0], m_vertices[1], m_vertices[2], m_vertices[3]);
0183   handle_face(m_vertices[4], m_vertices[5], m_vertices[6], m_vertices[7]);
0184   handle_face(m_vertices[0], m_vertices[3], m_vertices[7], m_vertices[4]);
0185   handle_face(m_vertices[1], m_vertices[2], m_vertices[6], m_vertices[5]);
0186   handle_face(m_vertices[2], m_vertices[3], m_vertices[7], m_vertices[6]);
0187   handle_face(m_vertices[1], m_vertices[0], m_vertices[4], m_vertices[5]);
0188 }
0189
0190 std::vector<double> Acts::GenericCuboidVolumeBounds::values() const {
0191   std::vector<double> rvalues;
0192   rvalues.reserve(BoundValues::eSize);
0193   for (std::size_t iv = 0; iv < 8; ++iv) {
0194     for (std::size_t ic = 0; ic < 3; ++ic) {
0195       rvalues.push_back(m_vertices[iv][ic]);
0196     }
0197   }
0198   return rvalues;
0199 }
0200
0201 Acts::Volume::BoundingBox Acts::GenericCuboidVolumeBounds::boundingBox(
0202     const Acts::Transform3* trf, const Vector3& envelope,
0203     const Volume* entity) const {
0204   Vector3 vmin, vmax;
0205
0206   Transform3 transform = Transform3::Identity();
0207   if (trf != nullptr) {
0208     transform = *trf;
0209   }
0210
0211   vmin = transform * m_vertices[0];
0212   vmax = transform * m_vertices[0];
0213
0214   for (std::size_t i = 1; i < 8; i++) {
0215     Vector3 vtx = transform * m_vertices[i];
0216     vmin = vmin.cwiseMin(vtx);
0217     vmax = vmax.cwiseMax(vtx);
0218   }
0219
0220   return {entity, vmin - envelope, vmax + envelope};
0221 }
0222
0223 void Acts::GenericCuboidVolumeBounds::draw(IVisualization3D& helper,
0224                                            const Transform3& transform) const {
0225   auto draw_face = [&](const auto& a, const auto& b, const auto& c,
0226                        const auto& d) {
0227     helper.face(std::vector<Vector3>(
0228         {transform * a, transform * b, transform * c, transform * d}));
0229   };
0230
0231   draw_face(m_vertices[0], m_vertices[1], m_vertices[2], m_vertices[3]);
0232   draw_face(m_vertices[4], m_vertices[5], m_vertices[6], m_vertices[7]);
0233   draw_face(m_vertices[0], m_vertices[3], m_vertices[7], m_vertices[4]);
0234   draw_face(m_vertices[1], m_vertices[2], m_vertices[6], m_vertices[5]);
0235   draw_face(m_vertices[2], m_vertices[3], m_vertices[7], m_vertices[6]);
0236   draw_face(m_vertices[1], m_vertices[0], m_vertices[4], m_vertices[5]);
0237 }