Forgot to add this part of code to svn previously.

.gitattributes

@@ -3776,6 +3776,8 @@ Surface_mesher/examples/Surface_mesher/skip_vcproj_auto_generation -text
 Surface_mesher/find_debug_macros -text
 Surface_mesher/find_non_debug_macros -text
 Surface_mesher/test/Surface_mesher/combined_spheres.cin -text
+Surface_modeling/examples/Surface_modeling/Deform_mesh.cpp -text
+Surface_modeling/examples/Surface_modeling/Deform_mesh.h -text
 Surface_modeling/examples/Surface_modeling/example.cpp -text
 Surface_modeling/examples/Surface_modeling/k-ring_BGL.cpp -text
 Surface_reconstruction_points_3/demo/Surface_reconstruction_points_3/GlSplat/cmake/FindGLEW.cmake -text

Surface_modeling/examples/Surface_modeling/Deform_mesh.cpp

@@ -0,0 +1,134 @@
+#include "Deform_mesh.h"
+
+
+Deform_mesh::Deform_mesh(Polyhedron &P)
+{
+	polyhedron = P;
+}
+
+Deform_mesh::~Deform_mesh(void)
+{
+}
+
+void Deform_mesh::region_of_interest(Vertex_iterator begin, Vertex_iterator end, size_t k)
+{
+	roi.clear();
+	for (Vertex_iterator vit = begin; vit != end; vit ++)
+	{
+		Vertex_handle handle = vit;
+		roi.push_back(handle);
+	}
+	roi.push_back(end);
+
+	int idx_lv = 0;    // pointing the neighboring vertices on current level
+	int idx_lv_end;
+
+	for (size_t lv = 0; lv < k; lv++)
+	{
+		idx_lv_end = roi.size(); 
+		for ( ;idx_lv < idx_lv_end; idx_lv++ )
+		{
+			Vertex_handle vh = roi[idx_lv];
+			HV_circulator wc = vh->vertex_begin(), done(wc);
+			do {
+				Vertex_handle wh = wc->opposite()->vertex();
+				vector<Vertex_handle> ::iterator result = find(roi.begin(), roi.end(), wh);
+				if (result == roi.end())
+				{
+					roi.push_back(wh);
+				}
+				++wc;
+			}while(wc != done);
+		}
+	}
+}
+
+void Deform_mesh::handles(Vertex_iterator begin, Vertex_iterator end)
+{
+	hdl.clear();
+	for (Vertex_iterator vit = begin; vit != end; vit ++)
+	{
+		Vertex_handle handle = vit;
+		hdl.push_back(handle);
+	}
+	hdl.push_back(end);
+}
+
+void Deform_mesh::preprocess()
+{
+	CGAL_TRACE("Calls preprocess()\n");
+
+	double time_init = clock();
+
+	double duration_assembly = 0.0;
+	double duration_prefactor = 0.0;
+
+	// get #variables
+	unsigned int nb_variables = polyhedron.size_of_vertices();
+
+	CGAL_TRACE("  Creates matrix...\n");
+	// Assemble linear system A*X=B
+	Taucs_solver_traits<double>::Matrix A(nb_variables); // matrix is symmetric definite positive
+	Taucs_solver_traits<double>::Vector X(nb_variables), B(nb_variables);
+
+	assemble_laplacian(A, "uni");
+
+	duration_assembly = (clock() - time_init)/CLOCKS_PER_SEC;
+	CGAL_TRACE("  Creates matrix: done (%.2lf s)\n", duration_assembly);
+	CGAL_TRACE("  Pre-factorizing linear system...\n");
+
+	// Pre-factorizing the linear system A*X=B
+	time_init = clock();
+	double D;
+	if(!solver.pre_factor(A, D))
+		return;
+	duration_prefactor = (clock() - time_init)/CLOCKS_PER_SEC;
+
+	CGAL_TRACE("  Pre-factorizing linear system: done (%.2lf s)\n", duration_prefactor);
+
+}
+
+void Deform_mesh::assemble_laplacian(Taucs_solver_traits<double>::Matrix& A, string type)
+{
+	map<Vertex_handle, int> idx;
+	int index = 0;
+	for (Vertex_iterator vit = polyhedron.vertices_begin(); vit != polyhedron.vertices_end(); vit++)
+	{
+		idx[vit] = index;
+		index++;
+	}
+
+	for (Vertex_iterator vit = polyhedron.vertices_begin(); vit != polyhedron.vertices_end(); vit++)
+	{
+		Vertex_handle vi = vit;
+		double diagonal = 0;
+		int idx_i = idx[vi];
+		HV_circulator wc = vi->vertex_begin(), done(wc);
+		do {
+			Vertex_handle vj = wc->opposite()->vertex();
+			double wij = 1;
+			if (type == "cot")   // cotangent Laplacian weights
+			{
+				;
+			}
+			int idx_j = idx[vj];
+			A.set_coef(idx_i, idx_j, -wij, true);	// off-diagonal coefficient
+			diagonal += wij;
+			++wc;
+		}while(wc != done);
+		// diagonal coefficient
+		A.set_coef(idx_i, idx_i, diagonal, true);
+	}
+
+	// handle constraints
+	for (int i = 0; i < hdl.size(); i++)
+	{
+		int idx_i = idx[hdl[i]];
+		A.set_coef(idx_i, idx_i, 1.0);
+	}
+}
+
+void Deform_mesh::operator ()(Point p, Vector v)
+{
+	p = p - v;
+}

Surface_modeling/examples/Surface_modeling/Deform_mesh.h

@@ -0,0 +1,57 @@
+#pragma once
+
+#include <CGAL/trace.h>
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/Polyhedron_incremental_builder_3.h>
+#include <CGAL/Polyhedron_3.h>
+#include <CGAL/IO/Polyhedron_iostream.h>
+#include <CGAL/Taucs_solver_traits.h>
+#include <iostream>
+#include <string>
+#include <fstream>
+#include <ctime>
+
+typedef CGAL::Simple_cartesian<double>      Kernel;
+typedef Kernel::Vector_3                    Vector;
+typedef Kernel::Point_3                     Point;
+typedef CGAL::Polyhedron_3<Kernel>          Polyhedron;
+
+typedef Polyhedron::Vertex_const_handle                      Vertex_handle;
+typedef Polyhedron::Vertex_iterator							 Vertex_iterator;
+typedef Polyhedron::Halfedge_around_vertex_const_circulator  HV_circulator;
+
+using namespace std;
+using namespace CGAL;
+
+class Deform_mesh
+{
+private:
+	Polyhedron polyhedron;                // target mesh
+	vector<Vertex_handle> roi;
+	vector<Vertex_handle> hdl;
+	vector<Vertex_handle> dsplc;         // displacement of handles
+
+	Taucs_solver_traits<double> solver;
+
+public:
+	// The constructor gets the Polyhedron that we will model
+	Deform_mesh(Polyhedron &P);
+
+	// Release ressources
+	~Deform_mesh(void);
+
+	// The region of interest and the handles are a set of vertices
+	void region_of_interest(Vertex_iterator begin, Vertex_iterator end, size_t k);
+
+	void handles(Vertex_iterator begin, Vertex_iterator end);
+
+	// Before we can model we have to do some precomputation
+	void preprocess();
+
+	// Assemble Laplacian matrix A of linear system A*X=B 
+	void assemble_laplacian(Taucs_solver_traits<double>::Matrix& A, string type);
+
+	// The operator will be called in a real time loop from the GUI.
+	void operator()(Point p, Vector v);
+
+};