ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Frustum.ipp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Frustum.ipp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2018-2019 CERN for the benefit of the Acts project
4 //
5 // This Source Code Form is subject to the terms of the Mozilla Public
6 // License, v. 2.0. If a copy of the MPL was not distributed with this
7 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
8 
10 
11 template <typename value_t, size_t DIM, size_t SIDES>
12 template <size_t D, std::enable_if_t<D == 2, int>>
14  const VertexType& dir,
15  value_type opening_angle)
16  : m_origin(origin) {
17  using rotation_t = Eigen::Rotation2D<value_type>;
18 
19  static_assert(SIDES == 2, "2D frustum can only have 2 sides");
20  assert(opening_angle < M_PI);
21 
22  translation_t translation(origin);
24  Eigen::Rotation2D<value_type> rot(angle);
25 
26  value_type normal_angle = 0.5 * M_PI - 0.5 * opening_angle;
27  VertexType normal1 = rotation_t(normal_angle) * VertexType::UnitX();
28  VertexType normal2 = rotation_t(-normal_angle) * VertexType::UnitX();
29 
30  m_normals = {rot * VertexType::UnitX(), rot * normal1, rot * normal2};
31 }
32 
33 template <typename value_t, size_t DIM, size_t SIDES>
34 template <size_t D, std::enable_if_t<D == 3, int>>
36  const VertexType& dir,
37  value_type opening_angle)
38  : m_origin(origin) {
39  static_assert(SIDES > 2, "3D frustum must have 3 or more sides");
40  assert(opening_angle < M_PI);
41  using angle_axis_t = Eigen::AngleAxis<value_type>;
42 
43  const VertexType ldir = VertexType::UnitZ();
44  const VertexType lup = VertexType::UnitX();
45 
47  transform = (Eigen::Quaternion<value_type>().setFromTwoVectors(ldir, dir));
48 
49  m_normals[0] = ldir;
50 
51  const value_type phi_sep = 2 * M_PI / sides;
52  transform_type rot;
53  rot = angle_axis_t(phi_sep, ldir);
54 
55  value_type half_opening_angle = opening_angle / 2.;
56  auto calculate_normal =
57  [&ldir, &half_opening_angle](const VertexType& out) -> VertexType {
58  const VertexType tilt_axis = -1 * out.cross(ldir);
59  return (-1 * (angle_axis_t(half_opening_angle, tilt_axis) * out))
60  .normalized();
61  };
62 
63  VertexType current_outward = lup;
64  m_normals[1] = calculate_normal(current_outward);
65 
66  for (size_t i = 1; i < sides; i++) {
67  current_outward = rot * current_outward;
68  m_normals[i + 1] = calculate_normal(current_outward);
69  }
70 
71  for (auto& normal : m_normals) {
72  normal = transform * normal;
73  }
74 }
75 
76 template <typename value_t, size_t DIM, size_t SIDES>
77 template <size_t D, std::enable_if_t<D == 3, int>>
79  value_type far_distance) const {
80  static_assert(DIM == 3, "Drawing is only supported in 3D");
81 
82  // Iterate around normals, calculate cross with "far" plane
83  // to get intersection lines.
84  // Work in local reference frame of the frustum, and only convert to global
85  // right before drawing.
86  VertexType far_normal = m_normals[0]; // far has same normal as pseudo-near
87  VertexType far_center = m_normals[0] * far_distance;
88  std::array<std::pair<VertexType, VertexType>, SIDES> planeFarIXs;
89 
90  auto ixPlanePlane = [](const auto& n1, const auto& p1, const auto& n2,
91  const auto& p2) -> std::pair<VertexType, VertexType> {
92  const VertexType m = n1.cross(n2).normalized();
93  const double j = (n2.dot(p2 - p1)) / (n2.dot(n1.cross(m)));
94  const VertexType q = p1 + j * n1.cross(m);
95  return {m, q};
96  };
97 
98  auto ixLineLine = [](const auto& p1, const auto& d1, const auto& p2,
99  const auto& d2) -> VertexType {
100  return p1 + (((p2 - p1).cross(d2)).norm() / (d1.cross(d2)).norm()) * d1;
101  };
102 
103  // skip i=0 <=> pseudo-near
104  for (size_t i = 1; i < n_normals; i++) {
105  const auto ixLine =
106  ixPlanePlane(far_normal, far_center, m_normals[i], VertexType::Zero());
107  planeFarIXs.at(i - 1) = ixLine;
108  }
109 
110  std::array<VertexType, SIDES> points;
111 
112  for (size_t i = 0; i < std::size(planeFarIXs); i++) {
113  size_t j = (i + 1) % std::size(planeFarIXs);
114  const auto& l1 = planeFarIXs.at(i);
115  const auto& l2 = planeFarIXs.at(j);
116  const VertexType ix =
117  m_origin + ixLineLine(l1.second, l1.first, l2.second, l2.first);
118  points.at(i) = ix;
119  }
120 
121  for (size_t i = 0; i < std::size(points); i++) {
122  size_t j = (i + 1) % std::size(points);
123  helper.face(
124  std::vector<VertexType>({m_origin, points.at(i), points.at(j)}));
125  }
126 }
127 
128 template <typename value_t, size_t DIM, size_t SIDES>
129 template <size_t D, std::enable_if_t<D == 2, int>>
130 std::ostream& Acts::Frustum<value_t, DIM, SIDES>::svg(std::ostream& os,
131  value_type w,
132  value_type h,
133  value_type far_distance,
134  value_type unit) const {
135  static_assert(DIM == 2, "SVG is only supported in 2D");
136 
137  VertexType mid(w / 2., h / 2.);
138 
139  // set up transform for svg. +y is down, normally, and unit is pixels.
140  // We flip the y axis, and scale up by `unit`.
141  transform_type trf = transform_type::Identity();
142  trf.translate(mid);
143  trf = trf * Eigen::Scaling(VertexType(1, -1));
144  trf.scale(unit);
145 
146  std::array<std::string, 3> colors({"orange", "blue", "red"});
147 
148  auto draw_line = [&](const VertexType& left_, const VertexType& right_,
149  std::string color, size_t width) {
150  VertexType left = trf * left_;
151  VertexType right = trf * right_;
152  os << "<line ";
153 
154  os << "x1=\"" << left.x() << "\" ";
155  os << "y1=\"" << left.y() << "\" ";
156  os << "x2=\"" << right.x() << "\" ";
157  os << "y2=\"" << right.y() << "\" ";
158 
159  os << " stroke=\"" << color << "\" stroke-width=\"" << width << "\"/>\n";
160  };
161 
162  auto draw_point = [&](const VertexType& p_, std::string color, size_t r) {
163  VertexType p = trf * p_;
164  os << "<circle ";
165  os << "cx=\"" << p.x() << "\" cy=\"" << p.y() << "\" r=\"" << r << "\"";
166  os << " fill=\"" << color << "\"";
167  os << "/>\n";
168  };
169 
170  using vec3 = ActsVector<value_type, 3>;
171  auto ixLineLine = [](const VertexType& p1_2, const VertexType& d1_2,
172  const VertexType& p2_2,
173  const VertexType& d2_2) -> VertexType {
174  const vec3 p1(p1_2.x(), p1_2.y(), 0);
175  const vec3 p2(p2_2.x(), p2_2.y(), 0);
176  const vec3 d1(d1_2.x(), d1_2.y(), 0);
177  const vec3 d2(d2_2.x(), d2_2.y(), 0);
178 
179  vec3 num = (p2 - p1).cross(d2);
180  vec3 den = d1.cross(d2);
181 
182  value_type f = 1.;
183  value_type dot = num.normalized().dot(den.normalized());
184  if (std::abs(dot) - 1 < 1e-9 && dot < 0) {
185  f = -1.;
186  }
187 
188  const vec3 p = p1 + f * (num.norm() / den.norm()) * d1;
189  assert(std::abs(p.z()) < 1e-9);
190  return {p.x(), p.y()};
191  };
192 
193  const VertexType far_dir = {m_normals[0].y(), -m_normals[0].x()};
194  const VertexType far_point = m_normals[0] * far_distance;
195 
196  std::array<VertexType, 2> points;
197 
198  for (size_t i = 1; i < n_normals; i++) {
199  VertexType plane_dir(m_normals[i].y(), -m_normals[i].x());
200 
201  const VertexType ix = ixLineLine(far_point, far_dir, {0, 0}, plane_dir);
202  draw_point(m_origin + ix, "black", 3);
203  draw_line(m_origin, m_origin + ix, "black", 2);
204  points.at(i - 1) = ix;
205  }
206 
207  draw_line(m_origin + points.at(0), m_origin + points.at(1), "black", 2);
208 
209  draw_line(m_origin, m_origin + m_normals[0] * 2, colors[0], 3);
210 
211  draw_point({0, 0}, "yellow", 5);
212  draw_point(m_origin, "green", 5);
213 
214  return os;
215 }
216 
217 template <typename value_t, size_t DIM, size_t SIDES>
220  const transform_type& trf) const {
221  const auto& rot = trf.rotation();
222 
223  std::array<VertexType, n_normals> new_normals;
224  for (size_t i = 0; i < n_normals; i++) {
225  new_normals[i] = rot * m_normals[i];
226  }
227 
228  return Frustum<value_t, DIM, SIDES>(trf * m_origin, std::move(new_normals));
229 }