skbl correction bug constructor + contraction add garland heckbert example

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

Former-commit-id: 1789e903625df3c4e7c689fa4888bebd86e616eb

5 files changed, 573 insertions, 0 deletions

diff --git a/src/Contraction/example/CMakeLists.txt b/src/Contraction/example/CMakeLists.txt
index ba411bd4..6d65de73 100644
--- a/src/Contraction/example/CMakeLists.txt
+++ b/src/Contraction/example/CMakeLists.txt
@@ -1,7 +1,10 @@
 cmake_minimum_required(VERSION 2.6)
 project(GUDHIskbl)
 
+
 add_executable(RipsContraction Rips_contraction.cpp)
+add_executable(GarlandHeckbert Garland_heckbert.cpp)
+
 target_link_libraries(RipsContraction ${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)
diff --git a/src/Contraction/example/Garland_heckbert.cpp b/src/Contraction/example/Garland_heckbert.cpp
new file mode 100644
index 00000000..912199a4
--- /dev/null
+++ b/src/Contraction/example/Garland_heckbert.cpp
@@ -0,0 +1,115 @@
+/*
+ * Garland_heckbert.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_H_
+#define GARLAND_HECKBERT_H_
+
+
+#include <iostream>
+#include "Garland_heckbert/Point.h"
+#include "Garland_heckbert/Error_quadric.h"
+#include "Garland_heckbert/Garland_heckbert_policies.h"
+
+#include "gudhi/Edge_contraction.h"
+#include "gudhi/Skeleton_blocker.h"
+#include "gudhi/Off_reader.h"
+
+using namespace std;
+using namespace Gudhi;
+using namespace skbl;
+using namespace contraction;
+
+
+struct Geometry_trait{
+	typedef Point_d Point;
+};
+
+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 Skeleton_blocker_contractor<Complex> Complex_contractor;
+
+
+
+
+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;
+	}
+
+	Complex complex;
+
+	// load the points
+	Skeleton_blocker_off_reader<Complex> off_reader(argv[1],complex);
+	if(!off_reader.is_valid()){
+		std::cerr << "Unable to read file:"<<argv[1]<<std::endl;
+		return EXIT_FAILURE;
+	}
+
+	int num_contractions = atoi(argv[2]);
+
+	std::cout << "Initial complex has "<<
+			complex.num_vertices()<<" vertices, "<<
+			complex.num_blockers()<<" blockers, "<<
+			complex.num_edges()<<" edges and" <<
+			complex.num_triangles()<<" triangles.";
+
+	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)
+	);
+
+	std::cout<<"Contract "<<num_contractions<<" edges\n";
+	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);
+
+
+	return EXIT_SUCCESS;
+}
+
+
+
+#endif /* GARLAND_HECKBERT_H_ */
diff --git a/src/Contraction/example/Garland_heckbert/Error_quadric.h b/src/Contraction/example/Garland_heckbert/Error_quadric.h
new file mode 100755
index 00000000..7a30c35a
--- /dev/null
+++ b/src/Contraction/example/Garland_heckbert/Error_quadric.h
@@ -0,0 +1,164 @@
+/*

+ * Error_quadric.h

+ *

+ *  Created on: 24 janv. 2014

+ *      Author: dsalinas

+ */

+

+#ifndef ERROR_QUADRIC_H_

+#define ERROR_QUADRIC_H_

+

+#include <vector>

+#include <utility>

+#include <boost/optional/optional.hpp>

+

+

+template <typename Point> class Error_quadric{

+private :

+	double coeff[10];

+

+public :

+	Error_quadric(){

+		clear();

+	}

+

+	/**

+	 * Quadric corresponding to the L2 distance to the plane.

+	 *

+	 * According to the notation of Garland Heckbert, they

+	 * denote a quadric symetric matrix as :

+	 * Q = [ q11 q12 q13 q14]

+	 *     [ q12 q22 q23 q24]

+	 *     [ q13 q23 q33 q34]

+	 *     [ q14 q24 q34 q44]

+	 *

+	 * which is represented by a vector with 10 elts that

+	 * are denoted ci for clarity with :

+	 * Q = [ c0  c1  c2  c3 ]

+	 *     [ c1  c4  c5  c6 ]

+	 *     [ c2  c5  c7  c8 ]

+	 *     [ c3  c6  c8  c9 ]

+	 *

+	 * The constructor return the quadrics that represents

+	 * the squared distance to the plane defined by triangle p0,p1,p2

+	 * times the area of triangle p0,p1,p2.

+	 */

+	Error_quadric(const Point & p0,const Point & p1,const Point & p2){

+

+		Point normal(unit_normal(p0,p1,p2));

+		double a=normal[0];

+		double b=normal[1];

+		double c=normal[2];

+		double d= -a*p0[0]-b*p0[1]-c*p0[2];

+		coeff[0] = a*a ;

+		coeff[1] = a*b ;

+		coeff[2] = a*c ;

+		coeff[3] = a*d ;

+		coeff[4] = b*b ;

+		coeff[5] = b*c ;

+		coeff[6] = b*d ;

+		coeff[7] = c*c ;

+		coeff[8] = c*d ;

+		coeff[9] = d*d ;

+

+		double area_p0p1p2 = std::sqrt(squared_area(p0,p1,p2));

+		for(auto& x : coeff)

+			x*= area_p0p1p2;

+	}

+

+

+	inline double squared_area(const Point& p0,const Point& p1,const Point& p2) {

+		//if (x1,x2,x3) = p1-p0 and (y1,y2,y3) = p2-p0

+		//then the squared area is = (u^2+v^2+w^2)/4

+		//with:		u = x2 * y3 - x3 * y2;

+		//		    v = x3 * y1 - x1 * y3;

+		//	     	w = x1 * y2 - x2 * y1;

+		Point p0p1(p1-p0);

+		Point p0p2(p2-p0);

+		double A = p0p1[1] * p0p2[2] - p0p1[2] * p0p2[1];

+		double B = p0p1[2] * p0p2[0] - p0p1[0] * p0p2[2];

+		double C = p0p1[0] * p0p2[1] - p0p1[1] * p0p2[0];

+		return 1./4. * (A*A+B*B+C*C);

+	}

+

+

+	void clear(){

+		for(auto& x:coeff)

+			x=0;

+	}

+

+	Error_quadric& operator+=(const Error_quadric& other){

+		if(this!=&other)

+			for(int i = 0 ; i < 10; ++i)

+				coeff[i] += other.coeff[i];

+		return *this;

+	}

+

+	/**

+	 * @return The quadric quost defined by the scalar product v^T Q v where Q is the quadratic form of Garland/Heckbert

+	 */

+	inline double cost(const Point& point) const{

+		double cost =

+				coeff[0]*point.x()*point.x()+coeff[4]*point.y()*point.y()+coeff[7]*point.z()*point.z()

+				+2*(coeff[1]*point.x()*point.y()+coeff[5]*point.y()*point.z()+coeff[2]*point.z()*point.x())

+				+2*(coeff[3]*point.x()+coeff[6]*point.y()+coeff[8]*point.z())

+				+coeff[9];

+		if(cost<0) return 0;

+		else {

+			return cost;

+		}

+	}

+

+	inline double grad_determinant() const{

+		return

+				coeff[0] * coeff[4] * coeff[7]

+				                            -  coeff[0] * coeff[5] * coeff[5]

+				                                                           -  coeff[1] * coeff[1] * coeff[7]

+				                                                                                          +2*coeff[1] * coeff[5] * coeff[2]

+				                                                                                                                         -  coeff[4] * coeff[2] * coeff[2];

+	}

+

+	/**

+	 * Return the point such that it minimizes the gradient of the quadric.

+	 * 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(

+				(-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

+		);

+	}

+

+

+	/**

+	 * returns the point that minimizes the quadric.

+	 * It inverses the quadric if its determinant is higher that a given threshold .

+	 * If the determinant is lower than this value the returned value is uninitialized.

+	 */

+	boost::optional<Point> min_cost(double scale=1) const{

+		//		const double min_determinant = 1e-4 * scale*scale;

+		const double min_determinant = 1e-5;

+		boost::optional<Point> pt_res;

+		double det = grad_determinant();

+		if (std::abs(det)>min_determinant)

+			pt_res = solve_linear_gradient(det);

+		return pt_res;

+	}

+

+	friend std::ostream& operator<< (std::ostream& stream, const Error_quadric& quadric) {

+		stream << "\n[ "<<quadric.coeff[0]<<","<<quadric.coeff[1]<<","<<quadric.coeff[2]<<","<<quadric.coeff[3]<<";\n";

+		stream <<    " "<<quadric.coeff[1]<<","<<quadric.coeff[4]<<","<<quadric.coeff[5]<<","<<quadric.coeff[6]<<";\n";

+		stream <<    " "<<quadric.coeff[2]<<","<<quadric.coeff[5]<<","<<quadric.coeff[7]<<","<<quadric.coeff[8]<<";\n";

+		stream <<    " "<<quadric.coeff[3]<<","<<quadric.coeff[6]<<","<<quadric.coeff[8]<<","<<quadric.coeff[9]<<"]";

+		return stream;

+	}

+

+

+};

+

+

+

+

+#endif /* ERROR_QUADRIC_H_ */

+

diff --git a/src/Contraction/example/Garland_heckbert/Garland_heckbert_policies.h b/src/Contraction/example/Garland_heckbert/Garland_heckbert_policies.h
new file mode 100644
index 00000000..c792b5b3
--- /dev/null
+++ b/src/Contraction/example/Garland_heckbert/Garland_heckbert_policies.h
@@ -0,0 +1,115 @@
+/*
+ * 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
new file mode 100644
index 00000000..eb250347
--- /dev/null
+++ b/src/Contraction/example/Garland_heckbert/Point.h
@@ -0,0 +1,176 @@
+/*
+ * 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_ */