Benchmark for SNC_intersection

Intersections_benchmark.cpp

#define CGAL_NO_ASSERTIONS 1
#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <CGAL/Kernel/global_functions.h>
#include <CGAL/point_generators_3.h>
#include <CGAL/Timer.h>
#include <CGAL/Nef_polyhedron_3.h>
#include <random>
#include <algorithm>

using K = CGAL::Exact_predicates_exact_constructions_kernel;

using Ray_3 = CGAL::Ray_3<K>;
using Segment_3 = CGAL::Segment_3<K>;
using Point_3 = CGAL::Point_3<K>;
using Line_3 = CGAL::Line_3<K>;
using Vector_3 = CGAL::Vector_3<K>;
using Plane_3 = CGAL::Plane_3<K>;
using Object = CGAL::Object;
using Transform = CGAL::Aff_transformation_3<K>;

CGAL::SNC_intersection<CGAL::Nef_polyhedron_3<K>::SNC_structure> is;

bool ray_does_intersect_internally_old( const Ray_3& s1,
                                        const Segment_3& s2,
                                        Point_3& p) {
    return is.does_intersect_internally(s1,s2,p);
}

bool ray_does_intersect_internally( const Ray_3& s1,
                                    const Segment_3& s2,
                                    Point_3& p) {
    if(s1.has_on(s2.source()) || s1.has_on(s2.target()) ||
       s2.has_on(s1.source()))
        return false;
    Object o = intersection(s1, s2);
    return CGAL::assign(p,o);
}

bool segment_does_intersect_internally_old( const Segment_3& s1,
                                            const Segment_3& s2,
                                            Point_3& p) {
    return is.does_intersect_internally(s1,s2,p);
}

bool segment_does_intersect_internally( const Segment_3& s1,
                                        const Segment_3& s2,
                                        Point_3& p) {
    if(s1.has_on(s2.source()) || s1.has_on(s2.target()) ||
       s2.has_on(s1.source()) || s2.has_on(s1.target()))
        return false;
    Object o = intersection(s1, s2);
    return CGAL::assign(p,o);
}

template <typename Point>
Segment_3 random_axis_of_sphere(CGAL::Random_points_on_sphere_3<Point>& rs) {
    const Point_3 p1(*rs++),p2(CGAL::ORIGIN-(p1-CGAL::ORIGIN));
    return Segment_3(p1,p2);
}

template <typename Point>
Segment_3 random_chord_of_sphere(CGAL::Random_points_on_sphere_3<Point>& rs) {
    return Segment_3(*rs++,*rs++);
}

template <typename Function1, typename Function2>
void timing_test(std::vector<std::pair<Segment_3,Segment_3>> segments,
                 Function1 ray_test,
                 Function2 segment_test)
{
    CGAL::Timer ts;
    const size_t iterations=200;
    Point_3 p;
    ts.start();
    for(size_t n=0; n<iterations; ++n)
    {
        for(auto&& a: segments) {
            Ray_3 r(a.first.source(),a.first.to_vector());
            ray_test(r,a.second,p);
        }
    }
    ts.stop();
    std::cout << ts.time() << "  |  ";
    ts.start();
    for(size_t n=0; n<iterations; ++n)
    {
        for(auto&& a: segments) {
            segment_test(a.first,a.second,p);
        }
    }
    ts.stop();
    std::cout << ts.time() << std::endl;
}

int main()
{
    CGAL::Random_points_on_sphere_3<Point_3> rs;

    const size_t count=1000;
    std::vector<std::pair<Segment_3,Segment_3>> segments;
    while (segments.size()<count) {
        Segment_3 s1 = random_axis_of_sphere(rs);
        Segment_3 s2 = random_axis_of_sphere(rs);
        Object o = intersection(s1,s2);
        Point_3 p;
        if(CGAL::assign(p,o))
            segments.push_back(std::make_pair(s1,s2));
    }

    while (segments.size()<count*2) {
        Segment_3 s1 = random_chord_of_sphere(rs);
        Segment_3 s2 = random_chord_of_sphere(rs);
        Object o = intersection(s1,s2);
        Point_3 p; Segment_3 s3;
        if(!CGAL::assign(p,o) && !CGAL::assign(s3,o))
            segments.push_back(std::make_pair(s1,s2));
    }

    Transform t(CGAL::TRANSLATION, Vector_3(0,0,1));
    while (segments.size()<count*3) {
        Segment_3 s1 = random_chord_of_sphere(rs);
        Segment_3 s2 = s1.transform(t);
        Object o = intersection(s1,s2);
        Point_3 p; Segment_3 s3;
        if(!CGAL::assign(p,o) && !CGAL::assign(s3,o))
            segments.push_back(std::make_pair(s1,s2));
    }

    std::random_device rd;
    std::mt19937 g(rd());
    std::shuffle(segments.begin(), segments.end(), g);

    std::cout << "Ray      |  Segment" << std::endl;
    std::cout << std::setw(7);
    timing_test(segments,ray_does_intersect_internally_old,segment_does_intersect_internally_old);
    timing_test(segments,ray_does_intersect_internally,segment_does_intersect_internally);

}

Results:

Ray      |  Segment
7.16862  |  13.2177
14.0318  |  24.6894