Merge branch 'main' of github.com:acts-project/acts into tmp-distance…

…-to-boundary

Core/include/Acts/Surfaces/BoundaryTolerance.hpp

@@ -67,7 +67,12 @@ class BoundaryTolerance {
 
     AbsoluteBound() = default;
     AbsoluteBound(double tolerance0_, double tolerance1_)
-        : tolerance0(tolerance0_), tolerance1(tolerance1_) {}
+        : tolerance0(tolerance0_), tolerance1(tolerance1_) {
+      if (tolerance0 < 0 || tolerance1 < 0) {
+        throw std::invalid_argument(
+            "AbsoluteBound: Tolerance must be non-negative");
+      }
+    }
   };
 
   /// Absolute tolerance in Cartesian coordinates
@@ -77,7 +82,12 @@ class BoundaryTolerance {
 
     AbsoluteCartesian() = default;
     AbsoluteCartesian(double tolerance0_, double tolerance1_)
-        : tolerance0(tolerance0_), tolerance1(tolerance1_) {}
+        : tolerance0(tolerance0_), tolerance1(tolerance1_) {
+      if (tolerance0 < 0 || tolerance1 < 0) {
+        throw std::invalid_argument(
+            "AbsoluteCartesian: Tolerance must be non-negative");
+      }
+    }
   };
 
   /// Absolute tolerance in Euclidean distance
@@ -98,6 +108,12 @@ class BoundaryTolerance {
         : maxChi2(maxChi2_), weight(weight_) {}
   };
 
+  enum class ToleranceMode {
+    Extend,  // Extend the boundary
+    None,    // No tolerance
+    Shrink   // Shrink the boundary
+  };
+
   /// Underlying variant type
   using Variant = std::variant<Infinite, None, AbsoluteBound, AbsoluteCartesian,
                                AbsoluteEuclidean, Chi2Bound>;
@@ -132,7 +148,7 @@ class BoundaryTolerance {
   bool hasChi2Bound() const;
 
   /// Check if any tolerance is set.
-  bool hasTolerance() const;
+  ToleranceMode toleranceMode() const;
 
   /// Get the tolerance as absolute bound.
   AbsoluteBound asAbsoluteBound(bool isCartesian = false) const;

Core/include/Acts/Surfaces/SurfaceBounds.hpp

@@ -95,16 +95,21 @@ class SurfaceBounds {
   /// @return boolean indicator for the success of this operation
   virtual bool inside(const Vector2& lposition,
                       const BoundaryTolerance& boundaryTolerance) const {
+    using enum BoundaryTolerance::ToleranceMode;
+
     if (boundaryTolerance.isInfinite()) {
       return true;
     }
 
-    if (inside(lposition)) {
-      return true;
+    BoundaryTolerance::ToleranceMode mode = boundaryTolerance.toleranceMode();
+    bool strictlyInside = inside(lposition);
+
+    if (mode == None) {
+      return strictlyInside;
     }
 
-    if (!boundaryTolerance.hasTolerance()) {
-      return false;
+    if (mode == Extend && strictlyInside) {
+      return true;
     }
 
     std::optional<SquareMatrix2> jacobian;
@@ -120,6 +125,11 @@ class SurfaceBounds {
 
     Vector2 closest = closestPoint(lposition, metric);
     Vector2 distance = closest - lposition;
+
+    if (mode == Shrink) {
+      return boundaryTolerance.isTolerated(distance, jacobian) &&
+             strictlyInside;
+    }
     return boundaryTolerance.isTolerated(distance, jacobian);
   }
 

Core/src/Surfaces/AnnulusBounds.cpp

@@ -294,10 +294,9 @@ Vector2 AnnulusBounds::closestPoint(
   jacobianStripPCToModulePC(1, 0) = -(B * O_y - C * O_x) / A;
   jacobianStripPCToModulePC(1, 1) = r_strip * (B * O_x + C * O_y + r_strip) / A;
 
-  // Mahalanobis distance uses the coordinate weights (usually inverse
-  // covariance) as a metric
-  auto metricModulePC = jacobianStripPCToModulePC.transpose() *
-                        metric.value_or(SquareMatrix2::Identity()) *
+  // Mahalanobis distance uses inverse covariance as weights
+  const auto& weightStripPC = metric.value_or(SquareMatrix2::Identity());
+  auto weightModulePC = jacobianStripPCToModulePC.transpose() * weightStripPC *
                         jacobianStripPCToModulePC;
 
   double minDist = std::numeric_limits<double>::max();
@@ -309,35 +308,36 @@ Vector2 AnnulusBounds::closestPoint(
 
   // first: STRIP system. locpo is in STRIP PC already
   currentClosest = closestOnSegment(m_inLeftStripPC, m_outLeftStripPC,
-                                    locpo_rotated, metricModulePC);
-  currentDist = squaredNorm(locpo_rotated - currentClosest, metricModulePC);
+                                    locpo_rotated, weightStripPC);
+  currentDist = squaredNorm(locpo_rotated - currentClosest, weightStripPC);
   minDist = currentDist;
 
   currentClosest = closestOnSegment(m_inRightStripPC, m_outRightStripPC,
-                                    locpo_rotated, metricModulePC);
-  currentDist = squaredNorm(locpo_rotated - currentClosest, metricModulePC);
+                                    locpo_rotated, weightStripPC);
+  currentDist = squaredNorm(locpo_rotated - currentClosest, weightStripPC);
   if (currentDist < minDist) {
     minDist = currentDist;
   }
 
   // now: MODULE system. Need to transform locpo to MODULE PC
   //  transform is STRIP PC -> STRIP XY -> MODULE XY -> MODULE PC
-  Vector2 locpoStripXY(locpo_rotated[0] * std::cos(locpo_rotated[1]),
-                       locpo_rotated[0] * std::sin(locpo_rotated[1]));
+  Vector2 locpoStripXY(
+      locpo_rotated[eBoundLoc0] * std::cos(locpo_rotated[eBoundLoc1]),
+      locpo_rotated[eBoundLoc0] * std::sin(locpo_rotated[eBoundLoc1]));
   Vector2 locpoModulePC = stripXYToModulePC(locpoStripXY);
 
   // now check edges in MODULE PC (inner and outer circle) assuming Mahalanobis
   // distances are of same unit if covariance is correctly transformed
   currentClosest = closestOnSegment(m_inLeftModulePC, m_inRightModulePC,
-                                    locpoModulePC, metricModulePC);
-  currentDist = squaredNorm(locpoModulePC - currentClosest, metricModulePC);
+                                    locpoModulePC, weightModulePC);
+  currentDist = squaredNorm(locpoModulePC - currentClosest, weightModulePC);
   if (currentDist < minDist) {
     minDist = currentDist;
   }
 
   currentClosest = closestOnSegment(m_outLeftModulePC, m_outRightModulePC,
-                                    locpoModulePC, metricModulePC);
-  currentDist = squaredNorm(locpoModulePC - currentClosest, metricModulePC);
+                                    locpoModulePC, weightModulePC);
+  currentDist = squaredNorm(locpoModulePC - currentClosest, weightModulePC);
   if (currentDist < minDist) {
     minDist = currentDist;
   }

Core/src/Surfaces/BoundaryTolerance.cpp

@@ -60,38 +60,59 @@ bool BoundaryTolerance::hasChi2Bound() const {
   return holdsVariant<Chi2Bound>();
 }
 
-bool BoundaryTolerance::hasTolerance() const {
+BoundaryTolerance::ToleranceMode BoundaryTolerance::toleranceMode() const {
+  using enum ToleranceMode;
   if (isInfinite()) {
-    return true;
+    return Extend;
   }
 
   if (isNone()) {
-    return false;
+    return None;
   }
 
   if (const auto* absoluteBound = getVariantPtr<AbsoluteBound>();
       absoluteBound != nullptr) {
-    return absoluteBound->tolerance0 != 0. || absoluteBound->tolerance1 != 0.;
+    if (absoluteBound->tolerance0 == 0. && absoluteBound->tolerance1 == 0.) {
+      return None;
+    }
+
+    return Extend;
   }
 
   if (const auto* absoluteCartesian = getVariantPtr<AbsoluteCartesian>();
       absoluteCartesian != nullptr) {
-    return absoluteCartesian->tolerance0 != 0. ||
-           absoluteCartesian->tolerance1 != 0.;
+    if (absoluteCartesian->tolerance0 == 0. &&
+        absoluteCartesian->tolerance1 == 0.) {
+      return None;
+    }
+
+    return Extend;
   }
 
   if (const auto* absoluteEuclidean = getVariantPtr<AbsoluteEuclidean>();
       absoluteEuclidean != nullptr) {
-    return absoluteEuclidean->tolerance != 0.;
+    if (absoluteEuclidean->tolerance == 0.) {
+      return None;
+    } else if (absoluteEuclidean->tolerance > 0.) {
+      return Extend;
+    } else {
+      return Shrink;
+    }
   }
 
   if (const auto* chi2Bound = getVariantPtr<Chi2Bound>();
       chi2Bound != nullptr) {
-    return chi2Bound->maxChi2 != 0.;
+    if (chi2Bound->maxChi2 == 0.) {
+      return None;
+    } else if (chi2Bound->maxChi2 >= 0.) {
+      return Extend;
+    } else {
+      return Shrink;
+    }
   }
 
   assert(false && "Unsupported tolerance type");
-  return false;
+  return None;
 }
 
 BoundaryTolerance::AbsoluteBound BoundaryTolerance::asAbsoluteBound(
@@ -148,9 +169,13 @@ bool BoundaryTolerance::isTolerated(
 
   if (const auto* chi2Bound = getVariantPtr<Chi2Bound>();
       chi2Bound != nullptr) {
-    double chi2 = distance.transpose() * chi2Bound->weight * distance;
     // Mahalanobis distances mean is 2 in 2-dim. cut is 1-d sigma.
-    return chi2 <= 2 * chi2Bound->maxChi2;
+    double chi2 = distance.transpose() * chi2Bound->weight * distance;
+    if (chi2Bound->maxChi2 < 0) {
+      return chi2 > 2 * std::abs(chi2Bound->maxChi2);
+    } else {
+      return chi2 <= 2 * chi2Bound->maxChi2;
+    }
   }
 
   bool isCartesian = !jacobianOpt.has_value();
@@ -170,7 +195,11 @@ bool BoundaryTolerance::isTolerated(
 
   if (const auto* absoluteEuclidean = getVariantPtr<AbsoluteEuclidean>();
       absoluteEuclidean != nullptr) {
-    return cartesianDistance.norm() <= absoluteEuclidean->tolerance;
+    if (absoluteEuclidean->tolerance < 0) {
+      return cartesianDistance.norm() > std::abs(absoluteEuclidean->tolerance);
+    } else {
+      return cartesianDistance.norm() <= absoluteEuclidean->tolerance;
+    }
   }
 
   throw std::logic_error("Unsupported tolerance type");

Core/src/Surfaces/RectangleBounds.cpp

@@ -53,9 +53,15 @@ bool RectangleBounds::inside(const Vector2& lposition) const {
 Vector2 RectangleBounds::closestPoint(
     const Vector2& lposition,
     const std::optional<SquareMatrix2>& metric) const {
+  // If no metric is provided we can use a shortcut for the rectangle
+  if (!metric.has_value()) {
+    return detail::VerticesHelper::computeEuclideanClosestPointOnRectangle(
+        lposition, m_min, m_max);
+  }
+
+  // Otherwise we need to compute the closest point on the polygon
   std::array<Vector2, 4> vertices = {
       {m_min, {m_max[0], m_min[1]}, m_max, {m_min[0], m_max[1]}}};
-
   return detail::VerticesHelper::computeClosestPointOnPolygon(
       lposition, vertices, metric.value_or(SquareMatrix2::Identity()));
 }