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

Fix Do_intersect_3(Bbox_3, Segment_3) 

svn merge --reintegrate \^/branches/features/Intersection_3-fix_do_intersect_bbox_segment-lrineau
This commit is contained in:
Laurent Rineau 2012-06-25 08:37:32 +00:00
parent 57ffa4fe02 19ab08d6e9
commit eef769cf9a
5 changed files with 724 additions and 500 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
.gitignore vendored
View File

@ -168,6 +168,13 @@ Interpolation/demo/Interpolation/cgal_test_with_cmake
Intersections_3/test/Intersections_3/CMakeLists.txt
Intersections_3/test/Intersections_3/bbox_other_do_intersect_test
Intersections_3/test/Intersections_3/cgal_test_with_cmake
Intersections_3/test/Intersections_3/circle_other
Intersections_3/test/Intersections_3/line_line
Intersections_3/test/Intersections_3/segment_segment
Intersections_3/test/Intersections_3/test_intersections_3
Intersections_3/test/Intersections_3/triangle_3_triangle_3_intersection
Intersections_3/test/Intersections_3/triangle_other
Intersections_3/test/Intersections_3/triangle_other_intersection_test
Jet_fitting_3/examples/Jet_fitting_3/*.exe
Jet_fitting_3/examples/Jet_fitting_3/*.sln
Jet_fitting_3/examples/Jet_fitting_3/*.vcproj

352
Filtered_kernel/include/CGAL/internal/Static_filters/Do_intersect_3.h
View File

@ -28,7 +28,10 @@
#include <CGAL/Profile_counter.h>
#include <CGAL/internal/Static_filters/Static_filter_error.h>
#include <CGAL/internal/Static_filters/tools.h>
#include <cmath>
#include <CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h>
// for CGAL::internal::do_intersect_bbox_segment_aux
#include <iostream>
// inspired from http://cag.csail.mit.edu/~amy/papers/box-jgt.pdf
@ -111,9 +114,6 @@ public:
const Point_3& q = s.target();
double px, py, pz, qx, qy, qz;
double bxmin = b.xmin(), bymin = b.ymin(), bzmin = b.zmin(),
bxmax = b.xmax(), bymax = b.ymax(), bzmax = b.zmax();
if (fit_in_double(get_approx(p).x(), px) && fit_in_double(get_approx(p).y(), py) &&
fit_in_double(get_approx(p).z(), pz) &&
fit_in_double(get_approx(q).x(), qx) && fit_in_double(get_approx(q).y(), qy) &&
@ -121,178 +121,18 @@ public:
{
CGAL_BRANCH_PROFILER_BRANCH_1(tmp);
// -----------------------------------
// treat x coord
// -----------------------------------
double dmin, dmax, tmin, tmax;
if ( qx >= px ) // this is input and needs no epsilon
{
if(px > bxmax) return false; // segment on the right of bbox
if(qx < bxmin) return false; // segment on the left of bbox
const Uncertain<result_type> ub =
do_intersect_bbox_segment_aux
<double,
true, // bounded at t=0
true, // bounded at t=1
true> // do use static filters
(px, py, pz,
qx, qy, qz,
b);
tmax = bxmax - px;
dmax = qx - px;
if ( bxmin < px ) // tmin < 0 means px is in the x-range of bbox
{
tmin = 0;
dmin = 1;
} else {
tmin = bxmin - px;
dmin = dmax;
}
}
else
{
if(qx > bxmax) return false; // segment on the right of bbox
if(px < bxmin) return false; // segment on the left of bbox
tmax = px - bxmin;
dmax = px - qx;
if ( px < bxmax ) // tmin < 0 means px is in the x-range of bbox
{
tmin = 0;
dmin = 1;
} else {
tmin = px - bxmax;
dmin = dmax;
}
}
double m = CGAL::abs(tmin), m2;
m2 = CGAL::abs(tmax); if(m2 > m) { m = m2; }
m2 = CGAL::abs(dmin); if(m2 > m) { m = m2; }
if(m < 7e-294) {
// underflow in the computation of 'error'
return Base::operator()(s,b);
}
const double EPS_1 = 3.55618e-15;
double error = EPS_1 * m;
switch(sign_with_error( tmax - dmax, error)) {
case POSITIVE:
tmax = 1;
dmax = 1;
break;
case NEGATIVE:
break;
default:
// ambiguity of the comparison tmax > dmin
// let's call the exact predicate
return Base::operator()(s,b);
}
// -----------------------------------
// treat y coord
// -----------------------------------
double d_, tmin_, tmax_;
if ( qy >= py ) // this is input and needs no epsilon
{
tmin_ = bymin - py;
tmax_ = bymax - py;
d_ = qy - py;
}
else
{
tmin_ = py - bymax;
tmax_ = py - bymin;
d_ = py - qy;
}
m2 = CGAL::abs(tmin_); if(m2 > m) { m = m2; }
m2 = CGAL::abs(tmax_); if(m2 > m) { m = m2; }
m2 = CGAL::abs(d_); if(m2 > m) { m = m2; }
if(m < 3e-147) {
// underflow in the computation of 'error'
return Base::operator()(s,b);
}
error = EPS_1 * m * m;
if(m > 1e153) { /* sqrt(max_double [hadamard]/2) */
// potential overflow on the computation of 'sign1' and 'sign2'
return Base::operator()(s,b);
}
Sign sign1 = sign_with_error( (d_*tmin) - (dmin*tmax_) , error);
Sign sign2 = sign_with_error( (dmax*tmin_) - (d_*tmax) , error);
if(sign1 == POSITIVE || sign2 == POSITIVE)
return false; // We are *sure* the segment is outside the box, on one
// side or the other.
if(sign1 == ZERO || sign2 == ZERO) {
return Base::operator()(s,b); // We are *unsure*: one *may be*
// positive.
}
// Here we are sure the two signs are negative. We can continue with
// the rest of the function...
// epsilons needed
switch(sign_with_error((dmin*tmin_) - (d_*tmin) , error)) {
case POSITIVE:
tmin = tmin_;
dmin = d_;
break;
case NEGATIVE:
break;
default: // uncertainty
return Base::operator()(s,b);
}
// epsilons needed
switch(sign_with_error((d_*tmax) - (dmax*tmax_) , error)) {
case POSITIVE:
tmax = tmax_;
dmax = d_;
case NEGATIVE:
break;
default: // uncertainty
return Base::operator()(s,b);
}
// -----------------------------------
// treat z coord
// -----------------------------------
if ( qz >= pz ) // this is input and needs no epsilon
{
tmin_ = bzmin - pz;
tmax_ = bzmax - pz;
d_ = qz - pz;
}
else
{
tmin_ = pz - bzmax;
tmax_ = pz - bzmin;
d_ = pz - qz;
}
m2 = CGAL::abs(tmin_); if(m2 > m) { m = m2; }
m2 = CGAL::abs(tmax_); if(m2 > m) { m = m2; }
m2 = CGAL::abs(d_); if(m2 > m) { m = m2; }
// m may have changed
error = EPS_1 * m * m;
if(m > 1e153) { /* sqrt(max_double [hadamard]/2) */
// potential overflow on the computation of 'sign1' and 'sign2'
return Base::operator()(s,b);
}
sign1 = sign_with_error( (dmin*tmax_) - (d_*tmin) , error);
sign2 = sign_with_error( (d_*tmax) - (dmax*tmin_) , error);
if(sign1 == NEGATIVE || sign2 == NEGATIVE) {
return false; // We are *sure* the segment is outside the box, on one
// side or the other.
}
if(sign1 == ZERO || sign2 == ZERO) {
return Base::operator()(s,b); // We are *unsure*: one *may be*
// negative.
}
return true; // We are *sure* the two signs are positive.
if(is_indeterminate(ub)) return Base::operator()(s,b);
else return ub.sup();
}
return Base::operator()(s,b);
}
@ -317,9 +157,6 @@ public:
const Point_3& q = r.second_point();
double px, py, pz, qx, qy, qz;
double bxmin = b.xmin(), bymin = b.ymin(), bzmin = b.zmin(),
bxmax = b.xmax(), bymax = b.ymax(), bzmax = b.zmax();
if (fit_in_double(get_approx(p).x(), px) && fit_in_double(get_approx(p).y(), py) &&
fit_in_double(get_approx(p).z(), pz) &&
fit_in_double(get_approx(q).x(), qx) && fit_in_double(get_approx(q).y(), qy) &&
@ -327,155 +164,18 @@ public:
{
CGAL_BRANCH_PROFILER_BRANCH_1(tmp);
// -----------------------------------
// treat x coord
// -----------------------------------
double dmin, dmax, tmin, tmax;
if ( qx >= px ) // this is input and needs no epsilon
{
if(px > bxmax) return false; // if tmax < 0, ray on the right of bbox
tmax = bxmax - px;
const Uncertain<result_type> ub =
do_intersect_bbox_segment_aux
<double,
true, // bounded at t=0
false,// not bounded at t=1
true> // do use static filters
(px, py, pz,
qx, qy, qz,
b);
dmax = qx - px;
if ( bxmin < px ) // tmin < 0 means px is in the x-range of bbox
{
tmin = 0;
dmin = 1;
} else {
if( px == qx ) return false; // if dmin == 0
tmin = bxmin - px;
dmin = dmax;
}
}
else
{
if(px < bxmin) return false; // if tmax < 0, ray on the left of bbox
tmax = px - bxmin;
dmax = px - qx;
if ( px < bxmax ) // tmin < 0 means px is in the x-range of bbox
{
tmin = 0;
dmin = 1;
} else {
if( px == qx ) return false; // if dmin == 0
tmin = px - bxmax;
dmin = dmax;
}
}
// -----------------------------------
// treat y coord
// -----------------------------------
double d_, tmin_, tmax_;
if ( qy >= py ) // this is input and needs no epsilon
{
tmin_ = bymin - py;
tmax_ = bymax - py;
d_ = qy - py;
}
else
{
tmin_ = py - bymax;
tmax_ = py - bymin;
d_ = py - qy;
}
double m = CGAL::abs(tmin), m2;
m2 = CGAL::abs(tmax); if(m2 > m) { m = m2; }
m2 = CGAL::abs(dmin); if(m2 > m) { m = m2; }
m2 = CGAL::abs(tmin_); if(m2 > m) { m = m2; }
m2 = CGAL::abs(tmax_); if(m2 > m) { m = m2; }
m2 = CGAL::abs(d_); if(m2 > m) { m = m2; }
if(m < 3e-147) {
// underflow in the computation of 'error'
return Base::operator()(r,b);
}
const double EPS_1 = 3.55618e-15;
double error = EPS_1 * m * m;
if(m > 1e153) { /* sqrt(max_double [hadamard]/2) */
// potential overflow on the computation of 'sign1' and 'sign2'
return Base::operator()(r,b);
}
Sign sign1 = sign_with_error( (d_*tmin) - (dmin*tmax_) , error);
Sign sign2 = sign_with_error( (dmax*tmin_) - (d_*tmax) , error);
if(sign1 == POSITIVE || sign2 == POSITIVE)
return false; // We are *sure* the ray is outside the box, on one
// side or the other.
if(sign1 == ZERO || sign2 == ZERO) {
return Base::operator()(r,b); // We are *unsure*: one *may be*
// positive.
}
// Here we are sure the two signs are negative. We can continue with
// the rest of the function...
// epsilons needed
switch(sign_with_error((dmin*tmin_) - (d_*tmin) , error)) {
case POSITIVE:
tmin = tmin_;
dmin = d_;
break;
case NEGATIVE:
break;
default: // uncertainty
return Base::operator()(r,b);
}
// epsilons needed
switch(sign_with_error((d_*tmax) - (dmax*tmax_) , error)) {
case POSITIVE:
tmax = tmax_;
dmax = d_;
case NEGATIVE:
break;
default: // uncertainty
return Base::operator()(r,b);
}
// -----------------------------------
// treat z coord
// -----------------------------------
if ( qz >= pz ) // this is input and needs no epsilon
{
tmin_ = bzmin - pz;
tmax_ = bzmax - pz;
d_ = qz - pz;
}
else
{
tmin_ = pz - bzmax;
tmax_ = pz - bzmin;
d_ = pz - qz;
}
m2 = CGAL::abs(tmin_); if(m2 > m) { m = m2; }
m2 = CGAL::abs(tmax_); if(m2 > m) { m = m2; }
m2 = CGAL::abs(d_); if(m2 > m) { m = m2; }
// m may have changed
error = EPS_1 * m * m;
if(m > 1e153) { /* sqrt(max_double [hadamard]/2) */
// potential overflow on the computation of 'sign1' and 'sign2'
return Base::operator()(r,b);
}
sign1 = sign_with_error( (dmin*tmax_) - (d_*tmin) , error);
sign2 = sign_with_error( (d_*tmax) - (dmax*tmin_) , error);
if(sign1 == NEGATIVE || sign2 == NEGATIVE) {
return false; // We are *sure* the ray is outside the box, on one
// side or the other.
}
if(sign1 == ZERO || sign2 == ZERO) {
return Base::operator()(r,b); // We are *unsure*: one *may be*
// negative.
}
return true; // We are *sure* the two signs are positive.
if(is_indeterminate(ub)) return Base::operator()(r,b);
else return ub.sup();
}
return Base::operator()(r,b);
}

107
Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Ray_3_do_intersect.h
View File

@ -25,105 +25,15 @@
#include <CGAL/Ray_3.h>
#include <CGAL/Bbox_3.h>
#include <CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h>
// for CGAL::internal::do_intersect_bbox_segment_aux
// inspired from http://cag.csail.mit.edu/~amy/papers/box-jgt.pdf
namespace CGAL {
namespace internal {
template <typename FT>
inline
bool
bbox_ray_do_intersect_aux(const FT& px, const FT& py, const FT& pz,
const FT& qx, const FT& qy, const FT& qz,
const FT& bxmin, const FT& bymin, const FT& bzmin,
const FT& bxmax, const FT& bymax, const FT& bzmax)
{
// (px, py, pz) is the source
// -----------------------------------
// treat x coord
// -----------------------------------
FT dmin, tmin, tmax;
if ( qx >= px )
{
tmin = bxmin - px;
tmax = bxmax - px;
dmin = qx - px;
if ( tmax < FT(0) )
return false;
}
else
{
tmin = px - bxmax;
tmax = px - bxmin;
dmin = px - qx;
if ( tmax < FT(0) )
return false;
}
FT dmax = dmin;
if ( tmin > FT(0) )
{
if ( dmin == FT(0) )
return false;
}
else
{
tmin = FT(0);
dmin = FT(1);
}
// -----------------------------------
// treat y coord
// -----------------------------------
FT d_, tmin_, tmax_;
if ( qy >= py )
{
tmin_ = bymin - py;
tmax_ = bymax - py;
d_ = qy - py;
}
else
{
tmin_ = py - bymax;
tmax_ = py - bymin;
d_ = py - qy;
}
if ( (dmin*tmax_) < (d_*tmin) || (dmax*tmin_) > (d_*tmax) )
return false;
if( (dmin*tmin_) > (d_*tmin) )
{
tmin = tmin_;
dmin = d_;
}
if( (dmax*tmax_) < (d_*tmax) )
{
tmax = tmax_;
dmax = d_;
}
// -----------------------------------
// treat z coord
// -----------------------------------
if ( qz >= pz )
{
tmin_ = bzmin - pz;
tmax_ = bzmax - pz;
d_ = qz - pz;
}
else
{
tmin_ = pz - bzmax;
tmax_ = pz - bzmin;
d_ = pz - qz;
}
return ( (dmin*tmax_) >= (d_*tmin) && (dmax*tmin_) <= (d_*tmax) );
}
template <class K>
bool do_intersect(const typename K::Ray_3& ray,
const CGAL::Bbox_3& bbox,
@ -135,11 +45,16 @@ namespace internal {
const Point_3& source = ray.source();
const Point_3& point_on_ray = ray.second_point();
return bbox_ray_do_intersect_aux(
return do_intersect_bbox_segment_aux
<FT,
true, // bounded at t=0
false, // not bounded at t=1
false> // do not use static filters
(
source.x(), source.y(), source.z(),
point_on_ray.x(), point_on_ray.y(), point_on_ray.z(),
FT(bbox.xmin()), FT(bbox.ymin()), FT(bbox.zmin()),
FT(bbox.xmax()), FT(bbox.ymax()), FT(bbox.zmax()) );
bbox
);
}
} // namespace internal

360
Intersections_3/include/CGAL/internal/Intersections_3/Bbox_3_Segment_3_do_intersect.h
View File

@ -1,5 +1,4 @@
// Copyright (c) 2008 ETH Zurich (Switzerland)
// Copyright (c) 2008-2009 INRIA Sophia-Antipolis (France)
// Copyright (c) 2012 GeometryFactory Sarl (France)
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org); you can redistribute it and/or
@ -17,7 +16,7 @@
// $Id$
//
//
// Author(s) : Camille Wormser, Jane Tournois, Pierre Alliez, Stephane Tayeb
// Author(s) : Laurent Rineau
#ifndef CGAL_INTERNAL_INTERSECTIONS_3_BBOX_3_SEGMENT_3_DO_INTERSECT_H
@ -25,104 +24,374 @@
#include <CGAL/Segment_3.h>
#include <CGAL/Bbox_3.h>
#include <CGAL/Kernel/Same_uncertainty.h>
#include <CGAL/assertions.h>
#include <boost/type_traits/is_same.hpp>
// inspired from http://cag.csail.mit.edu/~amy/papers/box-jgt.pdf
// This algorithm intersects the line with the x-, y-, and z-slabs of the
// bounding box, and computes the interval [t1, t2], in the
// parameterization of the line given by the segment (for t=0, that is the
// source of the segment, and for t=1 that is its target), where the line
// intersects the three slabs of the bounding box.
// For a segment, the intersection is non-empty iff
// [t1, t2] intersects [0, 1].
namespace CGAL {
namespace internal {
template <typename FT>
template <typename FT, bool bounded_0, bool use_static_filters = false>
struct Do_intersect_bbox_segment_aux_is_greater
{
typedef typename Same_uncertainty<bool, FT>::type result_type;
void register_new_input_values(const FT& t, const FT& d) {}
void compute_new_error_bound() {}
bool bound_overflow() { return false; }
bool value_might_underflow() { return false; }
static result_type uncertain() {
return true;
}
result_type operator()(const FT& a, const FT& b) const {
return a > b;
}
}; // end struct template Do_intersect_bbox_segment_aux_is_greater
template <typename FT, bool bounded_0>
class Do_intersect_bbox_segment_aux_is_greater<FT, bounded_0, true>
{
double error;
double tmax;
double dmax;
public:
CGAL_static_assertion((boost::is_same<FT, double>::value));
Do_intersect_bbox_segment_aux_is_greater() : error(0.), tmax(0.), dmax(0.) {}
void register_new_input_values(const double& t, const double& d) {
if(bounded_0) {
if(t > tmax) tmax = t;
if(d > dmax) dmax = d;
} else {
const double at = CGAL::abs(t);
const double ad = CGAL::abs(d);
if(at > tmax) tmax = at;
if(ad > dmax) dmax = ad;
}
}
void compute_new_error_bound() {
const double EPS = 8.8872057372592798e-16;
error = tmax * dmax * EPS;
}
bool bound_overflow() {
const double OVERF = 1e153;
return dmax > OVERF || tmax > OVERF;
}
bool value_might_underflow() {
const double UNDERF = 1e-146;
return dmax < UNDERF || tmax < UNDERF;
}
typedef Uncertain<bool> result_type;
static result_type uncertain() {
return result_type::indeterminate();
}
result_type operator()(const FT& a, const FT& b) const {
const FT x = a - b;
if(x > error) return true;
else if(x < -error) return false;
else return uncertain();
}
}; // end specialization Do_intersect_bbox_segment_aux_is_greater<FT, true>
template <typename FT,
bool bounded_0,
bool bounded_1,
bool use_static_filters>
inline
bool
typename Do_intersect_bbox_segment_aux_is_greater
<
FT,
bounded_0,
use_static_filters
>::result_type
do_intersect_bbox_segment_aux(
const FT& px, const FT& py, const FT& pz,
const FT& qx, const FT& qy, const FT& qz,
const FT& bxmin, const FT& bymin, const FT& bzmin,
const FT& bxmax, const FT& bymax, const FT& bzmax)
const Bbox_3& bbox)
{
const double& bxmin = bbox.xmin();
const double& bymin = bbox.ymin();
const double& bzmin = bbox.zmin();
const double& bxmax = bbox.xmax();
const double& bymax = bbox.ymax();
const double& bzmax = bbox.zmax();
// The following code encode t1 and t2 by:
// t1 = tmin/dmin
// t2 = tmax/dmax
// For the first lines, dmax==dmin and is not explicitly defined.
// -----------------------------------
// treat x coord
// -----------------------------------
FT dmin, tmin, tmax;
FT dmin, tmin, tmax, dmax;
if ( qx >= px )
{
tmin = bxmin - px;
if(bounded_0 && px > bxmax) return false; // segment on the right of bbox
if(bounded_1 && qx < bxmin) return false; // segment on the left of bbox
if(bounded_1 && bxmax > qx) {
tmax = 1;
dmax = 1;
} else {
tmax = bxmax - px;
dmax = qx - px;
}
tmin = bxmin - px;
dmin = qx - px;
}
else
{
tmin = px - bxmax;
if(bounded_1 && qx > bxmax) return false; // segment on the right of bbox
if(bounded_0 && px < bxmin) return false; // segment on the left of bbox
if(bounded_1 && bxmin < qx) {
tmax = 1;
dmax = 1;
} else {
tmax = px - bxmin;
dmax = px - qx;
}
tmin = px - bxmax;
dmin = px - qx;
}
if ( tmax < FT(0) || tmin > dmin )
return false;
if(bounded_0) tmin = (CGAL::max)(FT(0), tmin);
FT dmax = dmin;
if ( tmin < FT(0) )
// If the query is vertical for x, then check its x-coordinate is in
// the x-slab.
if( (px == qx) && // <=> (dmin == 0)
(! (bounded_0 && bounded_1) ) ) // do not check for a segment
{
tmin = FT(0);
dmin = FT(1);
if(px > bxmax || px < bxmin) return false;
// Note: for a segment the condition has already been tested by the two
// previous tests tmax<0 || tmin>dmin (with dmin==0).
}
if ( tmax > dmax )
{
tmax = FT(1);
dmax = FT(1);
// If dmin == 0, at this point, [t1, t2] == ]-inf, +inf[, or t1 or t2
// is a NaN. But the case with NaNs is treated as if the interval
// [t1, t2] was ]-inf, +inf[.
CGAL_assertion(dmin >= 0);
CGAL_assertion(dmax >= 0);
if(bounded_0) {
CGAL_assertion(tmin >= 0);
CGAL_assertion(tmax >= 0);
}
// -----------------------------------
// treat y coord
// -----------------------------------
FT d_, tmin_, tmax_;
FT dymin, tymin, tymax, dymax;
if ( qy >= py )
{
tmin_ = bymin - py;
tmax_ = bymax - py;
d_ = qy - py;
if(bounded_0 && py > bymax) return false; // segment on the right of bbox
if(bounded_1 && qy < bymin) return false; // segment on the left of bbox
if(bounded_1 && bymax > qy) {
tymax = 1;
dymax = 1;
} else {
tymax = bymax - py;
dymax = qy - py;
}
tymin = bymin - py;
dymin = qy - py;
}
else
{
tmin_ = py - bymax;
tmax_ = py - bymin;
d_ = py - qy;
if(bounded_1 && qy > bymax) return false; // segment on the right of bbox
if(bounded_0 && py < bymin) return false; // segment on the left of bbox
if(bounded_1 && bymin < qy) {
tymax = 1;
dymax = 1;
} else {
tymax = py - bymin;
dymax = py - qy;
}
if ( (dmin*tmax_) < (d_*tmin) || (dmax*tmin_) > (d_*tmax) )
return false;
tymin = py - bymax;
dymin = py - qy;
}
if( (dmin*tmin_) > (d_*tmin) )
if(bounded_0) tymin = (CGAL::max)(FT(0), tymin);
// If the query is vertical for y, then check its y-coordinate is in
// the y-slab.
if( (py == qy) && // <=> (dmin == 0)
(! (bounded_0 && bounded_1) ) ) // do not check for a segment
{
tmin = tmin_;
dmin = d_;
if(py > bymax || py < bymin) return false;
}
if( (dmax*tmax_) < (d_*tmax) )
{
tmax = tmax_;
dmax = d_;
// If dmin == 0, at this point, [t1, t2] == ]-inf, +inf[, or t1 or t2
// is a NaN. But the case with NaNs is treated as if the interval
// [t1, t2] was ]-inf, +inf[.
CGAL_assertion(dymin >= 0);
CGAL_assertion(dymax >= 0);
if(bounded_0) {
CGAL_assertion(tymin >= 0);
CGAL_assertion(tymax >= 0);
}
// -----------------------------------
// treat z coord
// -----------------------------------
FT dzmin, tzmin, tzmax, dzmax;
if ( qz >= pz )
{
tmin_ = bzmin - pz;
tmax_ = bzmax - pz;
d_ = qz - pz;
if(bounded_0 && pz > bzmax) return false; // segment on the right of bbox
if(bounded_1 && qz < bzmin) return false; // segment on the left of bbox
if(bounded_1 && bzmax > qz) {
tzmax = 1;
dzmax = 1;
} else {
tzmax = bzmax - pz;
dzmax = qz - pz;
}
tzmin = bzmin - pz;
dzmin = qz - pz;
}
else
{
tmin_ = pz - bzmax;
tmax_ = pz - bzmin;
d_ = pz - qz;
if(bounded_1 && qz > bzmax) return false; // segment on the right of bbox
if(bounded_0 && pz < bzmin) return false; // segment on the left of bbox
if(bounded_1 && bzmin < qz) {
tzmax = 1;
dzmax = 1;
} else {
tzmax = pz - bzmin;
dzmax = pz - qz;
}
return ( (dmin*tmax_) >= (d_*tmin) && (dmax*tmin_) <= (d_*tmax) );
tzmin = pz - bzmax;
dzmin = pz - qz;
}
if(bounded_0) tzmin = (CGAL::max)(FT(0), tzmin);
// If the query is vertical for z, then check its z-coordinate is in
// the z-slab.
if( (pz == qz) && // <=> (dmin == 0)
(! (bounded_0 && bounded_1) ) ) // do not check for a segment
{
if(pz > bzmax || pz < bzmin) return false;
}
// If dmin == 0, at this point, [t1, t2] == ]-inf, +inf[, or t1 or t2
// is a NaN. But the case with NaNs is treated as if the interval
// [t1, t2] was ]-inf, +inf[.
CGAL_assertion(dzmin >= 0);
CGAL_assertion(dzmax >= 0);
if(bounded_0) {
CGAL_assertion(tzmin >= 0);
CGAL_assertion(tzmax >= 0);
}
typedef Do_intersect_bbox_segment_aux_is_greater
<FT,
bounded_0,
use_static_filters
> Is_greater;
typedef typename Is_greater::result_type Is_greater_value;
Is_greater is_greater;
is_greater.register_new_input_values(tmin, dmin);
is_greater.register_new_input_values(tymin, dymin);
is_greater.register_new_input_values(tmax, dmax);
is_greater.register_new_input_values(tymax, dymax);
is_greater.compute_new_error_bound();
if(is_greater.bound_overflow() || is_greater.value_might_underflow())
return Is_greater::uncertain();
// If t1 > tymax/dymax || tymin/dymin > t2, return false.
if( py != qy && px != qx ) { // dmin > 0, dymax >0, dmax > 0, dymin > 0
const Is_greater_value b1 = is_greater(dymax* tmin, dmin*tymax);
if(possibly(b1)) return !b1; // if(is_greater) return false; // or uncertain
const Is_greater_value b2 = is_greater( dmax*tymin, dymin* tmax);
if(possibly(b2)) return !b2;
}
Is_greater_value b = Is_greater_value();
// If tymin/dymin > t1, set t1 = tymin/dymin.
if( (px == qx) || // <=> (dmin == 0)
( (py != qy) && // <=> (dymin > 0)
certainly(b = is_greater( dmin*tymin, dymin* tmin)) ) )
{
tmin = tymin;
dmin = dymin;
}
if(is_indeterminate(b)) return b; // Note that the default-constructed
// Is_greater_value cannot be
// indeterminate.
// If tymax/dymax < t2, set t2 = tymax/dymax.
if( (px == qx) || // <=> (dmax > 0)
( (py != qy) && // <=> dymax > 0
certainly(b = is_greater(dymax* tmax, dmax*tymax)) ) )
{
tmax = tymax;
dmax = dymax;
}
if(is_indeterminate(b)) return b;
CGAL_assertion(dmin >= 0);
CGAL_assertion(dmax >= 0);
// If t1 > tzmax || tzmin > t2, return false.
if( (px != qx ||
py != qy ) &&
(pz != qz) ) // dmin > 0, dmax > 0, dzmax > 0, dzmin > 0
{
is_greater.register_new_input_values(tzmin, dzmin);
is_greater.register_new_input_values(tzmax, dzmax);
is_greater.compute_new_error_bound();
if(is_greater.bound_overflow() || is_greater.value_might_underflow())
return Is_greater::uncertain();
const Is_greater_value b1 = is_greater(dzmax* tmin, dmin*tzmax);
if(possibly(b1)) return !b1; // if(is_greater) return false; // or uncertain
const Is_greater_value b2 = is_greater( dmax*tzmin, dzmin* tmax);
if(possibly(b2)) return !b2; // if(is_greater) return false; // or uncertain
}
return true;
}
template <class K>
@ -136,11 +405,10 @@ namespace internal {
const Point_3& source = segment.source();
const Point_3& target = segment.target();
return do_intersect_bbox_segment_aux(
return do_intersect_bbox_segment_aux<FT, true, true, false>(
source.x(), source.y(), source.z(),
target.x(), target.y(), target.z(),
FT(bbox.xmin()), FT(bbox.ymin()), FT(bbox.zmin()),
FT(bbox.xmax()), FT(bbox.ymax()), FT(bbox.zmax()) );
bbox);
}
template <class K>

392
Intersections_3/test/Intersections_3/bbox_other_do_intersect_test.cpp
View File

@ -37,6 +37,7 @@
#include <iomanip>
#include <boost/math/special_functions/next.hpp> // for nextafter
double random_in(const double a,
@ -61,17 +62,211 @@ template <class T>
bool test_aux(const T& t,
const std::string& name,
const CGAL::Bbox_3& bbox,
bool expected)
bool expected, bool exact_predicates = false)
{
bool b = CGAL::do_intersect(t,bbox);
if ( b != expected )
std::cout << "ERROR: do_intersect(" << name
std::cerr << "ERROR: do_intersect(" << name
<< ") did not answer the expected result !" << std::endl;
return (b == expected);
}
template <class K,
class FT>
bool test_case(const FT& px, const FT& py, const FT& pz,
const FT& qx, const FT& qy, const FT& qz,
FT bxmin, FT bymin, FT bzmin,
FT bxmax, FT bymax, FT bzmax,
const bool expected,
const bool exact_k = false,
const bool exactness_issue = false,
const bool change_signs = true,
const bool swap_coords = true,
const bool opposite_seg = true,
const bool translate = true,
const bool scale = true)
{
bool b = true;
if(change_signs) {
b &= test_case<K>( px, py, pz,
qx, qy, qz,
bxmin, bymin, bzmin,
bxmax, bymax, bzmax, expected, exact_k, exactness_issue, false);
b &= test_case<K>(-px, py, pz,
-qx, qy, qz,
-bxmin, bymin, bzmin,
-bxmax, bymax, bzmax, expected, exact_k, exactness_issue, false);
b &= test_case<K>( px, -py, pz,
qx, -qy, qz,
bxmin, -bymin, bzmin,
bxmax, -bymax, bzmax, expected, exact_k, exactness_issue, false);
b &= test_case<K>( px, py, -pz,
qx, qy, -qz,
bxmin, bymin, -bzmin,
bxmax, bymax, -bzmax, expected, exact_k, exactness_issue, false);
b &= test_case<K>(-px, -py, pz,
-qx, -qy, qz,
-bxmin, -bymin, bzmin,
-bxmax, -bymax, bzmax, expected, exact_k, exactness_issue, false);
b &= test_case<K>( px, -py, -pz,
qx, -qy, -qz,
bxmin, -bymin, -bzmin,
bxmax, -bymax, -bzmax, expected, exact_k, exactness_issue, false);
b &= test_case<K>(-px, py, -pz,
-qx, qy, -qz,
-bxmin, bymin, -bzmin,
-bxmax, bymax, -bzmax, expected, exact_k, exactness_issue, false);
b &= test_case<K>(-px, -py, -pz,
-qx, -qy, -qz,
-bxmin, -bymin, -bzmin,
-bxmax, -bymax, -bzmax, expected, exact_k, exactness_issue, false);
} else if(swap_coords) {
// xyz
b &= test_case<K>( px, py, pz,
qx, qy, qz,
bxmin, bymin, bzmin,
bxmax, bymax, bzmax, expected, exact_k, exactness_issue,
false, false);
// xzy
b &= test_case<K>( px, pz, py,
qx, qz, qy,
bxmin, bzmin, bymin,
bxmax, bzmax, bymax, expected, exact_k, exactness_issue,
false, false);
// yxz
b &= test_case<K>( py, px, pz,
qy, qx, qz,
bymin, bxmin, bzmin,
bymax, bxmax, bzmax, expected, exact_k, exactness_issue,
false, false);
// zxy
b &= test_case<K>( pz, px, py,
qz, qx, qy,
bzmin, bxmin, bymin,
bzmax, bxmax, bymax, expected, exact_k, exactness_issue,
false, false);
// yzx
b &= test_case<K>( py, pz, px,
qy, qz, qx,
bymin, bzmin, bxmin,
bymax, bzmax, bxmax, expected, exact_k, exactness_issue,
false, false);
// zyx
b &= test_case<K>( pz, py, px,
qz, qy, qx,
bzmin, bymin, bxmin,
bzmax, bymax, bxmax, expected, exact_k, exactness_issue,
false, false);
} else if(opposite_seg) {
b &= test_case<K>(px, py, pz,
qx, qy, qz,
bxmin, bymin, bzmin,
bxmax, bymax, bzmax, expected, exact_k, exactness_issue,
false, false, false);
b &= test_case<K>(qx, qy, qz,
px, py, pz,
bxmin, bymin, bzmin,
bxmax, bymax, bzmax, expected, exact_k, exactness_issue,
false, false, false);
} else if(translate) {
b &= test_case<K>(px, py, pz,
qx, qy, qz,
bxmin, bymin, bzmin,
bxmax, bymax, bzmax, expected, exact_k, exactness_issue,
false, false, false, false);
FT dx = 10, dy = 20, dz = 30;
b &= test_case<K>(dx + px, dy + py, dz + pz,
dx + qx, dy + qy, dz + qz,
dx + bxmin, dy + bymin, dz + bzmin,
dx + bxmax, dy + bymax, dz + bzmax, expected, exact_k, exactness_issue,
false, false, false, false);
dx = (1 >> 10), dy = dx, dz = dx;
b &= test_case<K>(dx + px, dy + py, dz + pz,
dx + qx, dy + qy, dz + qz,
dx + bxmin, dy + bymin, dz + bzmin,
dx + bxmax, dy + bymax, dz + bzmax, expected, exact_k, exactness_issue,
false, false, false, false);
dx = -(1 >> 10), dy = dx, dz = dx;
b &= test_case<K>(dx + px, dy + py, dz + pz,
dx + qx, dy + qy, dz + qz,
dx + bxmin, dy + bymin, dz + bzmin,
dx + bxmax, dy + bymax, dz + bzmax, expected, exact_k, exactness_issue,
false, false, false, false);
} else if(scale) {
b &= test_case<K>(px, py, pz,
qx, qy, qz,
bxmin, bymin, bzmin,
bxmax, bymax, bzmax, expected, exact_k, exactness_issue,
false, false, false, false, false);
FT delta = 9;
b &= test_case<K>(delta * qx, delta * qy, delta * qz,
delta * px, delta * py, delta * pz,
delta * bxmin, delta * bymin, delta * bzmin,
delta * bxmax, delta * bymax, delta * bzmax, expected, exact_k, exactness_issue,
false, false, false, false, false);
delta = (1 << 10);
b &= test_case<K>(delta * qx, delta * qy, delta * qz,
delta * px, delta * py, delta * pz,
delta * bxmin, delta * bymin, delta * bzmin,
delta * bxmax, delta * bymax, delta * bzmax, expected, exact_k, exactness_issue,
false, false, false, false, false);
delta = (1 << 10);
delta = 1/delta;
b &= test_case<K>(delta * qx, delta * qy, delta * qz,
delta * px, delta * py, delta * pz,
delta * bxmin, delta * bymin, delta * bzmin,
delta * bxmax, delta * bymax, delta * bzmax, expected, exact_k, exactness_issue,
false, false, false, false, false);
delta = 7;
delta /= 2;
b &= test_case<K>(delta * qx, delta * qy, delta * qz,
delta * px, delta * py, delta * pz,
delta * bxmin, delta * bymin, delta * bzmin,
delta * bxmax, delta * bymax, delta * bzmax, expected, exact_k, exactness_issue,
false, false, false, false, false);
delta = 1;
delta /= 8;
b &= test_case<K>(delta * qx, delta * qy, delta * qz,
delta * px, delta * py, delta * pz,
delta * bxmin, delta * bymin, delta * bzmin,
delta * bxmax, delta * bymax, delta * bzmax, expected, exact_k, exactness_issue,
false, false, false, false, false);
} else {
using CGAL::do_intersect;
using CGAL::Bbox_3;
typedef typename K::Point_3 Point_3;
typedef typename K::Segment_3 Segment_3;
if(bxmin > bxmax) std::swap(bxmin, bxmax);
if(bymin > bymax) std::swap(bymin, bymax);
if(bzmin > bzmax) std::swap(bzmin, bzmax);
if(do_intersect(Bbox_3(bxmin, bymin, bzmin,
bxmax, bymax, bzmax),
Segment_3(Point_3(px, py, pz),
Point_3(qx, qy, qz))) != expected)
{
if(!exactness_issue || exact_k) {
b = false;
CGAL::set_pretty_mode(std::cerr);
std::cerr.precision(17);
std::cerr << "Wrong result for do_intersect("
<< Bbox_3(bxmin, bymin, bzmin,
bxmax, bymax, bzmax)
<< ",\n"
<< " "
<< Segment_3(Point_3(px, py, pz),
Point_3(qx, qy, qz))
<< ")\n"
<< " it should have been " << std::boolalpha << expected
<< std::endl;
}
}
}
return b;
}
template <class T>
void speed(const std::string& name)
{
@ -100,7 +295,7 @@ void speed(const std::string& name)
CGAL::Timer timer;
timer.start();
std::size_t success = 0;
while ( timer.time() < 0.1 )
while ( timer.time() < 5. )
{
for ( typename std::vector<T>::iterator it = segment_vector.begin();
it != segment_vector.end() ; ++it )
@ -132,7 +327,50 @@ void test_speed()
}
template <class K>
bool test()
bool intensive_test(bool exact_predicates = true)
{
bool b = true;
// Test vertical segments
for(double x = 4.; x <= 7.; x+=1.)
for(double ymin = 0.; ymin <= 7.; ymin+=1.)
for(double ymax = ymin; ymax <= 7.; ymax+=1.)
{
const bool expected =
x >= -5. && x<= 5. &&
ymin <= 5. && ymax >= -5;
b &= test_case<K>(x, ymin, 0.,
x, ymax, 0.,
-5., -5., -5., 5., 5., 5., expected, exact_predicates);
}
// Test slanted segments
for(double x = -7.; x <= 6.; x+=1.)
for(double y = -1.; y <= 6.; y+=1.)
{
const bool expected =
x >= -6. && x <= 5. &&
y >= -6. && y <= 5. &&
y <= x + 10. && y >= x - 10.;
b &= test_case<K>(x, y, 0.,
x + 1., y + 1., 0.,
-5., -5., -5., 5., 5., 5., expected, exact_predicates);
}
for(double x = -9.; x <= 6.; x+=1.)
for(double y = -3.; y <= 6.; y+=1.)
{
const bool expected =
x >= -8. && x <= 5. &&
y >= -7. && y <= 5. &&
3 * y <= 2 * x + 25. && 3 * y >= 2 * x - 25.;
b &= test_case<K>(x, y, 0.,
x + 3., y + 2., 0.,
-5., -5., -5., 5., 5., 5., expected, exact_predicates);
}
return b;
}
template <class K>
bool test(bool exact_kernel = false)
{
// types
typedef typename K::FT FT;
@ -190,6 +428,46 @@ bool test()
b &= test_aux(sEC,"sEC",bbox,false);
CGAL::Bbox_3 bbox_elem_1834(2, 1.81818, 0.166666,
2.18182, 2.18182, 0.333333);
Point source_1834(2, 2, 1);
Point target_1834(2, 2, 0.75);
Segment segment_query_1834( source_1834, target_1834 );
Point source_1834a(2, 2, 0.5);
Point target_1834a(2, 2, 0.75);
Segment segment_query_1834a( source_1834a, target_1834a );
b &= test_aux(segment_query_1834,
"segment_query_1834", bbox_elem_1834, false);
b &= test_aux(segment_query_1834.opposite(),
"segment_query_1834.opposite()", bbox_elem_1834, false);
b &= test_aux(segment_query_1834a,
"segment_query_1834a", bbox_elem_1834, false);
b &= test_aux(segment_query_1834a.opposite(),
"segment_query_1834a.opposite()", bbox_elem_1834, false);
b &= test_case<K>(2., 2., 1.,
2., 2., 0.75,
1.81818, 2., 0.,
2., 2.18182, 0.333333, false);
b &= test_case<K>(2., 2., 0.5,
2., 2., 0.75,
1.81818, 2., 0.,
2., 2.18182, 0.333333, false);
CGAL::Bbox_3 bbox_elem_1834b(1.81818, 2, 0,
2, 2.18182, 0.333333);
Segment segment_query_1834b(Point(2, 2, 1),
Point(2, 2, 0.75));
b &= test_aux(segment_query_1834b,
"segment_query_1834b", bbox_elem_1834b, false);
b &= test_aux(segment_query_1834b.opposite(),
"segment_query_1834b.opposite()", bbox_elem_1834b, false);
b &= test_case<K>(2., 2., 1.,
2., 2., 0.75,
1.81818, 2., 0.,
2., 2.18182, 0.333333, false);
Ray r12(p1,p2);
Ray r13(p1,p3);
Ray r14(p1,p4);
@ -309,21 +587,84 @@ bool test()
Line line3(seg3);
Line line4(seg4);
test_aux(seg2, "seg2", bbox2, false);
test_aux(seg3, "seg3", bbox3, true);
test_aux(seg4, "seg4", bbox4, false);
b &= test_aux(seg2, "seg2", bbox2, false);
b &= test_aux(seg3, "seg3", bbox3, true);
b &= test_aux(seg4, "seg4", bbox4, false);
test_aux(ray2, "ray2", bbox2, false);
test_aux(ray3, "ray3", bbox3, true);
test_aux(ray4, "ray4", bbox4, false);
b &= test_aux(ray2, "ray2", bbox2, false);
b &= test_aux(ray3, "ray3", bbox3, true);
b &= test_aux(ray4, "ray4", bbox4, false);
test_aux(line2, "line2", bbox2, false);
test_aux(line3, "line3", bbox3, true);
test_aux(line4, "line4", bbox4, false);
b &= test_aux(line2, "line2", bbox2, false);
b &= test_aux(line3, "line3", bbox3, true);
b &= test_aux(line4, "line4", bbox4, false);
// Use do_intersect(bbox,bbox)
CGAL::do_intersect(bbox2,bbox4);
b &= test_case<K>(1., 1., 0.,
1., 1., 1.,
0., 0., 0.,
1., 1., 1., true);
b &= test_case<K>(0.5, 0.5, -0.5,
0.5, 0.5, 0.5,
-0.5, -0.5, -0.5,
0.5, 0.5, 0.5, true);
float f = 0.5f;
double d = boost::math::nextafter(f, f+1);
double d2 = boost::math::nextafter(f, f-1);
b &= test_case<K>(d, 0.5, -0.5,
d, 0.5, 0.5,
-0.5, -0.5, -0.5,
0.5, 0.5, 0.5, false, exact_kernel, true,
false, false, false, false, false);
b &= test_case<K>(d2, 0.5, -0.5,
d, 0.5, 0.5,
-0.5, -0.5, -0.5,
0.5, 0.5, 0.5, true, exact_kernel, true,
false, false, false, false, false);
b &= test_case<K>(1., 1., 0.,
2., 2., 2.,
0., 0., 0.,
1., 1., 1., true);
b &= test_case<K>(1., 1., 1.,
1., 1., 1.,
0., 0., 0.,
1., 1., 1., true);
b &= test_case<K>(0.9, 0.9, 0.9,
0.9, 0.9, 0.9,
0., 0., 0.,
1., 1., 1., true);
b &= test_case<K>(0., 0., 0.,
0., 0., 0.,
0., 0., 0.,
1., 1., 1., true);
b &= test_case<K>(0.1, 0., 0.1,
0.1, 0., 0.1,
0., 0., 0.,
1., 1., 1., true);
b &= test_case<K>(0.1, -0.1, 0.1,
0.1, -0.1, 0.1,
0., 0., 0.,
1., 1., 1., false);
b &= test_case<K>(0.1, 0.1, 0.1,
0.1, 0.1, 0.1,
0., 0., 0.,
1., 1., 1., true);
b &= test_case<K>(1., 1., 1.1,
1., 1., 1.1,
0., 0., 0.,
1., 1., 1., false);
return b;
}
template <typename K>
bool test_kernel(bool exact_predicates = true, K k = K())
{
bool b = test<K>(exact_predicates) &&
intensive_test<K>(exact_predicates);
test_speed<K>();
return b;
}
@ -332,38 +673,31 @@ int main()
srand(0);
std::cout << std::setprecision(5);
bool b;
std::cout << "Testing with Simple_cartesian<float>..." << std::endl ;
bool b = test<CGAL::Simple_cartesian<float> >();
test_speed<CGAL::Simple_cartesian<float> >();
b = test_kernel<CGAL::Simple_cartesian<float> >(false);
std::cout << std::endl << "Testing with Simple_cartesian<double>..." << std::endl ;
b &= test<CGAL::Simple_cartesian<double> >();
test_speed<CGAL::Simple_cartesian<double> >();
b &= test_kernel<CGAL::Simple_cartesian<double> >(true);
std::cout << std::endl << "Testing with Simple_cartesian<Gmpq>..." << std::endl ;
b &= test<CGAL::Simple_cartesian<CGAL::Gmpq> >();
test_speed<CGAL::Simple_cartesian<CGAL::Gmpq> >();
b &= test_kernel<CGAL::Simple_cartesian<CGAL::Gmpq> >(true);
std::cout << std::endl << "Testing with Cartesian<float>..." << std::endl ;
b &= test<CGAL::Cartesian<float> >();
test_speed<CGAL::Cartesian<float> >();
b &= test_kernel<CGAL::Cartesian<float> >(false);
std::cout << std::endl << "Testing with Cartesian<double>..." << std::endl ;
b &= test<CGAL::Cartesian<double> >();
test_speed<CGAL::Cartesian<double> >();
b &= test_kernel<CGAL::Cartesian<double> >(true);
std::cout << std::endl << "Testing with Filtered_kernel<Simple_cartesian<double> > without static filters..." << std::endl ;
typedef CGAL::Filtered_kernel<CGAL::Simple_cartesian<double>, false> Fk_no_static;
b &= test<Fk_no_static>();
test_speed<Fk_no_static>();
b &= test_kernel<Fk_no_static>();
std::cout << std::endl << "Testing with Exact_predicates_inexact_constructions_kernel..." << std::endl ;
b &= test<CGAL::Exact_predicates_inexact_constructions_kernel>();
test_speed<CGAL::Exact_predicates_inexact_constructions_kernel>();
b &= test_kernel<CGAL::Exact_predicates_inexact_constructions_kernel>();
std::cout << std::endl << "Testing with Exact_predicates_exact_constructions_kernel..." << std::endl ;
b &= test<CGAL::Exact_predicates_exact_constructions_kernel>();
test_speed<CGAL::Exact_predicates_exact_constructions_kernel>();
b &= test_kernel<CGAL::Exact_predicates_exact_constructions_kernel>();
if ( b )
return EXIT_SUCCESS;