[Geometry] Bugfixes and add unit tests by fredroy · Pull Request #5987 · sofa-framework/sofa

fredroy · 2026-03-03T00:21:38Z

Fix various Edge and triangles functions

intersectionWith* (test against abs value to check both sides) and collinearity
fix 2d and 3d cases
isPointInTriangle: add checks if the point is in the triangle plane but also add a boolean to bypass this check (as it was before the "fix")
mat_LCP (used for proximity): the solver returned true even it fails to converge (max iterations reached)

and more unit tests for Sofa.Geometry module

[with-all-tests]

By submitting this pull request, I acknowledge that
I have read, understand, and agree SOFA Developer Certificate of Origin (DCO).

Reviewers will merge this pull-request only if

it builds with SUCCESS for all platforms on the CI.
it does not generate new warnings.
it does not generate new unit test failures.
it does not generate new scene test failures.
it does not break API compatibility.
it is more than 1 week old (or has fast-merge label).

fredroy · 2026-03-03T00:22:14Z

[ci-build][with-all-tests]

bakpaul · 2026-03-26T10:18:59Z

+            // In 3D, use signed areas via dot product with triangle normal
+            // to get correct sign for outside-triangle points
+            const auto N = sofa::type::cross(n1 - n0, n2 - n0);
+            const auto NdotN = sofa::type::dot(N, N);
+
+            if (NdotN < std::numeric_limits<T>::epsilon() * std::numeric_limits<T>::epsilon()) // triangle is flat
+            {
+                return sofa::type::Vec<3, T>(-1, -1, -1);
+            }
+
+            sofa::type::Vec<3, T> baryCoefs(type::NOINIT);
+            baryCoefs[0] = sofa::type::dot(N, sofa::type::cross(n2 - n1, p0 - n1)) / NdotN;
+            baryCoefs[1] = sofa::type::dot(N, sofa::type::cross(n0 - n2, p0 - n2)) / NdotN;
+            baryCoefs[2] = 1 - baryCoefs[0] - baryCoefs[1];
+
+            if (fabs(baryCoefs[2]) <= std::numeric_limits<T>::epsilon()){
+                baryCoefs[2] = 0;
+            }
+
+            return baryCoefs;


Ok, before this the barycentric mapping was a bit wrong when the point wasn't on the triangle. Because of 104, the value of the third barycentric coordinate was always underestimated because the other ones where overestimated.

@epernod this fix the cases when the point is not on the same plane as the triangle. But does it ever happen ? Anyway I think that the api should not expect anything and check because otherwise we would rely on the user to make the check before.

yes good catch I wrote all those functions as helper methods taking into account user know what they do but I agree to add all those isInPlane safeguards.

However here, in the Triangle::getBarycentricCoordinates method, if the point is not in a triangle, I think we should not do a projection of this point and send some valid values.
For me you should just call at start your new method: isPointOnPlane and if it is false you exit and return incorrect value like return sofa::type::Vec<3, T>(-1, -1, -1); for flat triangle.

you could also add a bool& isInPlane input parameter, that could be useful, but it means changing the signature of the method

you could also add a bool& isInPlane input parameter, that could be useful, but it means changing the signature of the method

I would just assume that if it returns false, the caller (if he wants) would investigate on the reason it is false (either by calling isOnPlane itself, or checking the resulting coeffs)

the problem here is getBarycentricCoordinates return the vec3 of barycoefs, not a bool. That's why I was speaking of having a additional bool& parameter.
Maybe changing the method signature to take the vec3 as input/output and returning a bool would be more efficient to avoid checking against -1 -1 -1 value when using the method.

The issue with what you propose Erik is that it will change the behavior of the API. Wasn't this previous "projection-like" mechanism useful in some cases where people used a barycentric mapping of points on a triangular mesh without making sure it was actually on the mesh ? Anyway, I expect great regression with this change. I would do it in another PR.

epernod

didn't check all the tests but the safeguard on checking if point is on a plane is indeed a good practice. However in this file as we request value regarding the triangle, if the point is out of the plane we should always return invalid value in my opinion

epernod · 2026-04-08T21:50:46Z

    >
-        static constexpr bool isPointInTriangle(const Node& p0, const Node& n0, const Node& n1, const Node& n2, sofa::type::Vec<3, T>& baryCoefs)
+    [[nodiscard]]
+    static constexpr bool isPointInTriangle(const Node& p0, const Node& n0, const Node& n1, const Node& n2, sofa::type::Vec<3, T>& baryCoefs, bool assumePointIsOnPlane = true)


shouldn't be the default value of assumePointIsOnPlane to false to always test it by default?

As it was not tested before, setting to false would potentially decrease the perfs automatically.

An other solution would be to set to false, and adding true to all existing call to isPointInTriangle(). But this would not cover the usage in external projects...
Last solution would be: no bool and two functions (isPointInTriangle and fastIsPointInTriangle or something like that)

yes maybe the last solution is better to avoid any missunderstanding.
a 2nd method isPointInTriangleNoPlaneCheck 😆

epernod · 2026-04-08T21:59:48Z

+    [[nodiscard]]
+    static constexpr bool isPointInTriangle(const Node& p0, const Node& n0, const Node& n1, const Node& n2, sofa::type::Vec<3, T>& baryCoefs, bool assumePointIsOnPlane = true)
    {
        baryCoefs = Triangle::getBarycentricCoordinates(p0, n0, n1, n2);


I would call this computation after the if/else to avoid some unnecessary calculation if point is not on plane. and init the baryCoef to same invalid value sofa::type::Vec<3, T>(-1, -1, -1);

aha you didn't see this one Claude!

Not sure about the best way if not on plane: setting the {-1, -1, -1 } as you suggest ; or do not touch the Vec3 as it will return false, and it is the caller's responsibility to manage this case

epernod · 2026-04-08T22:09:17Z

+            // In 3D, use signed areas via dot product with triangle normal
+            // to get correct sign for outside-triangle points
+            const auto N = sofa::type::cross(n1 - n0, n2 - n0);
+            const auto NdotN = sofa::type::dot(N, N);
+
+            if (NdotN < std::numeric_limits<T>::epsilon() * std::numeric_limits<T>::epsilon()) // triangle is flat
+            {
+                return sofa::type::Vec<3, T>(-1, -1, -1);
+            }
+
+            sofa::type::Vec<3, T> baryCoefs(type::NOINIT);
+            baryCoefs[0] = sofa::type::dot(N, sofa::type::cross(n2 - n1, p0 - n1)) / NdotN;
+            baryCoefs[1] = sofa::type::dot(N, sofa::type::cross(n0 - n2, p0 - n2)) / NdotN;
+            baryCoefs[2] = 1 - baryCoefs[0] - baryCoefs[1];
+
+            if (fabs(baryCoefs[2]) <= std::numeric_limits<T>::epsilon()){
+                baryCoefs[2] = 0;
+            }
+
+            return baryCoefs;


yes good catch I wrote all those functions as helper methods taking into account user know what they do but I agree to add all those isInPlane safeguards.

However here, in the Triangle::getBarycentricCoordinates method, if the point is not in a triangle, I think we should not do a projection of this point and send some valid values.
For me you should just call at start your new method: isPointOnPlane and if it is false you exit and return incorrect value like return sofa::type::Vec<3, T>(-1, -1, -1); for flat triangle.

epernod · 2026-04-08T22:11:10Z

+            // In 3D, use signed areas via dot product with triangle normal
+            // to get correct sign for outside-triangle points
+            const auto N = sofa::type::cross(n1 - n0, n2 - n0);
+            const auto NdotN = sofa::type::dot(N, N);
+
+            if (NdotN < std::numeric_limits<T>::epsilon() * std::numeric_limits<T>::epsilon()) // triangle is flat
+            {
+                return sofa::type::Vec<3, T>(-1, -1, -1);
+            }
+
+            sofa::type::Vec<3, T> baryCoefs(type::NOINIT);
+            baryCoefs[0] = sofa::type::dot(N, sofa::type::cross(n2 - n1, p0 - n1)) / NdotN;
+            baryCoefs[1] = sofa::type::dot(N, sofa::type::cross(n0 - n2, p0 - n2)) / NdotN;
+            baryCoefs[2] = 1 - baryCoefs[0] - baryCoefs[1];
+
+            if (fabs(baryCoefs[2]) <= std::numeric_limits<T>::epsilon()){
+                baryCoefs[2] = 0;
+            }
+
+            return baryCoefs;


you could also add a bool& isInPlane input parameter, that could be useful, but it means changing the signature of the method

epernod · 2026-04-08T22:19:17Z

+
+    // out of plan with assumePointIsInPlane=true (default): bary coords are all positive
+    // so the function returns true even though the point is not in the triangle plane
+    res = sofa::geometry::Triangle::isPointInTriangle(p0, a, b, c, bary);
+    EXPECT_TRUE(res);
+    EXPECT_NEAR(bary[0], 0.6, 1e-4);
+    EXPECT_NEAR(bary[1], 0.3, 1e-4);
+    EXPECT_NEAR(bary[2], 0.1, 1e-4);


for me this should not be allowed

epernod · 2026-04-08T22:20:23Z

+    EXPECT_EQ(bary[0], -1);
    EXPECT_EQ(bary[1], 1);
-    EXPECT_EQ(bary[2], -1);
+    EXPECT_EQ(bary[2], 1);


this change is strange but ok if it passes

if (denominator < EQUALITY_THRESHOLD) the denominator is a signed dot product. When negative (edge crosses plane from the "back"), this is always true, so the function rejects ~half of all valid intersections.

1-The 2D path computes alpha (parameter on edge AB) and checks 0 <= alpha <= 1, but never computes or checks beta (parameter on edge CD). Reports false intersections when the infinite line CD crosses segment AB but outside segment CD. 2-if (alphaDenom < std::numeric_limits<T>::epsilon()) same pattern as Bug 1. The cross product alphaDenom is signed. When negative (non-collinear edges in one orientation), incorrectly reports them as collinear.

solveLCP returns true when solver fails to converge (max iterations reached)

While isPointInTriangle still works (it only checks sign pattern), any code that uses the returned coordinates for interpolation, projection to the nearest point on the triangle, or distance computation gets wrong values. This is a public utility function likely consumed by many parts of the codebase.

…o bypass the plane check (as before)

…or for beta) Co-authored-by: Paul Baksic <30337881+bakpaul@users.noreply.github.com>

fredroy added pr: fix Fix a bug pr: status to review To notify reviewers to review this pull-request pr: test PR adding test(s) in SOFA pr: AI-aided Label notifying the reviewers that part or all of the PR has been generated with the help of an AI labels Mar 3, 2026

fredroy force-pushed the add_tests_geometry branch from ae273ee to 69cba4f Compare March 3, 2026 06:55

hugtalbot requested review from bakpaul and epernod March 12, 2026 09:33

bakpaul reviewed Mar 26, 2026

View reviewed changes

fredroy force-pushed the add_tests_geometry branch 2 times, most recently from 9cd1ebb to 5a274e4 Compare March 31, 2026 07:14

epernod requested changes Apr 8, 2026

View reviewed changes

fredroy and others added 11 commits April 9, 2026 08:46

fix test for intersectionWithPlane (if zero with tolerance)

3871364

if (denominator < EQUALITY_THRESHOLD) the denominator is a signed dot product. When negative (edge crosses plane from the "back"), this is always true, so the function rejects ~half of all valid intersections.

fix equality for 2d case of ispointonedge

a5a7c3d

Mat_solve_LCP: fix when reaching max nb loop

29f5ef4

solveLCP returns true when solver fails to converge (max iterations reached)

isPointInTriangle(3d) was not checking the plane, and add a boolean t…

88ccfd7

…o bypass the plane check (as before)

add more unit tests

f420a90

Update Sofa/framework/Geometry/src/sofa/geometry/Edge.h (fix sign err…

5451f63

…or for beta) Co-authored-by: Paul Baksic <30337881+bakpaul@users.noreply.github.com>

add isPointOnPlane function and add its unit tests

5bf1fbe

compute coeffs only if in plane; set -1,-1,-1 if not

645e987

update unit tests

e089be3

fredroy force-pushed the add_tests_geometry branch from 5a274e4 to e089be3 Compare April 9, 2026 00:09

epernod approved these changes Apr 9, 2026

View reviewed changes

bakpaul added pr: status ready Approved a pull-request, ready to be squashed and removed pr: status to review To notify reviewers to review this pull-request labels Apr 9, 2026

epernod merged commit 912335e into sofa-framework:master Apr 9, 2026
14 of 17 checks passed

Conversation

fredroy commented Mar 3, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

fredroy commented Mar 3, 2026

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

epernod left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

fredroy Apr 8, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

fredroy Apr 8, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

fredroy commented Mar 3, 2026 •

edited

Loading

fredroy Apr 8, 2026 •

edited

Loading

fredroy Apr 8, 2026 •

edited

Loading