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

15 files changed, 261 insertions, 230 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/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);
 	}
 
diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_off_io.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_off_io.h
index aaa682c3..c98b0b45 100644
--- a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_off_io.h
+++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_off_io.h
@@ -134,7 +134,7 @@ public:
 	Skeleton_blocker_off_reader(const std::string & name_file, Complex& read_complex, bool read_only_points = false,bool is_flag = false):valid_(false) {
 		std::ifstream stream(name_file);
 		if (stream.is_open()) {
-			if(is_flag){
+			if(is_flag || read_only_points){
 				Skeleton_blocker_off_flag_visitor_reader<Complex> off_visitor(read_complex, read_only_points);
 				Off_reader off_reader(stream);
 				valid_ = off_reader.read(off_visitor);
@@ -184,8 +184,8 @@ public:
 			for(auto v : save_complex.vertex_range()){
 				vertex_num[v]=current_vertex++;
 				const auto& pt(save_complex.point(v));
-				for(auto it = pt.begin();it!=pt.end();++it)
-					stream<<*it<<" ";
+				for(auto x : pt)
+					stream<<x<<" ";
 				stream<<std::endl;
 			}
 
diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_edges_iterators.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_edges_iterators.h
index feae3d57..0be6c74d 100644
--- a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_edges_iterators.h
+++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_edges_iterators.h
@@ -30,8 +30,8 @@ namespace Gudhi{
 namespace skbl {
 
 template<typename SkeletonBlockerComplex>
-class Complex_edge_around_vertex_iterator :
- public boost::iterator_facade < Complex_edge_around_vertex_iterator<SkeletonBlockerComplex>
+class Edge_around_vertex_iterator :
+ public boost::iterator_facade < Edge_around_vertex_iterator<SkeletonBlockerComplex>
 		, typename SkeletonBlockerComplex::Edge_handle const
 		, boost::forward_traversal_tag
 		, typename SkeletonBlockerComplex::Edge_handle const
@@ -54,10 +54,10 @@ private:
 
 public:
 
-	Complex_edge_around_vertex_iterator():complex(NULL){
+	Edge_around_vertex_iterator():complex(NULL){
 	}
 
-	Complex_edge_around_vertex_iterator(const Complex* complex_,Vertex_handle v_):
+	Edge_around_vertex_iterator(const Complex* complex_,Vertex_handle v_):
 		complex(complex_),
 		v(v_)
 	{
@@ -67,7 +67,7 @@ public:
 	/**
 	 * returns an iterator to the end
 	 */
-	Complex_edge_around_vertex_iterator(const Complex* complex_,Vertex_handle v_,int end):
+	Edge_around_vertex_iterator(const Complex* complex_,Vertex_handle v_,int end):
 		complex(complex_),
 		v(v_)
 	{
@@ -75,7 +75,7 @@ public:
 		set_end();
 	}
 
-	bool equal(const Complex_edge_around_vertex_iterator& other) const{
+	bool equal(const Edge_around_vertex_iterator& other) const{
 		return (complex== other.complex) && (v == other.v) && (current_ == other.current_) && (end_ == other.end_);
 	}
 
@@ -102,8 +102,8 @@ private:
  *
  */
 template<typename SkeletonBlockerComplex>
-class Complex_edge_iterator :
-public boost::iterator_facade < Complex_edge_iterator<SkeletonBlockerComplex>
+class Edge_iterator :
+public boost::iterator_facade < Edge_iterator<SkeletonBlockerComplex>
 , typename SkeletonBlockerComplex::Edge_handle const
 , boost::forward_traversal_tag
 , typename SkeletonBlockerComplex::Edge_handle const
@@ -120,10 +120,10 @@ public:
 	const Complex* complex;
 	std::pair<boost_edge_iterator,boost_edge_iterator> edge_iterator ;
 
-	Complex_edge_iterator():complex(NULL){
+	Edge_iterator():complex(NULL){
 	}
 
-	Complex_edge_iterator(const SkeletonBlockerComplex* complex_):
+	Edge_iterator(const SkeletonBlockerComplex* complex_):
 		complex(complex_),
 		edge_iterator(boost::edges(complex_->skeleton))
 	{
@@ -132,7 +132,7 @@ public:
 	/**
 	 * return an iterator to the end
 	 */
-	Complex_edge_iterator(const SkeletonBlockerComplex* complex_,int end):
+	Edge_iterator(const SkeletonBlockerComplex* complex_,int end):
 		complex(complex_),
 		edge_iterator(boost::edges(complex_->skeleton))
 	{
@@ -140,7 +140,7 @@ public:
 	}
 
 
-	bool equal(const Complex_edge_iterator& other) const{
+	bool equal(const Edge_iterator& other) const{
 		return (complex == other.complex) && (edge_iterator == other.edge_iterator);
 	}
 
diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_simplices_iterators.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_simplices_iterators.h
index 0b397f56..666ce430 100644
--- a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_simplices_iterators.h
+++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_simplices_iterators.h
@@ -67,7 +67,6 @@ class Simplex_around_vertex_iterator :
 	typedef typename Link::Vertex_handle Link_vertex_handle;
 	// Link_vertex_handle == Complex_Vertex_handle but this renaming helps avoiding confusion
 
-
 	typedef typename Gudhi::skbl::Trie<Simplex_handle> Trie;
 
 
@@ -97,8 +96,7 @@ public:
 		v(other.v),
 		link_v(other.link_v),
 		trie(other.trie),
-		nodes_to_be_seen(other.nodes_to_be_seen)
-	{
+		nodes_to_be_seen(other.nodes_to_be_seen){
 		if(!other.is_end()){
 		}
 	}
@@ -237,7 +235,7 @@ class Simplex_iterator :
 	typedef typename Complex::Edge_handle Edge_handle;
 	typedef typename Complex::Simplex_handle Simplex_handle;
 
-	typedef typename Complex::CVI CVI;
+	typedef typename Complex::Complex_vertex_iterator Complex_vertex_iterator;
 
 
 	typedef typename Link::Vertex_handle Link_vertex_handle;
@@ -245,7 +243,7 @@ class Simplex_iterator :
 private:
 
 	const Complex* complex_;
-	CVI current_vertex_;
+	Complex_vertex_iterator current_vertex_;
 
 	typedef Simplex_around_vertex_iterator<SkeletonBlockerComplex,Link> SAVI;
 	SAVI current_simplex_around_current_vertex_;
diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_triangles_iterators.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_triangles_iterators.h
index 1fdbdfc9..e137d1ea 100644
--- a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_triangles_iterators.h
+++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_triangles_iterators.h
@@ -49,7 +49,7 @@ private:
 	typedef typename LinkType::Vertex_handle Vertex_handle;
 	typedef typename LinkType::Root_vertex_handle Root_vertex_handle;
 	typedef typename LinkType::Simplex_handle Simplex_handle;
-	typedef Complex_edge_iterator<Complex> Complex_edge_iterator_;
+	typedef typename Complex::Complex_edge_iterator Complex_edge_iterator_;
 
 	const Complex* complex_;
 	Vertex_handle v_;
@@ -133,9 +133,10 @@ private:
 	typedef typename SkeletonBlockerComplex::Root_vertex_handle Root_vertex_handle;
 	typedef typename SkeletonBlockerComplex::Simplex_handle Simplex_handle;
 	typedef typename SkeletonBlockerComplex::Superior_triangle_around_vertex_iterator STAVI;
+	typedef typename SkeletonBlockerComplex::Complex_vertex_iterator Complex_vertex_iterator;
 
 	const SkeletonBlockerComplex* complex_;
-	Complex_vertex_iterator<SkeletonBlockerComplex> current_vertex_;
+	Complex_vertex_iterator current_vertex_;
 	STAVI current_triangle_;
 	bool is_end_;
 public:
@@ -176,7 +177,7 @@ public:
 
 	Triangle_iterator& operator=(const Triangle_iterator & other){
 		complex_ = other.complex_;
-		Complex_vertex_iterator<SkeletonBlockerComplex> current_vertex_;
+		Complex_vertex_iterator current_vertex_;
 		STAVI current_triangle_;
 		return *this;
 	}
diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_vertices_iterators.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_vertices_iterators.h
index 4c90ee51..a9d4e373 100644
--- a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_vertices_iterators.h
+++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_vertices_iterators.h
@@ -35,8 +35,8 @@ namespace skbl {
  *@remark Incrementation increases Vertex_handle.
  */
 template<typename SkeletonBlockerComplex>
-class Complex_vertex_iterator : public boost::iterator_facade
-< Complex_vertex_iterator <SkeletonBlockerComplex>
+class Vertex_iterator : public boost::iterator_facade
+< Vertex_iterator <SkeletonBlockerComplex>
   , typename SkeletonBlockerComplex::Vertex_handle const
   , boost::forward_traversal_tag
   , typename SkeletonBlockerComplex::Vertex_handle const
@@ -52,10 +52,10 @@ private:
 
 
 public:
-	Complex_vertex_iterator():complex(NULL){
+	Vertex_iterator():complex(NULL){
 	}
 
-	Complex_vertex_iterator(const SkeletonBlockerComplex* complex_):
+	Vertex_iterator(const SkeletonBlockerComplex* complex_):
 		complex(complex_),
 		vertexIterator(vertices(complex_->skeleton)){
 		if(!finished() && !is_active()) {
@@ -66,7 +66,7 @@ public:
 	/**
 	 * return an iterator to the end.
 	 */
-	Complex_vertex_iterator(const SkeletonBlockerComplex* complex_,int end):
+	Vertex_iterator(const SkeletonBlockerComplex* complex_,int end):
 		complex(complex_),vertexIterator(vertices(complex_->skeleton)){
 		vertexIterator.first = vertexIterator.second ;
 	}
@@ -77,11 +77,11 @@ public:
 		return(Vertex_handle(*(vertexIterator.first)));
 	}
 
-	bool equal(const Complex_vertex_iterator& other) const{
+	bool equal(const Vertex_iterator& other) const{
 		return vertexIterator == other.vertexIterator &&  complex == other.complex;
 	}
 
-	bool operator<(const Complex_vertex_iterator& other) const{
+	bool operator<(const Vertex_iterator& other) const{
 		return dereference()<other.dereference();
 	}
 
@@ -107,8 +107,8 @@ private:
 
 
 template<typename SkeletonBlockerComplex>
-class Complex_neighbors_vertices_iterator
-: public boost::iterator_facade < Complex_neighbors_vertices_iterator<SkeletonBlockerComplex>
+class Neighbors_vertices_iterator
+: public boost::iterator_facade < Neighbors_vertices_iterator<SkeletonBlockerComplex>
   , typename SkeletonBlockerComplex::Vertex_handle const
   , boost::forward_traversal_tag
   , typename SkeletonBlockerComplex::Vertex_handle const
@@ -133,10 +133,10 @@ public:
 	//	boost_adjacency_iterator ai, ai_end;
 	//	for (tie(ai, ai_end) = adjacent_vertices(v.vertex, skeleton); ai != ai_end; ++ai){
 
-	Complex_neighbors_vertices_iterator():complex(NULL){
+	Neighbors_vertices_iterator():complex(NULL){
 	}
 
-	Complex_neighbors_vertices_iterator(const Complex* complex_,Vertex_handle v_):
+	Neighbors_vertices_iterator(const Complex* complex_,Vertex_handle v_):
 		complex(complex_),
 		v(v_){
 		tie(current_,end_) = adjacent_vertices(v.vertex, complex->skeleton);
@@ -145,7 +145,7 @@ public:
 	/**
 	 * returns an iterator to the end
 	 */
-	Complex_neighbors_vertices_iterator(const Complex* complex_,Vertex_handle v_,int end):
+	Neighbors_vertices_iterator(const Complex* complex_,Vertex_handle v_,int end):
 		complex(complex_),
 		v(v_){
 		tie(current_,end_) = adjacent_vertices(v.vertex, complex->skeleton);
@@ -162,7 +162,7 @@ public:
 		return(Vertex_handle(*current_));
 	}
 
-	bool equal(const Complex_neighbors_vertices_iterator& other) const{
+	bool equal(const Neighbors_vertices_iterator& other) const{
 		return (complex== other.complex) && (v == other.v) && (current_ == other.current_) && (end_ == other.end_);
 	}
 
diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_complex.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_complex.h
index 893c65c6..5b2c1f27 100644
--- a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_complex.h
+++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_complex.h
@@ -63,10 +63,10 @@ namespace skbl {
  */
 template<class SkeletonBlockerDS>
 class Skeleton_blocker_complex {
-	template<class ComplexType> friend class Complex_vertex_iterator;
-	template<class ComplexType> friend class Complex_neighbors_vertices_iterator;
-	template<class ComplexType> friend class Complex_edge_iterator;
-	template<class ComplexType> friend class Complex_edge_around_vertex_iterator;
+	template<class ComplexType> friend class Vertex_iterator;
+	template<class ComplexType> friend class Neighbors_vertices_iterator;
+	template<class ComplexType> friend class Edge_iterator;
+	template<class ComplexType> friend class Edge_around_vertex_iterator;
 
 	template<class ComplexType> friend class Skeleton_blocker_link_complex;
 	template<class ComplexType> friend class Skeleton_blocker_link_superior;
@@ -587,7 +587,6 @@ public:
 		if (found) {
 			if (visitor)
 				visitor->on_remove_edge(a, b);
-			// if (heapCollapse.Contains(edge)) heapCollapse.Delete(edge);
 			boost::remove_edge(a.vertex, b.vertex, skeleton);
 			degree_[a.vertex]--;
 			degree_[b.vertex]--;
@@ -1237,35 +1236,34 @@ public:
 	/** @name Vertex iterators
 	 */
 	//@{
-	typedef Complex_vertex_iterator<Skeleton_blocker_complex> CVI;  // todo rename
+	typedef Vertex_iterator<Skeleton_blocker_complex> Complex_vertex_iterator;
 
 	//
-	// @brief Range over the vertices of the simplicial complex.
+	// Range over the vertices of the simplicial complex.
 	// Methods .begin() and .end() return a Complex_vertex_iterator.
 	//
-	typedef boost::iterator_range<
-			Complex_vertex_iterator<Skeleton_blocker_complex> > Complex_vertex_range;
+	typedef boost::iterator_range<Complex_vertex_iterator> Complex_vertex_range;
 
 	/**
 	 * @brief Returns a Complex_vertex_range over all vertices of the complex
 	 */
 	Complex_vertex_range vertex_range() const {
-		auto begin = Complex_vertex_iterator<Skeleton_blocker_complex>(this);
-		auto end = Complex_vertex_iterator<Skeleton_blocker_complex>(this, 0);
+		auto begin = Complex_vertex_iterator(this);
+		auto end = Complex_vertex_iterator(this, 0);
 		return Complex_vertex_range(begin, end);
 	}
 
-	typedef boost::iterator_range<
-			Complex_neighbors_vertices_iterator<Skeleton_blocker_complex> > Complex_neighbors_vertices_range;
+	typedef Neighbors_vertices_iterator<Skeleton_blocker_complex> Complex_neighbors_vertices_iterator;
+
+
+	typedef boost::iterator_range<Complex_neighbors_vertices_iterator> Complex_neighbors_vertices_range;
 
 	/**
 	 * @brief Returns a Complex_edge_range over all edges of the simplicial complex that passes trough v
 	 */
 	Complex_neighbors_vertices_range vertex_range(Vertex_handle v) const {
-		auto begin = Complex_neighbors_vertices_iterator<Skeleton_blocker_complex>(
-				this, v);
-		auto end = Complex_neighbors_vertices_iterator<Skeleton_blocker_complex>(
-				this, v, 0);
+		auto begin = Complex_neighbors_vertices_iterator(this, v);
+		auto end = Complex_neighbors_vertices_iterator(this, v, 0);
 		return Complex_neighbors_vertices_range(begin, end);
 	}
 
@@ -1275,28 +1273,33 @@ public:
 	 */
 	//@{
 
-	typedef boost::iterator_range<
-			Complex_edge_iterator<Skeleton_blocker_complex<SkeletonBlockerDS>>> Complex_edge_range;
+	typedef Edge_iterator<Skeleton_blocker_complex> Complex_edge_iterator;
+
+
+	typedef boost::iterator_range<Complex_edge_iterator> Complex_edge_range;
 
 	/**
 	 * @brief Returns a Complex_edge_range over all edges of the simplicial complex
 	 */
 	Complex_edge_range edge_range() const {
-		auto begin = Complex_edge_iterator<Skeleton_blocker_complex < SkeletonBlockerDS >> (this);
-		auto end = Complex_edge_iterator<Skeleton_blocker_complex < SkeletonBlockerDS >> (this, 0);
+		auto begin = Complex_edge_iterator(this);
+		auto end = Complex_edge_iterator(this, 0);
 		return Complex_edge_range(begin, end);
 	}
 
-	typedef boost::iterator_range <Complex_edge_around_vertex_iterator<Skeleton_blocker_complex<SkeletonBlockerDS>>>
-			Complex_edge_around_vertex_range;
+
+	typedef Edge_around_vertex_iterator<Skeleton_blocker_complex> Complex_edge_around_vertex_iterator;
+
+
+	typedef boost::iterator_range <Complex_edge_around_vertex_iterator>	Complex_edge_around_vertex_range;
 
 	/**
 	 * @brief Returns a Complex_edge_range over all edges of the simplicial complex that passes
 	 * through 'v'
 	 */
 	Complex_edge_around_vertex_range edge_range(Vertex_handle v) const {
-		auto begin = Complex_edge_around_vertex_iterator<Skeleton_blocker_complex < SkeletonBlockerDS >> (this, v);
-		auto end = Complex_edge_around_vertex_iterator<Skeleton_blocker_complex < SkeletonBlockerDS >> (this, v, 0);
+		auto begin = Complex_edge_around_vertex_iterator(this, v);
+		auto end = Complex_edge_around_vertex_iterator(this, v, 0);
 		return Complex_edge_around_vertex_range(begin, end);
 	}
 
@@ -1327,6 +1330,10 @@ public:
 
 	typedef boost::iterator_range<Triangle_iterator<Skeleton_blocker_complex> > Complex_triangle_range;
 
+
+	typedef Triangle_iterator<Skeleton_blocker_complex> Complex_triangle_iterator;
+
+
 	/**
 	 * @brief Range over triangles of the simplicial complex.
 	 * Methods .begin() and .end() return a Triangle_around_vertex_iterator.
diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_geometric_complex.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_geometric_complex.h
index bb58d0dc..ee56ce3b 100644
--- a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_geometric_complex.h
+++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_geometric_complex.h
@@ -92,7 +92,8 @@ class Skeleton_blocker_geometric_complex :
 		unsigned current = 0;
 		complex=make_complex_from_top_faces<Skeleton_blocker_geometric_complex>(simplex_begin,simplex_end,is_flag_complex);
 		for(auto point = points_begin; point != points_end; ++point)
-			complex.point(Vertex_handle(current++)) = Point(point->begin(),point->end());
+			//			complex.point(Vertex_handle(current++)) = Point(point->begin(),point->end());
+			complex.point(Vertex_handle(current++)) = Point(*point);
 		return complex;
 	}
 
@@ -185,6 +186,31 @@ class Skeleton_blocker_geometric_complex :
 		return link;
 	}
 
+	typedef Skeleton_blocker_link_complex<Skeleton_blocker_complex<SkeletonBlockerGeometricDS>> Abstract_link;
+
+	/**
+	 * Constructs the abstract link of v (without points coordinates).
+	 */
+	Abstract_link abstract_link(Vertex_handle v) const {
+		return Abstract_link(*this, Simplex_handle(v));
+	}
+
+	/**
+	 * Constructs the link of 'simplex' with points coordinates.
+	 */
+	Geometric_link abstract_link(const Simplex_handle& simplex) const {
+		return Abstract_link(*this, simplex);
+	}
+
+	/**
+	 * Constructs the link of 'simplex' with points coordinates.
+	 */
+	Geometric_link abstract_link(Edge_handle edge) const {
+		return Abstract_link(*this, edge);
+	}
+
+
+
 		private:
 
 	void add_points_to_link(Geometric_link& link) const {
@@ -193,6 +219,9 @@ class Skeleton_blocker_geometric_complex :
 			link.point(v) = (*this).point(v_root);
 		}
 	}
+
+
+
 };
 
 } // namespace skbl
diff --git a/src/Skeleton_blocker/test/TestGeometricComplex.cpp b/src/Skeleton_blocker/test/TestGeometricComplex.cpp
index 1a81b7f7..bd7af89b 100644
--- a/src/Skeleton_blocker/test/TestGeometricComplex.cpp
+++ b/src/Skeleton_blocker/test/TestGeometricComplex.cpp
@@ -40,7 +40,7 @@ struct Geometry_trait{
 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 Complex::Vertex_handle Vertex_handle;
 
 
 bool test_constructor1(){
@@ -71,12 +71,47 @@ bool test_constructor1(){
 	return (complex==same);
 }
 
+bool test_constructor2(){
+	Complex complex;
+	Skeleton_blocker_off_reader<Complex> off_reader("test2.off",complex);
+	if(!off_reader.is_valid()){
+		std::cerr << "Unable to read file"<<std::endl;
+		return false;
+	}
+	std::cout << "complex has "<<
+			complex.num_vertices()<<" vertices, "<<
+			complex.num_blockers()<<" blockers, "<<
+			complex.num_edges()<<" edges and" <<
+			complex.num_triangles()<<" triangles.";
+
+	if(complex.num_vertices()!=7 || complex.num_edges()!=12 || complex.num_triangles() !=6)
+		return false;
+
+	auto link_0 = complex.abstract_link(Vertex_handle(0));
+
+
+	std::cout<<"\n link(0):"<<link_0.to_string()<<endl;
+
+	auto link_geometric_0 = complex.link(Vertex_handle(0));
+
+	auto print_point = [&](Vertex_handle v){for(auto x : link_geometric_0.point(v)) std::cout <<x<<" "; std::cout<<std::endl;};
+
+	std::for_each(link_geometric_0.vertex_range().begin(),link_geometric_0.vertex_range().end(),print_point);
+
+//	for(auto v : link_geometric_0.vertex_range())
+//		std::cout<<"point("<<v<<"):"<<link_geometric_0.point(v)<<std::endl;
+
+	return link_0.num_vertices()==2;
+}
+
+
 
 
 int main (int argc, char *argv[])
 {
 	Tests tests_geometric_complex;
 	tests_geometric_complex.add("Test constructor 1",test_constructor1);
+	tests_geometric_complex.add("Test constructor 2",test_constructor2);
 
 	if(tests_geometric_complex.run())
 		return EXIT_SUCCESS;
diff --git a/src/Contraction/example/Garland_heckbert/Point.h b/src/common/include/gudhi/Point.h
index eb250347..4023445b 100644
--- a/src/Contraction/example/Garland_heckbert/Point.h
+++ b/src/common/include/gudhi/Point.h
@@ -32,14 +32,14 @@
 #include <vector>
 #include <cassert>
 #include <cstddef>
+#include <initializer_list>
 
 class Point_d{
 public:
-	Point_d():coords_(3,0){}
-	Point_d(size_t dim):coords_(dim,0){}
+	Point_d(size_t dim=3):coords_(dim,0){}
 	Point_d(const Point_d& other):coords_(other.coords_){}
-	Point_d(double x,double y,double z):coords_({x,y,z}){}
-
+	Point_d(const std::initializer_list<double>& list):coords_(list) {
+	}
 	template<typename CoordsIt>
 	Point_d(CoordsIt begin,CoordsIt end):coords_(begin,end){}
 
@@ -93,11 +93,11 @@ public:
 		return res;
 	}
 
-	double squared_dist(const Point_d& other) const{
-		assert(dimension()==other.dimension());
+	friend double squared_dist(const Point_d& p1,const Point_d& p2){
+		assert(p1.dimension()==p2.dimension());
 		double res = 0;
-		for(unsigned i = 0; i < coords_.size(); ++i)
-			res+= coords_[i]*coords_[i] + other[i]*other[i];
+		for(unsigned i = 0; i < p1.coords_.size(); ++i)
+			res+= (p1[i]-p2[i])*(p1[i]-p2[i]);
 		return res;
 	}