To calculate intersections between a ray and a mesh

Question

Problem description

I want to calculate all intersections between a ray and a surface mesh which was built from points cloud using the function CGAL::advancing_front_surface_reconstruction(). The cost of time matters very much for me. So I want to use the package named Fast Intersection and Distance Computation. But when I built AABBtree from a surface mesh by the function tree(), a bug occurs. The IDE's bugs list shows: '<function-style-cast>' can not transform to CGAL::AABB_face_graph_triangle_primitive<meshCGAL, CGAL::Default, CGAL::Tag_true, CGAL::tag_false> from "initializer list".
I use the kernel named CGAL::Exact_predicates_inexact_constructions_kernel when I reconstruct mesh from points cloud with CGAL::advancing_front_surface_reconstruction(), but a Simple_cartesian<double> kernel is required in the Fast Intersection and Distance Computation. Is that the main reason? How can I use the two program package together?

Code

The code is as follows:

#include<iostream>
#include<vector>
#include<string>
#include<array>
#include<fstream>
#include<algorithm>

#include <CGAL/Simple_cartesian.h>
#include<CGAL/AABB_tree.h>
#include<CGAL/AABB_traits.h>
#include<CGAL/AABB_segment_primitive.h>
#include<CGAL/Surface_mesh.h>
#include <CGAL/AABB_face_graph_triangle_primitive.h>

#include<CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include<CGAL/Advancing_front_surface_reconstruction.h>
#include<CGAL/Surface_mesh.h>
#include <CGAL/boost/graph/IO/polygon_mesh_io.h>
#include<CGAL/Polygon_mesh_processing/corefinement.h>
#include<CGAL/polygon_mesh_processing/remesh.h>


using std::cout;
using std::cin;
using std::endl;
using std::string;

namespace PMP = CGAL::Polygon_mesh_processing;

typedef std::array<std::size_t, 3> Facet;
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::Point_3 Point_3;
typedef CGAL::Surface_mesh <Point_3> Mesh;

typedef CGAL::Simple_cartesian<double> K;
typedef K::FT FT;
typedef K::Point_3 pointCGAL;
typedef K::Vector_3 vectorCGAL;
typedef K::Ray_3 rayCGAL;
typedef CGAL::Surface_mesh<pointCGAL> meshCGAL;
typedef boost::graph_traits<meshCGAL>::face_descriptor face_descriptor;
typedef boost::graph_traits<meshCGAL>::halfedge_descriptor halfedge_descriptor;

typedef CGAL::AABB_face_graph_triangle_primitive<meshCGAL> Primitive;
typedef CGAL::AABB_traits<K, Primitive> Traits;
typedef CGAL::AABB_tree<Traits> Tree;
typedef boost::optional<Tree::Intersection_and_primitive_id<rayCGAL>::Type> Ray_intersection;


struct Construct {
    Mesh& mesh;
    template <typename PointIterator>
    Construct(Mesh& mesh, PointIterator b, PointIterator e) :mesh(mesh) {
        for (; b != e; ++b) {
            boost::graph_traits<Mesh>::vertex_descriptor v;
            v = add_vertex(mesh);
            mesh.point(v) = *b;
        }
    }

    Construct& operator=(const Facet f) {
        typedef boost::graph_traits<Mesh>::vertex_descriptor vertex_descriptor;
        typedef boost::graph_traits<Mesh>::vertices_size_type size_type;
        mesh.add_face(vertex_descriptor(static_cast<size_type>(f[0])),
            vertex_descriptor(static_cast<size_type>(f[1])),
            vertex_descriptor(static_cast<size_type>(f[2])));
        return *this;
    }

    Construct& operator*() { return *this; }
    Construct& operator++() { return *this; }
    Construct& operator++(int) { return *this; }
};


int main() {

    //load sv
    Mesh sv;
    std::cout << "load sv49..." << std::endl;
    std::ifstream fin("./svPoints/sv49.txt");
    std::vector<Point_3> points;
    std::vector<Facet> facets;
    std::copy(std::istream_iterator<Point_3>(fin),
        std::istream_iterator<Point_3>(),
        std::back_inserter(points));
    fin.close();

    //convert sv to Surface_mesh, repair
    Construct construct(sv, points.begin(), points.end());
    CGAL::advancing_front_surface_reconstruction(points.begin(), points.end(), construct);//convert sv to surface_mesh
    CGAL::Polygon_mesh_processing::remove_isolated_vertices(sv);//filter out isolated vertices
    if (!CGAL::Polygon_mesh_processing::is_outward_oriented(sv))
        CGAL::Polygon_mesh_processing::reverse_face_orientations(sv);//make sure that the meshes are outer oriented
    CGAL::IO::write_polygon_mesh("sv49.off",sv);

    //build AABBtree of sv
    Tree tree(faces(sv).first, faces(sv).second, sv);
    vectorCGAL orientation(0, 0, 1);
    pointCGAL start = pointCGAL(0, 0, 0);
    rayCGAL ray_query(start, orientation);
    std::vector<Ray_intersection> intersections;
    tree.all_intersections(ray_query, std::back_inserter(intersections));
    if (!intersections.empty()) {
        for (auto intersect : intersections) {
            auto a = intersect->first;
            if (pointCGAL* p = boost::get<pointCGAL>(&a)) {
                cout << *p << endl;
            }
        }
    }
    return 0;
}

Runtime environment

IDE: Visual Studio 2017
CGAL version: 5.3
Solution Configuration: Debug x64

File

File sv49.txt used in the program can download from GitHub, here's link.

Additional problem
When I use kernel CGAL::Exact_predicates_inexact_constructions_kernel instead of kernel Simple_cartesian, the program is compiled without errors. But when I test all intersections between a ray starting at point(-100.0,-50.0,-100.0) in the direction of vector(0.0, 0.0, 1.0) and the mesh named sv in the program, I got 46 intersections. The result is ridiculous. Then I try the kernel CGAL::Exact_predicates_exact_constructions_kernel, the result is also inexact. The result shows all intersections are as follows:

-100 -50 -10
-100 -50 -10
-100 -50 -10
-100 -50 -9
-100 -50 -8
-100 -50 -8
-100 -50 -9
-100 -50 -3
-100 -50 -1
-100 -50 -7
-100 -50 -7
-100 -50 -6
-100 -50 -5
-100 -50 -6
-100 -50 -5
-100 -50 -4
-100 -50 -4
-100 -50 -3
-100 -50 -2
-100 -50 -2
-100 -50 0
-100 -50 -0
-100 -50 -1
-100 -50 -0
-100 -50 -0

True intersections should be (-100, -50, -0) and (-100, -50, -10). The repeated intersections are reasonable, but the others are impossible intersections.

Answer 1

After I paint the point (-100,-50,-0) in the mesh sv, I got the reason why the result which seems to be ridiculous come up. Technically, it's not a wrong answer mathematically, it's just not what I was expecting. The ray starting at this point and in the direction of (0, 0, 1) is a ray intersect with many facets on the side of the grid. The CGAL library does a good job of calculating intersections in this case, but I need to filter out the intersections.