Skip to content

implemented fast triangle-triangle intersection of devillers and guigue #7

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

Open
wants to merge 1 commit into
base: develop
Choose a base branch
from
Open

implemented fast triangle-triangle intersection of devillers and guigue #7

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
350 changes: 324 additions & 26 deletions src/typed-geometry/functions/objects/intersection.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -256,6 +256,58 @@ template <class ScalarT, class B>

template <class A, class B>
using try_closest_intersection_parameter = decltype(closest_intersection_parameter(std::declval<A const&>(), std::declval<B const&>()));

// circular permutation to the vertices of triangle ta such that ta.pos0 is the only vertex that lies on positive halfspace induced by tb
template <class ScalarT>
void rotate_devillers_triangle(tg::triangle<3, ScalarT>& ta, tg::triangle<3, ScalarT>& tb, tg::comp3& determinants, tg::comp3& determinants_t2)
Copy link
Contributor

Choose a reason for hiding this comment

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

comp3 -> comp<3, ScalarT>

Copy link
Contributor

Choose a reason for hiding this comment

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

remove tg:: in all places

{
// implementation of triangle permutation according to: https://hal.inria.fr/inria-00072100/document

ScalarT d01 = determinants[0] * determinants[1];
ScalarT d02 = determinants[0] * determinants[2];

if (d01 > 0.0f) // vertices 0 and 1 on one side
Copy link
Contributor

Choose a reason for hiding this comment

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

ScalarT(0)

{
ta = {ta.pos2, ta.pos0, ta.pos1};
determinants = {determinants[2], determinants[0], determinants[1]};
}
else if (d02 > 0.0f) // vertices 0 and 2 on one side
{
ta = {ta.pos1, ta.pos2, ta.pos0};
determinants = {determinants[1], determinants[2], determinants[0]};
}
else if (determinants[0] == 0.f)
{
if (determinants[1] * determinants[2] < 0.f || determinants[1] == 0.f) // vertices 1 and 2 on different sides of plane and vertex 0 on plane
{
ta = {ta.pos2, ta.pos0, ta.pos1};
determinants = {determinants[2], determinants[0], determinants[1]};
}
else if (determinants[2] == 0.f) // vertices 0 and 2 on the plane
{
ta = {ta.pos1, ta.pos2, ta.pos0};
determinants = {determinants[1], determinants[2], determinants[0]};
}
}

// swap operation to map triangle1.pos0 to positive halfspace induced by plane of triangle2
if (determinants[0] < 0.f)
{
tb = {tb.pos0, tb.pos2, tb.pos1};
determinants = {-determinants[0], -determinants[1], -determinants[2]};
determinants_t2 = {determinants_t2[0], determinants_t2[2], determinants_t2[1]};
}

else if (determinants[0] == 0.f && (determinants[1] * determinants[2] > 0))
{
if (determinants[1] > 0) // swap
{
tb = {tb.pos0, tb.pos2, tb.pos1};
determinants = {-determinants[0], -determinants[1], -determinants[2]};
determinants_t2 = {determinants_t2[0], determinants_t2[2], determinants_t2[1]};
}
}
}
}


Expand Down Expand Up @@ -2418,50 +2470,296 @@ template <class ScalarT>
}

template <class ScalarT>
[[nodiscard]] constexpr optional<segment<3, ScalarT>> intersection(triangle<3, ScalarT> const& t1, triangle<3, ScalarT> const& t2)
[[nodiscard]] constexpr bool intersects(tg::triangle<2, ScalarT> const& t1, tg::triangle<2, ScalarT> const& t2)
{
// early out: check with plane clamped by triangle t1
auto const plane_t1 = plane_of(t1);
// Implementation of: https://hal.inria.fr/inria-00072100/document

if (!intersects(t2, plane_t1))
return {};
auto const determin = [](tg::pos<2, ScalarT> pa, tg::pos<2, ScalarT> pb, tg::pos<2, ScalarT> pc) -> float
{
auto m = tg::mat<2, 2, ScalarT>::from_data_colwise({pa.x - pc.x, pb.x - pc.x, pa.y - pc.y, pb.y - pc.y});
return tg::determinant(m);
};

auto const counter_clock = [determin](tg::triangle<2, ScalarT>& tri_a) -> void
{
// if (tg::cross(vec<2, ScalarT>(tri_a.pos1 - tri_a.pos0), vec<2, ScalarT>(tri_a.pos2 - tri_a.pos0)) < 0)
if (determin(tri_a.pos0, tri_a.pos1, tri_a.pos2) < 0)
tri_a = tg::triangle<2, ScalarT>(tri_a.pos0, tri_a.pos2, tri_a.pos1); // swap
};

auto const rotate = [determin](tg::triangle<2, ScalarT>& tri_b, tg::triangle<2, ScalarT>& tri_a, tg::comp3& determinants_a) -> void
{
// rotate triangle vertices by one
tri_b = tg::triangle<2, ScalarT>(tri_b.pos2, tri_b.pos0, tri_b.pos1);

array<pos<3, ScalarT>, 2> insecs;
array<segment<3, ScalarT>, 3> segments_t1 = edges_of(t1);
array<segment<3, ScalarT>, 3> segments_t2 = edges_of(t2);
int insec_count = 0;
// recompute determinants
determinants_a = tg::comp3(determin(tri_b.pos0, tri_b.pos1, tri_a.pos0), determin(tri_b.pos1, tri_b.pos2, tri_a.pos0),
determin(tri_b.pos2, tri_b.pos0, tri_a.pos0));
};

tg::triangle<2, ScalarT> ta = t1;
tg::triangle<2, ScalarT> tb = t2;

// ensure counter-clockwise orientation
counter_clock(ta);
counter_clock(tb);

auto det_ta0 = tg::comp3(determin(tb.pos0, tb.pos1, ta.pos0), determin(tb.pos1, tb.pos2, ta.pos0), determin(tb.pos2, tb.pos0, ta.pos0));
auto det_01 = det_ta0[0] * det_ta0[1];
auto det_12 = det_ta0[1] * det_ta0[2];
auto det_02 = det_ta0[0] * det_ta0[2];


if (det_01 > 0.f && det_12 > 0.f) // ta.pos0 interior of tb
return true;

if (det_01 == 0 && det_12 == 0 && det_02 == 0) // ta.pos0 lies on vertex of tb
return true;

// check intersection of t1 segments with t2
for (auto const& s : segments_t1)
if (det_01 == 0 && ((det_ta0[1] > 0 && det_ta0[2] > 0) || (det_ta0[0] > 0 && det_ta0[2] > 0))) // ta.pos0 interior of edge of tb
return true;

if (det_12 == 0 && ((det_ta0[0] > 0 && det_ta0[1] > 0) || (det_ta0[0] > 0 && det_ta0[2] > 0))) // ta.pos0 interior of edge of tb
return true;

// Rotate to Region1 or Region2
while (!(det_ta0[0] > 0 && det_ta0[2] < 0))
rotate(tb, ta, det_ta0);

// decision tree
// R1
if (det_ta0[1] > 0)
{
auto insec = intersection(s, t2);
// I
if (determin(tb.pos2, tb.pos0, ta.pos1) >= 0)
{
// II.a
if (determin(tb.pos2, ta.pos0, ta.pos1) < 0)
return false;
// III.a
else if (determin(ta.pos0, tb.pos0, ta.pos1) < 0)
{
// IV.a
if (determin(ta.pos0, tb.pos0, ta.pos2) < 0)
return false;
// V
else if (determin(ta.pos1, ta.pos2, tb.pos0) < 0)
return false;
}
return true;
}
// II.b
else if (determin(tb.pos2, tb.pos0, ta.pos2) < 0)
return false;
// III.b
else if (determin(ta.pos1, ta.pos2, tb.pos2) < 0)
return false;
// IV.b
else if (determin(ta.pos0, tb.pos0, ta.pos2) < 0)
return false;

return true;
}

if (insec.has_value())
// R2
else if (determin(tb.pos2, tb.pos0, ta.pos1) >= 0)
{
// II.a
if (determin(tb.pos1, tb.pos2, ta.pos1) >= 0)
{
insecs[insec_count++] = insec.value();
// III.a
if (determin(ta.pos0, tb.pos0, ta.pos1) >= 0)
{
// IV.a
if (determin(ta.pos0, tb.pos1, ta.pos1) > 0)
return false;

return true;
}
// IV.b
else if (determin(ta.pos0, tb.pos0, ta.pos2) < 0)
return false;
// V.a
else if (determin(tb.pos2, tb.pos0, ta.pos2) < 0)
{
if (determin(ta.pos1, tb.pos0, ta.pos2) > 0)
return false;
}

return true;
}
// III.b
else if (determin(ta.pos0, tb.pos1, ta.pos1) > 0)
return false;
// IV.c
else if (determin(tb.pos1, tb.pos2, ta.pos2) < 0)
return false;
// V.b
else if (determin(ta.pos1, ta.pos2, tb.pos1) < 0)
return false;

if (insec_count >= 2)
return segment<3, ScalarT>{insecs[0], insecs[1]};
return true;
}
// II.b
else if (determin(tb.pos2, tb.pos0, ta.pos2) < 0)
return false;
// III.c
else if (determin(ta.pos1, ta.pos2, tb.pos2) >= 0)
{
// IV.d
if (determin(ta.pos2, ta.pos0, tb.pos0) < 0)
return false;
}
// IV.e
else if (determin(ta.pos1, ta.pos2, tb.pos1) < 0)
return false;
// V.c
else if (determin(tb.pos1, tb.pos2, ta.pos2) < 0)
return false;

return true;
}

template <class ScalarT>
[[nodiscard]] constexpr bool intersects(tg::triangle<3, ScalarT> const& t1, tg::triangle<3, ScalarT> const& t2)
{
// Implementation of https://hal.inria.fr/inria-00072100/document

auto const determin = [](tg::pos<3, ScalarT> pa, tg::pos<3, ScalarT> pb, tg::pos<3, ScalarT> pc, tg::pos<3, ScalarT> pd) -> float
{
auto m = tg::mat<3, 3, ScalarT>::from_data_colwise(
{pa.x - pd.x, pb.x - pd.x, pc.x - pd.x, pa.y - pd.y, pb.y - pd.y, pc.y - pd.y, pa.z - pd.z, pb.z - pd.z, pc.z - pd.z});
return tg::determinant(m);
};

auto det_t2_t1 = tg::comp3(determin(t2.pos0, t2.pos1, t2.pos2, t1.pos0), determin(t2.pos0, t2.pos1, t2.pos2, t1.pos1),
determin(t2.pos0, t2.pos1, t2.pos2, t1.pos2));

auto const dt2_01 = det_t2_t1[0] * det_t2_t1[1];
auto const dt2_02 = det_t2_t1[0] * det_t2_t1[2];

// a) no insec of t1 with plane of t2
if (dt2_01 > 0.f && dt2_02 > 0.f)
return false;

// b) coplanar -> all dets are 0
if (det_t2_t1[0] == det_t2_t1[1] && det_t2_t1[1] == det_t2_t1[2] && det_t2_t1[2] == 0.f)
{
auto n = normal_of(t1);
// find appropriate projection plane
auto proj_plane = dot(n, tg::dir<3, ScalarT>(0.f, 1.f, 0.f)) == 0.f ? tg::plane<3, ScalarT>({0.f, 0.f, 1.f}, tg::pos<3, ScalarT>::zero)
: tg::plane<3, ScalarT>({0.f, 1.f, 0.f}, tg::pos<3, ScalarT>::zero);
// project triangles onto plane
auto t1_2D = proj_plane.normal.z == 0
? tg::triangle<2, ScalarT>(xz(project(t1.pos0, proj_plane)), xz(project(t1.pos1, proj_plane)), xz(project(t1.pos2, proj_plane)))
: tg::triangle<2, ScalarT>(xy(project(t1.pos0, proj_plane)), xy(project(t1.pos1, proj_plane)), xy(project(t1.pos2, proj_plane)));
auto t2_2D = proj_plane.normal.z == 0
? tg::triangle<2, ScalarT>(xz(project(t2.pos0, proj_plane)), xz(project(t2.pos1, proj_plane)), xz(project(t2.pos2, proj_plane)))
: tg::triangle<2, ScalarT>(xy(project(t2.pos0, proj_plane)), xy(project(t2.pos1, proj_plane)), xy(project(t2.pos2, proj_plane)));

return intersects(t1_2D, t2_2D);
}

auto det_t1_t2 = tg::comp3(determin(t1.pos0, t1.pos1, t1.pos2, t2.pos0), determin(t1.pos0, t1.pos1, t1.pos2, t2.pos1),
determin(t1.pos0, t1.pos1, t1.pos2, t2.pos2));

auto const dt1_01 = det_t1_t2[0] * det_t1_t2[1];
auto const dt1_02 = det_t1_t2[0] * det_t1_t2[2];

// a) no insec of t2 with plane of t1
if (dt1_01 > 0.f && dt1_02 > 0.f)
return false;

tg::triangle<3, ScalarT> ta = t1;
tg::triangle<3, ScalarT> tb = t2;

// check intersection of t2 segments with t1
for (auto const& s : segments_t2)
// circular permutation
detail::rotate_devillers_triangle(ta, tb, det_t2_t1, det_t1_t2);
detail::rotate_devillers_triangle(tb, ta, det_t1_t2, det_t2_t1);

// decision tree
if (determin(ta.pos0, ta.pos1, tb.pos0, tb.pos1) > 0)
return false;

if (determin(ta.pos0, ta.pos2, tb.pos2, tb.pos0) > 0)
return false;

return true;
}

template <class ScalarT>
[[nodiscard]] constexpr optional<segment<3, ScalarT>> intersection(triangle<3, ScalarT> const& t1, triangle<3, ScalarT> const& t2)
{
// Implementation of https://hal.inria.fr/inria-00072100/document

auto const determin = [](tg::pos<3, ScalarT> pa, tg::pos<3, ScalarT> pb, tg::pos<3, ScalarT> pc, tg::pos<3, ScalarT> pd) -> float
{
auto insec = intersection(s, t1);
if (insec.has_value())
auto m = tg::mat<3, 3, ScalarT>::from_data_colwise(
{pa.x - pd.x, pb.x - pd.x, pc.x - pd.x, pa.y - pd.y, pb.y - pd.y, pc.y - pd.y, pa.z - pd.z, pb.z - pd.z, pc.z - pd.z});
return tg::determinant(m);
};

auto det_t2_t1 = tg::comp3(determin(t2.pos0, t2.pos1, t2.pos2, t1.pos0), determin(t2.pos0, t2.pos1, t2.pos2, t1.pos1),
determin(t2.pos0, t2.pos1, t2.pos2, t1.pos2));

auto const dt2_01 = det_t2_t1[0] * det_t2_t1[1];
auto const dt2_02 = det_t2_t1[0] * det_t2_t1[2];

// a) no insec of t1 with plane of t2
if (dt2_01 > 0.f && dt2_02 > 0.f)
return {};

// b) coplanar -> case not yet handled (arbitrary polygonal return type required)
if (det_t2_t1[0] == det_t2_t1[1] && det_t2_t1[1] == det_t2_t1[2] && det_t2_t1[2] == 0.f)
return {};

auto det_t1_t2 = tg::comp3(determin(t1.pos0, t1.pos1, t1.pos2, t2.pos0), determin(t1.pos0, t1.pos1, t1.pos2, t2.pos1),
determin(t1.pos0, t1.pos1, t1.pos2, t2.pos2));

auto const dt1_01 = det_t1_t2[0] * det_t1_t2[1];
auto const dt1_02 = det_t1_t2[0] * det_t1_t2[2];

// a) no insec of t2 with plane of t1
if (dt1_01 > 0.f && dt1_02 > 0.f)
return {};

tg::triangle<3, ScalarT> ta = t1;
tg::triangle<3, ScalarT> tb = t2;

// circular permutation for triangle orientation
detail::rotate_devillers_triangle(ta, tb, det_t2_t1, det_t1_t2);
detail::rotate_devillers_triangle(tb, ta, det_t1_t2, det_t2_t1);

// decision tree
tg::plane<3, ScalarT> p1 = plane_of(ta);
tg::plane<3, ScalarT> p2 = plane_of(tb);

// intersects tests
if (determin(ta.pos0, ta.pos1, tb.pos0, tb.pos1) > 0)
return {};

else if (determin(ta.pos0, ta.pos2, tb.pos2, tb.pos0) > 0)
return {};

// decision tree to find intersection segment
else if (determin(ta.pos0, ta.pos2, tb.pos1, tb.pos0) > 0)
{
if (determin(ta.pos0, ta.pos1, tb.pos2, tb.pos0) > 0)
{
insecs[insec_count++] = insec.value();
return tg::segment<3, ScalarT>{intersection(tg::inf_of(tg::segment<3, ScalarT>{ta.pos0, ta.pos2}), p2).first(),
intersection(tg::inf_of(tg::segment<3, ScalarT>{tb.pos0, tb.pos2}), p1).first()};
}

if (insec_count >= 2)
return segment<3, ScalarT>{insecs[0], insecs[1]};
return tg::segment<3, ScalarT>{intersection(tg::inf_of(tg::segment<3, ScalarT>{ta.pos0, ta.pos2}), p2).first(),
intersection(tg::inf_of(tg::segment<3, ScalarT>{tb.pos0, tb.pos1}), p1).first()};
}

if (insec_count == 1)
return segment<3, ScalarT>{insecs[0], insecs[0]};
else if (determin(ta.pos0, ta.pos1, tb.pos2, tb.pos0) > 0)
return tg::segment<3, ScalarT>{intersection(tg::inf_of(tg::segment<3, ScalarT>{tb.pos0, tb.pos1}), p1).first(),
intersection(tg::inf_of(tg::segment<3, ScalarT>{tb.pos0, tb.pos2}), p1).first()};

return {};
return tg::segment<3, ScalarT>{intersection(tg::inf_of(tg::segment<3, ScalarT>{tb.pos0, tb.pos1}), p1).first(),
intersection(tg::inf_of(tg::segment<3, ScalarT>{ta.pos0, ta.pos1}), p2).first()};
}

template <class ScalarT>
Expand Down
Loading