AABB Tree 官方文档链接:CGAL 5.5 - 3D Fast Intersection and Distance Computation (AABB Tree): User Manual
1 介绍 AABB树提供了一个静态的数据结构和算法,能够对有限3D几何对象集合进行高效的相交和距离查询。
相交查询可以是任何类型,前提是在traits类中实现了相应的交集谓词和构造函数。 距离查询仅限于点的查询。 AABB树的数据结构将几何数据的迭代器范围作为输入,然后将其转换为primitives(图元)。在这些primitives中,构造了一个轴对齐边界框(axis-aligned bounding boxes)(AABB)的层次结构,用于加速相交和距离查询。每个图元都能访问一个输入几何对象(datum)和该对象的参考id。例如,一个图元将3D triangle作为datum,多面体表面的face handle作为id。而通过AABB tree进行相交和距离查询时,返回值中就包含了相交对象/最近点和相交图元id/最近图元id。
左图为表面三角网格模型,右图为其构建的AABB树。
2 接口 相交:
代码语言:javascript复制AABB_tree::do_intersect()
AABB_tree::number_of_intersected_primitives()
AABB_tree::all_intersected_primitives()
AABB_tree::any_intersected_primitive()
AABB_tree::first_intersected_primitive()
以上函数不会构造相交对象,仅做测试。
代码语言:javascript复制AABB_tree::all_intersections()
AABB_tree::any_intersection()
AABB_tree::first_intersection()
以上函数会构造相交的对象。
距离:
代码语言:javascript复制AABB_tree::closest_point()
AABB_tree::closest_point_and_primitive()
AABB_tree::accelerate_distance_queries()
注意,在AABB Tree中应避免出现退化的图元,防止算法出错。
3 几个栗子 下面例子中,三维三角形集合以list的形式存储。AABB图元将三角形(triangle)作为datum(数据),list里的迭代器作为id。程序中实现了射线与三角形集合的相交查询,点与三角形集合的最近点查询和距离计算。
代码语言:javascript复制// Author(s) : Camille Wormser, Pierre Alliez
#include <iostream>
#include <list>
#include <CGAL/Simple_cartesian.h>
#include <CGAL/AABB_tree.h>
#include <CGAL/AABB_traits.h>
#include <CGAL/AABB_triangle_primitive.h>
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::Triangle_3 Triangle;
typedef std::list<Triangle>::iterator Iterator; //三角形list迭代器
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;
int main()
{
Point a(1.0, 0.0, 0.0);
Point b(0.0, 1.0, 0.0);
Point c(0.0, 0.0, 1.0);
Point d(0.0, 0.0, 0.0);
std::list<Triangle> triangles;
triangles.push_back(Triangle(a,b,c));
triangles.push_back(Triangle(a,b,d));
triangles.push_back(Triangle(a,d,c));
// constructs AABB tree
Tree tree(triangles.begin(),triangles.end());
// counts #intersections
Ray ray_query(a,b);
std::cout << tree.number_of_intersected_primitives(ray_query)
<< " intersections(s) with ray query" << std::endl;
// compute closest point and squared distance
Point point_query(2.0, 2.0, 2.0);
Point closest_point = tree.closest_point(point_query);
std::cerr << "closest point is: " << closest_point << std::endl;
FT sqd = tree.squared_distance(point_query);
std::cout << "squared distance: " << sqd << std::endl;
return EXIT_SUCCESS;
}
在下面这个例子中,将创建一个多面体三角面片的AABB树。其中,AABB图元将三角形面片句柄包装为id,对应的面片作为几何对象(datum)。
代码语言:javascript复制// Author(s) : Camille Wormser, Pierre Alliez
#include <iostream>
#include <list>
#include <CGAL/Simple_cartesian.h>
#include <CGAL/AABB_tree.h>
#include <CGAL/AABB_traits.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/AABB_face_graph_triangle_primitive.h>
typedef CGAL::Simple_cartesian<double> K;
typedef K::Point_3 Point;
typedef K::Plane_3 Plane;
typedef K::Vector_3 Vector;
typedef K::Segment_3 Segment;
typedef K::Ray_3 Ray;
typedef CGAL::Polyhedron_3<K> Polyhedron;
typedef CGAL::AABB_face_graph_triangle_primitive<Polyhedron> Primitive;
typedef CGAL::AABB_traits<K, Primitive> Traits;
typedef CGAL::AABB_tree<Traits> Tree;
typedef boost::optional< Tree::Intersection_and_primitive_id<Segment>::Type > Segment_intersection;
typedef boost::optional< Tree::Intersection_and_primitive_id<Plane>::Type > Plane_intersection;
typedef Tree::Primitive_id Primitive_id;
int main()
{
Point p(1.0, 0.0, 0.0);
Point q(0.0, 1.0, 0.0);
Point r(0.0, 0.0, 1.0);
Point s(0.0, 0.0, 0.0);
Polyhedron polyhedron;
polyhedron.make_tetrahedron(p, q, r, s);
// constructs AABB tree
Tree tree(faces(polyhedron).first, faces(polyhedron).second, polyhedron);
// constructs segment query
Point a(-0.2, 0.2, -0.2);
Point b(1.3, 0.2, 1.3);
Segment segment_query(a,b);
// tests intersections with segment query
if(tree.do_intersect(segment_query))
std::cout << "intersection(s)" << std::endl;
else
std::cout << "no intersection" << std::endl;
// computes #intersections with segment query
std::cout << tree.number_of_intersected_primitives(segment_query)
<< " intersection(s)" << std::endl;
// computes first encountered intersection with segment query
// (generally a point)
Segment_intersection intersection =
tree.any_intersection(segment_query);
if(intersection)
{
// gets intersection object
const Point* p = boost::get<Point>(&(intersection->first));
if(p)
std::cout << "intersection object is a point " << *p << std::endl;
}
// computes all intersections with segment query (as pairs object - primitive_id)
std::list<Segment_intersection> intersections;
tree.all_intersections(segment_query, std::back_inserter(intersections));
// computes all intersected primitives with segment query as primitive ids
std::list<Primitive_id> primitives;
tree.all_intersected_primitives(segment_query, std::back_inserter(primitives));
// constructs plane query
Vector vec(0.0,0.0,1.0);
Plane plane_query(a,vec);
// computes first encountered intersection with plane query
// (generally a segment)
Plane_intersection plane_intersection = tree.any_intersection(plane_query);
if(plane_intersection)
{
if(boost::get<Segment>(&(plane_intersection->first)))
std::cout << "intersection object is a segment" << std::endl;
}
return EXIT_SUCCESS;
}
下面例子先读取一个闭合的多面体表面,然后以每个face的重心为起始点和垂直于face向模型内部的方向作射线,进行一个ray shooting query。
代码语言:javascript复制#include <CGAL/Simple_cartesian.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/AABB_tree.h>
#include <CGAL/AABB_traits.h>
#include <CGAL/AABB_face_graph_triangle_primitive.h>
#include <CGAL/Polygon_mesh_processing/compute_normal.h>
#include <CGAL/Polygon_mesh_processing/orientation.h>
#include <iostream>
#include <fstream>
typedef CGAL::Simple_cartesian<double> K;
typedef K::FT FT;
typedef K::Point_3 Point;
typedef K::Vector_3 Vector;
typedef K::Ray_3 Ray;
typedef CGAL::Surface_mesh<Point> Mesh;
typedef boost::graph_traits<Mesh>::face_descriptor face_descriptor;
typedef boost::graph_traits<Mesh>::halfedge_descriptor halfedge_descriptor;
typedef CGAL::AABB_face_graph_triangle_primitive<Mesh> Primitive;
typedef CGAL::AABB_traits<K, Primitive> Traits;
typedef CGAL::AABB_tree<Traits> Tree;
typedef boost::optional<Tree::Intersection_and_primitive_id<Ray>::Type> Ray_intersection;
struct Skip
{
face_descriptor fd;
Skip(const face_descriptor fd)
: fd(fd)
{}
bool operator()(const face_descriptor& t) const
{ if(t == fd){
std::cerr << "ignore " << t <<std::endl;
};
return(t == fd);
}
};
int main(int argc, char* argv[])
{
const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/tetrahedron.off");
Mesh mesh;
if(!CGAL::IO::read_polygon_mesh(filename, mesh))
{
std::cerr << "Invalid input." << std::endl;
return 1;
}
Tree tree(faces(mesh).first, faces(mesh).second, mesh);
double d = CGAL::Polygon_mesh_processing::is_outward_oriented(mesh)?-1:1;
for(face_descriptor fd : faces(mesh))
{
halfedge_descriptor hd = halfedge(fd,mesh);
Point p = CGAL::centroid(mesh.point(source(hd,mesh)),
mesh.point(target(hd,mesh)),
mesh.point(target(next(hd,mesh),mesh)));
Vector v = CGAL::Polygon_mesh_processing::compute_face_normal(fd,mesh);
Ray ray(p,d * v);
Skip skip(fd);
Ray_intersection intersection = tree.first_intersection(ray, skip);
if(intersection)
{
if(boost::get<Point>(&(intersection->first))){
const Point* p = boost::get<Point>(&(intersection->first) );
std::cout << *p << std::endl;
}
}
}
std::cerr << "done" << std::endl;
return 0;
}
因为重心计算属于浮点运算,因此射线第一个击中的面可能是起点质心所在的面。为了避免此状况,这里需要给first_intersection()传入一个skip functor将此面忽略。
上个例子是计算的射线与mesh的相交,下面这个例子展示如何查询一个点到mesh的squared distance和closest point及其所在的triangle。
代码语言:javascript复制// Author(s) : Pierre Alliez
#include <iostream>
#include <CGAL/Simple_cartesian.h>
#include <CGAL/AABB_tree.h>
#include <CGAL/AABB_traits.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/AABB_face_graph_triangle_primitive.h>
typedef CGAL::Simple_cartesian<double> K;
typedef K::FT FT;
typedef K::Point_3 Point;
typedef K::Segment_3 Segment;
typedef CGAL::Polyhedron_3<K> Polyhedron;
typedef CGAL::AABB_face_graph_triangle_primitive<Polyhedron> Primitive;
typedef CGAL::AABB_traits<K, Primitive> Traits;
typedef CGAL::AABB_tree<Traits> Tree;
typedef Tree::Point_and_primitive_id Point_and_primitive_id;
int main()
{
Point p(1.0, 0.0, 0.0);
Point q(0.0, 1.0, 0.0);
Point r(0.0, 0.0, 1.0);
Point s(0.0, 0.0, 0.0);
Polyhedron polyhedron;
polyhedron.make_tetrahedron(p, q, r, s);
// constructs AABB tree and computes internal KD-tree
// data structure to accelerate distance queries
Tree tree(faces(polyhedron).first, faces(polyhedron).second, polyhedron);
// query point
Point query(0.0, 0.0, 3.0);
// computes squared distance from query
FT sqd = tree.squared_distance(query);
std::cout << "squared distance: " << sqd << std::endl;
// computes closest point
Point closest = tree.closest_point(query);
std::cout << "closest point: " << closest << std::endl;
// computes closest point and primitive id
Point_and_primitive_id pp = tree.closest_point_and_primitive(query);
Point closest_point = pp.first;
Polyhedron::Face_handle f = pp.second; // closest primitive id
std::cout << "closest point: " << closest_point << std::endl;
std::cout << "closest triangle: ( "
<< f->halfedge()->vertex()->point() << " , "
<< f->halfedge()->next()->vertex()->point() << " , "
<< f->halfedge()->next()->next()->vertex()->point()
<< " )" << std::endl;
return EXIT_SUCCESS;
}
最后一个例子,是对于AABB树中图元的增量插入。虽然AABB树是一个静态的数据结构,但是它允许插入primitives(图元)。
代码语言:javascript复制#include <iostream>
#include <CGAL/Simple_cartesian.h>
#include <CGAL/AABB_tree.h>
#include <CGAL/AABB_traits.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/AABB_face_graph_triangle_primitive.h>
typedef CGAL::Simple_cartesian<double> K;
typedef K::FT FT;
typedef K::Point_3 Point;
typedef K::Segment_3 Segment;
typedef CGAL::Polyhedron_3<K> Polyhedron;
typedef CGAL::AABB_face_graph_triangle_primitive<Polyhedron, CGAL::Default, CGAL::Tag_false> Primitive;
typedef CGAL::AABB_traits<K, Primitive> Traits;
typedef CGAL::AABB_tree<Traits> Tree;
typedef Tree::Point_and_primitive_id Point_and_primitive_id;
int main()
{
Point p(1.0, 0.0, 0.0);
Point q(0.0, 1.0, 0.0);
Point r(0.0, 0.0, 1.0);
Point s(0.0, 0.0, 0.0);
Polyhedron polyhedron1;
polyhedron1.make_tetrahedron(p, q, r, s);
Point p2(11.0, 0.0, 0.0);
Point q2(10.0, 1.0, 0.0);
Point r2(10.0, 0.0, 1.0);
Point s2(10.0, 0.0, 0.0);
Polyhedron polyhedron2;
polyhedron2.make_tetrahedron(p2, q2, r2, s2);
// constructs AABB tree and computes internal KD-tree
// data structure to accelerate distance queries
Tree tree(faces(polyhedron1).first, faces(polyhedron1).second, polyhedron1);
tree.insert(faces(polyhedron2).first, faces(polyhedron2).second, polyhedron2);
// query point
Point query(0.0, 0.0, 3.0);
// computes squared distance from query
FT sqd = tree.squared_distance(query);
std::cout << "squared distance: " << sqd << std::endl;
return EXIT_SUCCESS;
}
上面这个例子中,首先使用polyhedron1构建tree,然后使用insert()函数将polyhedron2的faces作为primitives插入到tree中。