Commit 501ba7aa authored by Jan Möbius's avatar Jan Möbius
Browse files

Directed BSP ray shooting including unit tests

Fixed undirected shooting function returning wrong nearest triangle (was always the one with the most negative distance)

git-svn-id: http://www.openflipper.org/svnrepo/OpenFlipper/branches/Free@15172 383ad7c9-94d9-4d36-a494-682f7c89f535
parent 52f52dd9
......@@ -154,3 +154,4 @@ target_link_libraries ( ACG OpenMeshCore
${GLEW_LIBRARY}
${GLUT_LIBRARIES}
${ADDITIONAL_LINK_LIBRARIES} )
......@@ -72,7 +72,7 @@ nearest(const Point& _p) const
data.ref = _p;
data.dist = std::numeric_limits<Scalar>::infinity();
if (this->root_ == 0)
throw std::runtime_error("It seems like the BSP hasn't been built, yet. Did you call build(...)?");
throw std::runtime_error("It seems like the BSP hasn't been built, yet. Did you call build(...)?");
_nearest(this->root_, data);
return NearestNeighbor(data.nearest, sqrt(data.dist));
}
......@@ -90,34 +90,34 @@ _nearest(Node* _node, NearestNeighborData& _data) const
if (!_node->left_child_)
{
Scalar dist;
for (HandleIter it=_node->begin(); it!=_node->end(); ++it)
{
dist = this->traits_.sqrdist(*it, _data.ref);
if (dist < _data.dist)
{
_data.dist = dist;
_data.nearest = *it;
_data.dist = dist;
_data.nearest = *it;
}
}
}
// non-terminal node
else
{
Scalar dist = _node->plane_.distance(_data.ref);
if (dist > 0.0)
{
_nearest(_node->left_child_, _data);
if (dist*dist < _data.dist)
_nearest(_node->right_child_, _data);
_nearest(_node->right_child_, _data);
}
else
{
_nearest(_node->right_child_, _data);
if (dist*dist < _data.dist)
_nearest(_node->left_child_, _data);
_nearest(_node->left_child_, _data);
}
}
}
......@@ -129,16 +129,34 @@ typename BSPImplT<BSPCore>::RayCollision
BSPImplT<BSPCore>::
raycollision(const Point& _p, const Point& _r) const
{
// Prepare the struct for returning the data
RayCollisionData data;
data.ref = _p;
data.dist = FLT_MAX;
data.ray = _r;
data.hit_vertices.clear();
_raycollision(this->root_, data);
_raycollision_non_directional(this->root_, data);
return RayCollision(data.nearest, data.dist, data.hit_vertices);
}
template <class BSPCore>
typename BSPImplT<BSPCore>::RayCollision
BSPImplT<BSPCore>::
directionalRaycollision(const Point& _p, const Point& _r) const {
// Prepare the struct for returning the data
RayCollisionData data;
data.ref = _p;
data.dist = FLT_MAX;
data.ray = _r;
data.hit_vertices.clear();
_raycollision_directional(this->root_, data);
return RayCollision(data.nearest, data.dist, data.hit_vertices);
}
//-----------------------------------------------------------------------------
......@@ -146,44 +164,94 @@ raycollision(const Point& _p, const Point& _r) const
template <class BSPCore>
void
BSPImplT<BSPCore>::
_raycollision(Node* _node, RayCollisionData& _data) const
_raycollision_non_directional(Node* _node, RayCollisionData& _data) const
{
// terminal node
if (!_node->left_child_)
// terminal node
if (!_node->left_child_)
{
Scalar dist;
Point v0, v1, v2;
Scalar u, v;
for (HandleIter it=_node->begin(); it!=_node->end(); ++it)
{
Scalar dist;
Point v0, v1, v2;
Scalar u, v;
this->traits_.points(*it, v0, v1, v2);
if (ACG::Geometry::triangleIntersection(_data.ref, _data.ray, v0, v1, v2, dist, u, v)) {
for (HandleIter it=_node->begin(); it!=_node->end(); ++it)
_data.hit_vertices.push_back(*it);
// face intersects with ray. But is it closer than any that we have found so far?
if ( fabs(dist) < _data.dist)
{
this->traits_.points(*it, v0, v1, v2);
if (ACG::Geometry::triangleIntersection(_data.ref, _data.ray, v0, v1, v2, dist, u, v)) {
_data.hit_vertices.push_back(*it);
//face intersects with ray. But is it closer than any that we have found so far?
if (dist < _data.dist)
{
_data.dist = dist;
_data.nearest = *it;
}
}
_data.dist = fabs(dist);
_data.nearest = *it;
}
}
}
}
// non-terminal node
else
{
Scalar tmin, tmax;
bool used = false;
if ( _node->left_child_ && ACG::Geometry::axisAlignedBBIntersection( _data.ref, _data.ray, _node->left_child_->bb_min, _node->left_child_->bb_max, tmin, tmax)) {
_raycollision_non_directional(_node->left_child_, _data);
}
if ( _node->right_child_ && ACG::Geometry::axisAlignedBBIntersection( _data.ref, _data.ray, _node->right_child_->bb_min, _node->right_child_->bb_max, tmin, tmax)) {
_raycollision_non_directional(_node->right_child_, _data);
}
}
}
//-----------------------------------------------------------------------------
template <class BSPCore>
void
BSPImplT<BSPCore>::
_raycollision_directional(Node* _node, RayCollisionData& _data) const
{
// terminal node
if (!_node->left_child_)
{
Scalar dist;
Point v0, v1, v2;
Scalar u, v;
for (HandleIter it=_node->begin(); it!=_node->end(); ++it)
{
Scalar tmin, tmax;
bool used = false;
if ( _node->left_child_ && ACG::Geometry::axisAlignedBBIntersection( _data.ref, _data.ray, _node->left_child_->bb_min, _node->left_child_->bb_max, tmin, tmax)) {
_raycollision(_node->left_child_, _data);
}
if ( _node->right_child_ && ACG::Geometry::axisAlignedBBIntersection( _data.ref, _data.ray, _node->right_child_->bb_min, _node->right_child_->bb_max, tmin, tmax)) {
_raycollision(_node->right_child_, _data);
this->traits_.points(*it, v0, v1, v2);
if (ACG::Geometry::triangleIntersection(_data.ref, _data.ray, v0, v1, v2, dist, u, v)) {
if ( dist < 0.0 )
continue;
_data.hit_vertices.push_back(*it);
// face intersects with ray. But is it closer than any that we have found so far?
// Note as we rely on the direction of the hit, so we will never get negative directions here
if ( dist < _data.dist)
{
_data.dist = fabs(dist);
_data.nearest = *it;
}
}
}
}
// non-terminal node
else
{
Scalar tmin, tmax;
bool used = false;
if ( _node->left_child_ && ACG::Geometry::axisAlignedBBIntersection( _data.ref, _data.ray, _node->left_child_->bb_min, _node->left_child_->bb_max, tmin, tmax)) {
_raycollision_directional(_node->left_child_, _data);
}
if ( _node->right_child_ && ACG::Geometry::axisAlignedBBIntersection( _data.ref, _data.ray, _node->right_child_->bb_min, _node->right_child_->bb_max, tmin, tmax)) {
_raycollision_directional(_node->right_child_, _data);
}
}
}
......
......@@ -97,8 +97,8 @@ public: //---------------------------------------------------------------------
{
RayCollision() {}
RayCollision(Handle _h, Scalar _d, std::vector<Handle> _v) : handle(_h), dist(_d), hit_vertices(_v) {}
Handle handle;
Scalar dist;
Handle handle;
Scalar dist;
std::vector<Handle> hit_vertices;
};
......@@ -107,14 +107,27 @@ public: //---------------------------------------------------------------------
/** \brief intersect mesh with ray
*
* This function shots a ray through the ray and collects all intersected triangles and
* the handle of the closest face
* This function shots a ray through the mesh and collects all intersected triangles and
* the handle of the closest face ( non-directional, so no matter of the ray direction, the
* closest face handle is returned in either direction)
*
* @param _p Start point of the ray
* @param _r Ray direction
* @return Collision information
*/
RayCollision raycollision (const Point& _p, const Point& _r) const;
/** \brief intersect mesh with ray
*
* This function shots a ray through the mesh and collects all intersected triangles and
* the handle of the closest face ( directional, so the ray direction is taken into account!)
*
* @param _p Start point of the ray
* @param _r Ray direction
* @return Collision information
*/
RayCollision directionalRaycollision (const Point& _p, const Point& _r) const;
private: //---------------------------------------------------------------------
......@@ -146,8 +159,19 @@ private: //---------------------------------------------------------------------
// Recursive part of nearest()
void _nearest(Node* _node, NearestNeighborData& _data) const;
//resursive part of raycollide()
void _raycollision(Node* _node, RayCollisionData& _data) const;
/** \brief recursive part of raycollision()
*
* @param _node The current node in the tree
* @param _data Data pointer, used to collect the collision information
*/
void _raycollision_non_directional(Node* _node, RayCollisionData& _data) const;
/** \brief recursive part of directionalRaycollision()
*
* @param _node The current node in the tree
* @param _data Data pointer, used to collect the collision information
*/
void _raycollision_directional(Node* _node, RayCollisionData& _data) const;
};
//=============================================================================
......
......@@ -105,77 +105,78 @@ struct TreeNode
Plane plane_;
Point bb_min, bb_max;
// /// This visualizes the bounding boxes
// void visualizeTree(PolyMesh *_object, int _max_depth)
// {
// if (_max_depth > 0 && (left_child_ || right_child_) )
// {
// if (left_child_)
// left_child_->visualizeTree(_object, _max_depth-1);
// if (right_child_)
// right_child_->visualizeTree(_object, _max_depth-1);
// }
// else
// {
// Point size_ = bb_max - bb_min;
//
// std::vector<VertexHandle> vhandle(8);
// vhandle[0] = _object->add_vertex(bb_min+Point(0.0,0.0,size_[2]));
// vhandle[1] = _object->add_vertex(bb_min+Point(size_[0],0.0,size_[2]));
// vhandle[2] = _object->add_vertex(bb_min+Point(size_[0],size_[1],size_[2]));
// vhandle[3] = _object->add_vertex(bb_min+Point(0.0,size_[1],size_[2]));
// vhandle[4] = _object->add_vertex(bb_min+Point(0.0,0.0,0.0));
// vhandle[5] = _object->add_vertex(bb_min+Point(size_[0],0.0,0.0));
// vhandle[6] = _object->add_vertex(bb_min+Point(size_[0],size_[1],0.0));
// vhandle[7] = _object->add_vertex(bb_min+Point(0.0,size_[1],0.0));
//
//
// // generate (quadrilateral) faces
// std::vector<VertexHandle> face_vhandles;
//
// face_vhandles.clear();
// face_vhandles.push_back(vhandle[0]);
// face_vhandles.push_back(vhandle[1]);
// face_vhandles.push_back(vhandle[2]);
// face_vhandles.push_back(vhandle[3]);
// _object->add_face(face_vhandles);
//
// face_vhandles.clear();
// face_vhandles.push_back(vhandle[7]);
// face_vhandles.push_back(vhandle[6]);
// face_vhandles.push_back(vhandle[5]);
// face_vhandles.push_back(vhandle[4]);
// _object->add_face(face_vhandles);
//
// face_vhandles.clear();
// face_vhandles.push_back(vhandle[1]);
// face_vhandles.push_back(vhandle[0]);
// face_vhandles.push_back(vhandle[4]);
// face_vhandles.push_back(vhandle[5]);
// _object->add_face(face_vhandles);
//
// face_vhandles.clear();
// face_vhandles.push_back(vhandle[2]);
// face_vhandles.push_back(vhandle[1]);
// face_vhandles.push_back(vhandle[5]);
// face_vhandles.push_back(vhandle[6]);
// _object->add_face(face_vhandles);
//
// face_vhandles.clear();
// face_vhandles.push_back(vhandle[3]);
// face_vhandles.push_back(vhandle[2]);
// face_vhandles.push_back(vhandle[6]);
// face_vhandles.push_back(vhandle[7]);
// _object->add_face(face_vhandles);
//
// face_vhandles.clear();
// face_vhandles.push_back(vhandle[0]);
// face_vhandles.push_back(vhandle[3]);
// face_vhandles.push_back(vhandle[7]);
// face_vhandles.push_back(vhandle[4]);
// _object->add_face(face_vhandles);
// }
// }
/// This visualizes the bounding boxes
template< typename MeshT >
void visualizeTree(MeshT *_object, int _max_depth)
{
if (_max_depth > 0 && (left_child_ || right_child_) )
{
if (left_child_)
left_child_->visualizeTree(_object, _max_depth-1);
if (right_child_)
right_child_->visualizeTree(_object, _max_depth-1);
}
else
{
Point size_ = bb_max - bb_min;
std::vector<VertexHandle> vhandle(8);
vhandle[0] = _object->add_vertex(bb_min+Point(0.0,0.0,size_[2]));
vhandle[1] = _object->add_vertex(bb_min+Point(size_[0],0.0,size_[2]));
vhandle[2] = _object->add_vertex(bb_min+Point(size_[0],size_[1],size_[2]));
vhandle[3] = _object->add_vertex(bb_min+Point(0.0,size_[1],size_[2]));
vhandle[4] = _object->add_vertex(bb_min+Point(0.0,0.0,0.0));
vhandle[5] = _object->add_vertex(bb_min+Point(size_[0],0.0,0.0));
vhandle[6] = _object->add_vertex(bb_min+Point(size_[0],size_[1],0.0));
vhandle[7] = _object->add_vertex(bb_min+Point(0.0,size_[1],0.0));
// generate (quadrilateral) faces
std::vector<VertexHandle> face_vhandles;
face_vhandles.clear();
face_vhandles.push_back(vhandle[0]);
face_vhandles.push_back(vhandle[1]);
face_vhandles.push_back(vhandle[2]);
face_vhandles.push_back(vhandle[3]);
_object->add_face(face_vhandles);
face_vhandles.clear();
face_vhandles.push_back(vhandle[7]);
face_vhandles.push_back(vhandle[6]);
face_vhandles.push_back(vhandle[5]);
face_vhandles.push_back(vhandle[4]);
_object->add_face(face_vhandles);
face_vhandles.clear();
face_vhandles.push_back(vhandle[1]);
face_vhandles.push_back(vhandle[0]);
face_vhandles.push_back(vhandle[4]);
face_vhandles.push_back(vhandle[5]);
_object->add_face(face_vhandles);
face_vhandles.clear();
face_vhandles.push_back(vhandle[2]);
face_vhandles.push_back(vhandle[1]);
face_vhandles.push_back(vhandle[5]);
face_vhandles.push_back(vhandle[6]);
_object->add_face(face_vhandles);
face_vhandles.clear();
face_vhandles.push_back(vhandle[3]);
face_vhandles.push_back(vhandle[2]);
face_vhandles.push_back(vhandle[6]);
face_vhandles.push_back(vhandle[7]);
_object->add_face(face_vhandles);
face_vhandles.clear();
face_vhandles.push_back(vhandle[0]);
face_vhandles.push_back(vhandle[3]);
face_vhandles.push_back(vhandle[7]);
face_vhandles.push_back(vhandle[4]);
_object->add_face(face_vhandles);
}
}
private:
/*
......
......@@ -163,3 +163,5 @@ _build(Node* _node,
}
//=============================================================================
......@@ -116,12 +116,18 @@ public: //---------------------------------------------------------------------
*/
void build(unsigned int _max_handles, unsigned int _max_depth);
// /** \brief Create a PolyMesh object that visualizes the bounding boxes of the BSP tree
// *
// * @param _object The output mesh which the tree will be written into
// * @param _max_depth The maximal depth that will be visualized
// */
//void visualizeTree(PolyMesh *_object, int _max_depth);
/** \brief Create a PolyMesh object that visualizes the bounding boxes of the BSP tree
*
* @param _object The output mesh which the tree will be written into
* @param _max_depth The maximal depth that will be visualized
*/
template <typename MeshT>
void
visualizeTree(MeshT *_object, int _max_depth)
{
root_->visualizeTree(_object, _max_depth-1);
_object->update_normals();
}
private:
/*
......
......@@ -58,9 +58,6 @@
#include "BSPTreeNode.hh"
#include "TriangleBSPCoreT.hh"
#include "BSPImplT.hh"
//#include "Distance.hh"
#include <ACG/Geometry/Algorithms.hh>
//== CLASS DEFINITION =========================================================
#include <list>
......@@ -80,114 +77,116 @@ class OpenMeshTriangleBSPTraits
{
public:
typedef typename Mesh::Scalar Scalar;
typedef typename Mesh::Point Point;
typedef typename Mesh::FaceHandle Handle;
typedef std::vector<Handle> Handles;
typedef typename Handles::iterator HandleIter;
typedef TreeNode<Mesh> Node;
OpenMeshTriangleBSPTraits(const Mesh& _mesh) : mesh_(_mesh) {}
/// Returns the points belonging to the face handle _h
inline void points(const Handle _h, Point& _p0, Point& _p1, Point& _p2) const
{
typename Mesh::CFVIter fv_it(mesh_.cfv_iter(_h));
_p0 = mesh_.point(fv_it);
++fv_it;
_p1 = mesh_.point(fv_it);
++fv_it;
_p2 = mesh_.point(fv_it);
}
Scalar sqrdist(const Handle _h, const Point& _p) const
{
Point p0, p1, p2, q;
points(_h, p0, p1, p2);
return ACG::Geometry::distPointTriangleSquared(_p, p0, p1, p2, q);
typedef typename Mesh::Scalar Scalar;
typedef typename Mesh::Point Point;
typedef typename Mesh::FaceHandle Handle;
typedef std::vector<Handle> Handles;
typedef typename Handles::iterator HandleIter;
typedef TreeNode<Mesh> Node;
OpenMeshTriangleBSPTraits(const Mesh& _mesh) : mesh_(_mesh) {}
/// Returns the points belonging to the face handle _h
inline void points(const Handle _h, Point& _p0, Point& _p1, Point& _p2) const
{
typename Mesh::CFVIter fv_it(mesh_.cfv_iter(_h));
_p0 = mesh_.point(fv_it);
++fv_it;
_p1 = mesh_.point(fv_it);
++fv_it;
_p2 = mesh_.point(fv_it);
}
Scalar sqrdist(const Handle _h, const Point& _p) const
{
Point p0, p1, p2, q;
points(_h, p0, p1, p2);
return ACG::Geometry::distPointTriangleSquared(_p, p0, p1, p2, q);
}
void calculateBoundingBox(Node* _node, Point& median, int& axis)
{
//determine splitting axis
HandleIter it_h;
Point p0, p1, p2, bb_min, bb_max;
bb_min.vectorize(std::numeric_limits<typename Point::value_type>::infinity());
bb_max.vectorize(-std::numeric_limits<typename Point::value_type>::infinity());
std::list<Point> vertices;
for (it_h = _node->begin(); it_h != _node->end(); ++it_h) {
points(*it_h, p0, p1, p2);
/*
bb_min.minimize(p0);
bb_min.minimize(p1);
bb_min.minimize(p2);
bb_max.maximize(p0);
bb_max.maximize(p1);
bb_max.maximize(p2);*/
vertices.push_back(p0);
vertices.push_back(p1);
vertices.push_back(p2);
}
void calculateBoundingBox(
Node* _node,
Point& median,
int& axis)
{
//determine splitting axis
HandleIter it_h;
Point p0, p1, p2, bb_min, bb_max;
bb_min.vectorize(std::numeric_limits<typename Point::value_type>::infinity());
bb_max.vectorize(-std::numeric_limits<typename Point::value_type>::infinity());
std::list<Point> vertices;
for (it_h=_node->begin(); it_h!=_node->end(); ++it_h)
{
points(*it_h, p0, p1, p2);
/*
bb_min.minimize(p0);
bb_min.minimize(p1);
bb_min.minimize(p2);
bb_max.maximize(p0);
bb_max.maximize(p1);
bb_max.maximize(p2);*/
vertices.push_back (p0);
vertices.push_back (p1);
vertices.push_back (p2);
}
bb_min = _node->bb_min;
bb_max = _node->bb_max;
// split longest side of bounding box
Point bb = bb_max - bb_min;
Scalar length = bb[0];
axis = 0;
if (bb[1] > length) length = bb[ (axis=1) ];
if (bb[2] > length) length = bb[ (axis=2) ];
//calculate the median value in axis-direction
switch (axis) {
case 0: vertices.sort (x_sort() ); break;
case 1: vertices.sort (y_sort() ); break;
case 2: vertices.sort (z_sort() ); break; }
vertices.unique(); ///todo: does this work with Points?!
int size = vertices.size();
typename std::list<Point>::iterator it_v;
it_v = vertices.begin();
std::advance(it_v, size/2);
median = *it_v;
bb_min = _node->bb_min;
bb_max = _node->bb_max;
// split longest side of bounding box