Beef up the testcase

AABB_tree/test/AABB_tree/aabb_test_ray_intersection.cpp

@@ -4,6 +4,7 @@
 
 #include <boost/functional/value_factory.hpp>
 #include <boost/timer/timer.hpp>
+#include <boost/array.hpp>
 
 #include <CGAL/assertions.h>
 #include <CGAL/algorithm.h>
@@ -14,6 +15,7 @@
 #include <CGAL/AABB_tree.h>
 #include <CGAL/AABB_traits.h>
 #include <CGAL/boost/graph/graph_traits_Polyhedron_3.h>
+#include <CGAL/IO/Polyhedron_iostream.h>
 #include <CGAL/AABB_face_graph_triangle_primitive.h>
 
 typedef CGAL::Epeck K;
@@ -33,6 +35,8 @@ FT point_on_ray_dist(const Ray& ray, const Point& point) {
   return i_ray.squared_length();
 }
 
+static long accum = 0;
+
 boost::optional<
   Tree::template Intersection_and_primitive_id<Ray>::Type
   >
@@ -41,6 +45,7 @@ min_intersection(const Tree& tree, const Ray& ray) {
   IntersectionVector all_intersections;
 
   tree.all_intersections(ray, std::back_inserter(all_intersections));
+  accum += all_intersections.size();
   Tree::FT min_distance = DBL_MAX;
   boost::optional<
     Tree::template Intersection_and_primitive_id<Ray>::Type
@@ -65,31 +70,52 @@ int main()
 {
   Polyhedron polyhedron;
   {
-    Point p(1.0, 0.0, 0.0);
+    // Point p(1.0, 0.0, 0.0);
-    Point q(0.0, 1.0, 0.0);
+    // Point q(0.0, 1.0, 0.0);
-    Point r(0.0, 0.0, 1.0);
+    // Point r(0.0, 0.0, 1.0);
-    Point s(0.0, 0.0, 0.0);
+    // Point s(0.0, 0.0, 0.0);
-    polyhedron.make_tetrahedron(p, q, r, s);
+    // polyhedron.make_tetrahedron(p, q, r, s);
   }
-  // Polyhedron polyhedron;
+
-  // std::ifstream in("../bunny00.off");
+  std::ifstream in("../bunny00.off");
-  // if(in)
+  if(in)
-  //   in >> polyhedron;
+    in >> polyhedron;
-  // else
+  else
-  //   std::cout << "error reading bunny" << std::endl;
+    std::cout << "error reading bunny" << std::endl;
 
   Tree tree(faces(polyhedron).first, faces(polyhedron).second, polyhedron);
+  Tree::Bounding_box bbox = tree.bbox();
+  Vector bbox_center((bbox.xmin() + bbox.xmax()) / 2, 
+                     (bbox.ymin() + bbox.ymax()) / 2, 
+                     (bbox.zmin() + bbox.zmax()) / 2);
+  boost::array<double, 3> extents;
+  extents[0] = bbox.xmax() - bbox.xmin();
+  extents[1] = bbox.ymax() - bbox.ymin();
+  extents[2] = bbox.zmax() - bbox.zmin();
+  double max_extent = *std::max_element(extents.begin(), extents.end());
+
+  std::cout << bbox << std::endl;
+  std::cout << bbox_center << std::endl;
+  std::cout << max_extent << std::endl;
   
   const int NB_RAYS = 1000;
   std::vector<Point> v1, v2;
   v1.reserve(NB_RAYS); v2.reserve(NB_RAYS);
 
-  const float r = 2.0;
+  const float r = max_extent / 2;
-  // Generate NB_RAYS*2 points that lie on a sphere of radius r.
+  // Generate NB_RAYS*2 points that lie on a sphere of radius r, centered around bbox_center
   CGAL::Random rand = CGAL::Random(23); // fix the seed to yield the same results each run
   CGAL::cpp11::copy_n(CGAL::Random_points_on_sphere_3<Point>(r, rand), NB_RAYS, std::back_inserter(v1));
   CGAL::cpp11::copy_n(CGAL::Random_points_on_sphere_3<Point>(r, rand), NB_RAYS, std::back_inserter(v2));
 
+  for(std::vector<Point>::iterator it = v1.begin(); it != v1.end(); ++it) {
+    *it = *it + bbox_center;
+  }
+
+  for(std::vector<Point>::iterator it = v2.begin(); it != v2.end(); ++it) {
+    *it = *it + bbox_center;
+  }
+  
   // Generate NB_RAYS using v1 as source and v2 as target.
   std::vector<Ray> rays;
   rays.reserve(NB_RAYS);
@@ -99,17 +125,22 @@ int main()
   primitives1.reserve(NB_RAYS); primitives2.reserve(NB_RAYS);
 
 
+  {
     for(std::vector<Ray>::iterator it = rays.begin(); it != rays.end(); ++it) {
       primitives1.push_back(min_intersection(tree, *it));
     }
+  }
 
   for(std::vector<Ray>::iterator it = rays.begin(); it != rays.end(); ++it) {
     primitives2.push_back(tree.ray_intersection(*it));
   }
-
   CGAL_assertion_msg(primitives1.size() == primitives2.size(), "Different amount of primitives intersected.");
   CGAL_assertion_msg(std::equal(primitives1.begin(), primitives1.end(), primitives2.begin()),
                      "Primitives mismatch.");
+  std::size_t c = primitives1.size() - std::count(primitives1.begin(), primitives1.end(), boost::none);
+  std::cout << "Intersected " << c << " primitives with " << NB_RAYS << " rays" << std::endl;
+  std::cout << "Primitive method had to sort " << accum/NB_RAYS 
+            << " intersections on average." << std::endl;
 
   return 0;
 }