AABB tree demo: added menu to benchmark memory against number of facets.

I had to add a recursive longest edge bisection algorithm to refine a mesh without changing its geometry.

AABB_tree/demo/AABB_tree/MainWindow.cpp

@@ -283,3 +283,24 @@ void MainWindow::on_actionClear_distance_function_triggered()
 {
 	m_pScene->clear_distance_function();
 }
+
+void MainWindow::on_actionRefine_bisection_triggered()
+{
+	bool ok;
+	const double max_len = 
+		QInputDialog::getDouble(NULL, "Max edge len",
+		"Max edge len:",0.1,0.001,100.0,9,&ok);
+	if(!ok)
+		return;
+
+	QApplication::setOverrideCursor(Qt::WaitCursor);
+	m_pScene->refine_bisection(max_len * max_len);
+	QApplication::restoreOverrideCursor();
+}
+
+void MainWindow::on_actionBench_memory_triggered()
+{
+	QApplication::setOverrideCursor(Qt::WaitCursor);
+	m_pScene->bench_memory();
+	QApplication::restoreOverrideCursor();
+}

AABB_tree/demo/AABB_tree/MainWindow.h

@@ -50,12 +50,14 @@ public:
 			void on_actionEdge_points_triggered();
 			void on_actionInside_points_triggered();
 			void on_actionBoundary_points_triggered();
+			void on_actionRefine_bisection_triggered();
 			void on_actionBoundary_segments_triggered();
 			void on_actionSigned_distance_function_to_facets_triggered();
 			void on_actionUnsigned_distance_function_to_edges_triggered();
 			void on_actionUnsigned_distance_function_to_facets_triggered();
 
 			// benchmark menu
+			void on_actionBench_memory_triggered();
 			void on_actionBench_distances_triggered();
 			void on_actionBench_intersections_triggered();
 

AABB_tree/demo/AABB_tree/MainWindow.ui

@@ -60,6 +60,12 @@
     <property name="title" >
      <string>Algorithms</string>
     </property>
+    <widget class="QMenu" name="menuRefine" >
+     <property name="title" >
+      <string>Refine</string>
+     </property>
+     <addaction name="actionRefine_bisection" />
+    </widget>
     <addaction name="actionEdge_points" />
     <addaction name="actionInside_points" />
     <addaction name="actionBoundary_points" />
@@ -68,11 +74,14 @@
     <addaction name="actionSigned_distance_function_to_facets" />
     <addaction name="actionUnsigned_distance_function_to_facets" />
     <addaction name="actionUnsigned_distance_function_to_edges" />
+    <addaction name="separator" />
+    <addaction name="menuRefine" />
    </widget>
    <widget class="QMenu" name="menuBenchmarks" >
     <property name="title" >
      <string>Benchmark</string>
     </property>
+    <addaction name="actionBench_memory" />
     <addaction name="actionBench_distances" />
     <addaction name="actionBench_intersections" />
    </widget>
@@ -186,6 +195,21 @@
     <string>Clear distance function</string>
    </property>
   </action>
+  <action name="actionRefine_bisection" >
+   <property name="text" >
+    <string>Longest edge bisection</string>
+   </property>
+  </action>
+  <action name="actionLoop_subdivision" >
+   <property name="text" >
+    <string>Loop subdivision</string>
+   </property>
+  </action>
+  <action name="actionBench_memory" >
+   <property name="text" >
+    <string>Memory</string>
+   </property>
+  </action>
  </widget>
  <customwidgets>
   <customwidget>

AABB_tree/demo/AABB_tree/Refiner.h

@@ -0,0 +1,197 @@
+#ifndef _REFINER_H_
+#define _REFINER_H_
+
+#include <CGAL/basic.h>
+#include <CGAL/Polyhedron_3.h>
+#include <queue>
+
+template <class	Kernel,	class	Polyhedron>
+class CEdge
+{
+public:
+	typedef typename Kernel::FT FT;
+	typedef typename Polyhedron::Halfedge_handle Halfedge_handle;
+
+private:
+	FT m_sqlen;
+	Halfedge_handle m_he;
+
+public:
+	// life cycle
+	CEdge(const Halfedge_handle& he)
+	{
+		m_sqlen = squared_len(he);
+		m_he = he;
+	}
+	CEdge(const CEdge& edge) 
+		: m_sqlen(edge.sqlen()),
+		m_he(edge.he())
+	{
+	}
+	~CEdge() {}
+
+public:
+	// squared length of an edge
+	static FT squared_len(Halfedge_handle he)
+	{
+		return CGAL::squared_distance(he->vertex()->point(),
+			he->opposite()->vertex()->point());
+	}
+
+public:
+	// accessors 
+	FT& sqlen() { return m_sqlen; }
+	const FT sqlen() const { return m_sqlen; }
+
+	Halfedge_handle he() { return m_he; }
+	const Halfedge_handle he() const { return m_he; }
+};
+
+// functor for priority queue
+template<class Edge>
+struct less // read more priority
+{	
+	bool operator()(const Edge& e1, 
+		const Edge& e2) const
+	{	
+		return e1.sqlen() < e2.sqlen();
+	}
+};
+
+template <class	Kernel,	class	Polyhedron>
+class Refiner
+{
+	// types
+	typedef typename Kernel::FT FT;
+	typedef typename CEdge<Kernel,Polyhedron> Edge;
+	typedef typename Polyhedron::Halfedge_handle Halfedge_handle;
+	typedef typename Polyhedron::Edge_iterator   Edge_iterator;
+	typedef typename std::priority_queue<typename Edge,
+		std::vector<typename Edge>,
+		less<typename Edge> > PQueue;
+	// data
+	PQueue m_queue;
+	Polyhedron* m_pMesh;
+
+public :
+	// life cycle
+	Refiner(Polyhedron* pMesh) 
+	{
+		m_pMesh = pMesh;
+	}
+	~Refiner() {}
+
+public :
+
+	void fill_queue(const FT& max_sqlen)
+	{
+		for(Edge_iterator he = m_pMesh->edges_begin();
+			he != m_pMesh->edges_end();
+			he++)
+			if(Edge::squared_len(he) > max_sqlen)
+				m_queue.push(Edge(he));
+	}
+
+	void fill_queue()
+	{
+		for(Edge_iterator he = m_pMesh->edges_begin();
+			he != m_pMesh->edges_end();
+			he++)
+			m_queue.push(Edge(he));
+	}
+
+	// run
+	void run_nb_splits(const unsigned int nb_splits)
+	{
+		// fill queue
+		fill_queue();
+
+		unsigned int nb = 0;
+		while(nb < nb_splits)
+		{
+			// get a copy of the candidate edge
+			Edge edge = m_queue.top(); 
+			m_queue.pop();
+
+			Halfedge_handle he = edge.he();
+			// update point
+			Halfedge_handle hnew = m_pMesh->split_edge(he);
+			hnew->vertex()->point() = CGAL::midpoint(he->vertex()->point(),
+				he->opposite()->vertex()->point());
+
+			// hit has been split into two edges
+			m_queue.push(Edge(hnew));
+			m_queue.push(Edge(he));
+
+			// split facet if possible
+			if(!hnew->is_border())
+			{
+				m_pMesh->split_facet(hnew,hnew->next()->next());
+				m_queue.push(Edge(hnew->next()));
+			}
+
+			// split facet if possible
+			if(!hnew->opposite()->is_border())
+			{
+				m_pMesh->split_facet(hnew->opposite()->next(),
+					hnew->opposite()->next()->next()->next());
+				m_queue.push(Edge(hnew->opposite()->prev()));
+			}
+
+			nb++;
+		} // end while
+	} // end run
+
+
+
+	// run
+	unsigned int operator()(const FT& max_sqlen)
+	{
+		// fill queue
+		fill_queue(max_sqlen);
+
+		unsigned int nb_split = 0;
+		while(!m_queue.empty())
+		{
+			// get a copy of the candidate edge
+			Edge edge = m_queue.top(); 
+			m_queue.pop();
+
+			Halfedge_handle he = edge.he();
+			FT sqlen = Edge::squared_len(he);
+			if(sqlen > max_sqlen)
+			{
+				// update point
+				Halfedge_handle hnew = m_pMesh->split_edge(he);
+				hnew->vertex()->point() = CGAL::midpoint(he->vertex()->point(),
+					he->opposite()->vertex()->point());
+
+				// hit has been split into two edges
+				m_queue.push(Edge(hnew));
+				m_queue.push(Edge(he));
+
+				// split facet if possible
+				if(!hnew->is_border())
+				{
+					m_pMesh->split_facet(hnew,hnew->next()->next());
+					m_queue.push(Edge(hnew->next()));
+				}
+
+				// split facet if possible
+				if(!hnew->opposite()->is_border())
+				{
+					m_pMesh->split_facet(hnew->opposite()->next(),
+						hnew->opposite()->next()->next()->next());
+					m_queue.push(Edge(hnew->opposite()->prev()));
+				}
+
+				nb_split++;
+			} // end if(sqlen > max_sqlen)
+		} // end while(!m_queue.empty())
+		return nb_split;
+	} // end run
+};
+
+#endif // _REFINER_H_
+
+

AABB_tree/demo/AABB_tree/Scene.cpp

@@ -9,6 +9,7 @@
 #include <QTime>
 #include <QInputDialog>
 
+#include "Refiner.h"
 #include "render_edges.h"
 #include <CGAL/IO/Polyhedron_iostream.h>
 
@@ -556,9 +557,6 @@ unsigned_distance : m_max_distance_function;
 	m_signed_distance_function = true;
 }
 
-
-
-
 void Scene::toggle_view_poyhedron()
 {
 	m_view_polyhedron = !m_view_polyhedron;
@@ -579,4 +577,10 @@ void Scene::toggle_view_distance_function()
 	m_view_distance_function = !m_view_distance_function;
 }
 
-
+void Scene::refine_bisection(const FT max_sqlen)
+{
+	std::cout << "Refine through recursive longest edge bisection...";
+	Refiner<Kernel,Polyhedron> refiner(m_pPolyhedron);
+	refiner(max_sqlen);
+	std::cout << "done (" << m_pPolyhedron->size_of_facets() << " facets)" << std::endl;
+}

AABB_tree/demo/AABB_tree/Scene.h

@@ -38,7 +38,6 @@ public:
 
 public:
 	void draw();
-	Polyhedron* polyhedron() const;
 
 	typedef CGAL::Bbox_3 Bbox;
 	Bbox bbox() { return Bbox(); }
@@ -66,6 +65,7 @@ public:
 	void generate_inside_points(const unsigned int nb_points);
 	void generate_boundary_points(const unsigned int nb_points);
 	void generate_boundary_segments(const unsigned int nb_slices);
+	void refine_bisection(const FT max_sqlen);
 
 	// distance functions 
 	void signed_distance_function();
@@ -96,6 +96,8 @@ public:
 		  NB_INTERSECTIONS,
 		  ALL_INTERSECTIONS,
 		  ALL_INTERSECTED_PRIMITIVES};
+	void bench_memory();
+	unsigned int nb_digits(const unsigned int value);
 	void benchmark_intersections(const int duration);
 	void bench_intersection_rays(Facet_tree& tree,const int function,const int duration);
 	void bench_intersection_lines(Facet_tree& tree,const int function,const int duration);

AABB_tree/demo/AABB_tree/benchmarks.cpp

@@ -1,5 +1,6 @@
 #include "Scene.h"
 #include <QInputDialog>
+#include <CGAL/Memory_sizer.h>
 
 void Scene::benchmark_intersections(const int duration)
 {
@@ -35,6 +36,41 @@ void Scene::benchmark_distances(const int duration)
 	bench_closest_point_and_primitive(tree,duration);
 }
 
+unsigned int Scene::nb_digits(unsigned int value)
+{
+	unsigned int nb_digits = 0;
+	while(value > 0)
+	{
+		nb_digits++;
+		value /= 10;
+	}
+	return nb_digits;
+}
+
+// bench memory against number of facets in the tree
+// the tree is reconstructed each time in the mesh 
+// refinement loop
+void Scene::bench_memory()
+{
+	std::cout << std::endl << "Benchmark memory" << std::endl;
+	std::cout << "#Facets MBytes" << std::endl;
+	while(m_pPolyhedron->size_of_facets() < 2000000)
+	{
+		// refine mesh at increasing speed
+		Refiner<Kernel,Polyhedron> refiner(m_pPolyhedron);
+		unsigned int digits = nb_digits(m_pPolyhedron->size_of_facets());
+		unsigned int nb_splits = 0.2 * std::pow(10.0,(double)digits - 1.0);
+		refiner.run_nb_splits(nb_splits);
+
+		// constructs tree and measure memory before then after
+		long before = CGAL::Memory_sizer().virtual_size();
+		Facet_tree tree(m_pPolyhedron->facets_begin(),m_pPolyhedron->facets_end());
+		long after = CGAL::Memory_sizer().virtual_size();
+		double memory = (double)(after - before) / 1048576.0; // in MBytes
+		std::cout << m_pPolyhedron->size_of_facets() << " " << memory << std::endl;
+	}
+}
+
 //////////////////////////////////////////////////////////////
 // BENCH INTERSECTIONS
 //////////////////////////////////////////////////////////////