Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

[Geometry] Add methods getBarycentricCoordinates and isPointInTriangle in Triangle class #4053

Merged
merged 13 commits into from
Sep 7, 2023
Merged

[Geometry] Add methods getBarycentricCoordinates and isPointInTriangle in Triangle class #4053

Show file tree
Hide file tree
Changes from 12 commits
Commits
Show all changes
13 commits
Select commit Hold shift + click to select a range
c4aba12
backup work on pointBarycoefs
epernod Jul 12, 2023
518e5ea
[Geometry] Update Triangle barycoef tests to handle negative area
epernod Jul 19, 2023
068784c
[Geometry] Update isPointInTriangle to return barycoefs vec3 abd upda…
epernod Jul 19, 2023
1fbb186
[Geometry] Update methods doc
epernod Aug 28, 2023
0e0e3c7
debug Triangle_test.cpp
epernod Sep 4, 2023
623c818
Update Triangle.h
epernod Sep 4, 2023
b81956b
Update Triangle.h
epernod Sep 4, 2023
cbee20c
Update Triangle.h
epernod Sep 4, 2023
7b2cb24
Update Triangle.h
epernod Sep 4, 2023
ca56b6e
Update Triangle.h
epernod Sep 4, 2023
6a8f82d
Update Triangle_test.cpp
epernod Sep 4, 2023
5c2e7b5
Update Triangle.h
epernod Sep 5, 2023
c6bb38b
rename method into getBarycentricCoordinates and add comments
epernod Sep 6, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
67 changes: 67 additions & 0 deletions Sofa/framework/Geometry/src/sofa/geometry/Triangle.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -71,6 +71,45 @@ struct Triangle
}
}


/**
* @brief Compute the barycentric coefficients of the input point according to the Triangle vertices: (n0, n1, n2)
* @tparam Node iterable container
* @tparam T scalar
* @param p0: position of the input point to compute the coefficients
* @param n0, n1, n2: nodes of the triangle
* @return sofa::type::Vec<3, T> barycentric coefficients
*/
template<typename Node,
typename T = std::decay_t<decltype(*std::begin(std::declval<Node>()))>,
typename = std::enable_if_t<std::is_scalar_v<T>>
>
static constexpr auto pointBaryCoefs(const Node& p0, const Node& n0, const Node& n1, const Node& n2)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

pointBaryCoefs, the name looks... not really informative and uses abbreviations.
Maybe something like getBarycentricCoefficients() or even getBarycentricCoordinates() which seems more in par with the maths

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done and add more description on the method (based from wikipedia)

{
// Point can be written: p0 = a*n0 + b*n1 + c*n2
// with a = area(n1n2p0)/area(n0n1n2), b = area(n0n2p0)/area(n0n1n2) and c = area(n0n1p0)/area(n0n1n2)
const auto area = Triangle::area(n0, n1, n2);
if (fabs(area) < std::numeric_limits<T>::epsilon()) // triangle is flat
{
return sofa::type::Vec<3, T>(-1, -1, -1);
}

const auto A0 = Triangle::area(n1, n2, p0);
const auto A1 = Triangle::area(n0, p0, n2);

sofa::type::Vec<3, T> baryCoefs(type::NOINIT);
baryCoefs[0] = A0 / area;
baryCoefs[1] = A1 / area;
baryCoefs[2] = 1 - baryCoefs[0] - baryCoefs[1];

if (fabs(baryCoefs[2]) <= std::numeric_limits<T>::epsilon()){
baryCoefs[2] = 0;
}

return baryCoefs;
}


/**
* @brief Compute the normal of a triangle
* @remark triangle normal computation is only possible in 3D
Expand Down Expand Up @@ -109,6 +148,34 @@ struct Triangle

}


/**
* @brief Test if input point is inside Triangle (n0, n1, n2)
* @tparam Node iterable container
* @tparam T scalar
* @param p0: position of the point to test
* @param n0, n1, n2: nodes of the triangle
* @param output parameter: sofa::type::Vec<3, T> barycentric coefficients of the input point in Triangle
* @return bool result if point is inside Triangle.
*/
template<typename Node,
typename T = std::decay_t<decltype(*std::begin(std::declval<Node>()))>,
typename = std::enable_if_t<std::is_scalar_v<T>>
>
static constexpr bool isPointInTriangle(const Node& p0, const Node& n0, const Node& n1, const Node& n2, sofa::type::Vec<3, T>& baryCoefs)
{
baryCoefs = Triangle::pointBaryCoefs(p0, n0, n1, n2);

for (int i = 0; i < 3; ++i)
{
if (baryCoefs[i] < 0 || baryCoefs[i] > 1)
return false;
}

return true;
}


/**
* @brief Test if a ray intersects a triangle, and gives barycentric coordinates of the intersection if applicable
* @remark Implementation for 3D only
Expand Down
179 changes: 168 additions & 11 deletions Sofa/framework/Geometry/test/Triangle_test.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,6 +39,15 @@ class Geometry3DTriangle_test : public ::testing::Test
{
};

template <typename VecType>
class GeometryVec2DTriangle_test : public ::testing::Test
{
};
template <typename VecType>
class GeometryVec3DTriangle_test : public ::testing::Test
{
};

using Vec2Types = ::testing::Types <
std::array<float, 2>, sofa::type::fixed_array<float, 2>, sofa::type::Vec < 2, float >,
std::array<double, 2>, sofa::type::fixed_array<double, 2>, sofa::type::Vec < 2, double >>;
Expand All @@ -49,6 +58,13 @@ using Vec3Types = ::testing::Types <
std::array<double, 3>, sofa::type::fixed_array<double, 3>, sofa::type::Vec < 3, double >>;
TYPED_TEST_SUITE(Geometry3DTriangle_test, Vec3Types);

using SofaVec2Types = ::testing::Types <sofa::type::Vec < 2, float >, sofa::type::Vec < 2, double > >;
TYPED_TEST_SUITE(GeometryVec2DTriangle_test, SofaVec2Types);

using SofaVec3Types = ::testing::Types <sofa::type::Vec < 3, float >, sofa::type::Vec < 3, double > >;
TYPED_TEST_SUITE(GeometryVec3DTriangle_test, SofaVec3Types);



TYPED_TEST(Geometry2DTriangle_test, area)
{
Expand All @@ -58,7 +74,7 @@ TYPED_TEST(Geometry2DTriangle_test, area)

const auto testArea = sofa::geometry::Triangle::area(a, b, c);

EXPECT_FLOAT_EQ(testArea, 12.5f);
EXPECT_FLOAT_EQ(float(testArea), 12.5f);
}

TYPED_TEST(Geometry3DTriangle_test, area)
Expand All @@ -68,7 +84,7 @@ TYPED_TEST(Geometry3DTriangle_test, area)
const TypeParam c{ 2.f, -3.f, 4.f };

const auto testArea = sofa::geometry::Triangle::area(a, b, c);
EXPECT_FLOAT_EQ(testArea, 19.306734f);
EXPECT_FLOAT_EQ(float(testArea), 19.306734f);
}

TYPED_TEST(Geometry2DTriangle_test, flat_area)
Expand All @@ -78,7 +94,7 @@ TYPED_TEST(Geometry2DTriangle_test, flat_area)
const TypeParam c{ 0.f, 1.f };

const auto testArea = sofa::geometry::Triangle::area(a, b, c);
EXPECT_FLOAT_EQ(testArea, 0.f);
EXPECT_FLOAT_EQ(testArea, 0);
}

TYPED_TEST(Geometry3DTriangle_test, flat_area)
Expand All @@ -88,7 +104,7 @@ TYPED_TEST(Geometry3DTriangle_test, flat_area)
const TypeParam c{ 0.f, 1.f, 0.f };

const auto testArea = sofa::geometry::Triangle::area(a, b, c);
EXPECT_FLOAT_EQ(testArea, 0.f);
EXPECT_EQ(testArea, 0);
}


Expand All @@ -98,23 +114,164 @@ TYPED_TEST(Geometry3DTriangle_test, normal)
const TypeParam a{ 0.f, 0.f, 0.f };
const TypeParam b{ 0.f, 2.f, 0.f };
const TypeParam c{ 0.f, 0.f, 2.f };

auto normal = sofa::geometry::Triangle::normal(a, b, c);
EXPECT_FLOAT_EQ(normal[0], 4.f);
EXPECT_FLOAT_EQ(normal[1], 0.f);
EXPECT_FLOAT_EQ(normal[2], 0.f);
EXPECT_FLOAT_EQ(float(normal[0]), 4.f);
EXPECT_FLOAT_EQ(float(normal[1]), 0.f);
EXPECT_FLOAT_EQ(float(normal[2]), 0.f);

// flat triangle case
const TypeParam a2{ 0.f, 0.f, 0.f };
const TypeParam b2{ 0.f, 2.f, 0.f };
const TypeParam c2{ 0.f, 1.f, 0.f };

normal = sofa::geometry::Triangle::normal(a2, b2, c2);
EXPECT_FLOAT_EQ(normal[0], 0.f);
EXPECT_FLOAT_EQ(normal[1], 0.f);
EXPECT_FLOAT_EQ(normal[2], 0.f);
EXPECT_FLOAT_EQ(float(normal[0]), 0.f);
EXPECT_FLOAT_EQ(float(normal[1]), 0.f);
EXPECT_FLOAT_EQ(float(normal[2]), 0.f);
}


TYPED_TEST(GeometryVec2DTriangle_test, isPointInTriangle)
{
const TypeParam a{ 0., 0. };
const TypeParam b{ 2., 0. };
const TypeParam c{ 2., 2. };
sofa::type::Vec3 bary;

//// point inside
TypeParam p0{ 1.5, 0.5 };
auto res = sofa::geometry::Triangle::isPointInTriangle(p0, a, b, c, bary);
EXPECT_TRUE(res);
EXPECT_FLOAT_EQ(float(bary[0]), 0.25f);
EXPECT_FLOAT_EQ(float(bary[1]), 0.5f);
EXPECT_FLOAT_EQ(float(bary[2]), 0.25f);

// barycenter
p0 = { 4./3., 2./3. };
const auto bVal = float(1.f / 3.f);
res = sofa::geometry::Triangle::isPointInTriangle(p0, a, b, c, bary);
EXPECT_TRUE(res);
EXPECT_NEAR(bary[0], bVal, 1e-4);
EXPECT_NEAR(bary[1], bVal, 1e-4);
EXPECT_NEAR(bary[2], bVal, 1e-4);

// on edge
p0 = { 1., 0. };
res = sofa::geometry::Triangle::isPointInTriangle(p0, a, b, c, bary);
EXPECT_TRUE(res);
EXPECT_NEAR(bary[0], 0.5f, 1e-4);
EXPECT_NEAR(bary[1], 0.5f, 1e-4);
EXPECT_NEAR(bary[2], 0.0f, 1e-4);

// on corner
p0 = { 2., 0. };
res = sofa::geometry::Triangle::isPointInTriangle(p0, a, b, c, bary);
EXPECT_TRUE(res);
EXPECT_FLOAT_EQ(float(bary[0]), 0.0f);
EXPECT_FLOAT_EQ(float(bary[1]), 1.0f);
EXPECT_FLOAT_EQ(float(bary[2]), 0.0f);

// False case
p0 = { 4., 10. };
res = sofa::geometry::Triangle::isPointInTriangle(p0, a, b, c, bary);
EXPECT_FALSE(res);
EXPECT_NEAR(bary[0], -1.0f, 1e-4);
EXPECT_NEAR(bary[1], -3.0f, 1e-4);
EXPECT_NEAR(bary[2], 5.0f, 1e-4);


// Special cases
// flat triangle
const TypeParam d{ 1., 0. };
p0 = { 0.5, 0. };
res = sofa::geometry::Triangle::isPointInTriangle(p0, a, b, d, bary);
EXPECT_FALSE(res);
EXPECT_EQ(bary[0], -1);
EXPECT_EQ(bary[1], -1);
EXPECT_EQ(bary[2], -1);

//special False cases along edge
p0 = { 3., 3. };
res = sofa::geometry::Triangle::isPointInTriangle(p0, a, b, c, bary);
EXPECT_FALSE(res);
EXPECT_NEAR(bary[0], -0.5f, 1e-4);
EXPECT_NEAR(bary[1], 0.f, 1e-4);
EXPECT_NEAR(bary[2], 1.5f, 1e-4);
}


TYPED_TEST(GeometryVec3DTriangle_test, isPointInTriangle)
{
const TypeParam a{ 0., 0., 0. };
const TypeParam b{ 2., 0., 2. };
const TypeParam c{ 0., 2., 0. };

sofa::type::Vec3 bary;

// point inside
TypeParam p0{ 0.5, 0.5, 0.5 };
auto res = sofa::geometry::Triangle::isPointInTriangle(p0, a, b, c, bary);
EXPECT_TRUE(res);
EXPECT_NEAR(bary[0], 0.5f, 1e-4);
EXPECT_NEAR(bary[1], 0.25f, 1e-4);
EXPECT_NEAR(bary[2], 0.25f, 1e-4);

// barycenter
p0 = { 2. / 3., 2. / 3. , 2. / 3. };
const auto bVal = float(1.f / 3.f);
res = sofa::geometry::Triangle::isPointInTriangle(p0, a, b, c, bary);
EXPECT_TRUE(res);
EXPECT_NEAR(bary[0], bVal, 1e-4);
EXPECT_NEAR(bary[1], bVal, 1e-4);
EXPECT_NEAR(bary[2], bVal, 1e-4);

// on edge
p0 = { 1., 1., 1. };
res = sofa::geometry::Triangle::isPointInTriangle(p0, a, b, c, bary);
EXPECT_TRUE(res);
EXPECT_NEAR(bary[0], 0.f, 1e-4);
EXPECT_NEAR(bary[1], 0.5f, 1e-4);
EXPECT_NEAR(bary[2], 0.5f, 1e-4);

// on corner
p0 = { 0., 2., 0. };
res = sofa::geometry::Triangle::isPointInTriangle(p0, a, b, c, bary);
EXPECT_TRUE(res);
EXPECT_EQ(bary[0], 0);
EXPECT_EQ(bary[1], 0);
EXPECT_EQ(bary[2], 1);

// False cases
// out of plan
p0 = { 1., 0.2, 0.2 };
res = sofa::geometry::Triangle::isPointInTriangle(p0, a, b, c, bary);
EXPECT_FALSE(res);
EXPECT_NEAR(bary[0], 0.69282, 1e-4);
EXPECT_NEAR(bary[1], 0.360555, 1e-4);
EXPECT_NEAR(bary[2], -0.0533754, 1e-4);

// in plan but out of triangle
p0 = { 2., 2., 2. };
res = sofa::geometry::Triangle::isPointInTriangle(p0, a, b, c, bary);
EXPECT_FALSE(res);
EXPECT_EQ(bary[0], 1);
EXPECT_EQ(bary[1], 1);
EXPECT_EQ(bary[2], -1);


// Special cases
// flat triangle
const TypeParam d{ 1., 0., 1. };
p0 = { 0.5, 0., 0.5 };
res = sofa::geometry::Triangle::isPointInTriangle(p0, a, b, d, bary);
EXPECT_FALSE(res);
EXPECT_EQ(bary[0], -1);
EXPECT_EQ(bary[1], -1);
EXPECT_EQ(bary[2], -1);
}


TEST(GeometryTriangle_test, rayIntersectionVec3)
{
const sofa::type::Vec3 a{ 0., 3., 0. };
Expand Down
Loading