ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
BoundaryCheck.hpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file BoundaryCheck.hpp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2016-2020 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 
9 #pragma once
10 #include <cfloat>
11 #include <cmath>
12 #include <iterator>
13 #include <vector>
14 
18 
19 namespace Acts {
20 
37  public:
39  BoundaryCheck(bool check);
40 
47  BoundaryCheck(bool checkLocal0, bool checkLocal1, double tolerance0 = 0,
48  double tolerance1 = 0);
49 
54  BoundaryCheck(const ActsSymMatrixD<2>& localCovariance, double sigmaMax = 1);
55 
56  operator bool() const { return (m_type != Type::eNone); }
57  bool operator!() const { return !bool(*this); }
58 
69  template <typename Vector2DContainer>
70  bool isInside(const Vector2D& point, const Vector2DContainer& vertices) const;
71 
82  bool isInside(const Vector2D& point, const Vector2D& lowerLeft,
83  const Vector2D& upperRight) const;
84 
96  template <typename Vector2DContainer>
97  double distance(const Vector2D& point,
98  const Vector2DContainer& vertices) const;
99 
111  double distance(const Vector2D& point, const Vector2D& lowerLeft,
112  const Vector2D& upperRight) const;
113 
114  enum class Type {
115  eNone,
116  eAbsolute,
117  eChi2
118  };
119 
121  Type type() const;
122 
123  // Broadcast the tolerance
124  const Vector2D& tolerance() const;
125 
126  // Return the covariance matrix
128 
129  private:
134  BoundaryCheck transformed(const ActsMatrixD<2, 2>& jacobian) const;
135 
137  bool isTolerated(const Vector2D& delta) const;
138 
140  double squaredNorm(const Vector2D& x) const;
141 
143  template <typename Vector2DContainer>
145  const Vector2D& point, const Vector2DContainer& vertices) const;
146 
149  const Vector2D& point, const Vector2D& lowerLeft,
150  const Vector2D& upperRight) const;
151 
154 
158 
159  // To acces the m_type
160  friend class CylinderBounds;
161  friend class RectangleBounds;
162  // To be able to use `transformed`
163  friend class CylinderBounds;
164  friend class DiscTrapezoidBounds;
165  // EllipseBounds needs a custom implementation
166  friend class EllipseBounds;
167 };
168 
169 } // namespace Acts
170 
172  return m_type;
173 }
174 
176  return m_tolerance;
177 }
178 
180  return m_weight.inverse();
181 }
182 
184  : m_weight(ActsSymMatrixD<2>::Identity()),
185  m_tolerance(0, 0),
186  m_type(check ? Type::eAbsolute : Type::eNone) {}
187 
188 inline Acts::BoundaryCheck::BoundaryCheck(bool checkLocal0, bool checkLocal1,
189  double tolerance0, double tolerance1)
190  : m_weight(ActsSymMatrixD<2>::Identity()),
191  m_tolerance(checkLocal0 ? tolerance0 : DBL_MAX,
192  checkLocal1 ? tolerance1 : DBL_MAX),
193  m_type(Type::eAbsolute) {}
194 
196  const ActsSymMatrixD<2>& localCovariance, double sigmaMax)
197  : m_weight(localCovariance.inverse()),
198  m_tolerance(sigmaMax, 0),
199  m_type(Type::eChi2) {}
200 
202  const ActsMatrixD<2, 2>& jacobian) const {
203  BoundaryCheck bc = *this;
204  if (m_type == Type::eAbsolute) {
205  // project tolerances to the new system. depending on the jacobian we need
206  // to check both tolerances, even when the initial check does not.
207  bc.m_tolerance = (jacobian * m_tolerance).cwiseAbs();
208  } else /* Type::eChi2 */ {
209  bc.m_weight =
210  (jacobian * m_weight.inverse() * jacobian.transpose()).inverse();
211  }
212  return bc;
213 }
214 
215 template <typename Vector2DContainer>
217  const Vector2D& point, const Vector2DContainer& vertices) const {
218  if (m_type == Type::eNone) {
219  // The null boundary check always succeeds
220  return true;
221  } else if (detail::VerticesHelper::isInsidePolygon(point, vertices)) {
222  // If the point falls inside the polygon, the check always succeeds
223  return true;
224  } else if (m_tolerance == Vector2D(0., 0.)) {
225  // Outside of the polygon, since we've eliminated the case of an absence of
226  // check above, we know we'll always fail if the tolerance is zero.
227  //
228  // This allows us to avoid the expensive computeClosestPointOnPolygon
229  // computation in this simple case.
230  //
231  // TODO: When tolerance is not 0, we could also avoid this computation in
232  // some cases by testing against a bounding box of the polygon, padded
233  // on each side with our tolerance. Check if this optimization is
234  // worthwhile in some production workflows, and if so implement it.
235  return false;
236  } else {
237  // We are outside of the polygon, but there is a tolerance. Must find what
238  // the closest point on the polygon is and check if it's within tolerance.
239  auto closestPoint = computeClosestPointOnPolygon(point, vertices);
240  return isTolerated(closestPoint - point);
241  }
242 }
243 
245  const Vector2D& lowerLeft,
246  const Vector2D& upperRight) const {
247  if (detail::VerticesHelper::isInsideRectangle(point, lowerLeft, upperRight)) {
248  return true;
249  } else {
250  Vector2D closestPoint;
251 
252  if (m_type == Type::eNone || m_type == Type::eAbsolute) {
253  // absolute, can calculate directly
254  closestPoint =
255  computeEuclideanClosestPointOnRectangle(point, lowerLeft, upperRight);
256 
257  } else /* Type::eChi2 */ {
258  // need to calculate by projection and squarednorm
259  Vector2D vertices[] = {{lowerLeft[0], lowerLeft[1]},
260  {upperRight[0], lowerLeft[1]},
261  {upperRight[0], upperRight[1]},
262  {lowerLeft[0], upperRight[1]}};
263  closestPoint = computeClosestPointOnPolygon(point, vertices);
264  }
265 
266  return isTolerated(closestPoint - point);
267  }
268 }
269 
270 template <typename Vector2DContainer>
272  const Acts::Vector2D& point, const Vector2DContainer& vertices) const {
273  // TODO 2017-04-06 msmk: this should be calculable directly
274  double d = squaredNorm(point - computeClosestPointOnPolygon(point, vertices));
275  d = std::sqrt(d);
276  return detail::VerticesHelper::isInsidePolygon(point, vertices) ? -d : d;
277 }
278 
280  const Vector2D& lowerLeft,
281  const Vector2D& upperRight) const {
282  if (m_type == Type::eNone || m_type == Type::eAbsolute) {
283  // compute closest point on box
284  double d = (point - computeEuclideanClosestPointOnRectangle(
285  point, lowerLeft, upperRight))
286  .norm();
287  return detail::VerticesHelper::isInsideRectangle(point, lowerLeft,
288  upperRight)
289  ? -d
290  : d;
291 
292  } else /* Type::eChi2 */ {
293  Vector2D vertices[] = {{lowerLeft[0], lowerLeft[1]},
294  {upperRight[0], lowerLeft[1]},
295  {upperRight[0], upperRight[1]},
296  {lowerLeft[0], upperRight[1]}};
297  return distance(point, vertices);
298  }
299 }
300 
302  if (m_type == Type::eNone) {
303  return true;
304  } else if (m_type == Type::eAbsolute) {
305  return (std::abs(delta[0]) <= m_tolerance[0]) &&
306  (std::abs(delta[1]) <= m_tolerance[1]);
307  } else /* Type::eChi2 */ {
308  // Mahalanobis distances mean is 2 in 2-dim. cut is 1-d sigma.
309  return (squaredNorm(delta) < (2 * m_tolerance[0]));
310  }
311 }
312 
313 inline double Acts::BoundaryCheck::squaredNorm(const Vector2D& x) const {
314  return (x.transpose() * m_weight * x).value();
315 }
316 
317 template <typename Vector2DContainer>
319  const Acts::Vector2D& point, const Vector2DContainer& vertices) const {
320  // calculate the closest position on the segment between `ll0` and `ll1` to
321  // the point as measured by the metric induced by the weight matrix
322  auto closestOnSegment = [&](auto&& ll0, auto&& ll1) {
323  // normal vector and position of the closest point along the normal
324  auto n = ll1 - ll0;
325  auto weighted_n = m_weight * n;
326  auto f = n.dot(weighted_n);
327  auto u = std::isnormal(f)
328  ? (point - ll0).dot(weighted_n) / f
329  : 0.5; // ll0 and ll1 are so close it doesn't matter
330  // u must be in [0, 1] to still be on the polygon segment
331  return ll0 + std::clamp(u, 0.0, 1.0) * n;
332  };
333 
334  auto iv = std::begin(vertices);
335  Vector2D l0 = *iv;
336  Vector2D l1 = *(++iv);
337  Vector2D closest = closestOnSegment(l0, l1);
338  auto closestDist = squaredNorm(closest - point);
339  // Calculate the closest point on other connecting lines and compare distances
340  for (++iv; iv != std::end(vertices); ++iv) {
341  l0 = l1;
342  l1 = *iv;
343  Vector2D current = closestOnSegment(l0, l1);
344  auto currentDist = squaredNorm(current - point);
345  if (currentDist < closestDist) {
346  closest = current;
347  closestDist = currentDist;
348  }
349  }
350  // final edge from last vertex back to the first vertex
351  Vector2D last = closestOnSegment(l1, *std::begin(vertices));
352  if (squaredNorm(last - point) < closestDist) {
353  closest = last;
354  }
355  return closest;
356 }
357 
358 inline Acts::Vector2D
360  const Vector2D& point, const Vector2D& lowerLeft,
361  const Vector2D& upperRight) const {
362  /*
363  *
364  * | |
365  * IV | V | I
366  * | |
367  * ------------------------------
368  * | |
369  * | |
370  * VIII | INSIDE | VI
371  * | |
372  * | |
373  * ------------------------------
374  * | |
375  * III | VII | II
376  * | |
377  *
378  */
379 
380  double l0 = point[0], l1 = point[1];
381  double loc0Min = lowerLeft[0], loc0Max = upperRight[0];
382  double loc1Min = lowerLeft[1], loc1Max = upperRight[1];
383 
384  // check if inside
385  if (loc0Min <= l0 && l0 < loc0Max && loc1Min <= l1 && l1 < loc1Max) {
386  // INSIDE
387  double dist = std::abs(loc0Max - l0);
388  Vector2D cls(loc0Max, l1);
389 
390  double test = std::abs(loc0Min - l0);
391  if (test <= dist) {
392  dist = test;
393  cls = {loc0Min, l1};
394  }
395 
396  test = std::abs(loc1Max - l1);
397  if (test <= dist) {
398  dist = test;
399  cls = {l0, loc1Max};
400  }
401 
402  test = std::abs(loc1Min - l1);
403  if (test <= dist) {
404  return {l0, loc1Min};
405  }
406  return cls;
407  } else {
408  // OUTSIDE, check sectors
409  if (l0 > loc0Max) {
410  if (l1 > loc1Max) { // I
411  return {loc0Max, loc1Max};
412  } else if (l1 <= loc1Min) { // II
413  return {loc0Max, loc1Min};
414  } else { // VI
415  return {loc0Max, l1};
416  }
417  } else if (l0 < loc0Min) {
418  if (l1 > loc1Max) { // IV
419  return {loc0Min, loc1Max};
420  } else if (l1 <= loc1Min) { // III
421  return {loc0Min, loc1Min};
422  } else { // VIII
423  return {loc0Min, l1};
424  }
425  } else {
426  if (l1 > loc1Max) { // V
427  return {l0, loc1Max};
428  } else { // l1 <= loc1Min # VII
429  return {l0, loc1Min};
430  }
431  // third case not necessary, see INSIDE above
432  }
433  }
434 }