songsenand
/
cgal
mirror of https://github.com/CGAL/cgal
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #5680 from afabri/Intersections-fix_RT-GF 

Intersections_3: Fix do_intersect(sphere/triangle) if FT has no division
This commit is contained in:
Laurent Rineau 2021-05-27 19:12:35 +02:00
parent 908049d7f4 2820a9430e
commit 38812be5eb
7 changed files with 357 additions and 129 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

45
Distance_3/include/CGAL/squared_distance_3_0.h
View File

@ -26,6 +26,7 @@
#include <CGAL/Weighted_point_3.h>
#include <CGAL/Vector_3.h>
#include <CGAL/number_utils.h>
#include <CGAL/Rational_traits.h>
namespace CGAL {
@ -187,6 +188,21 @@ squared_distance(const typename K::Point_3 & pt1,
return k.compute_squared_distance_3_object()(pt1, pt2);
}
template <class K>
void
squared_distance_to_plane_RT(const typename K::Vector_3 & normal,
const typename K::Vector_3 & diff,
typename K::RT& num,
typename K::RT& den,
const K& k)
{
typedef typename K::RT RT;
RT dot, squared_length;
dot = wdot(normal, diff, k);
squared_length = wdot(normal, normal, k);
num = square(dot);
den = wmult((K*)0, squared_length, diff.hw(), diff.hw());
}
template <class K>
typename K::FT
@ -196,13 +212,24 @@ squared_distance_to_plane(const typename K::Vector_3 & normal,
{
typedef typename K::RT RT;
typedef typename K::FT FT;
RT dot, squared_length;
dot = wdot(normal, diff, k);
squared_length = wdot(normal, normal, k);
return FT(dot*dot) /
FT(wmult((K*)0, squared_length, diff.hw(), diff.hw()));
RT num, den;
squared_distance_to_plane_RT(normal, diff, num, den, k);
return Rational_traits<FT>().make_rational(num, den);
}
template <class K>
void
squared_distance_to_line_RT(const typename K::Vector_3 & dir,
const typename K::Vector_3 & diff,
typename K::RT& num,
typename K::RT& den,
const K& k)
{
typedef typename K::Vector_3 Vector_3;
Vector_3 wcr = wcross(dir, diff, k);
num = wdot(wcr, wcr, k);
den = wmult((K*)0, wdot(dir, dir, k), diff.hw(), diff.hw());
}
template <class K>
typename K::FT
@ -210,15 +237,13 @@ squared_distance_to_line(const typename K::Vector_3 & dir,
const typename K::Vector_3 & diff,
const K& k)
{
typedef typename K::Vector_3 Vector_3;
typedef typename K::RT RT;
typedef typename K::FT FT;
Vector_3 wcr = wcross(dir, diff, k);
return FT(wcr*wcr)/FT(wmult(
(K*)0, RT(wdot(dir, dir, k)), diff.hw(), diff.hw()));
RT num, den;
squared_distance_to_line_RT(dir, diff, num, den, k);
return Rational_traits<FT>().make_rational(num, den);
}
template <class K>
inline
bool

188
Distance_3/include/CGAL/squared_distance_3_1.h
View File

@ -31,17 +31,33 @@ namespace CGAL {
namespace internal {
template <class K>
typename K::FT
squared_distance(
void
squared_distance_RT(
const typename K::Point_3 &pt,
const typename K::Line_3 &line,
typename K::RT& num,
typename K::RT& den,
const K& k)
{
typedef typename K::Vector_3 Vector_3;
typename K::Construct_vector_3 construct_vector;
Vector_3 dir(line.direction().vector());
Vector_3 diff = construct_vector(line.point(), pt);
return internal::squared_distance_to_line(dir, diff, k);
return internal::squared_distance_to_line_RT(dir, diff, num, den, k);
}
template <class K>
typename K::FT
squared_distance(
const typename K::Point_3 &pt,
const typename K::Line_3 &line,
const K& k)
{
typedef typename K::RT RT;
typedef typename K::FT FT;
RT num, den;
squared_distance_RT(pt, line, num, den, k);
return Rational_traits<FT>().make_rational(num, den);
}
template <class K>
@ -56,6 +72,31 @@ squared_distance(
}
template <class K>
void
squared_distance_RT(
const typename K::Point_3 &pt,
const typename K::Ray_3 &ray,
typename K::RT& num,
typename K::RT& den,
const K& k)
{
typename K::Construct_vector_3 construct_vector;
typedef typename K::Vector_3 Vector_3;
typedef typename K::RT RT;
Vector_3 diff = construct_vector(ray.source(), pt);
const Vector_3 &dir = ray.direction().vector();
if (!is_acute_angle(dir,diff, k) )
{
num = wdot(diff, diff, k);
den = wmult((K*)0, RT(1), diff.hw(), diff.hw());
return;
}
squared_distance_to_line_RT(dir, diff, num, den, k);
}
template <class K>
typename K::FT
squared_distance(
@ -63,16 +104,25 @@ squared_distance(
const typename K::Ray_3 &ray,
const K& k)
{
typename K::Construct_vector_3 construct_vector;
// This duplicates code from the _RT functions, but it is a slowdown to do something like:
//
// RT num, den;
// squared_distance_RT(pt, ray, num, den, k);
// return Rational_traits<FT>().make_rational(num, den);
//
// See https://github.com/CGAL/cgal/pull/5680
typedef typename K::Vector_3 Vector_3;
Vector_3 diff = construct_vector(ray.source(), pt);
const Vector_3 &dir = ray.direction().vector();
if (!is_acute_angle(dir,diff, k) )
return (typename K::FT)(diff*diff);
return squared_distance_to_line(dir, diff, k);
}
typename K::Construct_vector_3 construct_vector = k.construct_vector_3_object();
Vector_3 diff = construct_vector(ray.source(), pt);
const Vector_3 &dir = ray.direction().vector();
if (!is_acute_angle(dir,diff, k) )
return (typename K::FT)(diff*diff);
return squared_distance_to_line(dir, diff, k);
}
template <class K>
inline
@ -85,75 +135,77 @@ squared_distance(
return squared_distance(pt, ray, k);
}
template <class K>
typename K::FT
squared_distance(
void
squared_distance_RT(
const typename K::Point_3 &pt,
const typename K::Segment_3 &seg,
const K& k,
const Homogeneous_tag)
typename K::RT& num,
typename K::RT& den,
const K& k)
{
typename K::Construct_vector_3 construct_vector;
typedef typename K::Vector_3 Vector_3;
typedef typename K::RT RT;
typedef typename K::FT FT;
// assert that the segment is valid (non zero length).
Vector_3 diff = construct_vector(seg.source(), pt);
Vector_3 segvec = construct_vector(seg.source(), seg.target());
RT d = wdot(diff,segvec, k);
if (d <= (RT)0)
return (FT(diff*diff));
RT e = wdot(segvec,segvec, k);
if ( (d * segvec.hw()) > (e * diff.hw()))
return squared_distance(pt, seg.target(), k);
typedef typename K::Vector_3 Vector_3;
typedef typename K::RT RT;
Vector_3 wcr = wcross(segvec, diff, k);
return FT(wcr*wcr)/FT(e * diff.hw() * diff.hw());
typename K::Construct_vector_3 construct_vector;
// assert that the segment is valid (non zero length).
const Vector_3 diff_s = construct_vector(seg.source(), pt);
const Vector_3 segvec = construct_vector(seg.source(), seg.target());
const RT d = wdot(diff_s, segvec, k);
if (d <= RT(0))
{
// this is squared_distance(pt, seg.source())
num = wdot(diff_s, diff_s, k);
den = wmult((K*)0, RT(1), diff_s.hw(), diff_s.hw());
return;
}
const RT e = wdot(segvec, segvec, k);
if (wmult((K*)0 ,d, segvec.hw()) > wmult((K*)0, e, diff_s.hw()))
{
// this is squared_distance(pt, seg.target())
const Vector_3 diff_t = construct_vector(seg.target(), pt);
num = wdot(diff_t, diff_t, k);
den = wmult((K*)0, RT(1), diff_t.hw(), diff_t.hw());
return;
}
// This is an expanded call to squared_distance_to_line_RT() to avoid recomputing 'e'
const Vector_3 wcr = wcross(segvec, diff_s, k);
num = wdot(wcr, wcr, k);
den = wmult((K*)0, e, diff_s.hw(), diff_s.hw());
}
template <class K>
typename K::FT
squared_distance(
const typename K::Point_3 &pt,
const typename K::Segment_3 &seg,
const K& k,
const Cartesian_tag&)
{
typename K::Construct_vector_3 construct_vector;
typedef typename K::Vector_3 Vector_3;
typedef typename K::RT RT;
typedef typename K::FT FT;
// assert that the segment is valid (non zero length).
Vector_3 diff = construct_vector(seg.source(), pt);
Vector_3 segvec = construct_vector(seg.source(), seg.target());
RT d = wdot(diff,segvec, k);
if (d <= (RT)0)
return (FT(diff*diff));
RT e = wdot(segvec,segvec, k);
if (d > e)
return squared_distance(pt, seg.target(), k);
Vector_3 wcr = wcross(segvec, diff, k);
return FT(wcr*wcr)/e;
}
template <class K>
inline
typename K::FT
squared_distance(
const typename K::Point_3 &pt,
const typename K::Segment_3 &seg,
const K& k)
{
typedef typename K::Kernel_tag Tag;
Tag tag;
return squared_distance(pt, seg, k, tag);
}
typedef typename K::Vector_3 Vector_3;
typedef typename K::RT RT;
typedef typename K::FT FT;
typename K::Construct_vector_3 construct_vector = k.construct_vector_3_object();
// assert that the segment is valid (non zero length).
Vector_3 diff = construct_vector(seg.source(), pt);
Vector_3 segvec = construct_vector(seg.source(), seg.target());
RT d = wdot(diff,segvec, k);
if (d <= RT(0))
return (FT(diff*diff));
RT e = wdot(segvec,segvec, k);
if (wmult((K*)0 ,d, segvec.hw()) > wmult((K*)0, e, diff.hw()))
return squared_distance(pt, seg.target(), k);
// This is an expanded call to squared_distance_to_line() to avoid recomputing 'e'
Vector_3 wcr = wcross(segvec, diff, k);
return FT(wcr*wcr) / wmult((K*)0, e, diff.hw(), diff.hw());
}
template <class K>
inline
@ -166,8 +218,6 @@ squared_distance(
return squared_distance(pt, seg, k);
}
template <class K>
typename K::Comparison_result
compare_distance_pssC3(
@ -187,7 +237,7 @@ compare_distance_pssC3(
Vector_3 diff = construct_vector(seg1.source(), pt);
Vector_3 segvec = construct_vector(seg1.source(), seg1.target());
RT d = wdot(diff,segvec, k);
if (d <= (RT)0){
if (d <= RT(0)){
d1 = (FT(diff*diff));
}else{
RT e = wdot(segvec,segvec, k);
@ -204,7 +254,7 @@ compare_distance_pssC3(
Vector_3 diff = construct_vector(seg2.source(), pt);
Vector_3 segvec = construct_vector(seg2.source(), seg2.target());
RT d = wdot(diff,segvec, k);
if (d <= (RT)0){
if (d <= RT(0)){
d2 = (FT(diff*diff));
}else{
RT e = wdot(segvec,segvec, k);
@ -240,7 +290,7 @@ compare_distance_ppsC3(
Vector_3 diff = construct_vector(seg.source(), pt);
Vector_3 segvec = construct_vector(seg.source(), seg.target());
RT d = wdot(diff,segvec, k);
if (d <= (RT)0){
if (d <= RT(0)){
d2 = (FT(diff*diff));
}else{
RT e = wdot(segvec,segvec, k);

141
Distance_3/include/CGAL/squared_distance_3_2.h
View File

@ -218,6 +218,91 @@ on_left_of_triangle_edge(const typename K::Point_3 & pt,
return result;
}
template <class K>
inline void
squared_distance_to_triangle_RT(
const typename K::Point_3 & pt,
const typename K::Point_3 & t0,
const typename K::Point_3 & t1,
const typename K::Point_3 & t2,
bool & inside,
typename K::RT& num,
typename K::RT& den,
const K& k)
{
typename K::Construct_vector_3 vector;
typedef typename K::Vector_3 Vector_3;
const Vector_3 e1 = vector(t0, t1);
const Vector_3 oe3 = vector(t0, t2);
const Vector_3 normal = wcross(e1, oe3, k);
if(normal == NULL_VECTOR) {
// The case normal==NULL_VECTOR covers the case when the triangle
// is colinear, or even more degenerate. In that case, we can
// simply take also the distance to the three segments.
squared_distance_RT(pt, typename K::Segment_3(t2, t0), num, den, k);
typename K::RT num2, den2;
squared_distance_RT(pt, typename K::Segment_3(t1, t2), num2, den2, k);
if(compare_quotients(num2,den2, num,den) == SMALLER) {
num = num2;
den = den2;
}
// should not be needed since at most 2 edges cover the full triangle in the degenerate case
squared_distance_RT(pt, typename K::Segment_3(t0, t1), num2, den2, k);
if(compare_quotients(num2,den2, num,den) == SMALLER){
num = num2;
den = den2;
}
return;
}
const bool b01 = on_left_of_triangle_edge(pt, normal, t0, t1, k);
if(!b01) {
squared_distance_RT(pt, typename K::Segment_3(t0, t1), num, den, k);
return;
}
const bool b12 = on_left_of_triangle_edge(pt, normal, t1, t2, k);
if(!b12) {
squared_distance_RT(pt, typename K::Segment_3(t1, t2), num, den, k);
return;
}
const bool b20 = on_left_of_triangle_edge(pt, normal, t2, t0, k);
if(!b20) {
squared_distance_RT(pt, typename K::Segment_3(t2, t0), num, den, k);
return;
}
// The projection of pt is inside the triangle
inside = true;
squared_distance_to_plane_RT(normal, vector(t0, pt), num, den, k);
}
template <class K>
void
squared_distance_RT(
const typename K::Point_3 & pt,
const typename K::Triangle_3 & t,
typename K::RT& num,
typename K::RT& den,
const K& k)
{
typename K::Construct_vertex_3 vertex;
bool inside = false;
squared_distance_to_triangle_RT(pt,
vertex(t, 0),
vertex(t, 1),
vertex(t, 2),
inside,
num,
den,
k);
}
template <class K>
inline typename K::FT
squared_distance_to_triangle(
@ -234,31 +319,37 @@ squared_distance_to_triangle(
const Vector_3 oe3 = vector(t0, t2);
const Vector_3 normal = wcross(e1, oe3, k);
if(normal != NULL_VECTOR
&& on_left_of_triangle_edge(pt, normal, t0, t1, k)
&& on_left_of_triangle_edge(pt, normal, t1, t2, k)
&& on_left_of_triangle_edge(pt, normal, t2, t0, k))
{
// the projection of pt is inside the triangle
inside = true;
return squared_distance_to_plane(normal, vector(t0, pt), k);
}
else {
// The case normal==NULL_VECTOR covers the case when the triangle
// is colinear, or even more degenerate. In that case, we can
// simply take also the distance to the three segments.
typename K::FT d1 = squared_distance(pt,
typename K::Segment_3(t2, t0),
k);
typename K::FT d2 = squared_distance(pt,
typename K::Segment_3(t1, t2),
k);
typename K::FT d3 = squared_distance(pt,
typename K::Segment_3(t0, t1),
k);
if(normal == NULL_VECTOR) {
// The case normal==NULL_VECTOR covers the case when the triangle
// is colinear, or even more degenerate. In that case, we can
// simply take also the distance to the three segments.
typename K::FT d1 = squared_distance(pt, typename K::Segment_3(t2, t0), k);
typename K::FT d2 = squared_distance(pt, typename K::Segment_3(t1, t2), k);
return (std::min)( (std::min)(d1, d2), d3);
}
// should not be needed since at most 2 edges cover the full triangle in the degenerate case
typename K::FT d3 = squared_distance(pt, typename K::Segment_3(t0, t1), k);
return (std::min)( (std::min)(d1, d2), d3);
}
const bool b01 = on_left_of_triangle_edge(pt, normal, t0, t1, k);
if(!b01) {
return internal::squared_distance(pt, typename K::Segment_3(t0, t1), k);
}
const bool b12 = on_left_of_triangle_edge(pt, normal, t1, t2, k);
if(!b12) {
return internal::squared_distance(pt, typename K::Segment_3(t1, t2), k);
}
const bool b20 = on_left_of_triangle_edge(pt, normal, t2, t0, k);
if(!b20) {
return internal::squared_distance(pt, typename K::Segment_3(t2, t0), k);
}
// The projection of pt is inside the triangle
inside = true;
return squared_distance_to_plane(normal, vector(t0, pt), k);
}
template <class K>
@ -268,7 +359,7 @@ squared_distance(
const typename K::Triangle_3 & t,
const K& k)
{
typename K::Construct_vertex_3 vertex;
typename K::Construct_vertex_3 vertex = k.construct_vertex_3_object();
bool inside = false;
return squared_distance_to_triangle(pt,
vertex(t, 0),

2
Distance_3/package_info/Distance_3/dependencies
View File

@ -2,7 +2,9 @@ Algebraic_foundations
Distance_2
Distance_3
Installation
Interval_support
Kernel_23
Modular_arithmetic
Number_types
Profiling_tools
STL_Extension

11
Distance_3/test/Distance_3/test_distance_3.cpp
View File

@ -112,6 +112,17 @@ struct Test {
{
std::cout << "Point - Triangle\n";
check_squared_distance (p(0, 1, 2), T(p(0, 0, 0), p( 2, 0, 0), p( 0, 2, 0)), 4);
T t(p(0,0,0), p(3,0,0), p(3,3,0));
check_squared_distance (p(-1, -1, 0), t, 2);
check_squared_distance (p(-1, 0, 0), t, 1);
check_squared_distance (p(0, 0, 0), t, 0);
check_squared_distance (p(1, 0, 0), t, 0);
check_squared_distance (p(4, 0, 0), t, 1);
check_squared_distance (p(1, -1, 0), t, 1);
check_squared_distance (p(1, 1, 1), t, 1);
check_squared_distance (p(1, 0, 1), t, 1);
check_squared_distance (p(0, 0, 1), t, 1);
}
void P_Tet()

74
Intersections_3/include/CGAL/Intersections_3/internal/Triangle_3_Sphere_3_do_intersect.h
View File

@ -15,7 +15,7 @@
#include <CGAL/squared_distance_3_2.h>
#include <CGAL/Intersection_traits_3.h>
#include <CGAL/Rational_traits.h>
namespace CGAL {
template <class K>
@ -36,9 +36,15 @@ inline
typename K::Boolean
do_intersect(const typename K::Sphere_3 &sp,
const typename K::Triangle_3 &tr,
const K & /* k */)
const K & k)
{
return squared_distance(sp.center(), tr) <= sp.squared_radius();
typedef typename K::RT RT;
RT num, den;
CGAL::internal::squared_distance_RT(sp.center(), tr, num, den, k);
return ! (compare_quotients<RT>(num, den,
Rational_traits<typename K::FT>().numerator(sp.squared_radius()),
Rational_traits<typename K::FT>().denominator(sp.squared_radius())) == LARGER);
}
template <class K>
@ -46,18 +52,9 @@ inline
typename K::Boolean
do_intersect(const typename K::Triangle_3 &tr,
const typename K::Sphere_3 &sp,
const K & /* k */)
const K & k)
{
return squared_distance(sp.center(), tr) <= sp.squared_radius();
}
template <class K>
inline
typename K::Boolean
do_intersect(const typename K::Sphere_3 &sp,
const typename K::Line_3 &lin,
const K & /* k */)
{
return squared_distance(sp.center(), lin) <= sp.squared_radius();
return do_intersect(sp, tr, k);
}
@ -66,9 +63,26 @@ inline
typename K::Boolean
do_intersect(const typename K::Line_3 &lin,
const typename K::Sphere_3 &sp,
const K & /* k */)
const K & k)
{
return squared_distance(sp.center(), lin) <= sp.squared_radius();
typedef typename K::RT RT;
RT num, den;
CGAL::internal::squared_distance_RT(sp.center(), lin, num, den, k);
return ! (compare_quotients<RT>(num, den,
Rational_traits<typename K::FT>().numerator(sp.squared_radius()),
Rational_traits<typename K::FT>().denominator(sp.squared_radius())) == LARGER);
}
template <class K>
inline
typename K::Boolean
do_intersect(const typename K::Sphere_3 &sp,
const typename K::Line_3 &lin,
const K & k)
{
return do_intersect(lin,sp,k);
}
@ -78,9 +92,15 @@ inline
typename K::Boolean
do_intersect(const typename K::Sphere_3 &sp,
const typename K::Ray_3 &ray,
const K & /* k */)
const K & k)
{
return squared_distance(sp.center(), ray) <= sp.squared_radius();
typedef typename K::RT RT;
RT num, den;
CGAL::internal::squared_distance_RT(sp.center(), ray, num, den, k);
return ! (compare_quotients<RT>(num, den,
Rational_traits<typename K::FT>().numerator(sp.squared_radius()),
Rational_traits<typename K::FT>().denominator(sp.squared_radius())) == LARGER);
}
@ -89,9 +109,9 @@ inline
typename K::Boolean
do_intersect(const typename K::Ray_3 &ray,
const typename K::Sphere_3 &sp,
const K & /* k */)
const K & k)
{
return squared_distance(sp.center(), ray) <= sp.squared_radius();
return do_intersect(sp,ray,k);
}
template <class K>
@ -99,9 +119,15 @@ inline
typename K::Boolean
do_intersect(const typename K::Sphere_3 &sp,
const typename K::Segment_3 &seg,
const K & /* k */)
const K & k)
{
return squared_distance(sp.center(), seg) <= sp.squared_radius();
typedef typename K::RT RT;
RT num, den;
CGAL::internal::squared_distance_RT(sp.center(), seg, num, den, k);
return ! (compare_quotients<RT>(num, den,
Rational_traits<typename K::FT>().numerator(sp.squared_radius()),
Rational_traits<typename K::FT>().denominator(sp.squared_radius())) == LARGER);
}
@ -110,9 +136,9 @@ inline
typename K::Boolean
do_intersect(const typename K::Segment_3 &seg,
const typename K::Sphere_3 &sp,
const K & /* k */)
const K & k)
{
return squared_distance(sp.center(), seg) <= sp.squared_radius();
return do_intersect(sp,seg,k);
}
} // namespace internal

25
Intersections_3/test/Intersections_3/test_intersections_3.cpp
View File

@ -161,13 +161,35 @@ struct Test {
void P_do_intersect()
{
P p(0,0,0), q(1,0,0), r(2,0,0), s(10,10,10);
P p(0,0,0), q(1,0,0), r(2,0,0), s(10,10,10), s2(10,10,0), s3(10,-10,-1), t(3,0,0);
Sph sph(p,1);
Cub cub(p,r);
assert(do_intersect(q,sph));
assert(do_intersect(sph,q));
assert(! do_intersect(s,cub));
assert(! do_intersect(cub,s));
assert(do_intersect(sph, Tr(p,q,s)));
assert(do_intersect(sph, Tr(r,q,s)));
assert(! do_intersect(sph, Tr(r,s,t)));
assert(! do_intersect(sph, Tr(s,s2,s3)));
assert(do_intersect(sph, S(p,r)));
assert(do_intersect(sph, S(q,r)));
assert(do_intersect(sph, S(s,q)));
assert(! do_intersect(sph, S(s,r)));
assert(do_intersect(sph, L(p,r)));
assert(do_intersect(sph, L(q,r)));
assert(do_intersect(sph, L(s,q)));
assert(! do_intersect(sph, L(s,r)));
assert(do_intersect(sph, R(p,r)));
assert(do_intersect(sph, R(q,r)));
assert(do_intersect(sph, R(s,q)));
assert(! do_intersect(sph, R(s,r)));
}
@ -1265,6 +1287,7 @@ struct Test {
Bbox_L();
Bbox_R();
Bbox_Tr();
}
};