Code optimized for x-axis

This commit is contained in:
Laurent Rineau 2012-03-20 16:56:13 +00:00
parent b0cfb5bc1f
commit 0358937f01
3 changed files with 157 additions and 95 deletions

View File

@ -139,8 +139,13 @@ namespace internal {
return do_intersect_bbox_segment_aux<FT, true, false>(
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.xmin(), bbox.ymin(), bbox.zmin(),
bbox.xmax(), bbox.ymax(), bbox.zmax() );
// const CGAL::cpp0x::array<FT, 6> rray = {source.x(), source.y(), source.z(),
// point_on_ray.x(), point_on_ray.y(), point_on_ray.z() };
// return do_intersect_bbox_segment_aux<FT, true, false>
// (rray,
// *reinterpret_cast<const CGAL::cpp0x::array<double, 6>*>(&*bbox.cartesian_begin()) );
}
} // namespace internal

View File

@ -41,16 +41,17 @@ namespace CGAL {
namespace internal {
template <typename FT,
template <typename FT,
bool bounded_0,
bool bounded_1>
// __attribute__ ((noinline))
inline
bool
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 double& bxmin, const double& bymin, const double& bzmin,
const double& bxmax, const double& bymax, const double& bzmax)
{
// The following code encode t1 and t2 by:
// t1 = tmin/dmin
@ -60,28 +61,56 @@ template <typename FT,
// -----------------------------------
// 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;
}
if(bounded_0 && bxmin < px) // tmin < 0 means px is in the x-range of bbox
{
tmin = 0;
dmin = 1;
} else {
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;
dmin = px - qx;
dmax = px - qx;
}
if ( bounded_0 && tmax < FT(0) ) // test t2 < 0, for a segment or a ray
return false;
if ( bounded_1 && tmin > dmin ) // test t1 > 1, for a segment
return false;
if(bounded_0 && px < bxmax) // tmin < 0 means px is in the x-range of bbox
{
tmin = 0;
dmin = 1;
} else {
tmin = px - bxmax;
dmin = px - qx;
}
}
// 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) ) &&
( CGAL::sign(tmin) * CGAL::sign(tmax) ) > 0 )
return false;
// Note: for a segment the condition sign(tmin)*sign(tmax) > 0 has
@ -92,22 +121,6 @@ template <typename FT,
// is a NaN. But the case with NaNs is treated as if the interval
// [t1, t2] was ]-inf, +inf[.
FT dmax = dmin;
// set t1=max(t1, 0), for a segment or a ray
if ( bounded_0 && tmin < FT(0) )
{
tmin = FT(0);
dmin = FT(1);
}
// set t2=min(t2, 1), for a segment
if ( bounded_1 && tmax > dmax )
{
tmax = FT(1);
dmax = FT(1);
}
CGAL_assertion(dmin >= 0);
CGAL_assertion(dmax >= 0);
@ -215,6 +228,40 @@ template <typename FT,
(!bounded_0 || tmax_ >= FT(0)) ); // t1 <= 1 && t2 >= 0
}
template <typename FT,
bool bounded_0,
bool bounded_1>
inline
bool
do_intersect_bbox_segment_aux(const CGAL::cpp0x::array<FT, 6> seg,
const CGAL::cpp0x::array<double, 6> box)
{
const FT& px = seg[0];
const FT& py = seg[1];
const FT& pz = seg[2];
const FT& qx = seg[3];
const FT& qy = seg[4];
const FT& qz = seg[5];
const double& bxmin = box[0];
const double& bymin = box[1];
const double& bzmin = box[2];
const double& bxmax = box[3];
const double& bymax = box[4];
const double& bzmax = box[5];
// for(int i = 0; i < 3; ++i) {
// const int sign = seg[3+i] > seg[i]; // (qx > px)?
// if(bounded_0 && seg[3*(1-sign) + i] > box[3+i]) return false; // segment on the right of bbox
// if(bounded_1 && seg[3*sign + i] < box[i]) return false; // segment on the left of bbox
// }
return do_intersect_bbox_segment_aux<FT, bounded_0, bounded_1>
(px, py, pz,
qy, qy, qz,
bxmin, bymin, bymax,
bxmax, bymax, bzmax);
}
template <class K>
bool do_intersect(const typename K::Segment_3& segment,
const CGAL::Bbox_3& bbox,
@ -229,8 +276,14 @@ template <typename FT,
return do_intersect_bbox_segment_aux<FT, true, true>(
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.xmin(), bbox.ymin(), bbox.zmin(),
bbox.xmax(), bbox.ymax(), bbox.zmax() );
// const CGAL::cpp0x::array<FT, 6> seg = {source.x(), source.y(), source.z(),
// target.x(), target.y(), target.z() };
// return do_intersect_bbox_segment_aux<FT, true, true>
// ( seg,
// *reinterpret_cast<const CGAL::cpp0x::array<double, 6>*>(&*bbox.cartesian_begin()) );
}
template <class K>

View File

@ -62,7 +62,7 @@ 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);
@ -80,6 +80,8 @@ bool test_case(const FT& px, const FT& py, const FT& pz,
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,
@ -91,146 +93,146 @@ bool test_case(const FT& px, const FT& py, const FT& pz,
b &= test_case<K>( px, py, pz,
qx, qy, qz,
bxmin, bymin, bzmin,
bxmax, bymax, bzmax, expected, false);
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, false);
-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, false);
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, false);
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, false);
-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, false);
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, false);
-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, false);
-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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
delta * bxmax, delta * bymax, delta * bzmax, expected, exact_k, exactness_issue,
false, false, false, false, false);
delta = ((unsigned int)-1 / 4);
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,
delta * bxmax, delta * bymax, delta * bzmax, expected, exact_k, exactness_issue,
false, false, false, false, false);
delta = ((unsigned int)-1 / 4);
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,
delta * bxmax, delta * bymax, delta * bzmax, expected, exact_k, exactness_issue,
false, false, false, false, false);
delta = 7;
delta /= 3;
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,
delta * bxmax, delta * bymax, delta * bzmax, expected, exact_k, exactness_issue,
false, false, false, false, false);
delta = 1;
delta /= 10;
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,
delta * bxmax, delta * bymax, delta * bzmax, expected, exact_k, exactness_issue,
false, false, false, false, false);
} else {
using CGAL::do_intersect;
@ -245,6 +247,7 @@ bool test_case(const FT& px, const FT& py, const FT& pz,
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);
@ -260,6 +263,7 @@ bool test_case(const FT& px, const FT& py, const FT& pz,
<< std::endl;
}
}
}
return b;
}
@ -291,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 )
@ -337,7 +341,7 @@ bool intensive_test(bool exact_predicates = true)
ymin <= 5. && ymax >= -5;
b &= test_case<K>(x, ymin, 0.,
x, ymax, 0.,
-5., -5., -5., 5., 5., 5., expected);
-5., -5., -5., 5., 5., 5., expected, exact_predicates);
}
// Test slanted segments
for(double x = -7.; x <= 6.; x+=1.)
@ -349,7 +353,7 @@ bool intensive_test(bool exact_predicates = true)
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);
-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.)
@ -360,13 +364,13 @@ bool intensive_test(bool exact_predicates = true)
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);
-5., -5., -5., 5., 5., 5., expected, exact_predicates);
}
return b;
}
template <class K>
bool test()
bool test(bool exact_kernel = false)
{
// types
typedef typename K::FT FT;
@ -612,12 +616,12 @@ bool test()
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,
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,
0.5, 0.5, 0.5, true, exact_kernel, true,
false, false, false, false, false);
b &= test_case<K>(1., 1., 0.,
@ -658,7 +662,7 @@ bool test()
template <typename K>
bool test_kernel(bool exact_predicates = true, K k = K())
{
bool b = test<K>() &&
bool b = test<K>(exact_predicates) &&
intensive_test<K>(exact_predicates);
test_speed<K>();
return b;
@ -673,16 +677,16 @@ int main()
bool b = test_kernel<CGAL::Simple_cartesian<float> >(false);
std::cout << std::endl << "Testing with Simple_cartesian<double>..." << std::endl ;
b &= test_kernel<CGAL::Simple_cartesian<double> >(false);
b &= test_kernel<CGAL::Simple_cartesian<double> >(true);
std::cout << std::endl << "Testing with Simple_cartesian<Gmpq>..." << std::endl ;
b &= test_kernel<CGAL::Simple_cartesian<CGAL::Gmpq> >(false);
b &= test_kernel<CGAL::Simple_cartesian<CGAL::Gmpq> >(true);
std::cout << std::endl << "Testing with Cartesian<float>..." << std::endl ;
b &= test_kernel<CGAL::Cartesian<float> >(false);
std::cout << std::endl << "Testing with Cartesian<double>..." << std::endl ;
b &= test_kernel<CGAL::Cartesian<double> >(false);
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;