File indexing completed on 2025-08-05 08:09:28
0001
0002
0003
0004
0005
0006
0007
0008
0009 #include "Acts/Utilities/VectorHelpers.hpp"
0010
0011 template <typename value_t, std::size_t DIM, std::size_t SIDES>
0012 template <std::size_t D, std::enable_if_t<D == 2, int>>
0013 Acts::Frustum<value_t, DIM, SIDES>::Frustum(const VertexType& origin,
0014 const VertexType& dir,
0015 value_type opening_angle)
0016 : m_origin(origin) {
0017 using rotation_t = Eigen::Rotation2D<value_type>;
0018
0019 static_assert(SIDES == 2, "2D frustum can only have 2 sides");
0020 assert(opening_angle < M_PI);
0021
0022 translation_t translation(origin);
0023 value_type angle = VectorHelpers::phi(dir);
0024 Eigen::Rotation2D<value_type> rot(angle);
0025
0026 value_type normal_angle = 0.5 * M_PI - 0.5 * opening_angle;
0027 VertexType normal1 = rotation_t(normal_angle) * VertexType::UnitX();
0028 VertexType normal2 = rotation_t(-normal_angle) * VertexType::UnitX();
0029
0030 m_normals = {rot * VertexType::UnitX(), rot * normal1, rot * normal2};
0031 }
0032
0033 template <typename value_t, std::size_t DIM, std::size_t SIDES>
0034 template <std::size_t D, std::enable_if_t<D == 3, int>>
0035 Acts::Frustum<value_t, DIM, SIDES>::Frustum(const VertexType& origin,
0036 const VertexType& dir,
0037 value_type opening_angle)
0038 : m_origin(origin) {
0039 static_assert(SIDES > 2, "3D frustum must have 3 or more sides");
0040 assert(opening_angle < M_PI);
0041 using angle_axis_t = Eigen::AngleAxis<value_type>;
0042
0043 const VertexType ldir = VertexType::UnitZ();
0044 const VertexType lup = VertexType::UnitX();
0045
0046 transform_type transform;
0047 transform = (Eigen::Quaternion<value_type>().setFromTwoVectors(ldir, dir));
0048
0049 m_normals[0] = ldir;
0050
0051 const value_type phi_sep = 2 * M_PI / sides;
0052 transform_type rot;
0053 rot = angle_axis_t(phi_sep, ldir);
0054
0055 value_type half_opening_angle = opening_angle / 2.;
0056 auto calculate_normal =
0057 [&ldir, &half_opening_angle](const VertexType& out) -> VertexType {
0058 const VertexType tilt_axis = -1 * out.cross(ldir);
0059 return (-1 * (angle_axis_t(half_opening_angle, tilt_axis) * out))
0060 .normalized();
0061 };
0062
0063 VertexType current_outward = lup;
0064 m_normals[1] = calculate_normal(current_outward);
0065
0066 for (std::size_t i = 1; i < sides; i++) {
0067 current_outward = rot * current_outward;
0068 m_normals[i + 1] = calculate_normal(current_outward);
0069 }
0070
0071 for (auto& normal : m_normals) {
0072 normal = transform * normal;
0073 }
0074 }
0075
0076 template <typename value_t, std::size_t DIM, std::size_t SIDES>
0077 template <std::size_t D, std::enable_if_t<D == 3, int>>
0078 void Acts::Frustum<value_t, DIM, SIDES>::draw(IVisualization3D& helper,
0079 value_type far_distance) const {
0080 static_assert(DIM == 3, "Drawing is only supported in 3D");
0081
0082
0083
0084
0085
0086 VertexType far_normal = m_normals[0];
0087 VertexType far_center = m_normals[0] * far_distance;
0088 std::array<std::pair<VertexType, VertexType>, SIDES> planeFarIXs;
0089
0090 auto ixPlanePlane = [](const auto& n1, const auto& p1, const auto& n2,
0091 const auto& p2) -> std::pair<VertexType, VertexType> {
0092 const VertexType m = n1.cross(n2).normalized();
0093 const double j = (n2.dot(p2 - p1)) / (n2.dot(n1.cross(m)));
0094 const VertexType q = p1 + j * n1.cross(m);
0095 return {m, q};
0096 };
0097
0098 auto ixLineLine = [](const auto& p1, const auto& d1, const auto& p2,
0099 const auto& d2) -> VertexType {
0100 return p1 + (((p2 - p1).cross(d2)).norm() / (d1.cross(d2)).norm()) * d1;
0101 };
0102
0103
0104 for (std::size_t i = 1; i < n_normals; i++) {
0105 const auto ixLine =
0106 ixPlanePlane(far_normal, far_center, m_normals[i], VertexType::Zero());
0107 planeFarIXs.at(i - 1) = ixLine;
0108 }
0109
0110 std::array<VertexType, SIDES> points;
0111
0112 for (std::size_t i = 0; i < std::size(planeFarIXs); i++) {
0113 std::size_t j = (i + 1) % std::size(planeFarIXs);
0114 const auto& l1 = planeFarIXs.at(i);
0115 const auto& l2 = planeFarIXs.at(j);
0116 const VertexType ix =
0117 m_origin + ixLineLine(l1.second, l1.first, l2.second, l2.first);
0118 points.at(i) = ix;
0119 }
0120
0121 for (std::size_t i = 0; i < std::size(points); i++) {
0122 std::size_t j = (i + 1) % std::size(points);
0123 helper.face(
0124 std::vector<VertexType>({m_origin, points.at(i), points.at(j)}));
0125 }
0126 }
0127
0128 template <typename value_t, std::size_t DIM, std::size_t SIDES>
0129 template <std::size_t D, std::enable_if_t<D == 2, int>>
0130 std::ostream& Acts::Frustum<value_t, DIM, SIDES>::svg(std::ostream& os,
0131 value_type w,
0132 value_type h,
0133 value_type far_distance,
0134 value_type unit) const {
0135 static_assert(DIM == 2, "SVG is only supported in 2D");
0136
0137 VertexType mid(w / 2., h / 2.);
0138
0139
0140
0141 transform_type trf = transform_type::Identity();
0142 trf.translate(mid);
0143 trf = trf * Eigen::Scaling(VertexType(1, -1));
0144 trf.scale(unit);
0145
0146 std::array<std::string, 3> colors({"orange", "blue", "red"});
0147
0148 auto draw_line = [&](const VertexType& left_, const VertexType& right_,
0149 const std::string& color, std::size_t width) {
0150 VertexType left = trf * left_;
0151 VertexType right = trf * right_;
0152 os << "<line ";
0153
0154 os << "x1=\"" << left.x() << "\" ";
0155 os << "y1=\"" << left.y() << "\" ";
0156 os << "x2=\"" << right.x() << "\" ";
0157 os << "y2=\"" << right.y() << "\" ";
0158
0159 os << " stroke=\"" << color << "\" stroke-width=\"" << width << "\"/>\n";
0160 };
0161
0162 auto draw_point = [&](const VertexType& p_, const std::string& color,
0163 std::size_t r) {
0164 VertexType p = trf * p_;
0165 os << "<circle ";
0166 os << "cx=\"" << p.x() << "\" cy=\"" << p.y() << "\" r=\"" << r << "\"";
0167 os << " fill=\"" << color << "\"";
0168 os << "/>\n";
0169 };
0170
0171 using vec3 = Eigen::Matrix<value_type, 3, 1>;
0172 auto ixLineLine = [](const VertexType& p1_2, const VertexType& d1_2,
0173 const VertexType& p2_2,
0174 const VertexType& d2_2) -> VertexType {
0175 const vec3 p1(p1_2.x(), p1_2.y(), 0);
0176 const vec3 p2(p2_2.x(), p2_2.y(), 0);
0177 const vec3 d1(d1_2.x(), d1_2.y(), 0);
0178 const vec3 d2(d2_2.x(), d2_2.y(), 0);
0179
0180 vec3 num = (p2 - p1).cross(d2);
0181 vec3 den = d1.cross(d2);
0182
0183 value_type f = 1.;
0184 value_type dot = num.normalized().dot(den.normalized());
0185 if (std::abs(dot) - 1 < 1e-9 && dot < 0) {
0186 f = -1.;
0187 }
0188
0189 const vec3 p = p1 + f * (num.norm() / den.norm()) * d1;
0190 assert(std::abs(p.z()) < 1e-9);
0191 return {p.x(), p.y()};
0192 };
0193
0194 const VertexType far_dir = {m_normals[0].y(), -m_normals[0].x()};
0195 const VertexType far_point = m_normals[0] * far_distance;
0196
0197 std::array<VertexType, 2> points;
0198
0199 for (std::size_t i = 1; i < n_normals; i++) {
0200 VertexType plane_dir(m_normals[i].y(), -m_normals[i].x());
0201
0202 const VertexType ix = ixLineLine(far_point, far_dir, {0, 0}, plane_dir);
0203 draw_point(m_origin + ix, "black", 3);
0204 draw_line(m_origin, m_origin + ix, "black", 2);
0205 points.at(i - 1) = ix;
0206 }
0207
0208 draw_line(m_origin + points.at(0), m_origin + points.at(1), "black", 2);
0209
0210 draw_line(m_origin, m_origin + m_normals[0] * 2, colors[0], 3);
0211
0212 draw_point({0, 0}, "yellow", 5);
0213 draw_point(m_origin, "green", 5);
0214
0215 return os;
0216 }
0217
0218 template <typename value_t, std::size_t DIM, std::size_t SIDES>
0219 Acts::Frustum<value_t, DIM, SIDES>
0220 Acts::Frustum<value_t, DIM, SIDES>::transformed(
0221 const transform_type& trf) const {
0222 const auto& rot = trf.rotation();
0223
0224 std::array<VertexType, n_normals> new_normals;
0225 for (std::size_t i = 0; i < n_normals; i++) {
0226 new_normals[i] = rot * m_normals[i];
0227 }
0228
0229 return Frustum<value_t, DIM, SIDES>(trf * m_origin, std::move(new_normals));
0230 }