skbl renaming, new link methods for abstract link of geometrical complex

git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/trunk@468 636b058d-ea47-450e-bf9e-a15bfbe3eedb

Former-commit-id: cb85e1ae86635857fbaaa0781526cb3eeaa9a50a

7 files changed, 120 insertions, 335 deletions

diff --git a/src/Contraction/example/CMakeLists.txt b/src/Contraction/example/CMakeLists.txt
index 6d65de73..4889b82f 100644
--- a/src/Contraction/example/CMakeLists.txt
+++ b/src/Contraction/example/CMakeLists.txt
@@ -6,6 +6,8 @@ add_executable(RipsContraction Rips_contraction.cpp)
 add_executable(GarlandHeckbert Garland_heckbert.cpp)
 
 target_link_libraries(RipsContraction ${Boost_TIMER_LIBRARY} ${Boost_SYSTEM_LIBRARY})	
+target_link_libraries(GarlandHeckbert ${Boost_TIMER_LIBRARY} ${Boost_SYSTEM_LIBRARY})	
+
 
 add_test(RipsContraction.sphere.0.2 ${CMAKE_CURRENT_BINARY_DIR}/RipsContraction ${CMAKE_SOURCE_DIR}/data/points/sphere3D_2646.off 0.2)
 add_test(RipsContraction.S0310000 ${CMAKE_CURRENT_BINARY_DIR}/RipsContraction ${CMAKE_SOURCE_DIR}/data/points/SO3_10000.off 0.3)
diff --git a/src/Contraction/example/Garland_heckbert.cpp b/src/Contraction/example/Garland_heckbert.cpp
index 912199a4..5b178ff9 100644
--- a/src/Contraction/example/Garland_heckbert.cpp
+++ b/src/Contraction/example/Garland_heckbert.cpp
@@ -28,16 +28,15 @@
 #ifndef GARLAND_HECKBERT_H_
 #define GARLAND_HECKBERT_H_
 
-
+#include <boost/timer/timer.hpp>
 #include <iostream>
-#include "Garland_heckbert/Point.h"
-#include "Garland_heckbert/Error_quadric.h"
-#include "Garland_heckbert/Garland_heckbert_policies.h"
-
+#include "gudhi/Point.h"
 #include "gudhi/Edge_contraction.h"
 #include "gudhi/Skeleton_blocker.h"
 #include "gudhi/Off_reader.h"
 
+#include "Garland_heckbert/Error_quadric.h"
+
 using namespace std;
 using namespace Gudhi;
 using namespace skbl;
@@ -48,28 +47,109 @@ struct Geometry_trait{
 	typedef Point_d Point;
 };
 
+/**
+ * The vertex stored in the complex contains a quadric.
+ */
 struct Garland_heckbert_traits : public Skeleton_blocker_simple_geometric_traits<Geometry_trait> {
 
 public:
-	// the vertex stored in the complex contains a quadric
 	struct Garland_heckbert_vertex : public Simple_geometric_vertex{
 		Error_quadric<Geometry_trait::Point> quadric;
 	};
 	typedef Garland_heckbert_vertex Graph_vertex;
-
 };
 
 typedef Skeleton_blocker_geometric_complex< Garland_heckbert_traits > Complex;
-typedef Edge_profile<Complex> Profile;
+typedef Edge_profile<Complex> EdgeProfile;
 typedef Skeleton_blocker_contractor<Complex> Complex_contractor;
 
 
+/**
+ * How the new vertex is placed after an edge collapse : here it is placed at
+ * the point minimizing the cost of the quadric.
+ */
+class GH_placement : public Gudhi::contraction::Placement_policy<EdgeProfile>{
+	Complex& complex_;
+public:
+	typedef typename Gudhi::contraction::Placement_policy<EdgeProfile>::Placement_type Placement_type;
+
+	GH_placement(Complex& complex):complex_(complex){}
+
+	Placement_type operator()(const EdgeProfile& profile) const override{
+		auto sum_quad(profile.v0().quadric);
+		sum_quad += profile.v1().quadric;
+
+		boost::optional<Point> min_quadric_pt(sum_quad.min_cost());
+		if (min_quadric_pt)
+			return Placement_type(*min_quadric_pt);
+		else
+			return profile.p0();
+	}
+};
+
+/**
+ * How much cost an edge collapse : here the costs is given by a quadric
+ * which expresses a squared distances with triangles planes.
+ */
+class GH_cost : public Gudhi::contraction::Cost_policy<EdgeProfile>{
+	Complex& complex_;
+public:
+
+	typedef typename Gudhi::contraction::Cost_policy<EdgeProfile>::Cost_type Cost_type;
+
+	GH_cost(Complex& complex):complex_(complex){}
+
+	Cost_type operator()( EdgeProfile const& profile, boost::optional<Point> const& new_point ) const override {
+		Cost_type res;
+		if (new_point){
+			auto sum_quad(profile.v0().quadric);
+			sum_quad += profile.v1().quadric;
+			res = sum_quad.cost(*new_point);
+		}
+		return res;
+	}
+};
+
+/**
+ * Visitor that is called at several moment.
+ * Here we initializes the quadrics of every vertex at the on_started call back
+ * and we update them when contracting an edge (the quadric become the sum of both quadrics).
+ */
+class GH_visitor: public Gudhi::contraction::Contraction_visitor<EdgeProfile> {
+	Complex& complex_;
+public:
+	GH_visitor(Complex& complex):complex_(complex){}
+
+	//Compute quadrics for every vertex v
+	//The quadric of v consists in the sum of quadric
+	//of every triangles passing through v weighted by its area
+	void on_started(Complex & complex) override{
+		for(auto v : complex.vertex_range()){
+			auto & quadric_v(complex[v].quadric);
+			for(auto t : complex.triangle_range(v)){
+				auto t_it = t.begin();
+				const auto& p0(complex.point(*t_it++));
+				const auto& p1(complex.point(*t_it++));
+				const auto& p2(complex.point(*t_it++));
+				quadric_v+=Error_quadric<Point>(p0,p1,p2);
+			}
+		}
+	}
+
+	/**
+	 * @brief Called when an edge is about to be contracted and replaced by a vertex whose position is *placement.
+	 */
+	void on_contracting(EdgeProfile const &profile, boost::optional< Point > placement)
+	override{
+		profile.v0().quadric += profile.v1().quadric;
+	}
+};
 
 
 int main(int argc, char *argv[]){
-	if (argc!=3){
-		std::cerr << "Usage "<<argv[0]<<" ../../../data/meshes/test.off N to load the file ../../data/test.off and contract N edges.\n";
-		return -1;
+	if (argc!=4){
+		std::cerr << "Usage "<<argv[0]<<" input.off output.off N to load the file input.off, contract N edges and save the result to output.off.\n";
+		return EXIT_FAILURE;
 	}
 
 	Complex complex;
@@ -81,31 +161,31 @@ int main(int argc, char *argv[]){
 		return EXIT_FAILURE;
 	}
 
-	int num_contractions = atoi(argv[2]);
+	std::cout << "Load complex with "<<complex.num_vertices()<<" vertices"<<std::endl;
 
-	std::cout << "Initial complex has "<<
-			complex.num_vertices()<<" vertices, "<<
-			complex.num_blockers()<<" blockers, "<<
-			complex.num_edges()<<" edges and" <<
-			complex.num_triangles()<<" triangles.";
 
+	int num_contractions = atoi(argv[3]);
+
+	boost::timer::auto_cpu_timer t;
+
+	// constructs the contractor object with Garland Heckbert policies.
 	Complex_contractor contractor(complex,
-			new GH_cost<Profile>(complex),
-			new GH_placement<Profile>(complex),
-			contraction::make_link_valid_contraction<Profile>(),
-			new GH_visitor<Profile>(complex)
+			new GH_cost(complex),
+			new GH_placement(complex),
+			contraction::make_link_valid_contraction<EdgeProfile>(),
+			new GH_visitor(complex)
 	);
 
-	std::cout<<"Contract "<<num_contractions<<" edges\n";
+	std::cout<<"Contract "<<num_contractions<<" edges"<<std::endl;
 	contractor.contract_edges(num_contractions);
 
 	std::cout << "Final complex has "<<
 			complex.num_vertices()<<" vertices, "<<
 			complex.num_edges()<<" edges and" <<
-			complex.num_triangles()<<" triangles.";
-
-	Skeleton_blocker_off_writer<Complex> off_writer("reduced.off",complex);
+			complex.num_triangles()<<" triangles."<<std::endl;
 
+	//write simplified complex
+	Skeleton_blocker_off_writer<Complex> off_writer(argv[2],complex);
 
 	return EXIT_SUCCESS;
 }
diff --git a/src/Contraction/example/Garland_heckbert/Error_quadric.h b/src/Contraction/example/Garland_heckbert/Error_quadric.h
index 7a30c35a..725a3a56 100755
--- a/src/Contraction/example/Garland_heckbert/Error_quadric.h
+++ b/src/Contraction/example/Garland_heckbert/Error_quadric.h
@@ -123,11 +123,11 @@ public :
 	 * Det must be passed with the determinant value of the gradient (should be non zero).

 	 */

 	inline Point solve_linear_gradient(double det = grad_determinant()) const{

-		return Point(

+		return Point({

 				(-coeff[1]*coeff[5]*coeff[8]+coeff[1]*coeff[7]*coeff[6]+coeff[2]*coeff[8]*coeff[4]-coeff[2]*coeff[5]*coeff[6]-coeff[3]*coeff[4]*coeff[7]+coeff[3]*coeff[5]*coeff[5])/ det,

 				(coeff[0]*coeff[5]*coeff[8]-coeff[0]*coeff[7]*coeff[6]-coeff[5]*coeff[2]*coeff[3]-coeff[1]*coeff[2]*coeff[8]+coeff[6]*coeff[2]*coeff[2]+coeff[1]*coeff[3]*coeff[7])/det,

 				(-coeff[8]*coeff[0]*coeff[4]+coeff[8]*coeff[1]*coeff[1]+coeff[2]*coeff[3]*coeff[4]+coeff[5]*coeff[0]*coeff[6]-coeff[5]*coeff[1]*coeff[3]-coeff[1]*coeff[2]*coeff[6])/det

-		);

+		});

 	}

 

 

diff --git a/src/Contraction/example/Garland_heckbert/Garland_heckbert_policies.h b/src/Contraction/example/Garland_heckbert/Garland_heckbert_policies.h
deleted file mode 100644
index c792b5b3..00000000
--- a/src/Contraction/example/Garland_heckbert/Garland_heckbert_policies.h
+++ /dev/null
@@ -1,115 +0,0 @@
-/*
- * Garland_heckbert_policies.h
- *  Created on: Feb 10, 2015
- * This file is part of the Gudhi Library. The Gudhi library 
- *    (Geometric Understanding in Higher Dimensions) is a generic C++ 
- *    library for computational topology.
- *
- *    Author(s):       David Salinas
- *
- *    Copyright (C) 2014  INRIA Sophia Antipolis-Méditerranée (France)
- *
- *    This program is free software: 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 3 of the License, or
- *    (at your option) any later version.
- *
- *    This program 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 General Public License for more details.
- *
- *    You should have received a copy of the GNU General Public License
- *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
- * 
- */
-
-
-#ifndef GARLAND_HECKBERT_POLICIES_H_
-#define GARLAND_HECKBERT_POLICIES_H_
-
-#include "gudhi/Edge_contraction.h"
-#include "Error_quadric.h"
-
-template<typename EdgeProfile>
-class GH_visitor: public Gudhi::contraction::Contraction_visitor<EdgeProfile> {
-	typedef typename EdgeProfile::Complex Complex;
-	typedef typename Complex::Point Point;
-	Complex& complex_;
-public:
-	GH_visitor(Complex& complex):complex_(complex){}
-
-
-	void on_started(Complex & complex) override{
-		//Compute quadrics for every vertex v
-		//The quadric of v consists in the sum of quadric
-		//of every triangles passing through v weighted by its area
-		for(auto v : complex.vertex_range()){
-			auto & quadric_v(complex[v].quadric);
-			for(auto t : complex.triangle_range(v)){
-				auto t_it = t.begin();
-				const auto& p0(complex.point(*t_it++));
-				const auto& p1(complex.point(*t_it++));
-				const auto& p2(complex.point(*t_it++));
-				quadric_v+=Error_quadric<Point>(p0,p1,p2);
-			}
-		}
-	}
-
-	/**
-	 * @brief Called when an edge is about to be contracted and replaced by a vertex whose position is *placement.
-	 */
-	void on_contracting(EdgeProfile const &profile, boost::optional< Point > placement)
-	override{
-		profile.v0().quadric += profile.v1().quadric;
-	}
-};
-
-
-template<typename EdgeProfile>
-class GH_placement : public Gudhi::contraction::Placement_policy<EdgeProfile>{
-	typedef typename EdgeProfile::Complex Complex;
-	Complex& complex_;
-public:
-	typedef typename EdgeProfile::Point Point;
-	typedef typename Gudhi::contraction::Placement_policy<EdgeProfile>::Placement_type Placement_type;
-
-	GH_placement(Complex& complex):complex_(complex){}
-
-	Placement_type operator()(const EdgeProfile& profile) const override{
-		auto sum_quad(profile.v0().quadric);
-		sum_quad += profile.v1().quadric;
-
-		boost::optional<Point> min_quadric_pt(sum_quad.min_cost());
-		if (min_quadric_pt)
-			return Placement_type(*min_quadric_pt);
-		else
-			return profile.p0();
-	}
-};
-
-template<typename EdgeProfile>
-class GH_cost : public Gudhi::contraction::Cost_policy<EdgeProfile>{
-	typedef typename EdgeProfile::Complex Complex;
-	Complex& complex_;
-public:
-
-	typedef typename Gudhi::contraction::Cost_policy<EdgeProfile>::Cost_type Cost_type;
-	typedef typename Complex::Point Point;
-
-	GH_cost(Complex& complex):complex_(complex){}
-
-	Cost_type operator()( EdgeProfile const& profile, boost::optional<Point> const& new_point ) const override {
-		Cost_type res;
-		if (new_point){
-			auto sum_quad(profile.v0().quadric);
-			sum_quad += profile.v1().quadric;
-			res = sum_quad.cost(*new_point);
-		}
-		return res;
-	}
-};
-
-
-
-#endif /* GARLAND_HECKBERT_POLICIES_H_ */
diff --git a/src/Contraction/example/Garland_heckbert/Point.h b/src/Contraction/example/Garland_heckbert/Point.h
deleted file mode 100644
index eb250347..00000000
--- a/src/Contraction/example/Garland_heckbert/Point.h
+++ /dev/null
@@ -1,176 +0,0 @@
-/*
- * Basic_geometry.h
- *  Created on: Feb 10, 2015
- * This file is part of the Gudhi Library. The Gudhi library 
- *    (Geometric Understanding in Higher Dimensions) is a generic C++ 
- *    library for computational topology.
- *
- *    Author(s):       David Salinas
- *
- *    Copyright (C) 2014  INRIA Sophia Antipolis-Méditerranée (France)
- *
- *    This program is free software: 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 3 of the License, or
- *    (at your option) any later version.
- *
- *    This program 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 General Public License for more details.
- *
- *    You should have received a copy of the GNU General Public License
- *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
- * 
- */
-
-
-#ifndef BASIC_GEOMETRY_H_
-#define BASIC_GEOMETRY_H_
-
-#include <cmath>
-#include <vector>
-#include <cassert>
-#include <cstddef>
-
-class Point_d{
-public:
-	Point_d():coords_(3,0){}
-	Point_d(size_t dim):coords_(dim,0){}
-	Point_d(const Point_d& other):coords_(other.coords_){}
-	Point_d(double x,double y,double z):coords_({x,y,z}){}
-
-	template<typename CoordsIt>
-	Point_d(CoordsIt begin,CoordsIt end):coords_(begin,end){}
-
-	size_t dimension() const{
-		return coords_.size();
-	}
-
-	double x() const{
-		return coords_[0];
-	}
-
-	double y() const{
-		return coords_[1];
-	}
-
-	double z() const{
-		return coords_[2];
-	}
-
-	double& x(){
-		return coords_[0];
-	}
-
-	double& y(){
-		return coords_[1];
-	}
-
-	double& z(){
-		return coords_[2];
-	}
-
-	std::vector<double>::const_iterator begin() const{
-		return coords_.begin();
-	}
-
-	std::vector<double>::const_iterator end() const{
-		return coords_.end();
-	}
-
-	double& operator[](unsigned i){
-		return coords_[i];
-	}
-	const double& operator[](unsigned i) const{
-		return coords_[i];
-	}
-
-	double squared_norm() const{
-		double res = 0;
-		for(auto x : coords_)
-			res+= x*x;
-		return res;
-	}
-
-	double squared_dist(const Point_d& other) const{
-		assert(dimension()==other.dimension());
-		double res = 0;
-		for(unsigned i = 0; i < coords_.size(); ++i)
-			res+= coords_[i]*coords_[i] + other[i]*other[i];
-		return res;
-	}
-
-	/**
-	 * dot product
-	 */
-	double operator*(const Point_d& other) const{
-		assert(dimension()==other.dimension());
-		double res = 0;
-		for(unsigned i = 0; i < coords_.size(); ++i)
-			res+= coords_[i]*other[i];
-		return res;
-	}
-
-	/**
-	 * only if points have dimension 3
-	 */
-	Point_d cross_product(const Point_d& other){
-		assert(dimension()==3 && other.dimension()==3);
-		Point_d res(3);
-		res[0] = (*this)[1] * other[2] - (*this)[2] * other[1];
-		res[1] = (*this)[2] * other[0] - (*this)[0] * other[2];
-		res[2] = (*this)[0] * other[1] - (*this)[1] * other[0];
-		return res;
-	}
-
-	Point_d operator+(const Point_d& other) const{
-		assert(dimension()==other.dimension());
-		Point_d res(dimension());
-		for(unsigned i = 0; i < coords_.size(); ++i)
-			res[i] = (*this)[i] + other[i];
-		return res;
-	}
-
-	Point_d operator*(double lambda) const{
-		Point_d res(dimension());
-		for(unsigned i = 0; i < coords_.size(); ++i)
-			res[i] = (*this)[i] * lambda;
-		return res;
-	}
-
-	Point_d operator/(double lambda) const{
-		Point_d res(dimension());
-		for(unsigned i = 0; i < coords_.size(); ++i)
-			res[i] = (*this)[i] / lambda;
-		return res;
-	}
-
-	Point_d operator-(const Point_d& other) const{
-		assert(dimension()==other.dimension());
-		Point_d res(dimension());
-		for(unsigned i = 0; i < coords_.size(); ++i)
-			res[i] = (*this)[i] - other[i];
-		return res;
-	}
-
-	friend Point_d unit_normal(const Point_d& p1,const Point_d& p2,const Point_d& p3){
-		assert(p1.dimension()==3);
-		assert(p2.dimension()==3);
-		assert(p3.dimension()==3);
-		Point_d p1p2 = p2 - p1;
-		Point_d p1p3 = p3 - p1;
-		Point_d res(p1p2.cross_product(p1p3));
-		return res / std::sqrt(res.squared_norm());
-	}
-
-
-private:
-	std::vector<double> coords_;
-};
-
-
-
-
-
-#endif /* BASIC_GEOMETRY_H_ */
diff --git a/src/Contraction/example/Rips_contraction.cpp b/src/Contraction/example/Rips_contraction.cpp
index 8d9b1ca2..bd0a8b8c 100644
--- a/src/Contraction/example/Rips_contraction.cpp
+++ b/src/Contraction/example/Rips_contraction.cpp
@@ -24,42 +24,34 @@
 #include "gudhi/Edge_contraction.h"
 #include "gudhi/Skeleton_blocker.h"
 #include "gudhi/Off_reader.h"
-
+#include "gudhi/Point.h"
 
 using namespace std;
 using namespace Gudhi;
 using namespace skbl;
 using namespace contraction;
 
+
 struct Geometry_trait{
-	typedef std::vector<double> Point;
+	typedef Point_d Point;
 };
 
+
 typedef Geometry_trait::Point Point;
 typedef Skeleton_blocker_simple_geometric_traits<Geometry_trait> Complex_geometric_traits;
 typedef Skeleton_blocker_geometric_complex< Complex_geometric_traits > Complex;
 typedef Edge_profile<Complex> Profile;
 typedef Skeleton_blocker_contractor<Complex> Complex_contractor;
 
-template<typename Point>
-double eucl_distance(const Point& a,const Point& b){
-	double res = 0;
-	auto a_coord = a.begin();
-	auto b_coord = b.begin();
-	for(; a_coord != a.end(); a_coord++, b_coord++){
-		res += (*a_coord - *b_coord) * (*a_coord - *b_coord);
-	}
-	return sqrt(res);
-}
-
 template<typename ComplexType>
 void build_rips(ComplexType& complex, double offset){
 	if (offset<=0) return;
 	auto vertices = complex.vertex_range();
 	for (auto p = vertices.begin(); p != vertices.end(); ++p)
-		for (auto q = p; ++q != vertices.end(); /**/)
-			if (eucl_distance(complex.point(*p), complex.point(*q)) < 2 * offset)
+		for (auto q = p; ++q != vertices.end(); /**/){
+			if ( squared_dist(complex.point(*p), complex.point(*q)) < 4 * offset * offset)
 				complex.add_edge(*p,*q);
+		}
 }
 
 int main (int argc, char *argv[])
@@ -71,13 +63,14 @@ int main (int argc, char *argv[])
 
 	Complex complex;
 
-	// load the points
+	// load only the points
 	Skeleton_blocker_off_reader<Complex> off_reader(argv[1],complex,true);
 	if(!off_reader.is_valid()){
 		std::cerr << "Unable to read file:"<<argv[1]<<std::endl;
 		return EXIT_FAILURE;
 	}
-	std::cout << "Build the Rips complex"<<std::endl;
+
+	std::cout << "Build the Rips complex with "<<complex.num_vertices()<<" vertices"<<std::endl;
 
 	build_rips(complex,atof(argv[2]));
 
diff --git a/src/Contraction/include/gudhi/Skeleton_blocker_contractor.h b/src/Contraction/include/gudhi/Skeleton_blocker_contractor.h
index ba2138a1..dcc05c22 100644
--- a/src/Contraction/include/gudhi/Skeleton_blocker_contractor.h
+++ b/src/Contraction/include/gudhi/Skeleton_blocker_contractor.h
@@ -338,6 +338,7 @@ private:
 	}
 
 	bool is_contraction_valid( Profile const& profile, Placement_type placement ) const{
+		if(!valid_contraction_policy_) return true;
 		return (*valid_contraction_policy_)(profile,placement);
 	}