diff --git a/Detectors/CADSupport/src/BoundedSurface.h b/Detectors/CADSupport/src/BoundedSurface.h index 20c21716c0483..5ebb6679cef58 100644 --- a/Detectors/CADSupport/src/BoundedSurface.h +++ b/Detectors/CADSupport/src/BoundedSurface.h @@ -2249,6 +2249,97 @@ inline std::vector sampleCurveWireByU(const CurveWire& wire, int segmentsP return samples; } +/// Ear-clipping of a simple parametric wire that clips the ear with the shortest 3D diagonal first, so a strip is not fanned. +inline std::vector> triangulateWireByShortestChord(const SurfaceWire& wire, + const std::vector& points) +{ + std::vector remainingIndices; + remainingIndices.reserve(wire.vertices.size()); + if (wire.signedArea() >= 0.) { + for (size_t vertexIndex = 0; vertexIndex < wire.vertices.size(); ++vertexIndex) { + remainingIndices.push_back(static_cast(vertexIndex)); + } + } else { + for (size_t reverseIndex = wire.vertices.size(); reverseIndex > 0; --reverseIndex) { + remainingIndices.push_back(static_cast(reverseIndex - 1)); + } + } + + // the 3D diagonal of the ear at a position, or infinity when it is not an ear (same tests as triangulateSimpleWire) + const double notAnEar = std::numeric_limits::infinity(); + auto earChord = [&](size_t indexPosition) { + const int previousIndex = remainingIndices[(indexPosition + remainingIndices.size() - 1) % remainingIndices.size()]; + const int currentIndex = remainingIndices[indexPosition]; + const int nextIndex = remainingIndices[(indexPosition + 1) % remainingIndices.size()]; + const auto& previousVertex = wire.vertices[previousIndex]; + const auto& currentVertex = wire.vertices[currentIndex]; + const auto& nextVertex = wire.vertices[nextIndex]; + if (cross2D(currentVertex - previousVertex, nextVertex - currentVertex) <= kTolerance) { + return notAnEar; + } + for (int candidateIndex : remainingIndices) { + if (candidateIndex == previousIndex || candidateIndex == currentIndex || candidateIndex == nextIndex) { + continue; + } + if (pointInTriangle(wire.vertices[candidateIndex], previousVertex, currentVertex, nextVertex)) { + return notAnEar; + } + } + return norm(points[nextIndex] - points[previousIndex]); + }; + + std::vector chords(remainingIndices.size()); + for (size_t indexPosition = 0; indexPosition < remainingIndices.size(); ++indexPosition) { + chords[indexPosition] = earChord(indexPosition); + } + + std::vector> triangles; + while (remainingIndices.size() > 3) { + // the shortest diagonal wins; ties go to the lowest position + size_t bestPosition = remainingIndices.size(); + double bestChord = notAnEar; + for (size_t indexPosition = 0; indexPosition < remainingIndices.size(); ++indexPosition) { + if (chords[indexPosition] < bestChord) { + bestChord = chords[indexPosition]; + bestPosition = indexPosition; + } + } + if (bestPosition == remainingIndices.size()) { + return triangulateSimpleWire(wire); // near-collinear runs can block every ear in this order + } + + const size_t count = remainingIndices.size(); + const int previousIndex = remainingIndices[(bestPosition + count - 1) % count]; + const int clippedIndex = remainingIndices[bestPosition]; + const int nextIndex = remainingIndices[(bestPosition + 1) % count]; + triangles.push_back({previousIndex, clippedIndex, nextIndex}); + remainingIndices.erase(remainingIndices.begin() + bestPosition); + chords.erase(chords.begin() + bestPosition); + + // the two neighbours change shape; any other non-ear can only change if the clipped vertex lay in its triangle + const size_t newCount = remainingIndices.size(); + const size_t beforePosition = (bestPosition + newCount - 1) % newCount; + const size_t afterPosition = bestPosition % newCount; + for (size_t indexPosition = 0; indexPosition < newCount; ++indexPosition) { + if (indexPosition == beforePosition || indexPosition == afterPosition) { + chords[indexPosition] = earChord(indexPosition); + } else if (chords[indexPosition] == notAnEar) { + const auto& previousVertex = wire.vertices[remainingIndices[(indexPosition + newCount - 1) % newCount]]; + const auto& currentVertex = wire.vertices[remainingIndices[indexPosition]]; + const auto& nextVertex = wire.vertices[remainingIndices[(indexPosition + 1) % newCount]]; + if (pointInTriangle(wire.vertices[clippedIndex], previousVertex, currentVertex, nextVertex)) { + chords[indexPosition] = earChord(indexPosition); + } + } + } + } + + if (remainingIndices.size() == 3) { + triangles.push_back({remainingIndices[0], remainingIndices[1], remainingIndices[2]}); + } + return triangles; +} + /// Append the display triangulation of a wire-trimmed quadric patch: the sampled outer loop, ear-clipped; holes are omitted. template void appendCurveTrimMesh(const CurveWire& outerWire, const MapUV& mapUV, std::vector& vertices, @@ -2263,7 +2354,8 @@ void appendCurveTrimMesh(const CurveWire& outerWire, const MapUV& mapUV, std::ve for (const auto& sample : sampledWire.vertices) { vertices.push_back(mapUV(sample.uCoord, sample.vCoord)); } - for (const auto& triangle : triangulateSimpleWire(sampledWire)) { + const std::vector points(vertices.begin() + firstVertexIndex, vertices.end()); + for (const auto& triangle : triangulateWireByShortestChord(sampledWire, points)) { triangles.push_back( {firstVertexIndex + triangle[0], firstVertexIndex + triangle[1], firstVertexIndex + triangle[2]}); } diff --git a/Detectors/CADSupport/test/testBVHSurfaceSolid.cxx b/Detectors/CADSupport/test/testBVHSurfaceSolid.cxx index 589aea8cc5aeb..f2a4a20e7757e 100644 --- a/Detectors/CADSupport/test/testBVHSurfaceSolid.cxx +++ b/Detectors/CADSupport/test/testBVHSurfaceSolid.cxx @@ -2186,6 +2186,64 @@ BOOST_AUTO_TEST_CASE(WireTrimmedConeMatchesCone) checkClose(solid.Capacity(), reference.Capacity(), 1.e-9); } +namespace +{ +// Check that a patch's display mesh keeps its vertex count and stays within the chord sag of radius on the patch; +// boundarySag adds the sag of a trim side that the wire sampling leaves as one chord. +void checkDisplayMeshFollowsPatch(const surf::BoundedSurface& patch, double radius, size_t expectedVertices, + double boundarySag = 0.) +{ + std::vector vertices; + std::vector> triangles; + patch.appendDisplayMesh(vertices, triangles); + BOOST_CHECK_EQUAL(vertices.size(), expectedVertices); + BOOST_CHECK_EQUAL(triangles.size(), expectedVertices - 2); + const double maxSag = 1.05 * (radius * (1. - std::cos(surf::kPi / surf::kArcSamples)) + boundarySag); + double worstSag = 0.; + for (const auto& triangle : triangles) { + const surf::Vec3& first = vertices[triangle[0]]; + const surf::Vec3& second = vertices[triangle[1]]; + const surf::Vec3& third = vertices[triangle[2]]; + const double longestSq = std::max({normSq(second - first), normSq(third - second), normSq(first - third)}); + BOOST_CHECK_GT(norm(cross(second - first, third - first)), 1.e-6 * longestSq); + for (const surf::Vec3& probe : {(first + second + third) * (1. / 3.), (first + second) * 0.5, + (second + third) * 0.5, (third + first) * 0.5}) { + worstSag = std::max(worstSag, std::sqrt(patch.distanceSqToPatch(probe))); + } + } + BOOST_CHECK_LE(worstSag, maxSag); +} +} // namespace + +// The display triangles of a wire-trimmed quadric follow the patch instead of fanning across it (O2-7284). +BOOST_AUTO_TEST_CASE(WireTrimmedDisplayMeshFollowsPatch) +{ + std::string error; + constexpr double phiMin = 1.6; + constexpr double phiMax = 4.683; + + // half-cylinder of the ALICE 3 tracker shell: 12 chords on each long side plus the two corners + surf::CylindricalBoundedSurface halfCylinder; + BOOST_REQUIRE(halfCylinder.initialize({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, 10., -20., 20., phiMin, + phiMax - phiMin, false, paramRectWireCurves(phiMin, phiMax, -20., 20.), {}, + error)); + checkDisplayMeshFollowsPatch(halfCylinder, 10., 26); + + // quarter-torus: 6 chords on each ring side; each tube side is sampled as a single chord + surf::TorusBoundedSurface quarterTorus; + BOOST_REQUIRE(quarterTorus.initialize({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, 10., 0.3, 0., surf::kHalfPi, 0., + surf::kPi / 3., false, + paramRectWireCurves(0., surf::kHalfPi, 0., surf::kPi / 3.), {}, error)); + checkDisplayMeshFollowsPatch(quarterTorus, 10.3, 14, 0.3 * (1. - std::cos(surf::kPi / 6.))); + + // half-cone + surf::ConicalBoundedSurface halfCone; + BOOST_REQUIRE(halfCone.initialize({0., 0., 0.}, {0., 0., 1.}, {1., 0., 0.}, 8., 12., -10., 10., phiMin, + phiMax - phiMin, false, paramRectWireCurves(phiMin, phiMax, -10., 10.), {}, + error)); + checkDisplayMeshFollowsPatch(halfCone, 12., 26); +} + // Green's theorem for wire-trimmed quadrics. // // The integrand does not depend on the second parameter for a cylinder, and depends on the first