| [eigen] Bounding Volume Hierarchies |
[ Thread Index |
Date Index
| More lists.tuxfamily.org/eigen Archives
]
- To: eigen@xxxxxxxxxxxxxxxxxxx
- Subject: [eigen] Bounding Volume Hierarchies
- From: Ilya Baran <baran37@xxxxxxxxx>
- Date: Wed, 4 Mar 2009 19:01:40 -0500
- Dkim-signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=gamma; h=domainkey-signature:mime-version:received:date:message-id:subject :from:to:content-type; bh=xN6z7qwUSDizFa/JFftQLsEeBNxn4QfyDjTIkIRQY1A=; b=qvvlAgcjDnfE2XP5AVyPUyGQW9Pj/XNyhWaXTlt0dCJ7K6P6mdaTQgk+9PhntP2gRx bRHhNwgJ/XAMd0XQOvA7X8syOdmkCf531B5MLtBs1x6hTp753T1H9hzIV5Wt5St7MV1+ GH2HTsPZh/CoqDV0Z+bjnQBtfz/WaNILpbes8=
- Domainkey-signature: a=rsa-sha1; c=nofws; d=gmail.com; s=gamma; h=mime-version:date:message-id:subject:from:to:content-type; b=PusJCF/EHOS+e6VfkPVcihm/jq8FrsfugUsZvHTRWjdPRKgb0LmSs+U/PJiC1S1dOo WX5dZgodYXKbtinv4bgkRJDN6VVRndLGAJ4PM0fr9Tim2IEdXjPj6gCDVjbBneVoubpG 5EkAb+qOfq0B8v1iZw/tcE/p4L1/KpFqc4dtE=
Hi,
Since Eigen has a Geometry module, I think it would be useful to have
a bounding volume hierarchy (BVH) module, since they are so useful in
accelerating tons of different kinds of geometric computations. I'm
attaching a patch (to unsupported) that provides what I think is a
pretty generic approach (criticism is, of course, welcomed)--it
decouples traversal algorithms from both the query and the hierarchy.
It also provides a simple BVH implementation based on AlignedBox. The
documentation in the patch describes the features in detail, but as
this is my first experience with CMake and doxygen, I may have gotten
some paths/settings wrong. In particular, I can't seem to get the
module documentation to display properly (at one point it did, but now
on even on a fresh checkout without my patch it doesn't). I'm also
not sure how well I followed Eigen's coding conventions.
To get this working, the patch also fixes a memory alignment bug in
AlignedBox (no pun intended :) and adds union and intersection
operators to it as well as a couple of other things. However, the
current AlignedBox design is slightly flawed: the intersection of two
disjoint boxes may be empty, but unioning it with another box can
change the box--checking for emptiness has to be done somewhere. I
think the best way is to check at construction/modification time
whether the box is empty and call setNull, if so. That has drawbacks,
though, such as making the non-const min() and max() either impossible
or requiring wrapper trickery.
Thanks,
-Ilya
Index: test/geo_alignedbox.cpp
===================================================================
--- test/geo_alignedbox.cpp (revision 935303)
+++ test/geo_alignedbox.cpp (working copy)
@@ -63,6 +63,13 @@
VERIFY_IS_APPROX(hp1f.template cast<Scalar>(),b0);
AlignedBox<Scalar,Dim> hp1d = b0.template cast<Scalar>();
VERIFY_IS_APPROX(hp1d.template cast<Scalar>(),b0);
+
+ // alignment -- make sure there are no alignment assertions
+ BoxType *bp0 = new BoxType(dim);
+ BoxType *bp1 = new BoxType(dim);
+ bp0->extend(*bp1);
+ delete bp0;
+ delete bp1;
}
void test_geo_alignedbox()
Index: Eigen/src/Geometry/AlignedBox.h
===================================================================
--- Eigen/src/Geometry/AlignedBox.h (revision 935303)
+++ Eigen/src/Geometry/AlignedBox.h (working copy)
@@ -41,7 +41,7 @@
class AlignedBox
{
public:
-EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_AmbientDim==Dynamic ? Dynamic : _AmbientDim+1)
+EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_AmbientDim)
enum { AmbientDimAtCompileTime = _AmbientDim };
typedef _Scalar Scalar;
typedef typename NumTraits<Scalar>::Real RealScalar;
@@ -85,6 +85,9 @@
/** \returns a non const reference to the maximal corner */
inline VectorType& max() { return m_max; }
+ /** \returns the center of the box */
+ inline VectorType center() const { return (m_min + m_max) / 2; }
+
/** \returns true if the point \a p is inside the box \c *this. */
inline bool contains(const VectorType& p) const
{ return (m_min.cwise()<=p).all() && (p.cwise()<=m_max).all(); }
@@ -105,6 +108,14 @@
inline AlignedBox& clamp(const AlignedBox& b)
{ m_min = m_min.cwise().max(b.m_min); m_max = m_max.cwise().min(b.m_max); return *this; }
+ /** Returns an AlignedBox that is the intersection of \a b and \c *this */
+ inline AlignedBox operator&(const AlignedBox &b) const
+ { return AlignedBox(m_min.cwise().max(b.m_min), m_max.cwise().min(b.m_max)); }
+
+ /** Returns an AlignedBox that is the union of \a b and \c *this */
+ inline AlignedBox operator|(const AlignedBox &b) const
+ { return AlignedBox(m_min.cwise().min(b.m_min), m_max.cwise().max(b.m_max)); }
+
/** Translate \c *this by the vector \a t and returns a reference to \c *this. */
inline AlignedBox& translate(const VectorType& t)
{ m_min += t; m_max += t; return *this; }
@@ -115,6 +126,12 @@
*/
inline Scalar squaredExteriorDistance(const VectorType& p) const;
+ /** \returns the squared distance between the boxes \a b and \c *this,
+ * and zero if the boxes intersect.
+ * \sa exteriorDistance()
+ */
+ inline Scalar squaredExteriorDistance(const AlignedBox& b) const;
+
/** \returns the distance between the point \a p and the box \c *this,
* and zero if \a p is inside the box.
* \sa squaredExteriorDistance()
@@ -122,6 +139,13 @@
inline Scalar exteriorDistance(const VectorType& p) const
{ return ei_sqrt(squaredExteriorDistance(p)); }
+ /** \returns the distance between the boxes \a b and \c *this,
+ * and zero if the boxes intersect.
+ * \sa squaredExteriorDistance()
+ */
+ inline Scalar exteriorDistance(const AlignedBox& b) const
+ { return ei_sqrt(squaredExteriorDistance(b)); }
+
/** \returns \c *this with scalar type casted to \a NewScalarType
*
* Note that if \a NewScalarType is equal to the current scalar type of \c *this
@@ -170,4 +194,19 @@
return dist2;
}
+template<typename Scalar,int AmbiantDim>
+inline Scalar AlignedBox<Scalar,AmbiantDim>::squaredExteriorDistance(const AlignedBox& b) const
+{
+ Scalar dist2 = 0.;
+ Scalar aux;
+ for (int k=0; k<dim(); ++k)
+ {
+ if ((aux = (b.m_min[k]-m_max[k]))>0.)
+ dist2 += aux*aux;
+ else if ( (aux = (m_min[k]-b.m_max[k]))>0. )
+ dist2 += aux*aux;
+ }
+ return dist2;
+}
+
#endif // EIGEN_ALIGNEDBOX_H
Index: unsupported/test/BVH.cpp
===================================================================
--- unsupported/test/BVH.cpp (revision 0)
+++ unsupported/test/BVH.cpp (revision 0)
@@ -0,0 +1,221 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra. Eigen itself is part of the KDE project.
+//
+// Copyright (C) 2009 Ilya Baran <ibaran@xxxxxxx>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#include <Eigen/StdVector>
+#include "main.h"
+#include <unsupported/Eigen/BVH>
+
+inline double SQR(double x) { return x * x; }
+
+template<int Dim>
+struct Ball
+{
+EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(double, Dim)
+
+ typedef Matrix<double, Dim, 1> VectorType;
+
+ Ball() {}
+ Ball(const VectorType &c, double r) : center(c), radius(r) {}
+
+ VectorType center;
+ double radius;
+};
+
+template<int Dim> AlignedBox<double, Dim> ei_bounding_box(const Matrix<double, Dim, 1> &v) { return AlignedBox<double, Dim>(v); }
+template<int Dim> AlignedBox<double, Dim> ei_bounding_box(const Ball<Dim> &b)
+{ return AlignedBox<double, Dim>(b.center.cwise() - b.radius, b.center.cwise() + b.radius); }
+
+template<int Dim>
+struct BallPointStuff //this class provides functions to be both an intersector and a minimizer, both for a ball and a point and for two trees
+{
+ typedef double Scalar;
+ typedef Matrix<double, Dim, 1> VectorType;
+ typedef Ball<Dim> BallType;
+ typedef AlignedBox<double, Dim> BoxType;
+
+ BallPointStuff() : calls(0), count(0) {}
+ BallPointStuff(const VectorType &inP) : p(inP), calls(0), count(0) {}
+
+
+ bool intersectVolume(const BoxType &r) { ++calls; return r.contains(p); }
+ bool intersectObject(const BallType &b) {
+ ++calls;
+ if((b.center - p).squaredNorm() < SQR(b.radius))
+ ++count;
+ return false; //continue
+ }
+
+ bool intersectVolumeVolume(const BoxType &r1, const BoxType &r2) { ++calls; return !(r1 & r2).isNull(); }
+ bool intersectVolumeObject(const BoxType &r, const BallType &b) { ++calls; return r.squaredExteriorDistance(b.center) < SQR(b.radius); }
+ bool intersectObjectVolume(const BallType &b, const BoxType &r) { ++calls; return r.squaredExteriorDistance(b.center) < SQR(b.radius); }
+ bool intersectObjectObject(const BallType &b1, const BallType &b2){
+ ++calls;
+ if((b1.center - b2.center).norm() < b1.radius + b2.radius)
+ ++count;
+ return false;
+ }
+ bool intersectVolumeObject(const BoxType &r, const VectorType &v) { ++calls; return r.contains(v); }
+ bool intersectObjectObject(const BallType &b, const VectorType &v){
+ ++calls;
+ if((b.center - v).squaredNorm() < SQR(b.radius))
+ ++count;
+ return false;
+ }
+
+ double minimumOnVolume(const BoxType &r) { ++calls; return r.squaredExteriorDistance(p); }
+ double minimumOnObject(const BallType &b) { ++calls; return std::max(0., (b.center - p).squaredNorm() - SQR(b.radius)); }
+ double minimumOnVolumeVolume(const BoxType &r1, const BoxType &r2) { ++calls; return r1.squaredExteriorDistance(r2); }
+ double minimumOnVolumeObject(const BoxType &r, const BallType &b) { ++calls; return SQR(std::max(0., r.exteriorDistance(b.center) - b.radius)); }
+ double minimumOnObjectVolume(const BallType &b, const BoxType &r) { ++calls; return SQR(std::max(0., r.exteriorDistance(b.center) - b.radius)); }
+ double minimumOnObjectObject(const BallType &b1, const BallType &b2){ ++calls; return SQR(std::max(0., (b1.center - b2.center).norm() - b1.radius - b2.radius)); }
+ double minimumOnVolumeObject(const BoxType &r, const VectorType &v) { ++calls; return r.squaredExteriorDistance(v); }
+ double minimumOnObjectObject(const BallType &b, const VectorType &v){ ++calls; return SQR(std::max(0., (b.center - v).norm() - b.radius)); }
+
+ VectorType p;
+ int calls;
+ int count;
+};
+
+template<int Dim>
+struct TreeTest
+{
+ typedef Matrix<double, Dim, 1> VectorType;
+ typedef Ball<Dim> BallType;
+ typedef AlignedBox<double, Dim> BoxType;
+
+ void testIntersect1()
+ {
+ std::vector<BallType> b;
+ for(int i = 0; i < 500; ++i) {
+ b.push_back(BallType(VectorType::Random(), 0.5 * ei_random(0., 1.)));
+ }
+ KdBVH<double, Dim, BallType> tree(b.begin(), b.end());
+
+ VectorType pt = VectorType::Random();
+ BallPointStuff<Dim> i1(pt), i2(pt);
+
+ for(int i = 0; i < (int)b.size(); ++i)
+ i1.intersectObject(b[i]);
+
+ BVIntersect(tree, i2);
+
+ VERIFY(i1.count == i2.count);
+ }
+
+ void testMinimize1()
+ {
+ std::vector<BallType> b;
+ for(int i = 0; i < 500; ++i) {
+ b.push_back(BallType(VectorType::Random(), 0.01 * ei_random(0., 1.)));
+ }
+ KdBVH<double, Dim, BallType> tree(b.begin(), b.end());
+
+ VectorType pt = VectorType::Random();
+ BallPointStuff<Dim> i1(pt), i2(pt);
+
+ double m1 = std::numeric_limits<double>::max(), m2 = m1;
+
+ for(int i = 0; i < (int)b.size(); ++i)
+ m1 = std::min(m1, i1.minimumOnObject(b[i]));
+
+ m2 = BVMinimize(tree, i2);
+
+ VERIFY_IS_APPROX(m1, m2);
+ }
+
+ void testIntersect2()
+ {
+ std::vector<BallType> b;
+ std::vector<VectorType> v;
+
+ for(int i = 0; i < 50; ++i) {
+ b.push_back(BallType(VectorType::Random(), 0.5 * ei_random(0., 1.)));
+ for(int j = 0; j < 3; ++j)
+ v.push_back(VectorType::Random());
+ }
+
+ KdBVH<double, Dim, BallType> tree(b.begin(), b.end());
+ KdBVH<double, Dim, VectorType> vTree(v.begin(), v.end());
+
+ BallPointStuff<Dim> i1, i2;
+
+ for(int i = 0; i < (int)b.size(); ++i)
+ for(int j = 0; j < (int)v.size(); ++j)
+ i1.intersectObjectObject(b[i], v[j]);
+
+ BVIntersect(tree, vTree, i2);
+
+ VERIFY(i1.count == i2.count);
+ }
+
+ void testMinimize2()
+ {
+ std::vector<BallType> b;
+ std::vector<VectorType> v;
+
+ for(int i = 0; i < 50; ++i) {
+ b.push_back(BallType(VectorType::Random(), 1e-7 + 1e-6 * ei_random(0., 1.)));
+ for(int j = 0; j < 3; ++j)
+ v.push_back(VectorType::Random());
+ }
+
+ KdBVH<double, Dim, BallType> tree(b.begin(), b.end());
+ KdBVH<double, Dim, VectorType> vTree(v.begin(), v.end());
+
+ BallPointStuff<Dim> i1, i2;
+
+ double m1 = std::numeric_limits<double>::max(), m2 = m1;
+
+ for(int i = 0; i < (int)b.size(); ++i)
+ for(int j = 0; j < (int)v.size(); ++j)
+ m1 = std::min(m1, i1.minimumOnObjectObject(b[i], v[j]));
+
+ m2 = BVMinimize(tree, vTree, i2);
+
+ VERIFY_IS_APPROX(m1, m2);
+ }
+};
+
+void test_BVH()
+{
+ for(int i = 0; i < g_repeat; i++) {
+ TreeTest<2> test2;
+ CALL_SUBTEST(test2.testIntersect1());
+ CALL_SUBTEST(test2.testMinimize1());
+ CALL_SUBTEST(test2.testIntersect2());
+ CALL_SUBTEST(test2.testMinimize2());
+
+ TreeTest<3> test3;
+ CALL_SUBTEST(test3.testIntersect1());
+ CALL_SUBTEST(test3.testMinimize1());
+ CALL_SUBTEST(test3.testIntersect2());
+ CALL_SUBTEST(test3.testMinimize2());
+
+ TreeTest<4> test4;
+ CALL_SUBTEST(test4.testIntersect1());
+ CALL_SUBTEST(test4.testMinimize1());
+ CALL_SUBTEST(test4.testIntersect2());
+ CALL_SUBTEST(test4.testMinimize2());
+ }
+}
Index: unsupported/test/CMakeLists.txt
===================================================================
--- unsupported/test/CMakeLists.txt (revision 935303)
+++ unsupported/test/CMakeLists.txt (working copy)
@@ -13,4 +13,6 @@
ei_add_test(forward_adolc " " ${ADOLC_LIBRARIES})
else(ADOLC_FOUND)
ei_add_property(EIGEN_MISSING_BACKENDS "Adolc")
-endif(ADOLC_FOUND)
\ No newline at end of file
+endif(ADOLC_FOUND)
+
+ei_add_test(BVH)
Index: unsupported/doc/Doxyfile.in
===================================================================
--- unsupported/doc/Doxyfile.in (revision 935303)
+++ unsupported/doc/Doxyfile.in (working copy)
@@ -626,7 +626,11 @@
EXAMPLE_PATH = "${Eigen_SOURCE_DIR}/doc/snippets" \
"${Eigen_BINARY_DIR}/doc/snippets" \
"${Eigen_SOURCE_DIR}/doc/examples" \
- "${Eigen_BINARY_DIR}/doc/examples"
+ "${Eigen_BINARY_DIR}/doc/examples" \
+ "${Eigen_SOURCE_DIR}/unsupported/doc/snippets" \
+ "${Eigen_BINARY_DIR}/unsupported/doc/snippets" \
+ "${Eigen_SOURCE_DIR}/unsupported/doc/examples" \
+ "${Eigen_BINARY_DIR}/unsupported/doc/examples"
# If the value of the EXAMPLE_PATH tag contains directories, you can use the
# EXAMPLE_PATTERNS tag to specify one or more wildcard pattern (like *.cpp
Index: unsupported/doc/CMakeLists.txt
===================================================================
--- unsupported/doc/CMakeLists.txt (revision 0)
+++ unsupported/doc/CMakeLists.txt (revision 0)
@@ -0,0 +1,2 @@
+
+add_subdirectory(examples)
Index: unsupported/doc/examples/BVH_Example.cpp
===================================================================
--- unsupported/doc/examples/BVH_Example.cpp (revision 0)
+++ unsupported/doc/examples/BVH_Example.cpp (revision 0)
@@ -0,0 +1,45 @@
+#include <unsupported/Eigen/BVH>
+
+using namespace Eigen;
+typedef AlignedBox<double, 2> Box2d;
+
+Box2d ei_bounding_box(const Vector2d &v) { return Box2d(v, v); } //compute the bounding box of a single point
+
+struct PointPointMinimizer //how to compute squared distances between points and rectangles
+{
+ PointPointMinimizer() : calls(0) {}
+ typedef double Scalar;
+
+ double minimumOnVolumeVolume(const Box2d &r1, const Box2d &r2) { ++calls; return r1.squaredExteriorDistance(r2); }
+ double minimumOnVolumeObject(const Box2d &r, const Vector2d &v) { ++calls; return r.squaredExteriorDistance(v); }
+ double minimumOnObjectVolume(const Vector2d &v, const Box2d &r) { ++calls; return r.squaredExteriorDistance(v); }
+ double minimumOnObjectObject(const Vector2d &v1, const Vector2d &v2) { ++calls; return (v1 - v2).squaredNorm(); }
+
+ int calls;
+};
+
+int main()
+{
+ std::vector<Vector2d> redPoints, bluePoints;
+ for(int i = 0; i < 100; ++i) { //initialize random set of red points and blue points
+ redPoints.push_back(Vector2d::Random());
+ bluePoints.push_back(Vector2d::Random());
+ }
+
+ PointPointMinimizer minimizer;
+ double minDistSq = std::numeric_limits<double>::max();
+
+ //brute force to find closest red-blue pair
+ for(int i = 0; i < (int)redPoints.size(); ++i)
+ for(int j = 0; j < (int)bluePoints.size(); ++j)
+ minDistSq = std::min(minDistSq, minimizer.minimumOnObjectObject(redPoints[i], bluePoints[j]));
+ std::cout << "Brute force distance = " << sqrt(minDistSq) << ", calls = " << minimizer.calls << std::endl;
+
+ //using BVH to find closest red-blue pair
+ minimizer.calls = 0;
+ KdBVH<double, 2, Vector2d> redTree(redPoints.begin(), redPoints.end()), blueTree(bluePoints.begin(), bluePoints.end()); //construct the trees
+ minDistSq = BVMinimize(redTree, blueTree, minimizer); //actual BVH minimization call
+ std::cout << "BVH distance = " << sqrt(minDistSq) << ", calls = " << minimizer.calls << std::endl;
+
+ return 0;
+}
Index: unsupported/doc/examples/CMakeLists.txt
===================================================================
--- unsupported/doc/examples/CMakeLists.txt (revision 0)
+++ unsupported/doc/examples/CMakeLists.txt (revision 0)
@@ -0,0 +1,17 @@
+FILE(GLOB examples_SRCS "*.cpp")
+
+ADD_CUSTOM_TARGET(unsupported_examples)
+
+FOREACH(example_src ${examples_SRCS})
+GET_FILENAME_COMPONENT(example ${example_src} NAME_WE)
+ADD_EXECUTABLE(${example} ${example_src})
+GET_TARGET_PROPERTY(example_executable
+ ${example} LOCATION)
+ADD_CUSTOM_COMMAND(
+ TARGET ${example}
+ POST_BUILD
+ COMMAND ${example_executable}
+ ARGS >${CMAKE_CURRENT_BINARY_DIR}/${example}.out
+)
+ADD_DEPENDENCIES(unsupported_examples ${example})
+ENDFOREACH(example_src)
Index: unsupported/Eigen/BVH
===================================================================
--- unsupported/Eigen/BVH (revision 0)
+++ unsupported/Eigen/BVH (revision 0)
@@ -0,0 +1,112 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra. Eigen itself is part of the KDE project.
+//
+// Copyright (C) 2009 Ilya Baran <ibaran@xxxxxxx>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_BVH_MODULE_H
+#define EIGEN_BVH_MODULE_H
+
+#include <Eigen/Geometry>
+#include <Eigen/StdVector>
+#include <algorithm>
+#include <queue>
+
+namespace Eigen {
+
+/** \ingroup Unsupported_modules
+ * \defgroup BVH_Module BVH module
+ * \brief This module provides generic bounding volume hierarchy algorithms
+ * and reference tree implementations.
+ *
+ *
+ * \code
+ * #include <unsupported/Eigen/BVH>
+ * \endcode
+ *
+ * A bounding volume hierarchy (BVH) can accelerate many geometric queries. This module provides a generic implementation
+ * of the two basic algorithms over a BVH: intersection of a query object against all objects in the hierarchy and minimization
+ * of a function over the objects in the hierarchy. It also provides intersection and minimization over a cartesian product of
+ * two BVH's. A BVH accelerates intersection by using the fact that if a query object does not intersect a volume, then it cannot
+ * intersect any object contained in that volume. Similarly, a BVH accelerates minimization because the minimum of a function
+ * over a volume is no greater than the minimum of a function over any object contained in it.
+ *
+ * Some sample queries that can be written in terms of intersection are:
+ * - Determine all points where a ray intersects a triangle mesh
+ * - Given a set of points, determine which are contained in a query sphere
+ * - Given a set of spheres, determine which contain the query point
+ * - Given a set of spheres, determine if any is completely contained in a query box (not an intersection query,
+ but can still be accelerated by pruning all spheres that do not intersect the query box)
+ * - Given a set of points, count how many pairs are \f$d\pm\epsilon\f$ apart (done by looking at the cartesian product of the set
+ * of points with itself)
+ *
+ * Some sample queries that can be written in terms of function minimization over a set of objects are:
+ * - Find the intersection between a ray and a triangle mesh closest to the ray origin (function is infinite off the ray)
+ * - Given a polyline and a query point, determine the closest point on the polyline to the query
+ * - Find the diameter of a point cloud (done by looking at the cartesian product and using negative distance as the function)
+ * - Determine how far two meshes are from colliding (this is also a cartesian product query)
+ *
+ * This implementation decouples the basic algorithms both from the type of hierarchy (and the types of the bounding volumes) and
+ * from the particulars of the query. To enable abstraction from the BVH, the BVH is required to implement a generic mechanism
+ * for traversal. To abstract from the query, the query is responsible for keeping track of results.
+ *
+ * To be used in the algorithms, a hierarchy must implement the following traversal mechanism (see KdBVH for a sample implementation): \code
+ typedef Volume //the type of bounding volume
+ typedef Object //the type of object in the hierarchy
+ typedef Index //a reference to a node in the hierarchy--typically an int or a pointer
+ typedef VolumeIterator //an iterator type over node children--returns Index
+ typedef ObjectIterator //an iterator over object (leaf) children--returns const Object &
+ Index getRootIndex() const //returns the index of the hierarchy root
+ const Volume &getVolume(Index index) const //returns the bounding volume of the node at given index
+ void getChildren(Index index, VolumeIterator &outVBegin, VolumeIterator &outVEnd,
+ ObjectIterator &outOBegin, ObjectIterator &outOEnd) const
+ //getChildren takes a node index and makes [outVBegin, outVEnd) range over its node children
+ //and [outOBegin, outOEnd) range over its object children
+ \endcode
+ *
+ * To use the hierarchy, call BVIntersect or BVMinimize, passing it a BVH (or two, for cartesian product) and a minimizer or intersector.
+ * For an intersection query on a single BVH, the intersector encapsulates the query and must provide two functions:
+ * \code
+ bool intersectVolume(const Volume &volume) //returns true if the query intersects the volume
+ bool intersectObject(const Object &object) //returns true if the intersection search should terminate immediately
+ \endcode
+ * The guarantee that BVIntersect provides is that intersectObject will be called on every object whose bounding volume
+ * intersects the query (but possibly on other objects too) unless the search is terminated prematurely. It is the
+ * responsibility of the intersectObject function to keep track of the results in whatever manner is appropriate.
+ * The cartesian product intersection and the BVMinimize queries are similar--see their individual documentation.
+ *
+ * \addexample BVH_Example \label How to use a BVH to find the closest pair between two point sets
+ *
+ * The following is a simple but complete example for how to use the BVH to accelerate the search for a closest red-blue point pair:
+ * \include BVH_Example.cpp
+ * Output: \verbinclude BVH_Example.out
+
+ */
+//@{
+
+#include "src/BVH/BVAlgorithms.h"
+#include "src/BVH/KdBVH.h"
+
+//@}
+
+}
+
+#endif // EIGEN_BVH_MODULE_H
Index: unsupported/Eigen/src/BVH/KdBVH.h
===================================================================
--- unsupported/Eigen/src/BVH/KdBVH.h (revision 0)
+++ unsupported/Eigen/src/BVH/KdBVH.h (revision 0)
@@ -0,0 +1,225 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra. Eigen itself is part of the KDE project.
+//
+// Copyright (C) 2009 Ilya Baran <ibaran@xxxxxxx>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef KDBVH_H_INCLUDED
+#define KDBVH_H_INCLUDED
+
+//internal pair class for the BVH--used instead of std::pair because of alignment
+template<typename Scalar, int Dim>
+struct ei_vector_int_pair
+{
+EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(Scalar, Dim)
+ typedef Matrix<Scalar, Dim, 1> VectorType;
+
+ ei_vector_int_pair(const VectorType &v, int i) : first(v), second(i) {}
+
+ VectorType first;
+ int second;
+};
+
+//these templates help the tree initializer get the bounding boxes either from a provided
+//iterator range or using ei_bounding_box in a unified way
+template<typename Object, typename Volume, typename BoxIter>
+struct ei_get_boxes_helper {
+ void operator()(const std::vector<Object> &objects, BoxIter boxBegin, BoxIter boxEnd, std::vector<Volume> &outBoxes)
+ {
+ outBoxes.insert(outBoxes.end(), boxBegin, boxEnd);
+ ei_assert(outBoxes.size() == objects.size());
+ }
+};
+
+template<typename Object, typename Volume>
+struct ei_get_boxes_helper<Object, Volume, int> {
+ void operator()(const std::vector<Object> &objects, int, int, std::vector<Volume> &outBoxes)
+ {
+ outBoxes.reserve(objects.size());
+ for(int i = 0; i < (int)objects.size(); ++i)
+ outBoxes.push_back(ei_bounding_box(objects[i]));
+ }
+};
+
+
+/** \ingroup BVTree_module
+ * \brief A simple bounding volume hierarchy based on AlignedBox
+ *
+ * \param _Scalar The underlying scalar type of the bounding boxes
+ * \param _Dim The dimension of the space in which the hierarchy lives
+ * \param _Object The object type that lives in the hierarchy. It must have value semantics. Either ei_bounding_box(_Object) must
+ * be defined and return an AlignedBox<_Scalar, _Dim> or bounding boxes must be provided to the tree initializer.
+ *
+ * This class provides a simple (as opposed to optimized) implementation of a bounding volume hierarchy analogous to a Kd-tree.
+ * Given a sequence of objects, it computes their bounding boxes, constructs a Kd-tree of their centers
+ * and builds a BVH with the structure of that Kd-tree. When the elements of the tree are too expensive to be copied around,
+ * it is useful for _Object to be a pointer.
+ */
+template<typename _Scalar, int _Dim, typename _Object> class KdBVH
+{
+public:
+ enum { Dim = _Dim };
+ typedef _Object Object;
+ typedef _Scalar Scalar;
+ typedef AlignedBox<Scalar, Dim> Volume;
+ typedef int Index;
+ typedef const int *VolumeIterator; //the iterators are just pointers into the tree's vectors
+ typedef const Object *ObjectIterator;
+
+ KdBVH() {}
+
+ /** Given an iterator range over \a Object references, constructs the BVH. Requires that ei_bounding_box(Object) return a Volume. */
+ template<typename Iter> KdBVH(Iter begin, Iter end) { init(begin, end, 0, 0); } //int is recognized by init as not being an iterator type
+
+ /** Given an iterator range over \a Object references and an iterator range over their bounding boxes, constructs the BVH */
+ template<typename OIter, typename BIter> KdBVH(OIter begin, OIter end, BIter boxBegin, BIter boxEnd) { init(begin, end, boxBegin, boxEnd); }
+
+ /** Given an iterator range over \a Object references, constructs the BVH, overwriting whatever is in there currently.
+ * Requires that ei_bounding_box(Object) return a Volume. */
+ template<typename Iter> void init(Iter begin, Iter end) { init(begin, end, 0, 0); }
+
+ /** Given an iterator range over \a Object references and an iterator range over their bounding boxes,
+ * constructs the BVH, overwriting whatever is in there currently. */
+ template<typename OIter, typename BIter> void init(OIter begin, OIter end, BIter boxBegin, BIter boxEnd)
+ {
+ objects.clear();
+ boxes.clear();
+ children.clear();
+
+ objects.insert(objects.end(), begin, end);
+ int n = objects.size();
+
+ if(n < 2)
+ return; //if we have at most one object, we don't need any internal nodes
+
+ std::vector<Volume> objBoxes;
+ std::vector<VIPair> objCenters;
+
+ ei_get_boxes_helper<Object, Volume, BIter>()(objects, boxBegin, boxEnd, objBoxes); //compute the bounding boxes depending on BIter type
+
+ objCenters.reserve(n);
+ boxes.reserve(n - 1);
+ children.reserve(2 * n - 2);
+
+ for(int i = 0; i < n; ++i)
+ objCenters.push_back(VIPair(objBoxes[i].center(), i));
+
+ build(objCenters, 0, n, objBoxes, 0); //the recursive part of the algorithm
+
+ std::vector<Object> tmp(n);
+ tmp.swap(objects);
+ for(int i = 0; i < n; ++i)
+ objects[i] = tmp[objCenters[i].second];
+ }
+
+ /** \returns the index of the root of the hierarchy */
+ inline Index getRootIndex() const { return (int)boxes.size() - 1; }
+
+ /** Given an \a index of a node, on exit, \a outVBegin and \a outVEnd range over the indices of the volume children of the node
+ * and \a outOBegin and \a outOEnd range over the object children of the node */
+ EIGEN_STRONG_INLINE void getChildren(Index index, VolumeIterator &outVBegin, VolumeIterator &outVEnd,
+ ObjectIterator &outOBegin, ObjectIterator &outOEnd) const
+ { //inlining this function should open lots of optimization opportunities to the compiler
+ if(index < 0) {
+ outVBegin = outVEnd;
+ if(!objects.empty())
+ outOBegin = &(objects[0]);
+ outOEnd = outOBegin + objects.size(); //output all objects--necessary when the tree has only one object
+ return;
+ }
+
+ int numBoxes = boxes.size();
+
+ int idx = index * 2;
+ if(children[idx + 1] < numBoxes) { //second index is always bigger
+ outVBegin = &(children[idx]);
+ outVEnd = outVBegin + 2;
+ outOBegin = outOEnd;
+ }
+ else if(children[idx] >= numBoxes) { //if both children are objects
+ outVBegin = outVEnd;
+ outOBegin = &(objects[children[idx] - numBoxes]);
+ outOEnd = outOBegin + 2;
+ } else { //if the first child is a volume and the second is an object
+ outVBegin = &(children[idx]);
+ outVEnd = outVBegin + 1;
+ outOBegin = &(objects[children[idx + 1] - numBoxes]);
+ outOEnd = outOBegin + 1;
+ }
+ }
+
+ /** \returns the bounding box of the node at \a index */
+ inline const Volume &getVolume(Index index) const
+ {
+ return boxes[index];
+ }
+
+private:
+ typedef ei_vector_int_pair<Scalar, Dim> VIPair;
+ typedef Matrix<Scalar, Dim, 1> VectorType;
+ struct VectorComparator //compares vectors, or, more specificall, VIPairs along a particular dimension
+ {
+ VectorComparator(int inDim) : dim(inDim) {}
+ inline bool operator()(const VIPair &v1, const VIPair &v2) const { return v1.first[dim] < v2.first[dim]; }
+ int dim;
+ };
+
+ //Build the part of the tree between objects[from] and objects[to] (not including objects[to]).
+ //This routine partitions the objCenters in [from, to) along the dimension dim, recursively constructs
+ //the two halves, and adds their parent node. TODO: a cache-friendlier layout
+ void build(std::vector<VIPair> &objCenters, int from, int to, const std::vector<Volume> &objBoxes, int dim)
+ {
+ ei_assert(to - from > 1);
+ if(to - from == 2) {
+ boxes.push_back(objBoxes[objCenters[from].second] | objBoxes[objCenters[from + 1].second]);
+ children.push_back(from + (int)objects.size() - 1); //there are objects.size() - 1 tree nodes
+ children.push_back(from + (int)objects.size());
+ }
+ else if(to - from == 3) {
+ int mid = from + 2;
+ std::nth_element(objCenters.begin() + from, objCenters.begin() + mid,
+ objCenters.begin() + to, VectorComparator(dim)); //partition
+ build(objCenters, from, mid, objBoxes, (dim + 1) % Dim);
+ int idx1 = (int)boxes.size() - 1;
+ boxes.push_back(boxes[idx1] | objBoxes[objCenters[mid].second]);
+ children.push_back(idx1);
+ children.push_back(mid + (int)objects.size() - 1);
+ }
+ else {
+ int mid = from + (to - from) / 2;
+ nth_element(objCenters.begin() + from, objCenters.begin() + mid,
+ objCenters.begin() + to, VectorComparator(dim)); //partition
+ build(objCenters, from, mid, objBoxes, (dim + 1) % Dim);
+ int idx1 = (int)boxes.size() - 1;
+ build(objCenters, mid, to, objBoxes, (dim + 1) % Dim);
+ int idx2 = (int)boxes.size() - 1;
+ boxes.push_back(boxes[idx1] | boxes[idx2]);
+ children.push_back(idx1);
+ children.push_back(idx2);
+ }
+ }
+
+ std::vector<int> children; //children of x are children[2x] and children[2x+1], indices bigger than boxes.size() index into objects.
+ std::vector<Volume> boxes;
+ std::vector<Object> objects;
+};
+
+#endif //KDBVH_H_INCLUDED
Index: unsupported/Eigen/src/BVH/BVAlgorithms.h
===================================================================
--- unsupported/Eigen/src/BVH/BVAlgorithms.h (revision 0)
+++ unsupported/Eigen/src/BVH/BVAlgorithms.h (revision 0)
@@ -0,0 +1,276 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra. Eigen itself is part of the KDE project.
+//
+// Copyright (C) 2009 Ilya Baran <ibaran@xxxxxxx>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_BVALGORITHMS_H
+#define EIGEN_BVALGORITHMS_H
+
+/** Given a BVH, runs the query encapsulated by \a intersector.
+ * The Intersector type must provide the following members: \code
+ bool intersectVolume(const BVH::Volume &volume) //returns true if volume intersects the query
+ bool intersectObject(const BVH::Object &object) //returns true if the search should terminate immediately
+ \endcode
+ */
+template<typename BVH, typename Intersector>
+void BVIntersect(const BVH &tree, Intersector &intersector)
+{
+ ei_intersect_helper(tree, intersector, tree.getRootIndex());
+}
+
+#ifndef EIGEN_PARSED_BY_DOXYGEN
+template<typename BVH, typename Intersector>
+bool ei_intersect_helper(const BVH &tree, Intersector &intersector, typename BVH::Index root)
+{
+ typedef typename BVH::Index Index;
+
+ typename BVH::VolumeIterator vBegin, vEnd;
+ typename BVH::ObjectIterator oBegin, oEnd;
+
+ std::vector<Index> todo(1, root);
+
+ while(!todo.empty()) {
+ tree.getChildren(todo.back(), vBegin, vEnd, oBegin, oEnd);
+ todo.pop_back();
+
+ for(; vBegin != vEnd; ++vBegin) //go through child volumes
+ if(intersector.intersectVolume(tree.getVolume(*vBegin)))
+ todo.push_back(*vBegin);
+
+ for(; oBegin != oEnd; ++oBegin) //go through child objects
+ if(intersector.intersectObject(*oBegin))
+ return true; //intersector said to stop query
+ }
+ return false;
+}
+#endif //not EIGEN_PARSED_BY_DOXYGEN
+
+template<typename Volume1, typename Object1, typename Object2, typename Intersector>
+struct ei_intersector_helper1
+{
+ ei_intersector_helper1(const Object2 &inStored, Intersector &in) : stored(inStored), intersector(in) {}
+ bool intersectVolume(const Volume1 &vol) { return intersector.intersectVolumeObject(vol, stored); }
+ bool intersectObject(const Object1 &obj) { return intersector.intersectObjectObject(obj, stored); }
+ Object2 stored;
+ Intersector &intersector;
+};
+
+template<typename Volume2, typename Object2, typename Object1, typename Intersector>
+struct ei_intersector_helper2
+{
+ ei_intersector_helper2(const Object1 &inStored, Intersector &in) : stored(inStored), intersector(in) {}
+ bool intersectVolume(const Volume2 &vol) { return intersector.intersectObjectVolume(stored, vol); }
+ bool intersectObject(const Object2 &obj) { return intersector.intersectObjectObject(stored, obj); }
+ Object1 stored;
+ Intersector &intersector;
+};
+
+/** Given two BVH's, runs the query on their cartesian product encapsulated by \a intersector.
+ * The Intersector type must provide the following members: \code
+ bool intersectVolumeVolume(const BVH1::Volume &v1, const BVH2::Volume &v2) //returns true if product of volumes intersects the query
+ bool intersectVolumeObject(const BVH1::Volume &v1, const BVH2::Object &o2) //returns true if the volume-object product intersects the query
+ bool intersectObjectVolume(const BVH1::Object &o1, const BVH2::Volume &v2) //returns true if the volume-object product intersects the query
+ bool intersectObjectObject(const BVH1::Object &o1, const BVH2::Object &o2) //returns true if the search should terminate immediately
+ \endcode
+ */
+template<typename BVH1, typename BVH2, typename Intersector>
+void BVIntersect(const BVH1 &tree1, const BVH2 &tree2, Intersector &intersector) //TODO: tandem descent when it makes sense
+{
+ typedef typename BVH1::Index Index1;
+ typedef typename BVH2::Index Index2;
+ typedef ei_intersector_helper1<typename BVH1::Volume, typename BVH1::Object, typename BVH2::Object, Intersector> Helper1;
+ typedef ei_intersector_helper2<typename BVH2::Volume, typename BVH2::Object, typename BVH1::Object, Intersector> Helper2;
+
+ typename BVH1::VolumeIterator vBegin1, vEnd1;
+ typename BVH1::ObjectIterator oBegin1, oEnd1;
+ typename BVH2::VolumeIterator vBegin2, vEnd2, vCur2;
+ typename BVH2::ObjectIterator oBegin2, oEnd2, oCur2;
+
+ std::vector<std::pair<Index1, Index2> > todo(1, std::make_pair(tree1.getRootIndex(), tree2.getRootIndex()));
+
+ while(!todo.empty()) {
+ tree1.getChildren(todo.back().first, vBegin1, vEnd1, oBegin1, oEnd1);
+ tree2.getChildren(todo.back().second, vBegin2, vEnd2, oBegin2, oEnd2);
+ todo.pop_back();
+
+ for(; vBegin1 != vEnd1; ++vBegin1) { //go through child volumes of first tree
+ const typename BVH1::Volume &vol1 = tree1.getVolume(*vBegin1);
+ for(vCur2 = vBegin2; vCur2 != vEnd2; ++vCur2) { //go through child volumes of second tree
+ if(intersector.intersectVolumeVolume(vol1, tree2.getVolume(*vCur2)))
+ todo.push_back(std::make_pair(*vBegin1, *vCur2));
+ }
+
+ for(oCur2 = oBegin2; oCur2 != oEnd2; ++oCur2) {//go through child objects of second tree
+ Helper1 helper(*oCur2, intersector);
+ if(ei_intersect_helper(tree1, helper, *vBegin1))
+ return; //intersector said to stop query
+ }
+ }
+
+ for(; oBegin1 != oEnd1; ++oBegin1) { //go through child objects of first tree
+ for(vCur2 = vBegin2; vCur2 != vEnd2; ++vCur2) { //go through child volumes of second tree
+ Helper2 helper(*oBegin1, intersector);
+ if(ei_intersect_helper(tree2, helper, *vCur2))
+ return; //intersector said to stop query
+ }
+
+ for(oCur2 = oBegin2; oCur2 != oEnd2; ++oCur2) {//go through child objects of second tree
+ if(intersector.intersectObjectObject(*oBegin1, *oCur2))
+ return; //intersector said to stop query
+ }
+ }
+ }
+}
+
+/** Given a BVH, runs the query encapsulated by \a minimizer.
+ * \returns the minimum value.
+ * The Minimizer type must provide the following members: \code
+ typedef Scalar //the numeric type of what is being minimized--not necessarily the Scalar type of the BVH (if it has one)
+ Scalar minimumOnVolume(const BVH::Volume &volume)
+ Scalar minimumOnObject(const BVH::Object &object)
+ \endcode
+ */
+template<typename BVH, typename Minimizer>
+typename Minimizer::Scalar BVMinimize(const BVH &tree, Minimizer &minimizer)
+{
+ return ei_minimize_helper(tree, minimizer, tree.getRootIndex(), std::numeric_limits<typename Minimizer::Scalar>::max());
+}
+
+#ifndef EIGEN_PARSED_BY_DOXYGEN
+template<typename BVH, typename Minimizer>
+typename Minimizer::Scalar ei_minimize_helper(const BVH &tree, Minimizer &minimizer, typename BVH::Index root, typename Minimizer::Scalar minimum)
+{
+ typedef typename Minimizer::Scalar Scalar;
+ typedef typename BVH::Index Index;
+ typedef std::pair<Scalar, Index> QueueElement; //first element is priority
+
+ typename BVH::VolumeIterator vBegin, vEnd;
+ typename BVH::ObjectIterator oBegin, oEnd;
+ std::priority_queue<QueueElement, std::vector<QueueElement>, std::greater<QueueElement> > todo; //smallest is at the top
+
+ todo.push(std::make_pair(Scalar(), root));
+
+ while(!todo.empty()) {
+ tree.getChildren(todo.top().second, vBegin, vEnd, oBegin, oEnd);
+ todo.pop();
+
+ for(; oBegin != oEnd; ++oBegin) //go through child objects
+ minimum = std::min(minimum, minimizer.minimumOnObject(*oBegin));
+
+ for(; vBegin != vEnd; ++vBegin) { //go through child volumes
+ Scalar val = minimizer.minimumOnVolume(tree.getVolume(*vBegin));
+ if(val < minimum)
+ todo.push(std::make_pair(val, *vBegin));
+ }
+ }
+
+ return minimum;
+}
+#endif //not EIGEN_PARSED_BY_DOXYGEN
+
+
+template<typename Volume1, typename Object1, typename Object2, typename Minimizer>
+struct ei_minimizer_helper1
+{
+ typedef typename Minimizer::Scalar Scalar;
+ ei_minimizer_helper1(const Object2 &inStored, Minimizer &m) : stored(inStored), minimizer(m) {}
+ Scalar minimumOnVolume(const Volume1 &vol) { return minimizer.minimumOnVolumeObject(vol, stored); }
+ Scalar minimumOnObject(const Object1 &obj) { return minimizer.minimumOnObjectObject(obj, stored); }
+ Object2 stored;
+ Minimizer &minimizer;
+};
+
+template<typename Volume2, typename Object2, typename Object1, typename Minimizer>
+struct ei_minimizer_helper2
+{
+ typedef typename Minimizer::Scalar Scalar;
+ ei_minimizer_helper2(const Object1 &inStored, Minimizer &m) : stored(inStored), minimizer(m) {}
+ Scalar minimumOnVolume(const Volume2 &vol) { return minimizer.minimumOnObjectVolume(stored, vol); }
+ Scalar minimumOnObject(const Object2 &obj) { return minimizer.minimumOnObjectObject(stored, obj); }
+ Object1 stored;
+ Minimizer &minimizer;
+};
+
+/** Given two BVH's, runs the query on their cartesian product encapsulated by \a minimizer.
+ * \returns the minimum value.
+ * The Minimizer type must provide the following members: \code
+ typedef Scalar //the numeric type of what is being minimized--not necessarily the Scalar type of the BVH (if it has one)
+ Scalar minimumOnVolumeVolume(const BVH1::Volume &v1, const BVH2::Volume &v2)
+ Scalar minimumOnVolumeObject(const BVH1::Volume &v1, const BVH2::Object &o2)
+ Scalar minimumOnObjectVolume(const BVH1::Object &o1, const BVH2::Volume &v2)
+ Scalar minimumOnObjectObject(const BVH1::Object &o1, const BVH2::Object &o2)
+ \endcode
+ */
+template<typename BVH1, typename BVH2, typename Minimizer>
+typename Minimizer::Scalar BVMinimize(const BVH1 &tree1, const BVH2 &tree2, Minimizer &minimizer)
+{
+ typedef typename Minimizer::Scalar Scalar;
+ typedef typename BVH1::Index Index1;
+ typedef typename BVH2::Index Index2;
+ typedef ei_minimizer_helper1<typename BVH1::Volume, typename BVH1::Object, typename BVH2::Object, Minimizer> Helper1;
+ typedef ei_minimizer_helper2<typename BVH2::Volume, typename BVH2::Object, typename BVH1::Object, Minimizer> Helper2;
+ typedef std::pair<Scalar, std::pair<Index1, Index2> > QueueElement; //first element is priority
+
+ typename BVH1::VolumeIterator vBegin1, vEnd1;
+ typename BVH1::ObjectIterator oBegin1, oEnd1;
+ typename BVH2::VolumeIterator vBegin2, vEnd2, vCur2;
+ typename BVH2::ObjectIterator oBegin2, oEnd2, oCur2;
+ std::priority_queue<QueueElement, std::vector<QueueElement>, std::greater<QueueElement> > todo; //smallest is at the top
+
+ Scalar minimum = std::numeric_limits<Scalar>::max();
+ todo.push(std::make_pair(Scalar(), std::make_pair(tree1.getRootIndex(), tree2.getRootIndex())));
+
+ while(!todo.empty()) {
+ tree1.getChildren(todo.top().second.first, vBegin1, vEnd1, oBegin1, oEnd1);
+ tree2.getChildren(todo.top().second.second, vBegin2, vEnd2, oBegin2, oEnd2);
+ todo.pop();
+
+ for(; oBegin1 != oEnd1; ++oBegin1) { //go through child objects of first tree
+ for(oCur2 = oBegin2; oCur2 != oEnd2; ++oCur2) {//go through child objects of second tree
+ minimum = std::min(minimum, minimizer.minimumOnObjectObject(*oBegin1, *oCur2));
+ }
+
+ for(vCur2 = vBegin2; vCur2 != vEnd2; ++vCur2) { //go through child volumes of second tree
+ Helper2 helper(*oBegin1, minimizer);
+ minimum = std::min(minimum, ei_minimize_helper(tree2, helper, *vCur2, minimum));
+ }
+ }
+
+ for(; vBegin1 != vEnd1; ++vBegin1) { //go through child volumes of first tree
+ const typename BVH1::Volume &vol1 = tree1.getVolume(*vBegin1);
+
+ for(oCur2 = oBegin2; oCur2 != oEnd2; ++oCur2) {//go through child objects of second tree
+ Helper1 helper(*oCur2, minimizer);
+ minimum = std::min(minimum, ei_minimize_helper(tree1, helper, *vBegin1, minimum));
+ }
+
+ for(vCur2 = vBegin2; vCur2 != vEnd2; ++vCur2) { //go through child volumes of second tree
+ Scalar val = minimizer.minimumOnVolumeVolume(vol1, tree2.getVolume(*vCur2));
+ if(val < minimum)
+ todo.push(std::make_pair(val, std::make_pair(*vBegin1, *vCur2)));
+ }
+ }
+ }
+ return minimum;
+}
+
+#endif // EIGEN_BVALGORITHMS_H
Index: unsupported/Eigen/src/BVH/CMakeLists.txt
===================================================================
--- unsupported/Eigen/src/BVH/CMakeLists.txt (revision 0)
+++ unsupported/Eigen/src/BVH/CMakeLists.txt (revision 0)
@@ -0,0 +1,6 @@
+FILE(GLOB Eigen_BVH_SRCS "*.h")
+
+INSTALL(FILES
+ ${Eigen_BVH_SRCS}
+ DESTINATION ${INCLUDE_INSTALL_DIR}/unsupported/Eigen/src/BVH COMPONENT Devel
+ )
Index: unsupported/Eigen/src/CMakeLists.txt
===================================================================
--- unsupported/Eigen/src/CMakeLists.txt (revision 935303)
+++ unsupported/Eigen/src/CMakeLists.txt (working copy)
@@ -1 +1,2 @@
ADD_SUBDIRECTORY(IterativeSolvers)
+ADD_SUBDIRECTORY(BVH)
Index: unsupported/Eigen/CMakeLists.txt
===================================================================
--- unsupported/Eigen/CMakeLists.txt (revision 935303)
+++ unsupported/Eigen/CMakeLists.txt (working copy)
@@ -1,4 +1,4 @@
-set(Eigen_HEADERS AdolcForward IterativeSolvers)
+set(Eigen_HEADERS AdolcForward BVH IterativeSolvers)
install(FILES
${Eigen_HEADERS}
Index: unsupported/CMakeLists.txt
===================================================================
--- unsupported/CMakeLists.txt (revision 935303)
+++ unsupported/CMakeLists.txt (working copy)
@@ -1,6 +1,8 @@
add_subdirectory(Eigen)
+add_subdirectory(doc)
+
if(EIGEN_BUILD_TESTS)
add_subdirectory(test)
endif(EIGEN_BUILD_TESTS)
Index: doc/unsupported_modules.dox
===================================================================
--- doc/unsupported_modules.dox (revision 935303)
+++ doc/unsupported_modules.dox (working copy)
@@ -16,4 +16,7 @@
/** \ingroup Unsupported_modules
* \defgroup IterativeSolvers_Module Iterative solvers module */
+/** \ingroup Unsupported_modules
+ * \defgroup BVH_Module BVH module */
+
} // namespace Eigen
\ No newline at end of file
Index: doc/CMakeLists.txt
===================================================================
--- doc/CMakeLists.txt (revision 935303)
+++ doc/CMakeLists.txt (working copy)
@@ -62,7 +62,7 @@
)
add_dependencies(doc-eigen all_snippets all_examples)
-add_dependencies(doc-unsupported doc-eigen)
+add_dependencies(doc-unsupported doc-eigen unsupported_examples)
# rerun doxygen to get eigen => unsupported cross references
add_custom_target(doc ALL COMMAND doxygen WORKING_DIRECTORY ${Eigen_BINARY_DIR}/doc)
add_dependencies(doc doc-eigen doc-unsupported)