CGAL users discussion list

[myociss <>]

*intersections.hpp:
*

#include <vector>
#include <list>
#include <array>
#include <iterator>
#include <functional>
#include <memory>
#include <Eigen/Dense>

#include <CGAL/Simple_cartesian.h>
#include <CGAL/AABB_tree.h>
#include <CGAL/AABB_traits.h>
#include <CGAL/AABB_triangle_primitive.h>


#ifndef INTERSECTIONS_HPP
#define INTERSECTIONS_HPP

using namespace std;
using namespace Eigen;

typedef CGAL::Simple_cartesian<double> K;

typedef K::FT FT;
//typedef K::Ray_3 Ray;
typedef K::Line_3 Line;
typedef K::Point_3 Point;
typedef K::Plane_3 Plane;
typedef K::Segment_3 Segment;
typedef K::Triangle_3 Triangle;
typedef list<Triangle>::iterator Iterator;
typedef CGAL::AABB_triangle_primitive<K, Iterator> Primitive;
typedef CGAL::AABB_traits<K, Primitive> AABB_triangle_traits;
typedef CGAL::AABB_tree<AABB_triangle_traits> Tree;


typedef Tree::Primitive_id Primitive_id;


class Intersection {
    bool exists;
    Vector3d point;
    Vector3d normal;
    vector<Vector3d> triangle;
  public:
    Intersection (bool _exists, Vector3d _point, Vector3d _normal,
vector<Vector3d> _triangle);
    bool Exists();
    Vector3d Point();
    Vector3d Normal();
    vector<Vector3d> Triangle();
};


class AABBTree {
    shared_ptr<Tree> tree;
  public:
    AABBTree (vector<array&lt;double, 3>> & _points, vector<array&lt;int,
3>> & _tris,
      Vector3d & p0, Vector3d & p1);
    Intersection getIntersection(Vector3d & p0, Vector3d & p1, 
      vector<Vector3d> & _lastIntersected);
};



*intersections.cpp
*

#include <vector>
#include <list>
#include <array>
#include <iterator>
#include <future>
#include <iostream>
#include <functional>
#include <thread>
#include <limits>
#include <Eigen/Dense>

#include <boost/optional.hpp>
#include <CGAL/Simple_cartesian.h>
#include <CGAL/AABB_tree.h>
#include <CGAL/AABB_traits.h>
#include <CGAL/AABB_triangle_primitive.h>

#include "intersections.hpp"

using namespace std;
using namespace Eigen;

typedef CGAL::Simple_cartesian<double> K;

typedef K::FT FT;
//typedef K::Ray_3 Ray;
typedef K::Line_3 Line;
typedef K::Point_3 Point;
typedef K::Direction_3 Direction;
typedef K::Segment_3 Segment;
typedef K::Triangle_3 Triangle;
typedef list<Triangle>::iterator Iterator;
typedef CGAL::AABB_triangle_primitive<K, Iterator> Primitive;
typedef CGAL::AABB_traits<K, Primitive> AABB_triangle_traits;
typedef CGAL::AABB_tree<AABB_triangle_traits> Tree;
typedef Tree::Primitive_id Primitive_id;
//typedef Tree::Intersection_and_primitive_id<Segment>
Intersection_and_primitive_id<Segment>
typedef
boost::optional<Tree::Intersection_and_primitive_id&lt;Segment>::Type>
Segment_intersection;




//AABBTree::AABBTree(vector<array&lt;array&lt;double, 3>, 3>> & _scene){
AABBTree::AABBTree(vector<array&lt;double, 3>> & _points,
vector<array&lt;int, 3>> & _tris, 
    Vector3d & p0, Vector3d & p1){



    vector<Point> points;
    list<Triangle> triangles;
    
    for (unsigned long int i=0; i < _points.size(); ++i){
        points.push_back(Point(_points[i][0], _points[i][1],
_points[i][2]));
    }

    for (unsigned long int i=0; i < _tris.size(); ++i){
        array<int, 3> _tri = _tris[i];
        triangles.push_back( Triangle(points[_tri[0]], points[_tri[1]],
points[_tri[2]]) );
    }


    tree = make_shared<Tree>(triangles.begin(), triangles.end());


    Triangle intersected;
    Point closestPoint;

    double closestDist = numeric_limits<double>::max();
    bool intersectionFound = false;


    Segment segment_query(Point(p0[0], p0[1], p0[2]), Point(p1[0], p1[1],
p1[2]));
    list<Segment_intersection> intersections;

*    //uncommenting this line causes the tree to return the correct
triangles when getIntersection is called later
*    *//cout << tree->do_intersect(segment_query) << endl;*
}


Intersection::Intersection(bool _exists, Vector3d _point, Vector3d _normal,
vector<Vector3d> _triangle){
    exists = _exists;
    point = _point;
    normal = _normal;
    triangle = _triangle;
}

bool Intersection::Exists(){
    return exists;
}

Vector3d Intersection::Point(){
    return point;
}
    
Vector3d Intersection::Normal(){
    return normal;
}

vector<Vector3d> Intersection::Triangle(){
    return triangle;
}


Intersection AABBTree::getIntersection(Vector3d & p0, Vector3d & p1, 
    vector<Vector3d> & _lastIntersected)
{

    Triangle lastIntersected;

    if(!_lastIntersected.empty()){
        Point l0(_lastIntersected[0][0], _lastIntersected[0][1],
_lastIntersected[0][2]);
            Point l1(_lastIntersected[1][0], _lastIntersected[1][1],
_lastIntersected[1][2]);
            Point l2(_lastIntersected[2][0], _lastIntersected[2][1],
_lastIntersected[2][2]);

        lastIntersected = Triangle(l0, l1, l2);
    }


    Triangle intersected;
    Point closestPoint;

    double closestDist = numeric_limits<double>::max();
    bool intersectionFound = false;


    Segment segment_query(Point(p0[0], p0[1], p0[2]), Point(p1[0], p1[1],
p1[2]));
    list<Segment_intersection> intersections;
    tree->all_intersections(segment_query, back_inserter(intersections));

    for (auto it = intersections.begin(); it != intersections.end(); ++it){

        Primitive pr = (*it)->second;
        Triangle tr = pr.datum();


        cout << "triangle points:" << endl;

        cout << tr.vertex(0) << endl;
        cout << tr.vertex(1) << endl;
        cout << tr.vertex(2) << endl;

        cout << "reference point" << endl;
        cout << pr.reference_point() << endl;


        if(boost::get<Point>(& (*it)->first )){
            intersectionFound = true;
            Point p = boost::get<Point>( (*it)->first );
            double squaredDist = pow(p.x() - p0[0], 2) + pow(p.y() - p0[1],
2) + pow(p.z() - p0[2], 2);

            if(squaredDist < closestDist){
                closestDist = squaredDist;
                intersected = tr;
                closestPoint = p;

            }
            }
    }

    cout << "closest point: " << endl;
    cout << closestPoint << endl;

    //Intersection foo;

    if(intersectionFound){
        //double pt[3] = {closestPoint.x(), closestPoint.y(),
closestPoint.z()};
        Vector3d pt(closestPoint.x(), closestPoint.y(), closestPoint.z());
        Plane triPlane = intersected.supporting_plane();
        Vector3d triNormal(triPlane.orthogonal_vector().x(),
triPlane.orthogonal_vector().y(), 
            triPlane.orthogonal_vector().z());

        if(triNormal.dot(p0 - pt) < 0){
            Plane reversed = triPlane.opposite();
            triNormal = Vector3d(reversed.orthogonal_vector().x(), 
                reversed.orthogonal_vector().y(),
reversed.orthogonal_vector().z());
        }
        triNormal.normalize();
        //double normal[3] = {triNormal[0], triNormal[1], triNormal[2]};

        Point v0 = intersected.vertex(0);
        Point v1 = intersected.vertex(1);
        Point v2 = intersected.vertex(2);

        /*cout << "..." << endl;
        cout << v0 << endl;
        cout << v1 << endl;
        cout << v2 << endl;*/

        vector<Vector3d> face;

        face.push_back(Vector3d(v0.x(), v0.y(), v0.z()));
        face.push_back(Vector3d(v1.x(), v1.y(), v1.z()));
        face.push_back(Vector3d(v2.x(), v2.y(), v2.z()));

        Intersection intersection(true, pt, triNormal, face);
        return intersection;
    } else {
        Vector3d pt(0.0,0.0,0.0);
        Vector3d normal(0.0,0.0,0.0);

        vector<Vector3d> face;
        
        face.push_back(Vector3d(0.0, 0.0, 0.0));
        face.push_back(Vector3d(0.0, 0.0, 0.0));
        face.push_back(Vector3d(0.0, 0.0, 0.0));
        
        Intersection intersection(false, pt, normal, face);
        return intersection;
    }

    //return foo;
}





--
Sent from: http://cgal-discuss.949826.n4.nabble.com/