Initial commit + started libigl-styling ShapeUp

Former-commit-id: d315c8e0466bba10f1b89e60fd0a2d30827e0673

include/igl/DiscreteShellsTraits.h

@@ -0,0 +1,315 @@
+// This file is part of libhedra, a library for polyhedral mesh processing
+//
+// Copyright (C) 2016 Amir Vaxman <avaxman@gmail.com>
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef HEDRA_DISCRETE_SHELLS_TRAITS_H
+#define HEDRA_DISCRETE_SHELLS_TRAITS_H
+#include <igl/igl_inline.h>
+#include <igl/harmonic.h>
+#include <Eigen/Core>
+#include <string>
+#include <vector>
+#include <cstdio>
+#include <set>
+
+
+namespace hedra { namespace optimization {
+        
+        //this class is a traits class for optimization of discrete shells deformation by given positional constraints. It is an implementation on [Froehlich and Botsch 2012] for general polyhedral meshes, using a triangulation of them
+    
+        //the solution vector is assumed to be arranged as xyzxyzxyz... where each triplet is a coordinate of the free vertices.
+    
+        class DiscreteShellsTraits{
+        public:
+            
+            //Requisites of the traits class
+            Eigen::VectorXi JRows, JCols;  //rows and column indices for the jacobian matrix
+            Eigen::VectorXd JVals;         //values for the jacobian matrix.
+            int xSize;                  //size of the solution
+            Eigen::VectorXd EVec;          //energy vector
+            
+            //These are for the the full Jacobian matrix without removing the handles
+            Eigen::VectorXi fullJRows, fullJCols;
+            Eigen::VectorXd fullJVals;
+            Eigen::MatrixXi flapVertexIndices;  //vertices (i,j,k,l) of a flap on edge e=(k,i) where the triangles are f=(i,j,k) and g=(i,k,l)
+            Eigen::MatrixXi EV;
+            Eigen::VectorXi h;              //list of handles
+            Eigen::MatrixXd qh;             //#h by 3 positions
+            Eigen::VectorXi a2x;            //from the entire set of variables "a" to the free variables in the optimization "x".
+            Eigen::VectorXi colMap;         //raw map of a2x into the columns from fullJCols into JCols
+            Eigen::VectorXd origLengths;    //the original edge lengths.
+            Eigen::VectorXd origDihedrals;  //Original dihedral angles
+            Eigen::MatrixXd VOrig;          //original positions
+            Eigen::MatrixXi T;              //triangles
+            Eigen::VectorXd Wd, Wl;         //geometric weights for the lengths and the dihedral angles.
+            double bendCoeff, lengthCoeff;  //coefficients for the different terms
+            
+            //Eigen::VectorXd x0;                 //the initial solution to the optimization
+            Eigen::MatrixXd fullSolution;       //The final solution of the last optimization
+            
+            void init(const Eigen::MatrixXd& _VOrig,
+                      const Eigen::MatrixXi& _T,
+                      const Eigen::VectorXi& _h,
+                      const Eigen::MatrixXi& _EV,
+                      const Eigen::MatrixXi& ET,
+                      const Eigen::MatrixXi& ETi,
+                      const Eigen::VectorXi& innerEdges){
+                
+                using namespace std;
+                using namespace Eigen;
+                
+                std::set<pair<int, int> > edgeIndicesList;
+                
+                VOrig=_VOrig;
+                T=_T;
+                h=_h;
+                EV=_EV;
+                lengthCoeff=10.0;
+                bendCoeff=1.0;
+                
+                //lengths of edges and diagonals
+                flapVertexIndices.resize(innerEdges.size(),4);
+                
+                //dihedral angles
+                for (int i=0;i<innerEdges.size();i++){
+                    int f=ET(innerEdges(i),0);
+                    int g=ET(innerEdges(i),1);
+                 
+                    int vi=EV(innerEdges(i), 1);
+                    int vj=T(f,(ETi(innerEdges(i),0)+2)%3);
+                    int vk=EV(innerEdges(i),0);
+                    int vl=T(g,(ETi(innerEdges(i),1)+2)%3);
+                    flapVertexIndices.row(i)<<vi,vj,vk,vl;
+                 }
+                
+                //across edges
+                EVec.resize(EV.rows()+flapVertexIndices.rows());
+                origLengths.resize(EV.rows());
+                origDihedrals.resize(flapVertexIndices.rows());
+                Wl.resize(EV.rows());
+                Wd.resize(flapVertexIndices.rows());
+                
+                
+                for (int i=0;i<EV.rows();i++){
+                    origLengths(i)=(VOrig.row(EV(i,1))-VOrig.row(EV(i,0))).norm();
+                    Wl(i)=1.0/origLengths(i);
+                }
+                
+                for (int i=0;i<flapVertexIndices.rows();i++){
+                    RowVector3d eji=VOrig.row(flapVertexIndices(i,0))-VOrig.row(flapVertexIndices(i,1));
+                    RowVector3d ejk=VOrig.row(flapVertexIndices(i,2))-VOrig.row(flapVertexIndices(i,1));
+                    RowVector3d eli=VOrig.row(flapVertexIndices(i,0))-VOrig.row(flapVertexIndices(i,3));
+                    RowVector3d elk=VOrig.row(flapVertexIndices(i,2))-VOrig.row(flapVertexIndices(i,3));
+                    RowVector3d eki=VOrig.row(flapVertexIndices(i,0))-VOrig.row(flapVertexIndices(i,2));
+                    
+                    RowVector3d n1 = (ejk.cross(eji));
+                    RowVector3d n2 = (eli.cross(elk));
+                    double sign=((n1.cross(n2)).dot(eki) >= 0 ? 1.0 : -1.0);
+                    double sinHalf=sign*sqrt((1.0-n1.normalized().dot(n2.normalized()))/2.0);
+                    origDihedrals(i)=2.0*asin(sinHalf);
+                    double areaSum=(n1.norm()+n2.norm())/2.0;
+                    Wd(i)=eki.norm()/sqrt(areaSum);
+                }
+                
+                
+                //creating the Jacobian pattern
+                xSize=3*(VOrig.rows()-h.size());
+                fullJRows.resize(6*EV.rows()+12*flapVertexIndices.rows());
+                fullJCols.resize(6*EV.rows()+12*flapVertexIndices.rows());
+                fullJVals.resize(6*EV.rows()+12*flapVertexIndices.rows());
+                
+                
+                //Jacobian indices for edge lengths
+                for (int i=0;i<EV.rows();i++){
+                    fullJRows.segment(6*i,6).setConstant(i);
+                    fullJCols.segment(6*i,3)<<3*EV(i,0),3*EV(i,0)+1,3*EV(i,0)+2;
+                    fullJCols.segment(6*i+3,3)<<3*EV(i,1),3*EV(i,1)+1,3*EV(i,1)+2;
+                }
+                
+                //Jacobian indices for dihedral angles
+                for (int i=0;i<flapVertexIndices.rows();i++){
+                    fullJRows.segment(6*EV.rows()+12*i,12).setConstant(EV.rows()+i);
+                    for (int k=0;k<4;k++)
+                        fullJCols.segment(6*EV.rows()+12*i+3*k,3)<<3*flapVertexIndices(i,k),3*flapVertexIndices(i,k)+1,3*flapVertexIndices(i,k)+2;
+                }
+                
+                
+                a2x=Eigen::VectorXi::Zero(VOrig.rows());
+                int CurrIndex=0;
+                for (int i=0;i<h.size();i++)
+                    a2x(h(i))=-1;
+                
+                for (int i=0;i<VOrig.rows();i++)
+                    if (a2x(i)!=-1)
+                        a2x(i)=CurrIndex++;
+                
+                colMap.resize(3*VOrig.rows());
+                for (int i=0;i<VOrig.rows();i++)
+                    if (a2x(i)!=-1)
+                        colMap.segment(3*i,3)<<3*a2x(i),3*a2x(i)+1,3*a2x(i)+2;
+                    else
+                        colMap.segment(3*i,3)<<-1,-1,-1;
+                
+                //setting up the Jacobian rows and columns
+                int actualGradCounter=0;
+                for (int i=0;i<fullJCols.size();i++){
+                    if (colMap(fullJCols(i))!=-1)  //not a removed variable
+                        actualGradCounter++;
+                }
+                
+                JRows.resize(actualGradCounter);
+                JCols.resize(actualGradCounter);
+                JVals.resize(actualGradCounter);
+                
+                
+                actualGradCounter=0;
+                for (int i=0;i<fullJCols.size();i++){
+                    if (colMap(fullJCols(i))!=-1){  //not a removed variable
+                        JRows(actualGradCounter)=fullJRows(i);
+                        JCols(actualGradCounter++)=colMap(fullJCols(i));
+                    }
+                }
+            }
+            
+            //provide the initial solution to the solver
+            void initial_solution(Eigen::VectorXd& x0){
+                //using biharmonic deformation fields
+                Eigen::MatrixXd origHandlePoses(qh.rows(),3);
+                for (int i=0;i<qh.rows();i++)
+                    origHandlePoses.row(i)=VOrig.row(h(i));
+                Eigen::MatrixXd D;
+                Eigen::MatrixXd D_bc = qh - origHandlePoses;
+                igl::harmonic(VOrig,T,h,D_bc,2,D);
+                Eigen::MatrixXd V0 = VOrig+D;
+                for (int i=0;i<VOrig.rows();i++)
+                    if (a2x(i)!=-1)
+                        x0.segment(3*a2x(i),3)<<V0.row(i).transpose();
+                
+            }
+            
+            void pre_iteration(const Eigen::VectorXd& prevx){}
+            bool post_iteration(const Eigen::VectorXd& x){return false;  /*never stop after an iteration*/}
+            
+            
+            //updating the energy vector for a given current solution
+            void update_energy(const Eigen::VectorXd& x){
+                
+                using namespace std;
+                using namespace Eigen;
+                
+                MatrixXd fullx(xSize+h.size(),3);
+                for (int i=0;i<a2x.size();i++)
+                    if (a2x(i)!=-1)
+                        fullx.row(i)<<x.segment(3*a2x(i),3).transpose();
+                
+                for (int i=0;i<h.size();i++)
+                    fullx.row(h(i))=qh.row(i);
+                
+                fullJVals.setZero();
+                for (int i=0;i<EV.rows();i++)
+                    EVec(i)=((fullx.row(EV(i,1))-fullx.row(EV(i,0))).norm()-origLengths(i))*Wl(i)*lengthCoeff;
+                
+                
+                for (int i=0;i<flapVertexIndices.rows();i++){
+                    RowVector3d eji=fullx.row(flapVertexIndices(i,0))-fullx.row(flapVertexIndices(i,1));
+                    RowVector3d ejk=fullx.row(flapVertexIndices(i,2))-fullx.row(flapVertexIndices(i,1));
+                    RowVector3d eli=fullx.row(flapVertexIndices(i,0))-fullx.row(flapVertexIndices(i,3));
+                    RowVector3d elk=fullx.row(flapVertexIndices(i,2))-fullx.row(flapVertexIndices(i,3));
+                    RowVector3d eki=fullx.row(flapVertexIndices(i,0))-fullx.row(flapVertexIndices(i,2));
+                    
+                    RowVector3d n1 = (ejk.cross(eji));
+                    RowVector3d n2 = (eli.cross(elk));
+                    double sign=((n1.cross(n2)).dot(eki) >= 0 ? 1.0 : -1.0);
+                    double dotn1n2=1.0-n1.normalized().dot(n2.normalized());
+                    if (dotn1n2<0.0) dotn1n2=0.0;  //sanitizing
+                    double sinHalf=sign*sqrt(dotn1n2/2.0);
+                    if (sinHalf>1.0) sinHalf=1.0; if (sinHalf<-1.0) sinHalf=-1.0;
+                    EVec(EV.rows()+i)=(2.0*asin(sinHalf)-origDihedrals(i))*Wd(i)*bendCoeff;
+                }
+                
+                for (int i=0;i<EVec.size();i++)
+                    if (isnan(EVec(i)))
+                        cout<<"nan in EVec("<<i<<")"<<endl;
+            }
+            
+            
+            //update the jacobian values for a given current solution
+            void update_jacobian(const Eigen::VectorXd& x){
+                using namespace std;
+                using namespace Eigen;
+                
+                MatrixXd fullx(xSize+h.size(),3);
+                for (int i=0;i<a2x.size();i++)
+                    if (a2x(i)!=-1)
+                        fullx.row(i)<<x.segment(3*a2x(i),3).transpose();
+                
+                for (int i=0;i<h.size();i++)
+                    fullx.row(h(i))=qh.row(i);
+                
+                fullJVals.setZero();
+                for (int i=0;i<EV.rows();i++){
+                   
+                    RowVector3d normedEdgeVector=(fullx.row(EV(i,1))-fullx.row(EV(i,0))).normalized();
+                    fullJVals.segment(6*i,3)<<-normedEdgeVector.transpose()*Wl(i)*lengthCoeff;
+                    fullJVals.segment(6*i+3,3)<<normedEdgeVector.transpose()*Wl(i)*lengthCoeff;
+                }
+                
+                for (int i=0;i<flapVertexIndices.rows();i++){
+                    RowVector3d eji=fullx.row(flapVertexIndices(i,0))-fullx.row(flapVertexIndices(i,1));
+                    RowVector3d ejk=fullx.row(flapVertexIndices(i,2))-fullx.row(flapVertexIndices(i,1));
+                    RowVector3d eli=fullx.row(flapVertexIndices(i,0))-fullx.row(flapVertexIndices(i,3));
+                    RowVector3d elk=fullx.row(flapVertexIndices(i,2))-fullx.row(flapVertexIndices(i,3));
+                    RowVector3d eki=fullx.row(flapVertexIndices(i,0))-fullx.row(flapVertexIndices(i,2));
+                    
+                    RowVector3d n1 = (ejk.cross(eji));
+                    RowVector3d n2 = (eli.cross(elk));
+                    double sign=((n1.cross(n2)).dot(eki) >= 0 ? 1.0 : -1.0);
+                    double dotn1n2=1.0-n1.normalized().dot(n2.normalized());
+                    if (dotn1n2<0.0) dotn1n2=0.0;  //sanitizing
+                    double sinHalf=sign*sqrt(dotn1n2/2.0);
+                    if (sinHalf>1.0) sinHalf=1.0; if (sinHalf<-1.0) sinHalf=-1.0;
+                    
+                    fullJVals.segment(6*EV.rows()+12*i,3)<<(Wd(i)*((ejk.dot(-eki)/(n1.squaredNorm()*eki.norm()))*n1+(elk.dot(-eki)/(n2.squaredNorm()*eki.norm()))*n2)).transpose()*bendCoeff;
+                    fullJVals.segment(6*EV.rows()+12*i+3,3)<<(Wd(i)*(-eki.norm()/n1.squaredNorm())*n1).transpose()*bendCoeff;
+                    fullJVals.segment(6*EV.rows()+12*i+6,3)<<(Wd(i)*((eji.dot(eki)/(n1.squaredNorm()*eki.norm()))*n1+(eli.dot(eki)/(n2.squaredNorm()*eki.norm()))*n2)).transpose()*bendCoeff;
+                    fullJVals.segment(6*EV.rows()+12*i+9,3)<<(Wd(i)*(-eki.norm()/n2.squaredNorm())*n2).transpose()*bendCoeff;
+                    
+                }
+                
+
+                int actualGradCounter=0;
+                for (int i=0;i<fullJCols.size();i++)
+                    if (colMap(fullJCols(i))!=-1)  //not a removed variable
+                        JVals(actualGradCounter++)=fullJVals(i);
+                
+                for (int i=0;i<JVals.size();i++)
+                    if (isnan(JVals(i)))
+                        cout<<"nan in JVals("<<i<<")"<<endl;
+            }
+
+            
+            
+            bool post_optimization(const Eigen::VectorXd& x){
+                fullSolution.conservativeResize(a2x.size(),3);
+                for (int i=0;i<a2x.size();i++)
+                    if (a2x(i)!=-1)
+                        fullSolution.row(i)<<x.segment(3*a2x(i),3).transpose();
+                
+                for (int i=0;i<h.size();i++)
+                    fullSolution.row(h(i))=qh.row(i);
+                
+                return true;  //stop optimization after this
+            }
+            
+            DiscreteShellsTraits(){}
+            ~DiscreteShellsTraits(){}
+        };
+        
+        
+    } }
+
+
+#endif

include/igl/LMSolver.h

@@ -0,0 +1,350 @@
+// This file is part of libhedra, a library for polyhedral mesh processing
+//
+// Copyright (C) 2016 Amir Vaxman <avaxman@gmail.com>
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef HEDRA_LEVENBERG_MARQUADT_SOLVER_H
+#define HEDRA_LEVENBERG_MARQUADT_SOLVER_H
+#include <igl/igl_inline.h>
+#include <igl/sortrows.h>
+#include <igl/speye.h>
+#include <Eigen/Core>
+#include <string>
+#include <vector>
+#include <cstdio>
+#include <iostream>
+
+namespace hedra {
+    namespace optimization
+    {
+        
+        template<class LinearSolver, class SolverTraits>
+        class LMSolver{
+        public:
+            Eigen::VectorXd x;      //current solution; always updated
+            Eigen::VectorXd prevx;  //the solution of the previous iteration
+            Eigen::VectorXd x0;     //the initial solution to the system
+            Eigen::VectorXd d;             //the direction taken.
+            Eigen::VectorXd currEnergy;    //the current value of the energy
+            Eigen::VectorXd prevEnergy;    //the previous value of the energy
+            
+            Eigen::VectorXi HRows, HCols;  //(row,col) pairs for H=J^T*J matrix
+            Eigen::VectorXd HVals;      //values for H matrix
+            Eigen::MatrixXi S2D;        //single J to J^J indices
+
+            LinearSolver* LS;
+            SolverTraits* ST;
+            int maxIterations;
+            double xTolerance;
+            double fooTolerance;
+            
+            /*void TestMatrixOperations(){
+             
+                using namespace Eigen;
+                using namespace std;
+                int RowSize=1000;
+                int ColSize=1000;
+                int TriSize=4000;
+             
+                VectorXi Rows;
+                VectorXi Cols;
+                VectorXd Vals=VectorXd::Random(TriSize);
+             
+                VectorXd dRows=VectorXd::Random(TriSize);
+                VectorXd dCols=VectorXd::Random(TriSize);
+             
+                cout<<"dRows Range: "<<dRows.minCoeff()<<","<<dRows.maxCoeff()<<endl;
+             
+                Rows=((dRows+MatrixXd::Constant(dRows.size(),1,1.0))*500.0).cast<int>();
+                Cols=((dCols+MatrixXd::Constant(dRows.size(),1,1.0))*500.0).cast<int>();
+                VectorXi SortRows, stub;
+                VectorXi MRows, MCols;
+                VectorXd MVals;
+             
+             igl::sortrows(Rows, true, SortRows, stub);
+             
+             MatrixXi S2D;
+                
+            double miu=15.0;
+             
+             MatrixPattern(SortRows, Cols,  MRows, MCols, S2D);
+             MVals.resize(MRows.size());
+             MatrixValues(MRows, MCols, Vals, S2D, miu, MVals);
+ 
+             //testing multiplyadjoint
+             SparseMatrix<double> M(SortRows.maxCoeff()+1,Cols.maxCoeff()+1);
+             
+             //cout<<"Our I"<<I<<endl;
+             //cout<<"Our J"<<J<<endl;
+             //cout<<"Our S"<<S<<endl;
+             
+             
+             vector<Triplet<double> > Tris;
+             
+             for (int i=0;i<SortRows.size();i++)
+             Tris.push_back(Triplet<double>(SortRows(i), Cols(i), Vals(i)));
+             
+             M.setFromTriplets(Tris.begin(), Tris.end());
+             Tris.clear();
+                
+                
+                SparseMatrix<double> I;
+                igl::speye(M.cols(),M.cols(),I);
+             
+             SparseMatrix<double> MtM1=M.transpose()*M+miu*I;
+             SparseMatrix<double> MtM2(MtM1.rows(), MtM1.cols());
+             for (int i=0;i<MRows.size();i++)
+                 Tris.push_back(Triplet<double>(MRows(i), MCols(i), MVals(i)));
+             
+             MtM2.setFromTriplets(Tris.begin(), Tris.end());
+             
+             bool Discrepancy=false;
+             SparseMatrix<double> DiffMat=MtM1-MtM2;
+             for (int k=0; k<DiffMat.outerSize(); ++k){
+             for (SparseMatrix<double>::InnerIterator it(DiffMat,k); it; ++it){
+             if ((abs(it.value())>10e-6)&&(it.row()<=it.col())){
+             cout<<"Discrepancy at ("<<it.row()<<","<<it.col()<<"): "<<it.value()<<endl;
+             cout<<"MtM Our Values:"<<MtM2.coeffRef(it.row(),it.col())<<endl;
+             cout<<"MtM Evaluated:"<<MtM1.coeffRef(it.row(),it.col())<<endl;
+             Discrepancy=true;
+             }
+             }
+             }
+             if (!Discrepancy) cout<<"Matrices are equal!"<<endl;
+             }*/
+
+            
+            //Input: pattern of matrix M by (iI,iJ) representation
+            //Output: pattern of matrix M^T*M by (oI, oJ) representation
+            //        map between values in the input to values in the output (Single2Double). The map is aggregating values from future iS to oS
+            //prerequisite: iI are sorted by rows (not necessary columns)
+            void MatrixPattern(const Eigen::VectorXi& iI,
+                               const Eigen::VectorXi& iJ,
+                               Eigen::VectorXi& oI,
+                               Eigen::VectorXi& oJ,
+                               Eigen::MatrixXi& S2D)
+            {
+                int CurrTri=0;
+                using namespace Eigen;
+                std::vector<int> oIlist;
+                std::vector<int> oJlist;
+                std::vector<std::pair<int, int> > S2Dlist;
+                do{
+                    int CurrRow=iI(CurrTri);
+                    int NumCurrTris=0;
+                    while ((CurrTri+NumCurrTris<iI.size())&&(iI(CurrTri+NumCurrTris)==CurrRow))
+                        NumCurrTris++;
+                    
+                    for (int i=CurrTri;i<CurrTri+NumCurrTris;i++){
+                        for (int j=CurrTri;j<CurrTri+NumCurrTris;j++){
+                            if (iJ(j)>=iJ(i)){
+                                oIlist.push_back(iJ(i));
+                                oJlist.push_back(iJ(j));
+                                S2Dlist.push_back(std::pair<int,int>(i,j));
+                            }
+                        }
+                    }
+                    CurrTri+=NumCurrTris;
+                }while (CurrTri!=iI.size());
+                
+                
+                //triplets for miu
+                for (int i=0;i<iJ.maxCoeff()+1;i++){
+                    oIlist.push_back(i);
+                    oJlist.push_back(i);
+                }
+                
+                oI.resize(oIlist.size());
+                oJ.resize(oJlist.size());
+                S2D.resize(S2Dlist.size(),2);
+                
+                for (int i=0;i<oIlist.size();i++){
+                    oI(i)=oIlist[i];
+                    oJ(i)=oJlist[i];
+                }
+                for (int i=0;i<S2Dlist.size();i++)
+                    S2D.row(i)<<S2Dlist[i].first, S2Dlist[i].second;
+                
+            }
+            
+            //returns the values of M^T*M+miu*I by multiplication and aggregating from Single2double list.
+            //prerequisite - oS is allocated
+            void MatrixValues(const Eigen::VectorXi& oI,
+                              const Eigen::VectorXi& oJ,
+                              const Eigen::VectorXd& iS,
+                              const Eigen::MatrixXi& S2D,
+                              const double miu,
+                              Eigen::VectorXd& oS)
+            {
+                for (int i=0;i<S2D.rows();i++)
+                    oS(i)=iS(S2D(i,0))*iS(S2D(i,1));
+                
+                //adding miu*I
+                for (int i=S2D.rows();i<oI.rows();i++)
+                    oS(i)=miu;
+                
+            }
+            
+            //returns M^t*ivec by (I,J,S) representation
+            void MultiplyAdjointVector(const Eigen::VectorXi& iI,
+                                       const Eigen::VectorXi& iJ,
+                                       const Eigen::VectorXd& iS,
+                                       const Eigen::VectorXd& iVec,
+                                       Eigen::VectorXd& oVec)
+            {
+                oVec.setZero();
+                for (int i=0;i<iI.size();i++)
+                    oVec(iJ(i))+=iS(i)*iVec(iI(i));
+            }
+            
+            
+        public:
+            
+            LMSolver(){};
+            
+            void init(LinearSolver* _LS,
+                      SolverTraits* _ST,
+                      int _maxIterations=100,
+                      double _xTolerance=10e-9,
+                      double _fooTolerance=10e-9){
+                
+                LS=_LS;
+                ST=_ST;
+                maxIterations=_maxIterations;
+                xTolerance=_xTolerance;
+                fooTolerance=_fooTolerance;
+                //analysing pattern
+                MatrixPattern(ST->JRows, ST->JCols,HRows,HCols,S2D);
+                HVals.resize(HRows.size());
+                
+                LS->analyze(HRows,HCols, true);
+                
+                d.resize(ST->xSize);
+                x.resize(ST->xSize);
+                x0.resize(ST->xSize);
+                prevx.resize(ST->xSize);
+                currEnergy.resize(ST->EVec.size());
+                prevEnergy.resize(ST->EVec.size());
+                
+                //TestMatrixOperations();
+            }
+            
+            
+            bool solve(const bool verbose) {
+                
+                using namespace Eigen;
+                using namespace std;
+                ST->initial_solution(x0);
+                prevx<<x0;
+                int currIter=0;
+                bool stop=false;
+                double currError, prevError;
+                VectorXd rhs(ST->xSize);
+                VectorXd direction;
+                if (verbose)
+                    cout<<"******Beginning Optimization******"<<endl;
+                
+                double tau=10e-3;
+                
+                //estimating initial miu
+                double miu=0.0;
+                ST->update_jacobian(prevx);
+                MatrixValues(HRows, HCols, ST->JVals, S2D, miu, HVals);
+                for (int i=0;i<HRows.size();i++)
+                    if (HRows(i)==HCols(i))  //on the diagonal
+                        miu=(miu < HVals(i) ? HVals(i) : miu);
+                miu*=tau;
+                double initmiu=miu;
+                cout<<"initial miu: "<<miu<<endl;
+                double beta=2.0;
+                double nu=beta;
+                double gamma=3.0;
+                do{
+                    currIter=0;
+                    stop=false;
+                    do{
+                        ST->pre_iteration(prevx);
+                        ST->update_energy(prevx);
+                        ST->update_jacobian(prevx);
+                        if (verbose)
+                            cout<<"Initial Energy for Iteration "<<currIter<<": "<<ST->EVec.template squaredNorm()<<endl;
+                        MatrixValues(HRows, HCols, ST->JVals, S2D,  miu, HVals);
+                        MultiplyAdjointVector(ST->JRows, ST->JCols, ST->JVals, -ST->EVec, rhs);
+                        
+                        double firstOrderOptimality=rhs.template lpNorm<Infinity>();
+                        if (verbose)
+                            cout<<"firstOrderOptimality: "<<firstOrderOptimality<<endl;
+                        
+                        if (firstOrderOptimality<fooTolerance){
+                            x=prevx;
+                            if (verbose){
+                                cout<<"First-order optimality has been reached"<<endl;
+                                break;
+                            }
+                        }
+                        
+                        //solving to get the GN direction
+                        if(!LS->factorize(HVals, true)) {
+                            // decomposition failed
+                            cout<<"Solver Failed to factorize! "<<endl;
+                            return false;
+                        }
+                        
+                        MatrixXd mRhs=rhs;
+                        MatrixXd mDirection;
+                        LS->solve(mRhs,mDirection);
+                        direction=mDirection.col(0);
+                        if (verbose)
+                            cout<<"direction magnitude: "<<direction.norm()<<endl;
+                        if (direction.norm() < xTolerance * prevx.norm()){
+                            x=prevx;
+                            if (verbose)
+                                cout<<"Stopping since direction magnitude small."<<endl;
+                            break;
+                        }
+                        VectorXd tryx=prevx+direction;
+                        ST->update_energy(prevx);
+                        double prevE=ST->EVec.squaredNorm();
+                        ST->update_energy(tryx);
+                        double currE=ST->EVec.squaredNorm();
+                        
+                        double rho=(prevE-currE)/(direction.dot(miu*direction+rhs));
+                        if (rho>0){
+                            x=tryx;
+                            //if (verbose){
+                                //cout<<"Energy: "<<currE<<endl;
+                            //    cout<<"1.0-(beta-1.0)*pow(2.0*rho-1.0,3): "<<1.0-(beta-1.0)*pow(2.0*rho-1.0,3)<<endl;
+                            //}
+                            miu*=(1.0/gamma > 1.0-(beta-1.0)*pow(2.0*rho-1.0,3) ? 1.0/gamma : 1.0-(beta-1.0)*pow(2.0*rho-1.0,3));
+                            nu=beta;
+                            //if (verbose)
+                            //    cout<<"rho, miu, nu: "<<rho<<","<<miu<<","<<nu<<endl;
+                        } else {
+                            x=prevx;
+                            miu = miu*nu;
+                            nu=2*nu;
+                            //if (verbose)
+                            //    cout<<"rho, miu, nu: "<<rho<<","<<miu<<","<<nu<<endl;
+                        }
+                                               
+                        //The SolverTraits can order the optimization to stop by giving "true" of to continue by giving "false"
+                        if (ST->post_iteration(x)){
+                            if (verbose)
+                                cout<<"ST->Post_iteration() gave a stop"<<endl;
+                            break;
+                        }
+                        currIter++;
+                        prevx=x;
+                    }while (currIter<=maxIterations);
+                }while (!ST->post_optimization(x));
+                return true;
+            }
+        };
+        
+    }
+}
+
+
+#endif

include/igl/SLSolver.h

@@ -0,0 +1,184 @@
+// This file is part of libhedra, a library for polyhedral mesh processing
+//
+// Copyright (C) 2016 Amir Vaxman <avaxman@gmail.com>
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef HEDRA_SINGLE_LINEAR_SOLVER_H
+#define HEDRA_SINGLE_LINEAR_SOLVER_H
+#include <igl/igl_inline.h>
+#include <igl/sortrows.h>
+#include <igl/speye.h>
+#include <Eigen/Core>
+#include <string>
+#include <vector>
+#include <cstdio>
+#include <iostream>
+
+namespace hedra {
+    namespace optimization
+    {
+        
+        template<class LinearSolver, class SolverTraits>
+        class SLSolver{
+        public:
+            Eigen::VectorXd x;      //current solution; always updated
+            Eigen::VectorXd prevx;  //the solution of the previous iteration
+            Eigen::VectorXd x0;     //the initial solution to the system
+            Eigen::VectorXd d;             //the direction taken.
+            Eigen::VectorXd currEnergy;    //the current value of the energy
+            Eigen::VectorXd prevEnergy;    //the previous value of the energy
+            
+            Eigen::VectorXi HRows, HCols;  //(row,col) pairs for H=J^T*J matrix
+            Eigen::VectorXd HVals;      //values for H matrix
+            Eigen::MatrixXi S2D;        //single J to J^J indices
+
+            LinearSolver* LS;
+            SolverTraits* ST;
+            
+
+            //Input: pattern of matrix M by (iI,iJ) representation
+            //Output: pattern of matrix M^T*M by (oI, oJ) representation
+            //        map between values in the input to values in the output (Single2Double). The map is aggregating values from future iS to oS
+            //prerequisite: iI are sorted by rows (not necessary columns)
+            void MatrixPattern(const Eigen::VectorXi& iI,
+                               const Eigen::VectorXi& iJ,
+                               Eigen::VectorXi& oI,
+                               Eigen::VectorXi& oJ,
+                               Eigen::MatrixXi& S2D)
+            {
+                int CurrTri=0;
+                using namespace Eigen;
+                std::vector<int> oIlist;
+                std::vector<int> oJlist;
+                std::vector<std::pair<int, int> > S2Dlist;
+                do{
+                    int CurrRow=iI(CurrTri);
+                    int NumCurrTris=0;
+                    while ((CurrTri+NumCurrTris<iI.size())&&(iI(CurrTri+NumCurrTris)==CurrRow))
+                        NumCurrTris++;
+                    
+                    for (int i=CurrTri;i<CurrTri+NumCurrTris;i++){
+                        for (int j=CurrTri;j<CurrTri+NumCurrTris;j++){
+                            if (iJ(j)>=iJ(i)){
+                                oIlist.push_back(iJ(i));
+                                oJlist.push_back(iJ(j));
+                                S2Dlist.push_back(std::pair<int,int>(i,j));
+                            }
+                        }
+                    }
+                    CurrTri+=NumCurrTris;
+                }while (CurrTri!=iI.size());
+            
+                
+                oI.resize(oIlist.size());
+                oJ.resize(oJlist.size());
+                S2D.resize(S2Dlist.size(),2);
+                
+                for (int i=0;i<oIlist.size();i++){
+                    oI(i)=oIlist[i];
+                    oJ(i)=oJlist[i];
+                }
+                for (int i=0;i<S2Dlist.size();i++)
+                    S2D.row(i)<<S2Dlist[i].first, S2Dlist[i].second;
+                
+            }
+            
+            //returns the values of M^T*M+miu*I by multiplication and aggregating from Single2double list.
+            //prerequisite - oS is allocated
+            void MatrixValues(const Eigen::VectorXi& oI,
+                              const Eigen::VectorXi& oJ,
+                              const Eigen::VectorXd& iS,
+                              const Eigen::MatrixXi& S2D,
+                              Eigen::VectorXd& oS)
+            {
+                for (int i=0;i<S2D.rows();i++)
+                    oS(i)=iS(S2D(i,0))*iS(S2D(i,1));
+                
+            }
+            
+            //returns M^t*ivec by (I,J,S) representation
+            void MultiplyAdjointVector(const Eigen::VectorXi& iI,
+                                       const Eigen::VectorXi& iJ,
+                                       const Eigen::VectorXd& iS,
+                                       const Eigen::VectorXd& iVec,
+                                       Eigen::VectorXd& oVec)
+            {
+                oVec.setZero();
+                for (int i=0;i<iI.size();i++)
+                    oVec(iJ(i))+=iS(i)*iVec(iI(i));
+            }
+            
+            
+        public:
+            
+            SLSolver(){};
+            
+            void init(LinearSolver* _LS,
+                      SolverTraits* _ST){
+                
+                LS=_LS;
+                ST=_ST;
+                //analysing pattern
+                MatrixPattern(ST->JRows, ST->JCols,HRows,HCols,S2D);
+                HVals.resize(HRows.size());
+                
+                LS->analyze(HRows,HCols);
+                
+                d.resize(ST->xSize);
+                x.resize(ST->xSize);
+                x0.resize(ST->xSize);
+                prevx.resize(ST->xSize);
+                currEnergy.resize(ST->EVec.size());
+                prevEnergy.resize(ST->EVec.size());
+                
+                //TestMatrixOperations();
+            }
+            
+            
+            bool solve(const bool verbose) {
+                
+                using namespace Eigen;
+                using namespace std;
+                ST->initial_solution(x0);
+                prevx<<x0;
+                int currIter=0;
+                bool stop=false;
+                double currError, prevError;
+                VectorXd rhs(ST->xSize);
+                VectorXd direction;
+                if (verbose)
+                    cout<<"******Beginning Optimization******"<<endl;
+
+                ST->update_jacobian(prevx);
+                do{
+                    ST->update_energy(prevx);
+                    ST->update_jacobian(prevx);
+                    ST->pre_iteration(prevx);
+                    MatrixValues(HRows, HCols, ST->JVals, S2D, HVals);
+                    MultiplyAdjointVector(ST->JRows, ST->JCols, ST->JVals, -ST->EVec, rhs);
+                    
+                    //solving to get the GN direction
+                    if(!LS->factorize(HVals)) {
+                        // decomposition failed
+                        cout<<"Solver Failed to factorize! "<<endl;
+                        return false;
+                    }
+                        
+                    LS->solve(rhs,direction);
+                    x=prevx+direction;
+                    ST->update_energy(x);
+                    ST->update_jacobian(x);
+                    ST->post_iteration(x);
+                    prevx=x;
+                }while (!ST->post_optimization(x));
+                return true;
+            }
+        };
+        
+    }
+}
+
+
+#endif

include/igl/check_traits.h

@@ -0,0 +1,88 @@
+// This file is part of libhedra, a library for polyhedral mesh processing
+//
+// Copyright (C) 2016 Amir Vaxman <avaxman@gmail.com>
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef HEDRA_CHECK_TRAITS_H
+#define HEDRA_CHECK_TRAITS_H
+#include <igl/igl_inline.h>
+#include <Eigen/Core>
+#include <string>
+#include <vector>
+#include <cstdio>
+
+namespace hedra {
+    namespace optimization
+    {
+        //This function checks the Jacobian of a traits class that is put for optimization by approximate finite differences, and reports the difference. It is important to use after coding, but it is not for the actual optimization process, since it is really brute-force and slow.
+        template<class SolverTraits>
+        void check_traits(SolverTraits& Traits){
+            using namespace Eigen;
+            using namespace std;
+            cout<<"WARNING: FE gradient checking, reducing performance!!!"<<endl;
+            cout<<"******************************************************"<<endl;
+            VectorXd CurrSolution(Traits.xSize);
+            CurrSolution<<VectorXd::Random(Traits.xSize, 1);
+            
+            cout<<"Solution Size: "<<Traits.xSize<<endl;
+            
+            Traits.update_energy(CurrSolution);
+            Traits.update_jacobian(CurrSolution);
+         
+            int MaxRow=Traits.JRows.maxCoeff()+1;
+            vector<Triplet<double> > GradTris;
+            
+            for (int i=0;i<Traits.JRows.size();i++)
+                GradTris.push_back(Triplet<double>(Traits.JRows(i), Traits.JCols(i), Traits.JVals(i)));
+            
+            
+            SparseMatrix<double> TraitGradient(MaxRow, CurrSolution.size());
+            TraitGradient.setFromTriplets(GradTris.begin(),GradTris.end());
+            
+            SparseMatrix<double> FEGradient(MaxRow, CurrSolution.size());
+            vector<Triplet<double> > FEGradientTris;
+            for (int i=0;i<CurrSolution.size();i++){
+                VectorXd vh(CurrSolution.size()); vh.setZero(); vh(i)=10e-5;
+                Traits.update_energy(CurrSolution+vh);
+                Traits.update_jacobian(CurrSolution+vh);
+                VectorXd EnergyPlus=Traits.EVec;
+                Traits.update_energy(CurrSolution-vh);
+                Traits.update_jacobian(CurrSolution-vh);
+                VectorXd EnergyMinus=Traits.EVec;
+                VectorXd CurrGradient=(EnergyPlus-EnergyMinus)/(2*10e-5);
+                //cout<<CurrGradient<<endl;
+                for (int j=0;j<CurrGradient.size();j++)
+                    if (abs(CurrGradient(j))>10e-7)
+                        FEGradientTris.push_back(Triplet<double>(j,i,CurrGradient(j)));
+            }
+            
+            FEGradient.setFromTriplets(FEGradientTris.begin(), FEGradientTris.end());
+            SparseMatrix<double> DiffMat=FEGradient-TraitGradient;
+            double maxcoeff=-32767.0;
+            int Maxi,Maxj;
+            for (int k=0; k<DiffMat.outerSize(); ++k)
+                for (SparseMatrix<double>::InnerIterator it(DiffMat,k); it; ++it){
+                    if (maxcoeff<abs(it.value())){
+                        maxcoeff=abs(it.value());
+                        Maxi=it.row();
+                        Maxj=it.col();
+                        
+                    }
+                    if (abs(it.value())>10e-7){
+                        cout<<"Gradient Discrepancy at: ("<<it.row()<<","<<it.col()<<") of "<<it.value()<<endl;
+                        cout<<"Our Value: "<<TraitGradient.coeffRef(it.row(), it.col())<<endl;
+                        cout<<"Computed Value: "<<FEGradient.coeffRef(it.row(), it.col())<<endl;
+                    }
+                }
+            cout<<"Maximum gradient difference: "<<maxcoeff<<endl;
+            cout<<"At Location: ("<<Maxi<<","<<Maxj<<")"<<endl;
+            
+            
+        }
+    }
+}
+
+
+#endif

include/igl/shapeup.cpp

@@ -0,0 +1,144 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2017 Amir Vaxman <avaxman@gmail.com>
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IGL_SHAPEUP_H
+#define IGL_SHAPEUP_H
+
+#include "shapeup.h"
+#include <igl/igl_inline.h>
+#include <igl/setdiff.h>
+#include <igl/cat.h>
+#include <Eigen/Core>
+#include <vector>
+
+
+//This file implements the following algorithm:
+
+//Boaziz et al.
+//Shape-Up: Shaping Discrete Geometry with Projections
+//Computer Graphics Forum (Proc. SGP) 31(5), 2012
+
+namespace igl
+{
+    
+    //This function does the precomputation of the necessary operators for the shape-up projection process.
+    
+
+    IGL_INLINE void shapeup_precomputation(const Eigen::MatrixXd& V,
+                                       const Eigen::VectorXi& D,
+                                       const Eigen::MatrixXi& F,
+                                       const Eigen::VectorXi& SD,
+                                       const Eigen::MatrixXi& S,
+                                       const Eigen::VectorXi& h,
+                                       const Eigen::VectorXd& w,
+                                       const double shapeCoeff,
+                                       const double closeCoeff,
+                                       struct ShapeupData& sudata)
+    {
+        using namespace Eigen;
+        //The integration solve is separable to x,y,z components
+        sudata.V=V; sudata.F=F; sudata.D=D; sudata.SD=SD; sudata.S=S; sudata.h=h; sudata.closeCoeff=closeCoeff; sudata.shapeCoeff=shapeCoeff;
+        sudata.Q.conservativeResize(SD.sum(), V.rows());  //Shape matrix (integration);
+        sudata.C.conservativeResize(h.rows(), V.rows());        //Closeness matrix for handles
+        
+        std::vector<Triplet<double> > QTriplets;
+        int currRow=0;
+        for (int i=0;i<S.rows();i++){
+            double avgCoeff=1.0/(double)SD(i);
+            
+            for (int j=0;j<SD(i);j++){
+                for (int k=0;k<SD(i);k++){
+                    if (j==k)
+                        QTriplets.push_back(Triplet<double>(currRow+j, S(i,k), (1.0-avgCoeff)));
+                    else
+                        QTriplets.push_back(Triplet<double>(currRow+j, S(i,k), (-avgCoeff)));
+                }
+            }
+            currRow+=SD(i);
+        }
+        
+        sudata.Q.setFromTriplets(QTriplets.begin(), QTriplets.end());
+        
+        std::vector<Triplet<double> > CTriplets;
+        for (int i=0;i<h.size();i++)
+            CTriplets.push_back(Triplet<double>(i,h(i), 1.0));
+        
+        sudata.C.setFromTriplets(CTriplets.begin(), CTriplets.end());
+        
+        igl::cat(1, sudata.Q, sudata.C, sudata.A);
+        sudata.At=sudata.A.transpose();  //to save up this expensive computation.
+        
+        //weight matrix
+        vector<Triplet<double> > WTriplets;
+        //std::cout<<"w: "<<w<<std::endl;
+        currRow=0;
+        for (int i=0;i<SD.rows();i++){
+            for (int j=0;j<SD(i);j++)
+                WTriplets.push_back(Triplet<double>(currRow+j,currRow+j,shapeCoeff*w(i)));
+            currRow+=SD(i);
+        }
+        for (int i=0;i<h.size();i++)
+            WTriplets.push_back(Triplet<double>(SD.sum()+i, SD.sum()+i, closeCoeff));
+        sudata.W.resize(SD.sum()+h.size(), SD.sum()+h.size());
+        sudata.W.setFromTriplets(WTriplets.begin(), WTriplets.end());
+        
+        sudata.E=sudata.At*sudata.W*sudata.A;
+        sudata.solver.compute(sudata.E);
+    }
+    
+    
+    
+    IGL_INLINE void shapeup_compute(void (*projection)(int , const hedra::ShapeupData&, const Eigen::MatrixXd& , Eigen::MatrixXd&),
+                                    const Eigen::MatrixXd& vh,
+                                    const struct ShapeupData& sudata,
+                                    Eigen::MatrixXd& currV,
+                                    const int maxIterations=50,
+                                    const double vTolerance=10e-6)
+    {
+        using namespace Eigen;
+        MatrixXd prevV=currV;
+        MatrixXd PV;
+        MatrixXd b(sudata.A.rows(),3);
+        b.block(sudata.Q.rows(), 0, sudata.h.rows(),3)=vh;  //this stays constant throughout the iterations
+        
+        //std::cout<<"vh: "<<vh<<std::endl;
+        //std::cout<<"V(h(0))"<<currV.row(sudata.h(0))<<std::endl;
+        PV.conservativeResize(sudata.SD.rows(), 3*sudata.SD.maxCoeff());
+        for (int i=0;i<maxIterations;i++){
+            //std::cout<<"A*prevV-b before local projection:"<<(sudata.W*(sudata.A*prevV-b)).squaredNorm()<<std::endl;
+            for (int j=0;j<sudata.SD.rows();j++)
+                projection(j, sudata, currV, PV);
+            //constructing the projection part of the right hand side
+            int currRow=0;
+            for (int i=0;i<sudata.S.rows();i++){
+                for (int j=0;j<sudata.SD(i);j++)
+                    b.row(currRow++)=PV.block(i, 3*j, 1,3);
+                //currRow+=sudata.SD(i);
+            }
+            //std::cout<<"A*prevV-b after local projection:"<<(sudata.W*(sudata.A*prevV-b)).squaredNorm()<<std::endl;
+            //std::cout<<"A*currV-b:"<<i<<(sudata.A*currV-b)<<std::endl;
+            currV=sudata.solver.solve(sudata.At*sudata.W*b);
+            //std::cout<<"b: "<<b<<std::endl;
+            //std::cout<<"A*cubbV-b after global solve:"<<
+            std::cout<<i<<","<<(sudata.W*(sudata.A*currV-b)).squaredNorm()<<std::endl;
+            //std::cout<<"b:"<<b.block(b.rows()-1, 0, 1, b.cols())<<std::endl;
+            //exit(0);
+            double currChange=(currV-prevV).lpNorm<Infinity>();
+            //std::cout<<"Iteration: "<<i<<", currChange: "<<currChange<<std::endl;
+            prevV=currV;
+            if (currChange<vTolerance)
+                break;
+            
+        }
+    }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#include "shapeup.cpp"
+#endif
+
+#endif

include/igl/shapeup.h

@@ -0,0 +1,108 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2017 Amir Vaxman <avaxman@gmail.com>
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IGL_SHAPEUP_H
+#define IGL_SHAPEUP_H
+
+#include <igl/igl_inline.h>
+#include <igl/setdiff.h>
+#include <igl/cat.h>
+#include <Eigen/Core>
+#include <vector>
+
+
+//This file implements the following algorithm:
+
+//Boaziz et al.
+//Shape-Up: Shaping Discrete Geometry with Projections
+//Computer Graphics Forum (Proc. SGP) 31(5), 2012
+
+namespace igl
+{
+    
+    struct ShapeupData{
+
+        //input data
+        Eigen::MatrixXd P;
+        Eigen::VectorXi SD;
+        Eigen::MatrixXi S;
+        Eigen::VectorXi b;
+        double shapeCoeff, closeCoeff;
+        std::function<bool(const MatrixXd&, const VectorXi&, const MatrixXi&,  Eigen::MatrixXd&)> local_projection;
+        
+        //Internally-used matrices
+        Eigen::SparseMatrix<double> A, Q, C, E, At, W;
+        
+         min_quad_with_fixed_data<double> solver_data;
+    };
+    
+    
+    
+    //This function precomputation the necessary matrices for the ShapeUp process, and prefactorizes them.
+    
+    //input:
+    //  P   #P by 3             point positions
+    //  SC  #Set by 1           cardinalities of sets in S
+    //  S   #Sets by max(SC)    independent sets where the local projection applies. Values beyond column SC(i) in row S(i,:) are "don't care"
+    //  b   #b by 1             boundary (fixed) vertices from P.
+    //  w   #Set by 1           weight for each set (used in the global step)
+    //  local_projection              function pointer taking (P,SC,S,projP),
+    // where the first three parameters are as defined, and "projP" is the output, as a #S by 3*max(SC) function in format xyzxyzxyz, and where it returns the projected points corresponding to each set in S in the same order.
+    //  shapecoeff,
+    //  closeCoeff,
+    //  smoothCoeff             energy coefficients as mentioned in the paper
+    
+    // Output:
+    //  sudata struct ShapeupData     the data necessary to solve the system in shapeup_solve
+
+    template <
+    typename DerivedP,
+    typename DerivedSC,
+    typename DerivedS,
+    typename Derivedb,
+    typename Derivedw>
+    IGL_INLINE void shapeup_precomputation(
+                                           const Eigen::PlainObjectBase<DerivedP>& P,
+                                           const Eigen::PlainObjectBase<DerivedSC>& SC,
+                                           const Eigen::PlainObjectBase<DerivedS>& S,
+                                           const Eigen::PlainObjectBase<Derivedb>& b,
+                                           const Eigen::PlainObjectBase<Derivedw>& w,
+                                           const std::function<bool(const Eigen::PlainObjectBase<DerivedP>&, const Eigen::PlainObjectBase<DerivedSX>&, const Eigen::PlainObjectBase<DerivedS>&,  Eigen::PlainObjectBase<Derivedb&)>& local_projection,
+                                           const double shapeCoeff,
+                                           const double closeCoeff,
+                                           const double smoothCoeff,
+                                           struct ShapeupData& sudata);
+    
+    //This function solve the shapeup project optimization. shapeup_precompute must be called before with the same sudata, or results are unpredictable
+    
+    //Input:
+    //bc                #b by 3 fixed point values corresonding to "b" in sudata
+    //NOTE: the input values in P0 don't need to correspond to prescribed values in bc; the iterations will project them automatically (by design).
+    //P0                #P by 3 initial solution (point positions)
+    //maxIterations     referring to number of local-global pairs.
+    //pTolerance        algorithm stops when max(|P_k-P_{k-1}|)<pTolerance.
+    //Output:
+    //P                 the solution to the problem, corresponding to P0.
+    template <
+    typename Derivedbc,
+    typename DerivedP>
+    IGL_INLINE void shapeup_compute(const Eigen::PlainObjectBase<Derivedbc>& bc,
+                                   const Eigen::PlainObjectBase<DerivedP>& P0,
+                                    const struct ShapeupData& sudata,
+                                    const int maxIterations=50,
+                                    const double pTolerance=10e-6
+                                    const Eigen::PlainObjectBase<DerivedP>& P)
+    {
+      
+    }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#include "shapeup.cpp"
+#endif
+
+#endif

tutorial/801_LevenbergMarquadt/CMakeLists.txt

@@ -0,0 +1,88 @@
+cmake_minimum_required(VERSION 2.8.12)
+project(levenberg-marquadt)
+
+set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${CMAKE_CURRENT_SOURCE_DIR}/cmake)
+
+find_package(LIBIGL QUIET)
+find_package(LIBHEDRA QUIET)
+
+if (NOT LIBIGL_FOUND)
+   message(FATAL_ERROR "libigl not found --- You can download it using: \n git clone --recursive https://github.com/libigl/libigl.git ${PROJECT_SOURCE_DIR}/../libigl")
+endif()
+
+if (NOT LIBHEDRA_FOUND)
+   message(FATAL_ERROR "libhedra not found --- You can download it in https://github.com/avaxman/libhedra.git")
+endif()
+
+# Compilation flags: adapt to your needs 
+if(MSVC)
+  # Enable parallel compilation
+  set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /MP /bigobj") 
+  set(CMAKE_RUNTIME_OUTPUT_DIRECTORY_DEBUG ${CMAKE_BINARY_DIR} )
+  set(CMAKE_RUNTIME_OUTPUT_DIRECTORY_RELEASE ${CMAKE_BINARY_DIR} )
+else()
+  # Libigl requires a modern C++ compiler that supports c++11
+  set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11") 
+  set(CMAKE_RUNTIME_OUTPUT_DIRECTORY "." )
+  set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wno-deprecated-declarations")
+endif()
+
+# libigl options: choose between header only and compiled static library
+# Header-only is preferred for small projects. For larger projects the static build
+# considerably reduces the compilation times
+option(LIBIGL_USE_STATIC_LIBRARY "Use LibIGL as static library" OFF)
+
+# add a customizable menu bar
+option(LIBIGL_WITH_NANOGUI     "Use Nanogui menu"   OFF)
+
+# libigl options: choose your dependencies (by default everything is OFF except opengl) 
+option(LIBIGL_WITH_VIEWER      "Use OpenGL viewer"  ON)
+option(LIBIGL_WITH_OPENGL      "Use OpenGL"         ON)
+option(LIBIGL_WITH_OPENGL_GLFW "Use GLFW"           ON)
+option(LIBIGL_WITH_BBW         "Use BBW"            OFF)
+option(LIBIGL_WITH_EMBREE      "Use Embree"         OFF)
+option(LIBIGL_WITH_PNG         "Use PNG"            OFF)
+option(LIBIGL_WITH_TETGEN      "Use Tetgen"         OFF)
+option(LIBIGL_WITH_TRIANGLE    "Use Triangle"       OFF)
+option(LIBIGL_WITH_XML         "Use XML"            OFF)
+option(LIBIGL_WITH_LIM         "Use LIM"            OFF)
+option(LIBIGL_WITH_COMISO      "Use CoMiso"         OFF)
+option(LIBIGL_WITH_MATLAB      "Use Matlab"         OFF) # This option is not supported yet
+option(LIBIGL_WITH_MOSEK       "Use MOSEK"          OFF) # This option is not supported yet
+option(LIBIGL_WITH_CGAL        "Use CGAL"           OFF)
+if(LIBIGL_WITH_CGAL) # Do not remove or move this block, the cgal build system fails without it
+  find_package(CGAL REQUIRED)
+  set(CGAL_DONT_OVERRIDE_CMAKE_FLAGS TRUE CACHE BOOL "CGAL's CMAKE Setup is super annoying ")
+  include(${CGAL_USE_FILE})
+endif()
+
+# Adding libigl: choose the path to your local copy libigl 
+# This is going to compile everything you requested 
+#message(FATAL_ERROR "${PROJECT_SOURCE_DIR}/../libigl/cmake")
+add_subdirectory("${LIBIGL_INCLUDE_DIR}/../shared/cmake" "libigl")
+
+# libigl information 
+message("libigl includes: ${LIBIGL_INCLUDE_DIRS}")
+message("libigl libraries: ${LIBIGL_LIBRARIES}")
+message("libigl extra sources: ${LIBIGL_EXTRA_SOURCES}")
+message("libigl extra libraries: ${LIBIGL_EXTRA_LIBRARIES}")
+message("libigl definitions: ${LIBIGL_DEFINITIONS}")
+
+message("libhedra includes: ${LIBHEDRA_INCLUDE_DIRS}")
+
+# Prepare the build environment
+include_directories(${LIBIGL_INCLUDE_DIRS})
+add_definitions(${LIBIGL_DEFINITIONS})
+
+include_directories(${LIBHEDRA_INCLUDE_DIRS})
+
+# Store location of data directory
+set(DATA_PATH ${CMAKE_CURRENT_SOURCE_DIR}/../data CACHE PATH "location of mesh data")
+add_definitions("-DDATA_PATH=\"${DATA_PATH}\"")
+
+include_directories(${CMAKE_CURRENT_SOURCE_DIR})
+
+# Add your project files
+FILE(GLOB SRCFILES *.cpp)
+add_executable(${PROJECT_NAME}_bin ${SRCFILES} ${LIBIGL_EXTRA_SOURCES})
+target_link_libraries(${PROJECT_NAME}_bin ${LIBIGL_LIBRARIES} ${LIBIGL_EXTRA_LIBRARIES})

tutorial/801_LevenbergMarquadt/cmake/FindLIBHEDRA.cmake

@@ -0,0 +1,24 @@
+# - Try to find the LIBHEDRA library
+# Once done this will define
+#
+#  LIBHEDRA_FOUND - system has LIBHEDRA
+#  LIBHEDRA_INCLUDE_DIR - **the** LIBHEDRA include directory
+#  LIBHEDRA_INCLUDE_DIRS - LIBHEDRA include directories
+#  LIBHEDRAL_SOURCES - the LIBHEDRA source files
+if(NOT LIBHEDRA_FOUND)
+message("hello")
+
+FIND_PATH(LIBHEDRA_INCLUDE_DIR hedra/polygonal_read_OFF.h
+   ${PROJECT_SOURCE_DIR}/../../include
+   ${PROJECT_SOURCE_DIR}/../include
+   ${PROJECT_SOURCE_DIR}/include
+   /usr/include
+   /usr/local/include
+)
+
+if(LIBHEDRA_INCLUDE_DIR)
+   set(LIBHEDRA_FOUND TRUE)
+   set(LIBHEDRA_INCLUDE_DIRS ${LIBHEDRA_INCLUDE_DIR})
+endif()
+
+endif()

tutorial/801_LevenbergMarquadt/cmake/FindLIBIGL.cmake

@@ -0,0 +1,35 @@
+# - Try to find the LIBIGL library
+# Once done this will define
+#
+#  LIBIGL_FOUND - system has LIBIGL
+#  LIBIGL_INCLUDE_DIR - **the** LIBIGL include directory
+#  LIBIGL_INCLUDE_DIRS - LIBIGL include directories
+#  LIBIGL_SOURCES - the LIBIGL source files
+if(NOT LIBIGL_FOUND)
+
+FIND_PATH(LIBIGL_INCLUDE_DIR igl/readOBJ.h
+   ${PROJECT_SOURCE_DIR}/../../include
+   ${PROJECT_SOURCE_DIR}/../include
+   ${PROJECT_SOURCE_DIR}/include
+   ${PROJECT_SOURCE_DIR}/../external/libigl/include
+   ${PROJECT_SOURCE_DIR}/../../external/libigl/include
+   $ENV{LIBIGL}/include
+   $ENV{LIBIGLROOT}/include
+   $ENV{LIBIGL_ROOT}/include
+   $ENV{LIBIGL_DIR}/include
+   $ENV{LIBIGL_DIR}/inc
+   /usr/include
+   /usr/local/include
+   /usr/local/igl/libigl/include
+)
+
+
+if(LIBIGL_INCLUDE_DIR)
+   set(LIBIGL_FOUND TRUE)
+   set(LIBIGL_INCLUDE_DIRS ${LIBIGL_INCLUDE_DIR}  ${LIBIGL_INCLUDE_DIR}/../external/Singular_Value_Decomposition)
+   #set(LIBIGL_SOURCES
+   #   ${LIBIGL_INCLUDE_DIR}/igl/viewer/Viewer.cpp
+   #)
+endif()
+
+endif()

tutorial/801_LevenbergMarquadt/main.cpp

@@ -0,0 +1,218 @@
+#include <igl/unproject_onto_mesh.h>
+#include <igl/viewer/Viewer.h>
+#include <igl/readDMAT.h>
+#include <igl/jet.h>
+#include <hedra/polygonal_read_OFF.h>
+#include <hedra/triangulate_mesh.h>
+#include <hedra/polygonal_edge_topology.h>
+#include <hedra/point_spheres.h>
+#include <hedra/scalar2RGB.h>
+#include <hedra/LMSolver.h>
+#include <hedra/EigenSolverWrapper.h>
+#include <hedra/DiscreteShellsTraits.h>
+#include <Eigen/SparseCholesky>
+#include <hedra/check_traits.h>
+
+
+std::vector<int> Handles;
+std::vector<Eigen::RowVector3d> HandlePoses;
+int CurrentHandle;
+Eigen::MatrixXd VOrig, V;
+Eigen::MatrixXi F, T;
+Eigen::VectorXi D, TF;
+Eigen::MatrixXi EV, EF, FE, EFi;
+Eigen::MatrixXd FEs;
+Eigen::VectorXi innerEdges;
+Eigen::Vector3d spans;
+bool Editing=false;
+bool ChoosingHandleMode=false;
+double CurrWinZ;
+typedef hedra::optimization::EigenSolverWrapper<Eigen::SimplicialLLT<Eigen::SparseMatrix<double> > > LinearSolver;
+hedra::optimization::DiscreteShellsTraits dst;
+LinearSolver esw;
+hedra::optimization::LMSolver<LinearSolver, hedra::optimization::DiscreteShellsTraits> lmSolver;
+
+
+
+bool UpdateCurrentView(igl::viewer::Viewer & viewer)
+{
+    using namespace Eigen;
+    using namespace std;
+    
+    MatrixXd sphereV;
+    MatrixXi sphereT;
+    MatrixXd sphereTC;
+    Eigen::MatrixXd bc(Handles.size(),V.cols());
+    for (int i=0;i<Handles.size();i++)
+        bc.row(i)=HandlePoses[i].transpose();
+    
+    
+    double sphereRadius=spans.sum()/200.0;
+    MatrixXd sphereGreens(Handles.size(),3);
+    sphereGreens.col(0).setZero();
+    sphereGreens.col(1).setOnes();
+    sphereGreens.col(2).setZero();
+
+    hedra::point_spheres(bc, sphereRadius, sphereGreens, 10, false, sphereV, sphereT, sphereTC);
+    
+    Eigen::MatrixXd bigV(V.rows()+sphereV.rows(),3);
+    Eigen::MatrixXi bigT(T.rows()+sphereT.rows(),3);
+    if (sphereV.rows()!=0){
+        bigV<<V, sphereV;
+        bigT<<T, sphereT+Eigen::MatrixXi::Constant(sphereT.rows(), sphereT.cols(), V.rows());
+    } else{
+        bigV<<V;
+        bigT<<T;
+    }
+    
+    viewer.core.show_lines=false;
+    Eigen::MatrixXd OrigEdgeColors(EV.rows(),3);
+    OrigEdgeColors.col(0)=Eigen::VectorXd::Constant(EV.rows(),0.0);
+    OrigEdgeColors.col(1)=Eigen::VectorXd::Constant(EV.rows(),0.0);
+    OrigEdgeColors.col(2)=Eigen::VectorXd::Constant(EV.rows(),0.0);
+    
+    viewer.data.clear();
+    viewer.data.set_mesh(bigV,bigT);
+    viewer.data.compute_normals();
+    viewer.data.set_edges(V,EV,OrigEdgeColors);
+    return true;
+}
+
+bool mouse_move(igl::viewer::Viewer& viewer, int mouse_x, int mouse_y)
+{
+    if (!Editing)
+        return false;
+    
+    double x = viewer.current_mouse_x;
+    double y = viewer.core.viewport(3) - viewer.current_mouse_y;
+    Eigen::RowVector3f NewPos=igl::unproject<float>(Eigen::Vector3f(x,y,CurrWinZ),
+                                                    viewer.core.view * viewer.core.model,
+                                                    viewer.core.proj,
+                                                    viewer.core.viewport);
+    
+    HandlePoses[HandlePoses.size()-1]=NewPos.cast<double>();
+    Eigen::RowVector3d Diff=HandlePoses[HandlePoses.size()-1]-VOrig.row(Handles[HandlePoses.size()-1]);
+    
+    Eigen::MatrixXd bc(Handles.size(),V.cols());
+    for (int i=0;i<Handles.size();i++)
+        bc.row(i)=HandlePoses[i].transpose();
+    
+    dst.qh=bc;
+    lmSolver.solve(true);
+    V=dst.fullSolution;
+    UpdateCurrentView(viewer);
+    return true;
+
+}
+
+
+bool mouse_up(igl::viewer::Viewer& viewer, int button, int modifier)
+{
+    if (((igl::viewer::Viewer::MouseButton)button==igl::viewer::Viewer::MouseButton::Left))
+        return false;
+    
+    Editing=false;
+ 
+    return true;
+}
+
+bool mouse_down(igl::viewer::Viewer& viewer, int button, int modifier)
+{
+    if (((igl::viewer::Viewer::MouseButton)button==igl::viewer::Viewer::MouseButton::Left))
+        return false;
+    int vid, fid;
+    Eigen::Vector3f bc;
+    double x = viewer.current_mouse_x;
+    double y = viewer.core.viewport(3) - viewer.current_mouse_y;
+    if (!ChoosingHandleMode){
+        Editing=true;
+        return false;
+    }
+    if(igl::unproject_onto_mesh(Eigen::Vector2f(x,y), viewer.core.view * viewer.core.model,
+                                viewer.core.proj, viewer.core.viewport, V, F, fid, bc))
+    {
+        //add the closest vertex to the handles
+        Eigen::MatrixXf::Index maxRow, maxCol;
+        bc.maxCoeff(&maxRow);
+        int CurrVertex=F(fid, maxRow);
+        bool Found=false;
+        for (int i=0;i<Handles.size();i++)
+            if (Handles[i]==CurrVertex){
+                CurrVertex=Handles[i];
+                Found=true;
+            }
+        
+        if (!Found){
+            Handles.push_back(CurrVertex);
+            HandlePoses.push_back(V.row(CurrVertex));
+        }
+    
+        Eigen::Vector3f WinCoords=igl::project<float>(V.row(CurrVertex).cast<float>(),
+                                               viewer.core.view * viewer.core.model,
+                                               viewer.core.proj,
+                                               viewer.core.viewport);
+
+        CurrWinZ=WinCoords(2);
+        std::cout<<"Choosing Vertex :"<<CurrVertex<<std::endl;
+
+        Eigen::VectorXi b(Handles.size());
+        for (int i=0;i<Handles.size();i++)
+            b(i)=Handles[i];
+        
+        dst.init(VOrig, T, b, EV, EF, EFi,innerEdges);
+        lmSolver.init(&esw, &dst, 250, 10e-6);
+        UpdateCurrentView(viewer);
+    }
+    return true;
+}
+
+bool key_up(igl::viewer::Viewer& viewer, unsigned char key, int modifiers)
+{
+    switch(key)
+    {
+            
+        case '1': ChoosingHandleMode=false;
+            break;
+    }
+    return false;
+}
+
+bool key_down(igl::viewer::Viewer& viewer, unsigned char key, int modifiers)
+{
+    switch(key)
+    {
+        case '1': ChoosingHandleMode=true;
+            break;
+    }
+    return false;
+}
+
+
+int main(int argc, char *argv[])
+{
+    
+    // Load a mesh in OFF format
+    using namespace std;
+    using namespace Eigen;
+    
+    hedra::polygonal_read_OFF(DATA_PATH "/moomoo.off", V, D, F);
+    hedra::triangulate_mesh(D, F, T, TF);
+    VectorXi DT=VectorXi::Constant(T.rows(),3);
+    hedra::polygonal_edge_topology(DT, T, EV, FE, EF,EFi,FEs,innerEdges);
+    
+    spans=V.colwise().maxCoeff()-V.colwise().minCoeff();
+    
+    VOrig=V;
+    igl::viewer::Viewer viewer;
+    viewer.callback_mouse_down = &mouse_down;
+    viewer.callback_mouse_move = &mouse_move;
+    viewer.callback_mouse_up=&mouse_up;
+    viewer.callback_key_down=&key_down;
+    viewer.callback_key_up=&key_up;
+    viewer.core.background_color<<0.75,0.75,0.75,1.0;
+    UpdateCurrentView(viewer);
+    viewer.launch();
+    
+    cout<<"press 1+right button to select new handles"<<endl;
+    cout<<"press the right button and drag the edit the mesh"<<endl;
+}

tutorial/802_ShapeUp/CMakeLists.txt

@@ -0,0 +1,89 @@
+cmake_minimum_required(VERSION 2.6) 
+project(shape-up)
+
+set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${CMAKE_CURRENT_SOURCE_DIR}/cmake)
+
+find_package(LIBIGL QUIET)
+find_package(LIBHEDRA QUIET)
+
+if (NOT LIBIGL_FOUND)
+   message(FATAL_ERROR "libigl not found --- You can download it using: \n git clone --recursive https://github.com/libigl/libigl.git ${PROJECT_SOURCE_DIR}/../libigl")
+endif()
+
+if (NOT LIBHEDRA_FOUND)
+   message(FATAL_ERROR "libhedra not found --- You can download it in https://github.com/avaxman/libhedra.git")
+endif()
+
+# Compilation flags: adapt to your needs 
+if(MSVC)
+  # Enable parallel compilation
+  set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /MP /bigobj") 
+  set(CMAKE_RUNTIME_OUTPUT_DIRECTORY_DEBUG ${CMAKE_BINARY_DIR} )
+  set(CMAKE_RUNTIME_OUTPUT_DIRECTORY_RELEASE ${CMAKE_BINARY_DIR} )
+else()
+  # Libigl requires a modern C++ compiler that supports c++11
+  set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11") 
+  set(CMAKE_RUNTIME_OUTPUT_DIRECTORY "." )
+endif()
+
+set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wno-deprecated-declarations")
+
+# libigl options: choose between header only and compiled static library
+# Header-only is preferred for small projects. For larger projects the static build
+# considerably reduces the compilation times
+option(LIBIGL_USE_STATIC_LIBRARY "Use LibIGL as static library" OFF)
+
+# add a customizable menu bar
+option(LIBIGL_WITH_NANOGUI     "Use Nanogui menu"   OFF)
+
+# libigl options: choose your dependencies (by default everything is OFF except opengl) 
+option(LIBIGL_WITH_VIEWER      "Use OpenGL viewer"  ON)
+option(LIBIGL_WITH_OPENGL      "Use OpenGL"         ON)
+option(LIBIGL_WITH_GLFW        "Use GLFW"           ON)
+option(LIBIGL_WITH_BBW         "Use BBW"            OFF)
+option(LIBIGL_WITH_EMBREE      "Use Embree"         OFF)
+option(LIBIGL_WITH_PNG         "Use PNG"            OFF)
+option(LIBIGL_WITH_TETGEN      "Use Tetgen"         OFF)
+option(LIBIGL_WITH_TRIANGLE    "Use Triangle"       OFF)
+option(LIBIGL_WITH_XML         "Use XML"            OFF)
+option(LIBIGL_WITH_LIM         "Use LIM"            OFF)
+option(LIBIGL_WITH_COMISO      "Use CoMiso"         OFF)
+option(LIBIGL_WITH_MATLAB      "Use Matlab"         OFF) # This option is not supported yet
+option(LIBIGL_WITH_MOSEK       "Use MOSEK"          OFF) # This option is not supported yet
+option(LIBIGL_WITH_CGAL        "Use CGAL"           OFF)
+if(LIBIGL_WITH_CGAL) # Do not remove or move this block, the cgal build system fails without it
+  find_package(CGAL REQUIRED)
+  set(CGAL_DONT_OVERRIDE_CMAKE_FLAGS TRUE CACHE BOOL "CGAL's CMAKE Setup is super annoying ")
+  include(${CGAL_USE_FILE})
+endif()
+
+# Adding libigl: choose the path to your local copy libigl 
+# This is going to compile everything you requested 
+#message(FATAL_ERROR "${PROJECT_SOURCE_DIR}/../libigl/cmake")
+add_subdirectory("${LIBIGL_INCLUDE_DIR}/../shared/cmake" "libigl")
+
+# libigl information 
+message("libigl includes: ${LIBIGL_INCLUDE_DIRS}")
+message("libigl libraries: ${LIBIGL_LIBRARIES}")
+message("libigl extra sources: ${LIBIGL_EXTRA_SOURCES}")
+message("libigl extra libraries: ${LIBIGL_EXTRA_LIBRARIES}")
+message("libigl definitions: ${LIBIGL_DEFINITIONS}")
+
+message("libhedra includes: ${LIBHEDRA_INCLUDE_DIRS}")
+
+# Prepare the build environment
+include_directories(${LIBIGL_INCLUDE_DIRS})
+add_definitions(${LIBIGL_DEFINITIONS})
+
+include_directories(${LIBHEDRA_INCLUDE_DIRS})
+
+# Store location of data directory
+set(DATA_PATH ${CMAKE_CURRENT_SOURCE_DIR}/../data CACHE PATH "location of mesh data")
+add_definitions("-DDATA_PATH=\"${DATA_PATH}\"")
+
+include_directories(${CMAKE_CURRENT_SOURCE_DIR})
+
+# Add your project files
+FILE(GLOB SRCFILES *.cpp)
+add_executable(${PROJECT_NAME}_bin ${SRCFILES} ${LIBIGL_EXTRA_SOURCES})
+target_link_libraries(${PROJECT_NAME}_bin ${LIBIGL_LIBRARIES} ${LIBIGL_EXTRA_LIBRARIES})

tutorial/802_ShapeUp/cmake/FindLIBHEDRA.cmake

@@ -0,0 +1,24 @@
+# - Try to find the LIBHEDRA library
+# Once done this will define
+#
+#  LIBHEDRA_FOUND - system has LIBHEDRA
+#  LIBHEDRA_INCLUDE_DIR - **the** LIBHEDRA include directory
+#  LIBHEDRA_INCLUDE_DIRS - LIBHEDRA include directories
+#  LIBHEDRAL_SOURCES - the LIBHEDRA source files
+if(NOT LIBHEDRA_FOUND)
+message("hello")
+
+FIND_PATH(LIBHEDRA_INCLUDE_DIR hedra/polygonal_read_OFF.h
+   ${PROJECT_SOURCE_DIR}/../../include
+   ${PROJECT_SOURCE_DIR}/../include
+   ${PROJECT_SOURCE_DIR}/include
+   /usr/include
+   /usr/local/include
+)
+
+if(LIBHEDRA_INCLUDE_DIR)
+   set(LIBHEDRA_FOUND TRUE)
+   set(LIBHEDRA_INCLUDE_DIRS ${LIBHEDRA_INCLUDE_DIR})
+endif()
+
+endif()

tutorial/802_ShapeUp/cmake/FindLIBIGL.cmake

@@ -0,0 +1,35 @@
+# - Try to find the LIBIGL library
+# Once done this will define
+#
+#  LIBIGL_FOUND - system has LIBIGL
+#  LIBIGL_INCLUDE_DIR - **the** LIBIGL include directory
+#  LIBIGL_INCLUDE_DIRS - LIBIGL include directories
+#  LIBIGL_SOURCES - the LIBIGL source files
+if(NOT LIBIGL_FOUND)
+
+FIND_PATH(LIBIGL_INCLUDE_DIR igl/readOBJ.h
+   ${PROJECT_SOURCE_DIR}/../../include
+   ${PROJECT_SOURCE_DIR}/../include
+   ${PROJECT_SOURCE_DIR}/include
+   ${PROJECT_SOURCE_DIR}/../external/libigl/include
+   ${PROJECT_SOURCE_DIR}/../../external/libigl/include
+   $ENV{LIBIGL}/include
+   $ENV{LIBIGLROOT}/include
+   $ENV{LIBIGL_ROOT}/include
+   $ENV{LIBIGL_DIR}/include
+   $ENV{LIBIGL_DIR}/inc
+   /usr/include
+   /usr/local/include
+   /usr/local/igl/libigl/include
+)
+
+
+if(LIBIGL_INCLUDE_DIR)
+   set(LIBIGL_FOUND TRUE)
+   set(LIBIGL_INCLUDE_DIRS ${LIBIGL_INCLUDE_DIR}  ${LIBIGL_INCLUDE_DIR}/../external/Singular_Value_Decomposition)
+   #set(LIBIGL_SOURCES
+   #   ${LIBIGL_INCLUDE_DIR}/igl/viewer/Viewer.cpp
+   #)
+endif()
+
+endif()

tutorial/802_ShapeUp/main.cpp

@@ -0,0 +1,245 @@
+#include <igl/unproject_onto_mesh.h>
+#include <igl/viewer/Viewer.h>
+#include <hedra/polygonal_read_OFF.h>
+#include <hedra/triangulate_mesh.h>
+#include <hedra/polygonal_edge_topology.h>
+#include <hedra/point_spheres.h>
+#include <hedra/shapeup.h>
+#include <hedra/planarity.h>
+#include <hedra/concyclity.h>
+#include <hedra/regularity.h>
+#include <hedra/scalar2RGB.h>
+
+enum ViewingMode{PLANARITY, CONCYCLITY, REGULARITY} ViewingMode;
+
+
+std::vector<int> Handles;
+std::vector<Eigen::RowVector3d> HandlePoses;
+int CurrentHandle;
+Eigen::MatrixXd VOrig, V, C;
+Eigen::MatrixXi F, T;
+Eigen::VectorXi D, TF;
+Eigen::MatrixXi EV, EF, FE, EFi;
+Eigen::MatrixXd FEs;
+Eigen::VectorXi innerEdges;
+Eigen::Vector3d spans;
+bool Editing=false;
+bool ChoosingHandleMode=false;
+double CurrWinZ;
+
+Eigen::VectorXd planarity, concyclity, regularity;
+
+hedra::ShapeupData sudata;
+Eigen::MatrixXi SFaceRegular;
+Eigen::VectorXi SDFaceRegular;
+
+
+
+bool UpdateCurrentView(igl::viewer::Viewer & viewer)
+{
+    using namespace Eigen;
+    using namespace std;
+    
+    hedra::planarity(V,D,F,planarity);
+    hedra::concyclity(V,D,F,concyclity);
+    hedra::regularity(V,D,F,regularity);
+    
+    MatrixXd sphereV;
+    MatrixXi sphereT;
+    MatrixXd sphereTC;
+    Eigen::MatrixXd bc(Handles.size(),V.cols());
+    for (int i=0;i<Handles.size();i++)
+        bc.row(i)=HandlePoses[i].transpose();
+    
+    switch(ViewingMode){
+        case PLANARITY: hedra::scalar2RGB(planarity, 0.0,1.0, C); break;
+        case CONCYCLITY:  hedra::scalar2RGB(concyclity, 0.0,5.0, C); break;
+        case REGULARITY:  hedra::scalar2RGB(regularity, 0.0,1.0, C); break;
+    }
+ 
+    double sphereRadius=spans.sum()/200.0;
+    MatrixXd sphereGreens(Handles.size(),3);
+    sphereGreens.col(0).setZero();
+    sphereGreens.col(1).setOnes();
+    sphereGreens.col(2).setZero();
+
+    hedra::point_spheres(bc, sphereRadius, sphereGreens, 10, false, sphereV, sphereT, sphereTC);
+    
+    Eigen::MatrixXd bigV(V.rows()+sphereV.rows(),3);
+    Eigen::MatrixXi bigT(T.rows()+sphereT.rows(),3);
+    Eigen::MatrixXd bigTC(C.rows()+sphereTC.rows(),3);
+    for (int i=0;i<T.rows();i++)
+        bigTC.row(i)=C.row(TF(i));
+    bigTC.block(T.rows(),0,sphereTC.rows(),3)=sphereTC;
+    if (sphereV.rows()!=0){
+        bigV<<V, sphereV;
+        bigT<<T, sphereT+Eigen::MatrixXi::Constant(sphereT.rows(), sphereT.cols(), V.rows());
+    } else{
+        bigV<<V;
+        bigT<<T;
+    }
+    
+    viewer.core.show_lines=false;
+    Eigen::MatrixXd OrigEdgeColors(EV.rows(),3);
+    OrigEdgeColors.col(0)=Eigen::VectorXd::Constant(EV.rows(),0.0);
+    OrigEdgeColors.col(1)=Eigen::VectorXd::Constant(EV.rows(),0.0);
+    OrigEdgeColors.col(2)=Eigen::VectorXd::Constant(EV.rows(),0.0);
+    
+    viewer.data.clear();
+    viewer.data.set_mesh(bigV,bigT);
+    viewer.data.set_colors(bigTC);
+    viewer.data.compute_normals();
+    viewer.data.set_edges(V,EV,OrigEdgeColors);
+    return true;
+}
+
+bool mouse_move(igl::viewer::Viewer& viewer, int mouse_x, int mouse_y)
+{
+    if (!Editing)
+        return false;
+    
+    double x = viewer.current_mouse_x;
+    double y = viewer.core.viewport(3) - viewer.current_mouse_y;
+    Eigen::RowVector3f NewPos=igl::unproject<float>(Eigen::Vector3f(x,y,CurrWinZ),
+                                                    viewer.core.view * viewer.core.model,
+                                                    viewer.core.proj,
+                                                    viewer.core.viewport);
+    
+    HandlePoses[HandlePoses.size()-1]=NewPos.cast<double>();
+    Eigen::RowVector3d Diff=HandlePoses[HandlePoses.size()-1]-VOrig.row(Handles[HandlePoses.size()-1]);
+    
+    Eigen::MatrixXd bc(Handles.size(),V.cols());
+    for (int i=0;i<Handles.size();i++)
+        bc.row(i)=HandlePoses[i].transpose();
+    
+   
+    UpdateCurrentView(viewer);
+    return true;
+
+}
+
+
+bool mouse_up(igl::viewer::Viewer& viewer, int button, int modifier)
+{
+    if (((igl::viewer::Viewer::MouseButton)button==igl::viewer::Viewer::MouseButton::Left))
+        return false;
+    
+    Editing=false;
+ 
+    return true;
+}
+
+bool mouse_down(igl::viewer::Viewer& viewer, int button, int modifier)
+{
+    if (((igl::viewer::Viewer::MouseButton)button==igl::viewer::Viewer::MouseButton::Left))
+        return false;
+    int vid, fid;
+    Eigen::Vector3f bc;
+    double x = viewer.current_mouse_x;
+    double y = viewer.core.viewport(3) - viewer.current_mouse_y;
+    if (!ChoosingHandleMode){
+        Editing=true;
+        return false;
+    }
+    if(igl::unproject_onto_mesh(Eigen::Vector2f(x,y), viewer.core.view * viewer.core.model,
+                                viewer.core.proj, viewer.core.viewport, V, F, fid, bc))
+    {
+        //add the closest vertex to the handles
+        Eigen::MatrixXf::Index maxRow, maxCol;
+        bc.maxCoeff(&maxRow);
+        int CurrVertex=F(fid, maxRow);
+        bool Found=false;
+        for (int i=0;i<Handles.size();i++)
+            if (Handles[i]==CurrVertex){
+                CurrVertex=Handles[i];
+                Found=true;
+            }
+        
+        if (!Found){
+            Handles.push_back(CurrVertex);
+            HandlePoses.push_back(V.row(CurrVertex));
+        }
+    
+        Eigen::Vector3f WinCoords=igl::project<float>(V.row(CurrVertex).cast<float>(),
+                                               viewer.core.view * viewer.core.model,
+                                               viewer.core.proj,
+                                               viewer.core.viewport);
+
+        CurrWinZ=WinCoords(2);
+        std::cout<<"Choosing Vertex :"<<CurrVertex<<std::endl;
+
+        Eigen::VectorXi b(Handles.size());
+        for (int i=0;i<Handles.size();i++)
+            b(i)=Handles[i];
+        
+        //hedra::shapeup_precompute(V,D, F,SDFaceRegular,SFaceRegular, b,1.0,100.0, sudata);
+
+        UpdateCurrentView(viewer);
+    }
+    return true;
+}
+
+bool key_up(igl::viewer::Viewer& viewer, unsigned char key, int modifiers)
+{
+    switch(key)
+    {
+            
+        case '1': ChoosingHandleMode=false;
+            break;
+        case '2': ViewingMode=PLANARITY; break;
+        case '3': ViewingMode=CONCYCLITY; break;
+        case '4': ViewingMode=REGULARITY; break;
+
+    }
+    UpdateCurrentView(viewer);
+    return false;
+}
+
+bool key_down(igl::viewer::Viewer& viewer, unsigned char key, int modifiers)
+{
+    switch(key)
+    {
+        case '1': ChoosingHandleMode=true;
+            break;
+    }
+    return false;
+}
+
+
+int main(int argc, char *argv[])
+{
+    
+    // Load a mesh in OFF format
+    using namespace std;
+    using namespace Eigen;
+    
+    std::cout<<R"(
+    1 choose handles (with right mouse button)
+    2 Show planarity between [0,1]
+    3 Show concyclity between [0,5]
+    4 Show regularity between [0,1]
+    )";
+    
+    hedra::polygonal_read_OFF(DATA_PATH "/intersection.off", V, D, F);
+    hedra::polygonal_edge_topology(D, F, EV, FE, EF,EFi,FEs,innerEdges);
+    hedra::triangulate_mesh(D, F, T, TF);
+
+    
+    
+    spans=V.colwise().maxCoeff()-V.colwise().minCoeff();
+    
+    VOrig=V;
+    ViewingMode=REGULARITY;
+    igl::viewer::Viewer viewer;
+    viewer.callback_mouse_down = &mouse_down;
+    viewer.callback_mouse_move = &mouse_move;
+    viewer.callback_mouse_up=&mouse_up;
+    viewer.callback_key_down=&key_down;
+    viewer.callback_key_up=&key_up;
+    viewer.core.background_color<<0.75,0.75,0.75,1.0;
+    UpdateCurrentView(viewer);
+    viewer.launch();
+    
+    cout<<"press 1+right button to select new handles"<<endl;
+    cout<<"press the right button and drag the edit the mesh"<<endl;
+}