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