From b4817ab0b6effc2179aaaf23d02a90ac60d12811 Mon Sep 17 00:00:00 2001 From: Paul Gessinger Date: Thu, 19 Dec 2024 11:54:49 +0100 Subject: [PATCH] add tolerance calculation including negative for rect, poly and annulus --- .../Acts/Surfaces/BoundaryTolerance.hpp | 14 +- .../Surfaces/detail/BoundaryCheckHelper.hpp | 28 +- Core/src/Surfaces/AnnulusBounds.cpp | 323 +++++++++++------- Core/src/Surfaces/BoundaryTolerance.cpp | 45 +-- .../Core/Surfaces/AnnulusBoundsTests.cpp | 155 +++++++++ .../Core/Surfaces/BoundaryToleranceTests.cpp | 167 ++++++--- 6 files changed, 517 insertions(+), 215 deletions(-) diff --git a/Core/include/Acts/Surfaces/BoundaryTolerance.hpp b/Core/include/Acts/Surfaces/BoundaryTolerance.hpp index fadd65f1a85..bc852d1e3ef 100644 --- a/Core/include/Acts/Surfaces/BoundaryTolerance.hpp +++ b/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 diff --git a/Core/include/Acts/Surfaces/detail/BoundaryCheckHelper.hpp b/Core/include/Acts/Surfaces/detail/BoundaryCheckHelper.hpp index 6df774e3590..a08204e1be0 100644 --- a/Core/include/Acts/Surfaces/detail/BoundaryCheckHelper.hpp +++ b/Core/include/Acts/Surfaces/detail/BoundaryCheckHelper.hpp @@ -69,15 +69,12 @@ inline bool insideAlignedBox(const Vector2& lowerLeft, point, vertices, metric); } - std::cout << "closestPoint: " << closestPoint.transpose() << std::endl; - std::cout << "point: " << point.transpose() << std::endl; Vector2 distance = closestPoint - point; - std::cout << "distance: " << distance.transpose() << std::endl; if (mode == Extend) { return tolerance.isTolerated(distance, jacobianOpt); } else { - return tolerance.isTolerated(-distance, jacobianOpt) && insideRectangle; + return tolerance.isTolerated(distance, jacobianOpt) && insideRectangle; } } @@ -93,23 +90,27 @@ inline bool insidePolygon(std::span vertices, const BoundaryTolerance& tolerance, const Vector2& point, const std::optional& jacobianOpt) { + using enum BoundaryTolerance::ToleranceMode; if (tolerance.isInfinite()) { // The null boundary check always succeeds return true; } - if (detail::VerticesHelper::isInsidePolygon(point, vertices)) { - // If the point falls inside the polygon, the check always succeeds - return true; - } + BoundaryTolerance::ToleranceMode mode = tolerance.toleranceMode(); + bool insidePolygon = detail::VerticesHelper::isInsidePolygon(point, vertices); - if (tolerance.toleranceMode() == BoundaryTolerance::ToleranceMode::None) { + if (mode == None) { + // If the point falls inside the polygon, the check always succeeds // 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. - return false; + return insidePolygon; + } + + if (mode == Extend && insidePolygon) { + return true; } // TODO: When tolerance is not 0, we could also avoid this computation in @@ -126,7 +127,12 @@ inline bool insidePolygon(std::span vertices, Vector2 distance = closestPoint - point; - return tolerance.isTolerated(distance, jacobianOpt); + if (mode == Extend) { + return tolerance.isTolerated(distance, jacobianOpt); + } else { + // @TODO: Check sign + return tolerance.isTolerated(-distance, jacobianOpt) && insidePolygon; + } } } // namespace Acts::detail diff --git a/Core/src/Surfaces/AnnulusBounds.cpp b/Core/src/Surfaces/AnnulusBounds.cpp index 2075e3cb59d..956a644ba5f 100644 --- a/Core/src/Surfaces/AnnulusBounds.cpp +++ b/Core/src/Surfaces/AnnulusBounds.cpp @@ -207,6 +207,7 @@ bool AnnulusBounds::inside(const Vector2& lposition, double tolR, bool AnnulusBounds::inside(const Vector2& lposition, const BoundaryTolerance& boundaryTolerance) const { + using enum BoundaryTolerance::ToleranceMode; if (boundaryTolerance.isInfinite()) { return true; } @@ -221,154 +222,218 @@ bool AnnulusBounds::inside(const Vector2& lposition, absoluteBound->tolerance1); } - // first check if inside. We don't need to look into the covariance if inside - if (inside(lposition, 0., 0.)) { - return true; - } + bool insideStrict = inside(lposition, 0., 0.); - if (!boundaryTolerance.hasChi2Bound()) { - throw std::logic_error("not implemented"); + BoundaryTolerance::ToleranceMode mode = boundaryTolerance.toleranceMode(); + if (mode == None) { + // first check if inside if we're in None tolerance mode. We don't need to + // look into the covariance if inside + return insideStrict; } - const auto& boundaryToleranceChi2 = boundaryTolerance.asChi2Bound(); + if (mode == Extend && insideStrict) { + return true; + } // locpo is PC in STRIP SYSTEM // we need to rotate the locpo Vector2 locpo_rotated = m_rotationStripPC * lposition; - // covariance is given in STRIP SYSTEM in PC we need to convert the covariance - // to the MODULE SYSTEM in PC via jacobian. The following transforms into - // STRIP XY, does the shift into MODULE XY, and then transforms into MODULE PC + // covariance is given in STRIP SYSTEM in PC we need to convert the + // covariance to the MODULE SYSTEM in PC via jacobian. The following + // transforms into STRIP XY, does the shift into MODULE XY, and then + // transforms into MODULE PC double dphi = get(eAveragePhi); double phi_strip = locpo_rotated[1]; double r_strip = locpo_rotated[0]; double O_x = m_shiftXY[0]; double O_y = m_shiftXY[1]; - // For a transformation from cartesian into polar coordinates - // - // [ _________ ] - // [ / 2 2 ] - // [ \/ x + y ] - // [ r' ] [ ] - // v = [ ] = [ / y \] - // [phi'] [2*atan|----------------|] - // [ | _________|] - // [ | / 2 2 |] - // [ \x + \/ x + y /] - // - // Where x, y are polar coordinates that can be rotated by dPhi - // - // [x] [O_x + r*cos(dPhi - phi)] - // [ ] = [ ] - // [y] [O_y - r*sin(dPhi - phi)] - // - // The general jacobian is: - // - // [d d ] - // [--(f_x) --(f_x)] - // [dx dy ] - // Jgen = [ ] - // [d d ] - // [--(f_y) --(f_y)] - // [dx dy ] - // - // which means in this case: - // - // [ d d ] - // [ ----------(rMod) ---------(rMod) ] - // [ dr_{strip} dphiStrip ] - // J = [ ] - // [ d d ] - // [----------(phiMod) ---------(phiMod)] - // [dr_{strip} dphiStrip ] - // - // Performing the derivative one gets: - // - // [B*O_x + C*O_y + rStrip rStrip*(B*O_y + O_x*sin(dPhi - phiStrip))] - // [---------------------- -----------------------------------------] - // [ ___ ___ ] - // [ \/ A \/ A ] - // J = [ ] - // [ -(B*O_y - C*O_x) rStrip*(B*O_x + C*O_y + rStrip) ] - // [ ----------------- ------------------------------- ] - // [ A A ] - // - // where - // 2 2 - // A = O_x + 2*O_x*rStrip*cos(dPhi - phiStrip) + O_y - // 2 - // - 2*O_y*rStrip*sin(dPhi - phiStrip) + rStrip - // B = cos(dPhi - phiStrip) - // C = -sin(dPhi - phiStrip) - - double cosDPhiPhiStrip = std::cos(dphi - phi_strip); - double sinDPhiPhiStrip = std::sin(dphi - phi_strip); - - double A = O_x * O_x + 2 * O_x * r_strip * cosDPhiPhiStrip + O_y * O_y - - 2 * O_y * r_strip * sinDPhiPhiStrip + r_strip * r_strip; - double sqrtA = std::sqrt(A); - - double B = cosDPhiPhiStrip; - double C = -sinDPhiPhiStrip; - SquareMatrix2 jacobianStripPCToModulePC; - jacobianStripPCToModulePC(0, 0) = (B * O_x + C * O_y + r_strip) / sqrtA; - jacobianStripPCToModulePC(0, 1) = - r_strip * (B * O_y + O_x * sinDPhiPhiStrip) / sqrtA; - 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 inverse covariance as weights - const auto& weightStripPC = boundaryToleranceChi2.weight; - auto weightModulePC = jacobianStripPCToModulePC.transpose() * weightStripPC * - jacobianStripPCToModulePC; - - double minDist = std::numeric_limits::max(); - - Vector2 currentClosest; - double currentDist = 0; - - // do projection in STRIP PC - - // first: STRIP system. locpo is in STRIP PC already - currentClosest = closestOnSegment(m_inLeftStripPC, m_outLeftStripPC, - locpo_rotated, weightStripPC); - currentDist = squaredNorm(locpo_rotated - currentClosest, weightStripPC); - minDist = currentDist; - - currentClosest = closestOnSegment(m_inRightStripPC, m_outRightStripPC, - locpo_rotated, weightStripPC); - currentDist = squaredNorm(locpo_rotated - currentClosest, weightStripPC); - if (currentDist < minDist) { + auto closestPointDistanceBound = [&](const SquareMatrix2& weight) { + // For a transformation from cartesian into polar coordinates + // + // [ _________ ] + // [ / 2 2 ] + // [ \/ x + y ] + // [ r' ] [ ] + // v = [ ] = [ / y \] + // [phi'] [2*atan|----------------|] + // [ | _________|] + // [ | / 2 2 |] + // [ \x + \/ x + y /] + // + // Where x, y are polar coordinates that can be rotated by dPhi + // + // [x] [O_x + r*cos(dPhi - phi)] + // [ ] = [ ] + // [y] [O_y - r*sin(dPhi - phi)] + // + // The general jacobian is: + // + // [d d ] + // [--(f_x) --(f_x)] + // [dx dy ] + // Jgen = [ ] + // [d d ] + // [--(f_y) --(f_y)] + // [dx dy ] + // + // which means in this case: + // + // [ d d ] + // [ ----------(rMod) ---------(rMod) ] + // [ dr_{strip} dphiStrip ] + // J = [ ] + // [ d d ] + // [----------(phiMod) ---------(phiMod)] + // [dr_{strip} dphiStrip ] + // + // Performing the derivative one gets: + // + // [B*O_x + C*O_y + rStrip rStrip*(B*O_y + O_x*sin(dPhi - phiStrip))] + // [---------------------- -----------------------------------------] + // [ ___ ___ ] + // [ \/ A \/ A ] + // J = [ ] + // [ -(B*O_y - C*O_x) rStrip*(B*O_x + C*O_y + rStrip) ] + // [ ----------------- ------------------------------- ] + // [ A A ] + // + // where + // 2 2 + // A = O_x + 2*O_x*rStrip*cos(dPhi - phiStrip) + O_y + // 2 + // - 2*O_y*rStrip*sin(dPhi - phiStrip) + rStrip + // B = cos(dPhi - phiStrip) + // C = -sin(dPhi - phiStrip) + + double cosDPhiPhiStrip = std::cos(dphi - phi_strip); + double sinDPhiPhiStrip = std::sin(dphi - phi_strip); + + double A = O_x * O_x + 2 * O_x * r_strip * cosDPhiPhiStrip + O_y * O_y - + 2 * O_y * r_strip * sinDPhiPhiStrip + r_strip * r_strip; + double sqrtA = std::sqrt(A); + + double B = cosDPhiPhiStrip; + double C = -sinDPhiPhiStrip; + SquareMatrix2 jacobianStripPCToModulePC; + jacobianStripPCToModulePC(0, 0) = (B * O_x + C * O_y + r_strip) / sqrtA; + jacobianStripPCToModulePC(0, 1) = + r_strip * (B * O_y + O_x * sinDPhiPhiStrip) / sqrtA; + 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 inverse covariance as weights + const auto& weightStripPC = weight; + auto weightModulePC = jacobianStripPCToModulePC.transpose() * + weightStripPC * jacobianStripPCToModulePC; + + double minDist = std::numeric_limits::max(); + Vector2 delta; + + Vector2 currentClosest; + double currentDist = 0; + + // do projection in STRIP PC + + // first: STRIP system. locpo is in STRIP PC already + currentClosest = closestOnSegment(m_inLeftStripPC, m_outLeftStripPC, + locpo_rotated, weightStripPC); + currentDist = squaredNorm(locpo_rotated - currentClosest, weightStripPC); minDist = currentDist; - } + delta = locpo_rotated - currentClosest; + currentClosest = closestOnSegment(m_inRightStripPC, m_outRightStripPC, + locpo_rotated, weightStripPC); + currentDist = squaredNorm(locpo_rotated - currentClosest, weightStripPC); + if (currentDist < minDist) { + minDist = currentDist; + delta = locpo_rotated - currentClosest; + } - // 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 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, weightModulePC); - currentDist = squaredNorm(locpoModulePC - currentClosest, weightModulePC); - 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[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, weightModulePC); + currentDist = squaredNorm(locpoModulePC - currentClosest, weightModulePC); + if (currentDist < minDist) { + minDist = currentDist; + delta = jacobianStripPCToModulePC.inverse() * + (currentClosest - locpoModulePC); + } - currentClosest = closestOnSegment(m_outLeftModulePC, m_outRightModulePC, - locpoModulePC, weightModulePC); - currentDist = squaredNorm(locpoModulePC - currentClosest, weightModulePC); - if (currentDist < minDist) { - minDist = currentDist; + currentClosest = closestOnSegment(m_outLeftModulePC, m_outRightModulePC, + locpoModulePC, weightModulePC); + currentDist = squaredNorm(locpoModulePC - currentClosest, weightModulePC); + if (currentDist < minDist) { + minDist = currentDist; + delta = jacobianStripPCToModulePC.inverse() * + (currentClosest - locpoModulePC); + } + + return std::tuple{delta, minDist}; + }; + + if (boundaryTolerance.hasChi2Bound()) { + const auto& boundaryToleranceChi2 = boundaryTolerance.asChi2Bound(); + + // Calculate minDist based on weight from the boundary tolerance object. + // That weight matrix is in STRIP PC + auto [delta, minDist] = + closestPointDistanceBound(boundaryToleranceChi2.weight); + + // compare resulting Mahalanobis distance to configured "number of sigmas" + // we square it b/c we never took the square root of the distance + if (mode == Extend) { + return minDist < + boundaryToleranceChi2.maxChi2 * boundaryToleranceChi2.maxChi2; + } else if (mode == Shrink) { + return minDist > boundaryToleranceChi2.maxChi2 * + boundaryToleranceChi2.maxChi2 && + insideStrict; + } + } else if (boundaryTolerance.hasAbsoluteEuclidean()) { + const auto& boundaryToleranceAbsoluteEuclidean = + boundaryTolerance.asAbsoluteEuclidean(); + + SquareMatrix2 jacobianPCToXY; + jacobianPCToXY(0, 0) = std::cos(phi_strip); + jacobianPCToXY(0, 1) = -r_strip * std::sin(phi_strip); + jacobianPCToXY(1, 0) = std::sin(phi_strip); + jacobianPCToXY(1, 1) = r_strip * std::cos(phi_strip); + + // This is J.T * J but we can also calculate it directly + SquareMatrix2 weightStripPC; + weightStripPC(0, 0) = 1; + weightStripPC(0, 1) = 0; + weightStripPC(1, 0) = 0; + weightStripPC(1, 1) = r_strip * r_strip; + + auto [delta, minDist] = closestPointDistanceBound(weightStripPC); + + Vector2 cartesianDistance = jacobianPCToXY * delta; + + if (mode == Extend) { + return cartesianDistance.squaredNorm() < + boundaryToleranceAbsoluteEuclidean.tolerance * + boundaryToleranceAbsoluteEuclidean.tolerance; + } else if (mode == Shrink) { + return cartesianDistance.squaredNorm() > + boundaryToleranceAbsoluteEuclidean.tolerance * + boundaryToleranceAbsoluteEuclidean.tolerance && + insideStrict; + } } - // compare resulting Mahalanobis distance to configured "number of sigmas" we - // square it b/c we never took the square root of the distance - return minDist < - boundaryToleranceChi2.maxChi2 * boundaryToleranceChi2.maxChi2; + throw std::logic_error("not implemented"); } Vector2 AnnulusBounds::stripXYToModulePC(const Vector2& vStripXY) const { diff --git a/Core/src/Surfaces/BoundaryTolerance.cpp b/Core/src/Surfaces/BoundaryTolerance.cpp index 1c21bc408de..1be955d0a48 100644 --- a/Core/src/Surfaces/BoundaryTolerance.cpp +++ b/Core/src/Surfaces/BoundaryTolerance.cpp @@ -128,7 +128,7 @@ BoundaryTolerance::ToleranceMode BoundaryTolerance::toleranceMode() const { chi2Bound != nullptr) { if (chi2Bound->maxChi2 == 0.) { return None; - } else if (chi2Bound->maxChi2 > 0.) { + } else if (chi2Bound->maxChi2 >= 0.) { return Extend; } else { return Shrink; @@ -187,28 +187,19 @@ bool BoundaryTolerance::isTolerated( if (const auto* absoluteBound = getVariantPtr(); absoluteBound != nullptr) { - bool tol0 = false; - if (absoluteBound->tolerance0 < 0) { - tol0 = std::abs(distance[0]) > std::abs(absoluteBound->tolerance0); - } else { - tol0 = std::abs(distance[0]) <= absoluteBound->tolerance0; - } - - bool tol1 = false; - if (absoluteBound->tolerance1 < 0) { - tol1 = std::abs(distance[1]) > std::abs(absoluteBound->tolerance1); - } else { - tol1 = std::abs(distance[1]) <= absoluteBound->tolerance1; - } - - return tol0 && tol1; + return std::abs(distance[0]) <= absoluteBound->tolerance0 && + std::abs(distance[1]) <= absoluteBound->tolerance1; } if (const auto* chi2Bound = getVariantPtr(); 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(); @@ -222,22 +213,8 @@ bool BoundaryTolerance::isTolerated( if (const auto* absoluteCartesian = getVariantPtr(); absoluteCartesian != nullptr) { - bool tol0 = false; - - if (absoluteCartesian->tolerance0 < 0) { - tol0 = cartesianDistance[0] > std::abs(absoluteCartesian->tolerance0); - } else { - tol0 = std::abs(cartesianDistance[0]) <= absoluteCartesian->tolerance0; - } - - bool tol1 = false; - if (absoluteCartesian->tolerance1 < 0) { - tol1 = cartesianDistance[1] > std::abs(absoluteCartesian->tolerance1); - } else { - tol1 = std::abs(cartesianDistance[1]) <= absoluteCartesian->tolerance1; - } - - return tol0 && tol1; + return std::abs(cartesianDistance[0]) <= absoluteCartesian->tolerance0 && + std::abs(cartesianDistance[1]) <= absoluteCartesian->tolerance1; } if (const auto* absoluteEuclidean = getVariantPtr(); diff --git a/Tests/UnitTests/Core/Surfaces/AnnulusBoundsTests.cpp b/Tests/UnitTests/Core/Surfaces/AnnulusBoundsTests.cpp index f2a060167c7..7857a5a3d33 100644 --- a/Tests/UnitTests/Core/Surfaces/AnnulusBoundsTests.cpp +++ b/Tests/UnitTests/Core/Surfaces/AnnulusBoundsTests.cpp @@ -141,6 +141,161 @@ BOOST_AUTO_TEST_CASE(AnnulusBoundsVertices) { BOOST_CHECK_EQUAL(vertices.size(), 14u); } +BOOST_AUTO_TEST_CASE(AnnulusBoundsNegativeTolerance) { + using namespace Acts::UnitLiterals; + AnnulusBounds aBounds(minRadius, maxRadius, minPhi, maxPhi, offset); + double phiAverage = (minPhi + maxPhi) / 2; + + auto check = [&](const BoundaryTolerance& tolerance, const Vector2& point) { + Vector2 pointAverage(point[0], phiAverage + point[1]); + + std::cout << "[" << pointAverage[0] << ", " << pointAverage[1] << "]," + << std::endl; + return aBounds.inside(pointAverage, tolerance); + }; + + double midRadius = (minRadius + maxRadius) / 2; + double hlPhi = (maxPhi - minPhi) / 2; + + { + auto tolerance = BoundaryTolerance::AbsoluteEuclidean(1); + + // Test points near radial boundaries + BOOST_CHECK(!check(tolerance, {minRadius - 1.5, 0})); + BOOST_CHECK(check(tolerance, {minRadius - 0.1, 0})); + BOOST_CHECK(check(tolerance, {minRadius + 0.4, 0})); + BOOST_CHECK(check(tolerance, {minRadius + 0.5, 0})); + BOOST_CHECK(check(tolerance, {minRadius + 1.5, 0})); + + BOOST_CHECK(check(tolerance, {maxRadius - 1.5, 0})); + BOOST_CHECK(check(tolerance, {maxRadius - 0.1, 0})); + BOOST_CHECK(check(tolerance, {maxRadius + 0.3, 0})); + BOOST_CHECK(check(tolerance, {maxRadius + 0.55, 0})); + BOOST_CHECK(check(tolerance, {maxRadius + 1.2, 0})); + BOOST_CHECK(!check(tolerance, {maxRadius + 1.7, 0})); + + // Check points near axial boundaries + BOOST_CHECK(!check(tolerance, {midRadius, -hlPhi * 1.5})); + BOOST_CHECK(check(tolerance, {midRadius, -hlPhi * 1.1})); + BOOST_CHECK(check(tolerance, {midRadius, -hlPhi * 0.8})); + BOOST_CHECK(check(tolerance, {midRadius, -hlPhi * 0.5})); + + BOOST_CHECK(check(tolerance, {midRadius, hlPhi * 0.5})); + BOOST_CHECK(check(tolerance, {midRadius, hlPhi * 0.8})); + BOOST_CHECK(check(tolerance, {midRadius, hlPhi * 1.1})); + BOOST_CHECK(!check(tolerance, {midRadius, hlPhi * 1.5})); + } + + { + auto tolerance = BoundaryTolerance::AbsoluteEuclidean(-1); + + // Test points near radial boundaries + BOOST_CHECK(!check(tolerance, {minRadius - 1.5, 0})); + BOOST_CHECK(!check(tolerance, {minRadius - 0.1, 0})); + BOOST_CHECK(!check(tolerance, {minRadius + 0.4, 0})); + BOOST_CHECK(!check(tolerance, {minRadius + 0.5, 0})); + BOOST_CHECK(check(tolerance, {minRadius + 1.5, 0})); + + BOOST_CHECK(check(tolerance, {maxRadius - 1.5, 0})); + BOOST_CHECK(!check(tolerance, {maxRadius - 0.1, 0})); + BOOST_CHECK(!check(tolerance, {maxRadius + 0.3, 0})); + BOOST_CHECK(!check(tolerance, {maxRadius + 0.55, 0})); + BOOST_CHECK(!check(tolerance, {maxRadius + 1.2, 0})); + BOOST_CHECK(!check(tolerance, {maxRadius + 1.7, 0})); + + // Check points near axial boundaries + BOOST_CHECK(!check(tolerance, {midRadius, -hlPhi * 1.5})); + BOOST_CHECK(!check(tolerance, {midRadius, -hlPhi * 1.1})); + BOOST_CHECK(!check(tolerance, {midRadius, -hlPhi * 0.8})); + BOOST_CHECK(check(tolerance, {midRadius, -hlPhi * 0.5})); + + BOOST_CHECK(check(tolerance, {midRadius, hlPhi * 0.5})); + BOOST_CHECK(!check(tolerance, {midRadius, hlPhi * 0.8})); + BOOST_CHECK(!check(tolerance, {midRadius, hlPhi * 1.1})); + BOOST_CHECK(!check(tolerance, {midRadius, hlPhi * 1.5})); + } + + { + auto tolerance = + BoundaryTolerance::Chi2Bound(SquareMatrix2::Identity(), 0.1); + + // Test points near radial boundaries + BOOST_CHECK(!check(tolerance, {minRadius - 1.5, 0})); + BOOST_CHECK(check(tolerance, {minRadius - 0.1, 0})); + BOOST_CHECK(check(tolerance, {minRadius + 0.4, 0})); + BOOST_CHECK(check(tolerance, {minRadius + 0.5, 0})); + BOOST_CHECK(check(tolerance, {minRadius + 1.5, 0})); + + BOOST_CHECK(check(tolerance, {maxRadius - 1.5, 0})); + BOOST_CHECK(check(tolerance, {maxRadius - 0.1, 0})); + BOOST_CHECK(check(tolerance, {maxRadius + 0.3, 0})); + BOOST_CHECK(check(tolerance, {maxRadius + 0.55, 0})); + BOOST_CHECK(!check(tolerance, {maxRadius + 1.2, 0})); + BOOST_CHECK(!check(tolerance, {maxRadius + 1.7, 0})); + + // Check points near axial boundaries + BOOST_CHECK(!check(tolerance, {midRadius, -hlPhi * 1.5})); + BOOST_CHECK(check(tolerance, {midRadius, -hlPhi * 1.1})); + BOOST_CHECK(check(tolerance, {midRadius, -hlPhi * 0.8})); + BOOST_CHECK(check(tolerance, {midRadius, -hlPhi * 0.5})); + + BOOST_CHECK(check(tolerance, {midRadius, hlPhi * 0.5})); + BOOST_CHECK(check(tolerance, {midRadius, hlPhi * 0.8})); + BOOST_CHECK(check(tolerance, {midRadius, hlPhi * 1.1})); + BOOST_CHECK(!check(tolerance, {midRadius, hlPhi * 1.5})); + // mc(tolerance); + // mc(BoundaryTolerance::Chi2Bound(SquareMatrix2::Identity(), 0.1)); + + // Test points near radial boundaries + // BOOST_CHECK(!check(tolerance, {minRadius - 0.2, 0})); + // BOOST_CHECK(!check(tolerance, {minRadius + 0.2, 0})); + // BOOST_CHECK(check(tolerance, {minRadius + 0.4, 0})); + // BOOST_CHECK(!check(tolerance, {maxRadius - 0.2, 0})); + // BOOST_CHECK(!check(tolerance, {maxRadius + 0.2, 0})); + // BOOST_CHECK(check(tolerance, {maxRadius - 0.4, 0})); + + // // Test points near angular boundaries + // BOOST_CHECK(!check(tolerance, {(minRadius + maxRadius) / 2, minPhi - + // 0.2})); BOOST_CHECK(!check(tolerance, {(minRadius + maxRadius) / 2, + // minPhi + 0.2})); BOOST_CHECK(check(tolerance, {(minRadius + maxRadius) / + // 2, minPhi + 0.4})); BOOST_CHECK(!check(tolerance, {(minRadius + + // maxRadius) / 2, maxPhi - 0.2})); BOOST_CHECK(!check(tolerance, + // {(minRadius + maxRadius) / 2, maxPhi + 0.2})); + // BOOST_CHECK(check(tolerance, {(minRadius + maxRadius) / 2, maxPhi - + // 0.4})); + } + + { + auto tolerance = + BoundaryTolerance::Chi2Bound(SquareMatrix2::Identity(), -0.1); + + // Test points near radial boundaries + BOOST_CHECK(!check(tolerance, {minRadius - 1.5, 0})); + BOOST_CHECK(!check(tolerance, {minRadius - 0.1, 0})); + BOOST_CHECK(!check(tolerance, {minRadius + 0.4, 0})); + BOOST_CHECK(check(tolerance, {minRadius + 0.5, 0})); + BOOST_CHECK(check(tolerance, {minRadius + 1.5, 0})); + + BOOST_CHECK(check(tolerance, {maxRadius - 1.5, 0})); + BOOST_CHECK(check(tolerance, {maxRadius - 0.1, 0})); + BOOST_CHECK(!check(tolerance, {maxRadius + 0.3, 0})); + BOOST_CHECK(!check(tolerance, {maxRadius + 0.55, 0})); + BOOST_CHECK(!check(tolerance, {maxRadius + 1.2, 0})); + BOOST_CHECK(!check(tolerance, {maxRadius + 1.7, 0})); + + // Check points near axial boundaries + BOOST_CHECK(!check(tolerance, {midRadius, -hlPhi * 1.5})); + BOOST_CHECK(!check(tolerance, {midRadius, -hlPhi * 1.1})); + BOOST_CHECK(!check(tolerance, {midRadius, -hlPhi * 0.8})); + BOOST_CHECK(check(tolerance, {midRadius, -hlPhi * 0.5})); + + BOOST_CHECK(check(tolerance, {midRadius, hlPhi * 0.5})); + BOOST_CHECK(!check(tolerance, {midRadius, hlPhi * 0.8})); + BOOST_CHECK(!check(tolerance, {midRadius, hlPhi * 1.1})); + BOOST_CHECK(!check(tolerance, {midRadius, hlPhi * 1.5})); + } +} + BOOST_AUTO_TEST_SUITE_END() } // namespace Acts::Test diff --git a/Tests/UnitTests/Core/Surfaces/BoundaryToleranceTests.cpp b/Tests/UnitTests/Core/Surfaces/BoundaryToleranceTests.cpp index 7d9d070985f..eb0b6159754 100644 --- a/Tests/UnitTests/Core/Surfaces/BoundaryToleranceTests.cpp +++ b/Tests/UnitTests/Core/Surfaces/BoundaryToleranceTests.cpp @@ -13,6 +13,7 @@ #include "Acts/Surfaces/BoundaryTolerance.hpp" #include "Acts/Surfaces/detail/BoundaryCheckHelper.hpp" #include "Acts/Tests/CommonHelpers/FloatComparisons.hpp" +#include "Acts/Utilities/Helpers.hpp" #include #include @@ -26,6 +27,64 @@ namespace Acts::Test { BOOST_AUTO_TEST_SUITE(Surfaces) +BOOST_AUTO_TEST_CASE(BoundaryToleranceConstructors) { + // Test AbsoluteBound constructor + { + // Valid positive tolerances + auto tolerance = BoundaryTolerance::AbsoluteBound(1.0, 2.0); + BOOST_CHECK_EQUAL(tolerance.tolerance0, 1.0); + BOOST_CHECK_EQUAL(tolerance.tolerance1, 2.0); + + // Negative tolerances should throw + BOOST_CHECK_THROW(BoundaryTolerance::AbsoluteBound(-1.0, 2.0), + std::invalid_argument); + BOOST_CHECK_THROW(BoundaryTolerance::AbsoluteBound(1.0, -2.0), + std::invalid_argument); + } + + // Test AbsoluteEuclidean constructor + { + // Valid positive tolerance + auto tolerance = BoundaryTolerance::AbsoluteEuclidean(1.0); + BOOST_CHECK_EQUAL(tolerance.tolerance, 1.0); + + // Valid negative tolerance + tolerance = BoundaryTolerance::AbsoluteEuclidean(-1.0); + BOOST_CHECK_EQUAL(tolerance.tolerance, -1.0); + } + + // Test AbsoluteCartesian constructor + { + // Valid positive tolerance + auto tolerance = BoundaryTolerance::AbsoluteCartesian(1.0, 2.0); + BOOST_CHECK_EQUAL(tolerance.tolerance0, 1.0); + BOOST_CHECK_EQUAL(tolerance.tolerance1, 2.0); + + // Negative tolerances should throw + BOOST_CHECK_THROW(BoundaryTolerance::AbsoluteCartesian(-1.0, 2.0), + std::invalid_argument); + BOOST_CHECK_THROW(BoundaryTolerance::AbsoluteCartesian(1.0, -2.0), + std::invalid_argument); + } + + // Test Chi2Bound constructor + { + SquareMatrix2 cov; + cov << 1, 0.5, 0.5, 2; + + // Valid positive chi2 bound + auto tolerance = BoundaryTolerance::Chi2Bound(cov, 3.0); + BOOST_CHECK_EQUAL(tolerance.maxChi2, 3.0); + + // Valid negative chi2 bound + tolerance = BoundaryTolerance::Chi2Bound(cov, -3.0); + BOOST_CHECK_EQUAL(tolerance.maxChi2, -3.0); + } + + // Test None constructor + BoundaryTolerance::None(); +} + // See: https://en.wikipedia.org/wiki/Bounding_volume // // Aligned box w/ simple check @@ -205,7 +264,7 @@ BOOST_AUTO_TEST_CASE(BoundaryCheckDifferentTolerances) { } } -BOOST_AUTO_TEST_CASE(BoundaryCheckNegativeTolerance) { +BOOST_AUTO_TEST_CASE(BoundaryCheckNegativeToleranceRect) { // Test points for boundary check with euclidean tolerance Vector2 ll(1, 1); Vector2 ur(3, 3); @@ -224,53 +283,83 @@ BOOST_AUTO_TEST_CASE(BoundaryCheckNegativeTolerance) { BOOST_CHECK(!check(tolerance, {2, 3.1})); BOOST_CHECK(!check(tolerance, {2, 2.8})); BOOST_CHECK(check(tolerance, {2, 2.7})); + + BOOST_CHECK(!check(tolerance, {0.8, 2})); + BOOST_CHECK(!check(tolerance, {1.2, 2})); + BOOST_CHECK(check(tolerance, {1.5, 2})); + BOOST_CHECK(!check(tolerance, {2, 0.8})); + BOOST_CHECK(!check(tolerance, {2, 1.2})); + BOOST_CHECK(check(tolerance, {2, 1.5})); } { - auto tolerance = BoundaryTolerance::AbsoluteBound(-0.25, -0.0); + auto tolerance = + BoundaryTolerance::Chi2Bound(SquareMatrix2::Identity(), -0.1); BOOST_CHECK(!check(tolerance, {2.8, 2})); BOOST_CHECK(!check(tolerance, {3.1, 2})); - BOOST_CHECK(check(tolerance, {2.7, 2})); + BOOST_CHECK(check(tolerance, {2.5, 2})); BOOST_CHECK(!check(tolerance, {2, 3.1})); - BOOST_CHECK(check(tolerance, {2, 2.8})); - BOOST_CHECK(check(tolerance, {2, 2.7})); + BOOST_CHECK(!check(tolerance, {2, 2.8})); + BOOST_CHECK(check(tolerance, {2, 2.5})); + + BOOST_CHECK(!check(tolerance, {0.8, 2})); + BOOST_CHECK(!check(tolerance, {1.4, 2})); + BOOST_CHECK(check(tolerance, {1.5, 2})); + BOOST_CHECK(!check(tolerance, {2, 0.8})); + BOOST_CHECK(!check(tolerance, {2, 1.4})); + BOOST_CHECK(check(tolerance, {2, 1.5})); + } +} + +BOOST_AUTO_TEST_CASE(BoundaryCheckNegativeToleranceTrap) { + Vector2 vertices[] = {{1.5, 1}, {2.5, 1}, {3, 3}, {1, 3}}; + + auto check = [&vertices](const BoundaryTolerance& tolerance, + const Vector2& point) { + return detail::insidePolygon(vertices, tolerance, point, std::nullopt); + }; + + { + auto tolerance = BoundaryTolerance::AbsoluteEuclidean(0.25); + // Axes + BOOST_CHECK(!check(tolerance, {3.1, 2})); + BOOST_CHECK(check(tolerance, {2.75, 2})); + BOOST_CHECK(check(tolerance, {2.5, 2})); + BOOST_CHECK(check(tolerance, {2.25, 2})); + BOOST_CHECK(check(tolerance, {2, 3.1})); + BOOST_CHECK(check(tolerance, {2, 2.75})); + BOOST_CHECK(check(tolerance, {2, 2.5})); + BOOST_CHECK(check(tolerance, {2, 2.25})); + BOOST_CHECK(check(tolerance, {2, 2})); + + // Corners + BOOST_CHECK(check(tolerance, {3.1, 3.2})); + BOOST_CHECK(check(tolerance, {0.9, 3.2})); + BOOST_CHECK(check(tolerance, {1.5, 0.8})); + BOOST_CHECK(check(tolerance, {2.5, 0.8})); } - // std::ofstream outFile("boundary_check_euclidean.csv"); - // outFile << "x,y,inone,ipos,ineg\n"; - // - // std::mt19937 gen(42); - // std::uniform_real_distribution<> dis(0.5, 3.5); - // - // for (int i = 0; i < 1000; i++) { - // double x = dis(gen); - // double y = dis(gen); - // Vector2 point(x, y); - // - // bool insideNone = - // detail::insideAlignedBox(ll, ur, none, point, std::nullopt); - // bool insidePos = detail::insideAlignedBox(ll, ur, pos, point, - // std::nullopt); bool insideNeg = detail::insideAlignedBox(ll, ur, neg, - // point, std::nullopt); - // - // outFile << x << "," << y << "," << insideNone << "," << insidePos << "," - // << insideNeg << "\n"; - // } - // - // outFile.close(); - // - // return; - // Verify some key points with known distances - // BOOST_CHECK( - // detail::insideAlignedBox(ll, ur, tolerance, {0, 0}, std::nullopt)); - // BOOST_CHECK(detail::insideAlignedBox(ll, ur, tolerance, {1.05, 0}, - // std::nullopt)); // Just within - // tolerance - // BOOST_CHECK(!detail::insideAlignedBox( - // ll, ur, tolerance, {1.2, 0}, std::nullopt)); // Just outside tolerance - // BOOST_CHECK( - // !detail::insideAlignedBox(ll, ur, tolerance, {2, 2}, std::nullopt)); + { + auto tolerance = BoundaryTolerance::AbsoluteEuclidean(-0.25); + // Axes + BOOST_CHECK(!check(tolerance, {3.0, 2})); + BOOST_CHECK(!check(tolerance, {2.5, 2})); + BOOST_CHECK(check(tolerance, {2.25, 2})); + BOOST_CHECK(!check(tolerance, {2, 3.1})); + BOOST_CHECK(!check(tolerance, {2, 2.9})); + BOOST_CHECK(check(tolerance, {2, 2.7})); + + // Corners + BOOST_CHECK(!check(tolerance, {2.7, 2.9})); + BOOST_CHECK(check(tolerance, {2.4, 2.6})); + BOOST_CHECK(!check(tolerance, {1.3, 2.9})); + BOOST_CHECK(check(tolerance, {1.6, 2.6})); + BOOST_CHECK(!check(tolerance, {2.4, 1.1})); + BOOST_CHECK(check(tolerance, {1.75, 1.4})); + BOOST_CHECK(!check(tolerance, {1.6, 1.1})); + BOOST_CHECK(check(tolerance, {2.25, 1.4})); + } } BOOST_AUTO_TEST_SUITE_END()