BoundaryCheck.hpp

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file BoundaryCheck.hpp

// This file is part of the Acts project.
//
// Copyright (C) 2016-2020 CERN for the benefit of the Acts project
//
// This Source Code Form is subject to the terms of the Mozilla Public
// License, v. 2.0. If a copy of the MPL was not distributed with this
// file, You can obtain one at http://mozilla.org/MPL/2.0/.

#pragma once
#include "Acts/Definitions/Algebra.hpp"
#include "Acts/Definitions/TrackParametrization.hpp"
#include "Acts/Surfaces/detail/VerticesHelper.hpp"

#include <cfloat>
#include <cmath>
#include <iterator>
#include <vector>

namespace Acts {

class BoundaryCheck {
 public:
  BoundaryCheck(bool check);

  BoundaryCheck(bool checkLocal0, bool checkLocal1, double tolerance0 = 0,
                double tolerance1 = 0);

  BoundaryCheck(const SquareMatrix2& localCovariance, double sigmaMax = 1);

  operator bool() const { return (m_type != Type::eNone); }
  bool operator!() const { return !bool(*this); }

  template <typename Vector2Container>
  bool isInside(const Vector2& point, const Vector2Container& vertices) const;

  bool isInside(const Vector2& point, const Vector2& lowerLeft,
                const Vector2& upperRight) const;

  template <typename Vector2Container>
  double distance(const Vector2& point, const Vector2Container& vertices) const;

  double distance(const Vector2& point, const Vector2& lowerLeft,
                  const Vector2& upperRight) const;

  enum class Type {
    eNone,      
    eAbsolute,  
    eChi2       
  };

  Type type() const;

  // Broadcast the tolerance
  const Vector2& tolerance() const;

  // Return the covariance matrix
  SquareMatrix2 covariance() const;

 private:
  BoundaryCheck transformed(const ActsMatrix<2, 2>& jacobian) const;

  bool isTolerated(const Vector2& delta) const;

  double squaredNorm(const Vector2& x) const;

  template <typename Vector2Container>
  Vector2 computeClosestPointOnPolygon(const Vector2& point,
                                       const Vector2Container& vertices) const;

  Vector2 computeEuclideanClosestPointOnRectangle(
      const Vector2& point, const Vector2& lowerLeft,
      const Vector2& upperRight) const;

  SquareMatrix2 m_weight;

  Vector2 m_tolerance;
  Type m_type;

  // To access the m_type
  friend class CylinderBounds;
  friend class RectangleBounds;
  // To be able to use `transformed`
  friend class CylinderBounds;
  friend class DiscTrapezoidBounds;
  // EllipseBounds needs a custom implementation
  friend class EllipseBounds;
};

}  // namespace Acts

inline Acts::BoundaryCheck::Type Acts::BoundaryCheck::type() const {
  return m_type;
}

inline const Acts::Vector2& Acts::BoundaryCheck::tolerance() const {
  return m_tolerance;
}

inline Acts::SquareMatrix2 Acts::BoundaryCheck::covariance() const {
  return m_weight.inverse();
}

inline Acts::BoundaryCheck::BoundaryCheck(bool check)
    : m_weight(SquareMatrix2::Identity()),
      m_tolerance(0, 0),
      m_type(check ? Type::eAbsolute : Type::eNone) {}

inline Acts::BoundaryCheck::BoundaryCheck(bool checkLocal0, bool checkLocal1,
                                          double tolerance0, double tolerance1)
    : m_weight(SquareMatrix2::Identity()),
      m_tolerance(checkLocal0 ? tolerance0 : DBL_MAX,
                  checkLocal1 ? tolerance1 : DBL_MAX),
      m_type(Type::eAbsolute) {}

inline Acts::BoundaryCheck::BoundaryCheck(const SquareMatrix2& localCovariance,
                                          double sigmaMax)
    : m_weight(localCovariance.inverse()),
      m_tolerance(sigmaMax, 0),
      m_type(Type::eChi2) {}

inline Acts::BoundaryCheck Acts::BoundaryCheck::transformed(
    const ActsMatrix<2, 2>& jacobian) const {
  BoundaryCheck bc = *this;
  if (m_type == Type::eAbsolute) {
    // project tolerances to the new system. depending on the jacobian we need
    // to check both tolerances, even when the initial check does not.
    bc.m_tolerance = (jacobian * m_tolerance).cwiseAbs();
  } else /* Type::eChi2 */ {
    bc.m_weight =
        (jacobian * m_weight.inverse() * jacobian.transpose()).inverse();
  }
  return bc;
}

template <typename Vector2Container>
inline bool Acts::BoundaryCheck::isInside(
    const Vector2& point, const Vector2Container& vertices) const {
  if (m_type == Type::eNone) {
    // The null boundary check always succeeds
    return true;
  } else if (detail::VerticesHelper::isInsidePolygon(point, vertices)) {
    // If the point falls inside the polygon, the check always succeeds
    return true;
  } else if (m_tolerance == Vector2(0., 0.)) {
    // Outside of the polygon, since we've eliminated the case of an absence of
    // check above, we know we'll always fail if the tolerance is zero.
    //
    // This allows us to avoid the expensive computeClosestPointOnPolygon
    // computation in this simple case.
    //
    // TODO: When tolerance is not 0, we could also avoid this computation in
    //       some cases by testing against a bounding box of the polygon, padded
    //       on each side with our tolerance. Check if this optimization is
    //       worthwhile in some production workflows, and if so implement it.
    return false;
  } else {
    // We are outside of the polygon, but there is a tolerance. Must find what
    // the closest point on the polygon is and check if it's within tolerance.
    auto closestPoint = computeClosestPointOnPolygon(point, vertices);
    return isTolerated(closestPoint - point);
  }
}

inline bool Acts::BoundaryCheck::isInside(const Vector2& point,
                                          const Vector2& lowerLeft,
                                          const Vector2& upperRight) const {
  if (detail::VerticesHelper::isInsideRectangle(point, lowerLeft, upperRight)) {
    return true;
  } else {
    Vector2 closestPoint;

    if (m_type == Type::eNone || m_type == Type::eAbsolute) {
      // absolute, can calculate directly
      closestPoint =
          computeEuclideanClosestPointOnRectangle(point, lowerLeft, upperRight);

    } else /* Type::eChi2 */ {
      // need to calculate by projection and squarednorm
      Vector2 vertices[] = {{lowerLeft[0], lowerLeft[1]},
                            {upperRight[0], lowerLeft[1]},
                            {upperRight[0], upperRight[1]},
                            {lowerLeft[0], upperRight[1]}};
      closestPoint = computeClosestPointOnPolygon(point, vertices);
    }

    return isTolerated(closestPoint - point);
  }
}

template <typename Vector2Container>
inline double Acts::BoundaryCheck::distance(
    const Acts::Vector2& point, const Vector2Container& vertices) const {
  // TODO 2017-04-06 msmk: this should be calculable directly
  double d = squaredNorm(point - computeClosestPointOnPolygon(point, vertices));
  d = std::sqrt(d);
  return detail::VerticesHelper::isInsidePolygon(point, vertices) ? -d : d;
}

inline double Acts::BoundaryCheck::distance(const Acts::Vector2& point,
                                            const Vector2& lowerLeft,
                                            const Vector2& upperRight) const {
  if (m_type == Type::eNone || m_type == Type::eAbsolute) {
    // compute closest point on box
    double d = (point - computeEuclideanClosestPointOnRectangle(
                            point, lowerLeft, upperRight))
                   .norm();
    return detail::VerticesHelper::isInsideRectangle(point, lowerLeft,
                                                     upperRight)
               ? -d
               : d;

  } else /* Type::eChi2 */ {
    Vector2 vertices[] = {{lowerLeft[0], lowerLeft[1]},
                          {upperRight[0], lowerLeft[1]},
                          {upperRight[0], upperRight[1]},
                          {lowerLeft[0], upperRight[1]}};
    return distance(point, vertices);
  }
}

inline bool Acts::BoundaryCheck::isTolerated(const Vector2& delta) const {
  if (m_type == Type::eNone) {
    return true;
  } else if (m_type == Type::eAbsolute) {
    return (std::abs(delta[0]) <= m_tolerance[0]) &&
           (std::abs(delta[1]) <= m_tolerance[1]);
  } else /* Type::eChi2 */ {
    // Mahalanobis distances mean is 2 in 2-dim. cut is 1-d sigma.
    return (squaredNorm(delta) < (2 * m_tolerance[0]));
  }
}

inline double Acts::BoundaryCheck::squaredNorm(const Vector2& x) const {
  return (x.transpose() * m_weight * x).value();
}

template <typename Vector2Container>
inline Acts::Vector2 Acts::BoundaryCheck::computeClosestPointOnPolygon(
    const Acts::Vector2& point, const Vector2Container& vertices) const {
  // calculate the closest position on the segment between `ll0` and `ll1` to
  // the point as measured by the metric induced by the weight matrix
  auto closestOnSegment = [&](auto&& ll0, auto&& ll1) {
    // normal vector and position of the closest point along the normal
    auto n = ll1 - ll0;
    auto weighted_n = m_weight * n;
    auto f = n.dot(weighted_n);
    auto u = std::isnormal(f)
                 ? (point - ll0).dot(weighted_n) / f
                 : 0.5;  // ll0 and ll1 are so close it doesn't matter
    // u must be in [0, 1] to still be on the polygon segment
    return ll0 + std::clamp(u, 0.0, 1.0) * n;
  };

  auto iv = std::begin(vertices);
  Vector2 l0 = *iv;
  Vector2 l1 = *(++iv);
  Vector2 closest = closestOnSegment(l0, l1);
  auto closestDist = squaredNorm(closest - point);
  // Calculate the closest point on other connecting lines and compare distances
  for (++iv; iv != std::end(vertices); ++iv) {
    l0 = l1;
    l1 = *iv;
    Vector2 current = closestOnSegment(l0, l1);
    auto currentDist = squaredNorm(current - point);
    if (currentDist < closestDist) {
      closest = current;
      closestDist = currentDist;
    }
  }
  // final edge from last vertex back to the first vertex
  Vector2 last = closestOnSegment(l1, *std::begin(vertices));
  if (squaredNorm(last - point) < closestDist) {
    closest = last;
  }
  return closest;
}

inline Acts::Vector2
Acts::BoundaryCheck::computeEuclideanClosestPointOnRectangle(
    const Vector2& point, const Vector2& lowerLeft,
    const Vector2& upperRight) const {
  /*
   *
   *        |                 |
   *   IV   |       V         | I
   *        |                 |
   *  ------------------------------
   *        |                 |
   *        |                 |
   *   VIII |     INSIDE      | VI
   *        |                 |
   *        |                 |
   *  ------------------------------
   *        |                 |
   *   III  |      VII        | II
   *        |                 |
   *
   */

  double l0 = point[0], l1 = point[1];
  double loc0Min = lowerLeft[0], loc0Max = upperRight[0];
  double loc1Min = lowerLeft[1], loc1Max = upperRight[1];

  // check if inside
  if (loc0Min <= l0 && l0 < loc0Max && loc1Min <= l1 && l1 < loc1Max) {
    // INSIDE
    double dist = std::abs(loc0Max - l0);
    Vector2 cls(loc0Max, l1);

    double test = std::abs(loc0Min - l0);
    if (test <= dist) {
      dist = test;
      cls = {loc0Min, l1};
    }

    test = std::abs(loc1Max - l1);
    if (test <= dist) {
      dist = test;
      cls = {l0, loc1Max};
    }

    test = std::abs(loc1Min - l1);
    if (test <= dist) {
      return {l0, loc1Min};
    }
    return cls;
  } else {
    // OUTSIDE, check sectors
    if (l0 > loc0Max) {
      if (l1 > loc1Max) {  // I
        return {loc0Max, loc1Max};
      } else if (l1 <= loc1Min) {  // II
        return {loc0Max, loc1Min};
      } else {  // VI
        return {loc0Max, l1};
      }
    } else if (l0 < loc0Min) {
      if (l1 > loc1Max) {  // IV
        return {loc0Min, loc1Max};
      } else if (l1 <= loc1Min) {  // III
        return {loc0Min, loc1Min};
      } else {  // VIII
        return {loc0Min, l1};
      }
    } else {
      if (l1 > loc1Max) {  // V
        return {l0, loc1Max};
      } else {  // l1 <= loc1Min # VII
        return {l0, loc1Min};
      }
      // third case not necessary, see INSIDE above
    }
  }
}