Implementation working.

Former-commit-id: f2556c22226f087342a80071d1b8e67a87c9aa03

include/igl/shapeup.cpp

@@ -20,15 +20,13 @@ namespace igl
     typename DerivedP,
     typename DerivedSC,
     typename DerivedS,
-    typename Derivedb,
     typename Derivedw>
     IGL_INLINE bool shapeup_precomputation(const Eigen::PlainObjectBase<DerivedP>& P,
                                            const Eigen::PlainObjectBase<DerivedSC>& SC,
                                            const Eigen::PlainObjectBase<DerivedS>& S,
                                            const Eigen::PlainObjectBase<DerivedS>& E,
-                                           const Eigen::PlainObjectBase<Derivedb>& b,
+                                           const Eigen::PlainObjectBase<DerivedSC>& b,
                                            const Eigen::PlainObjectBase<Derivedw>& w,
-                                           const std::function<bool(const Eigen::PlainObjectBase<DerivedP>&, const Eigen::PlainObjectBase<DerivedSC>&, const Eigen::PlainObjectBase<DerivedS>&,  Eigen::PlainObjectBase<DerivedP>&)>& local_projection,
                                            ShapeupData & sudata)
     {
         using namespace std;
@@ -37,24 +35,24 @@ namespace igl
         sudata.SC=SC;
         sudata.S=S;
         sudata.b=b;
-        //sudata.local_projection=local_projection;
+        typedef typename DerivedP::Scalar Scalar;
         
         sudata.DShape.conservativeResize(SC.sum(), P.rows());  //Shape matrix (integration);
         sudata.DClose.conservativeResize(b.rows(), P.rows());  //Closeness matrix for positional constraints
         sudata.DSmooth.conservativeResize(E.rows(), P.rows());  //smoothness matrix
         
         //Building shape matrix
-        std::vector<Triplet<double> > DShapeTriplets;
+        std::vector<Triplet<Scalar> > DShapeTriplets;
         int currRow=0;
         for (int i=0;i<S.rows();i++){
-            double avgCoeff=1.0/(double)SC(i);
+            double avgCoeff=1.0/(Scalar)SC(i);
             
             for (int j=0;j<SC(i);j++){
                 for (int k=0;k<SC(i);k++){
                     if (j==k)
-                        DShapeTriplets.push_back(Triplet<double>(currRow+j, S(i,k), (1.0-avgCoeff)));
+                        DShapeTriplets.push_back(Triplet<Scalar>(currRow+j, S(i,k), (1.0-avgCoeff)));
                     else
-                        DShapeTriplets.push_back(Triplet<double>(currRow+j, S(i,k), (-avgCoeff)));
+                        DShapeTriplets.push_back(Triplet<Scalar>(currRow+j, S(i,k), (-avgCoeff)));
                 }
             }
             currRow+=SC(i);
@@ -63,26 +61,25 @@ namespace igl
         sudata.DShape.setFromTriplets(DShapeTriplets.begin(), DShapeTriplets.end());
         
         //Building closeness matrix
-        std::vector<Triplet<double> > DCloseTriplets;
+        std::vector<Triplet<Scalar> > DCloseTriplets;
         for (int i=0;i<b.size();i++)
-            DCloseTriplets.push_back(Triplet<double>(i,b(i), 1.0));
+            DCloseTriplets.push_back(Triplet<Scalar>(i,b(i), 1.0));
         
         sudata.DClose.setFromTriplets(DCloseTriplets.begin(), DCloseTriplets.end());
 
 		//Building smoothness matrix
-		std::vector<Triplet<double> > DSmoothTriplets;
+		std::vector<Triplet<Scalar> > DSmoothTriplets;
 		for (int i = 0; i < E.rows(); i++) {
-			DSmoothTriplets.push_back(Triplet<double>(i, E(i, 0), -1));
-			DSmoothTriplets.push_back(Triplet<double>(i, E(i, 1), 1));
+			DSmoothTriplets.push_back(Triplet<Scalar>(i, E(i, 0), -1));
+			DSmoothTriplets.push_back(Triplet<Scalar>(i, E(i, 1), 1));
 		}
         
-        igl::cat(1, sudata.DShape, sudata.DClose, sudata.A);
-        //is this allowed? repeating A.
-        igl::cat(1, sudata.A, sudata.DSmooth, sudata.A);
-        //sudata.At=sudata.A.transpose();  //to save up this expensive computation.
+        SparseMatrix<Scalar> tempMat;
+        igl::cat(1, sudata.DShape, sudata.DClose, tempMat);
+        igl::cat(1, tempMat, sudata.DSmooth, sudata.A);
         
         //weight matrix
-        vector<Triplet<double> > WTriplets;
+        vector<Triplet<Scalar> > WTriplets;
         
         //one weight per set in S.
         currRow=0;
@@ -98,23 +95,25 @@ namespace igl
         for (int i=0;i<E.rows();i++)
             WTriplets.push_back(Triplet<double>(SC.sum()+b.size()+i, SC.sum()+b.size()+i, sudata.smoothCoeff));
         
-        
         sudata.W.conservativeResize(SC.sum()+b.size()+E.rows(), SC.sum()+b.size()+E.rows());
         sudata.W.setFromTriplets(WTriplets.begin(), WTriplets.end());
         
-        sudata.At=sudata.A.transpose();  //for efficieny, as we use the transpose a lot in the least squares
+        sudata.At=sudata.A.transpose();  //for efficieny, as we use the transpose a lot in the iteration
         sudata.Q=sudata.At*sudata.W*sudata.A;
-        
+
         return min_quad_with_fixed_precompute(sudata.Q,VectorXi(),SparseMatrix<double>(),true,sudata.solver_data);
     }
     
 
     template <
-    typename Derivedbc,
-    typename DerivedP>
-    IGL_INLINE bool shapeup_solve(const Eigen::PlainObjectBase<Derivedbc>& bc,
+    typename DerivedP,
+    typename DerivedSC,
+    typename DerivedS>
+    IGL_INLINE bool shapeup_solve(const Eigen::PlainObjectBase<DerivedP>& bc,
+                                  const std::function<bool(const Eigen::PlainObjectBase<DerivedP>&, const Eigen::PlainObjectBase<DerivedSC>&, const Eigen::PlainObjectBase<DerivedS>&,  Eigen::PlainObjectBase<DerivedP>&)>& local_projection,
                                   const Eigen::PlainObjectBase<DerivedP>& P0,
                                   const ShapeupData & sudata,
+                                  const bool quiet,
                                   Eigen::PlainObjectBase<DerivedP>& P)
     {
         using namespace Eigen;
@@ -125,30 +124,33 @@ namespace igl
 		MatrixXd rhs(sudata.A.rows(), 3); rhs.setZero();
         rhs.block(sudata.DShape.rows(), 0, sudata.b.rows(),3)=bc;  //this stays constant throughout the iterations
         
+        if (!quiet){
+            cout<<"Shapeup Iterations, "<<sudata.DShape.rows()<<" constraints, solution size "<<P0.size()<<endl;
+            cout<<"**********************************************************************************************"<<endl;
+        }
         projP.conservativeResize(sudata.SC.rows(), 3*sudata.SC.maxCoeff());
-        for (int i=0;i<sudata.maxIterations;i++){
+        for (int iter=0;iter<sudata.maxIterations;iter++){
+            
+           local_projection(currP, sudata.SC,sudata.S,projP);
             
-            for (int j=0;j<sudata.SC.rows();j++)
-                sudata.local_projection(currP, sudata.SC,sudata.S,projP);
             //constructing the projection part of the (DShape rows of the) right hand side
             int currRow=0;
-            for (int i=0;i<sudata.S.rows();i++){
-                for (int j=0;j<sudata.SC(i);j++){
+            for (int i=0;i<sudata.S.rows();i++)
+                for (int j=0;j<sudata.SC(i);j++)
 					rhs.row(currRow++)=projP.block(i, 3*j, 1,3);
-                }
-            }
-   
-            //the global solve is independent per dimension
+            
             Eigen::PlainObjectBase<DerivedP> lsrhs=-sudata.At*sudata.W*rhs;
-            min_quad_with_fixed_solve(sudata.solver_data, rhs,Eigen::PlainObjectBase<DerivedP>(),Eigen::PlainObjectBase<DerivedP>(), currP);
+            MatrixXd Y(0,3), Beq(0,3);
+            min_quad_with_fixed_solve(sudata.solver_data, lsrhs,Y,Beq,currP);
 
             double currChange=(currP-prevP).lpNorm<Infinity>();
-			cout << "currChange: " << currChange << endl;
+            if (!quiet)
+                cout << "Iteration "<<iter<<", integration Linf error: "<<currChange<< endl;
             prevP=currP;
             if (currChange<sudata.pTolerance)
                 break;
-            
         }
+        P=currP;
         return true;
     }
 }
@@ -158,7 +160,7 @@ namespace igl
 
 
 #ifdef IGL_STATIC_LIBRARY
-template bool igl::shapeup_precomputation<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&,  const std::function<bool(const Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, const Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, const Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&,  Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >& ) >& local_projection, igl::ShapeupData&);
+template bool igl::shapeup_precomputation< typename Eigen::Matrix<double, -1, -1, 0, -1, -1>, typename Eigen::Matrix<int, -1, 1, 0, -1, 1>, typename Eigen::Matrix<int, -1, -1, 0, -1, -1>, typename Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&,   igl::ShapeupData&);
 
-template bool igl::shapeup_solve<typename Eigen::Matrix<double, -1, -1, 0, -1, -1>, typename Eigen::Matrix<double, -1, -1, 0, -1, -1> >(const Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >& bc, const Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >& P0, const igl::ShapeupData & sudata, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >& P);
+template bool igl::shapeup_solve<typename Eigen::Matrix<double, -1, -1, 0, -1, -1>, typename Eigen::Matrix<int, -1, 1, 0, -1, 1>, typename Eigen::Matrix<int, -1, -1, 0, -1, -1> >(const Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >& bc, const std::function<bool(const Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, const Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, const Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&,  Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >& ) >& local_projection, const Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >& P0, const igl::ShapeupData & sudata, const bool quiet, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >& P);
 #endif

include/igl/shapeup.h

@@ -25,85 +25,79 @@
 
 namespace igl
 {
+  struct ShapeupData{
+    //input data
+    Eigen::MatrixXd P;
+    Eigen::VectorXi SC;
+    Eigen::MatrixXi S;
+    Eigen::VectorXi b;
+    int maxIterations; //referring to number of local-global pairs.
+    double pTolerance;   //algorithm stops when max(|P_k-P_{k-1}|)<pTolerance.
+    double shapeCoeff, closeCoeff, smoothCoeff;
+          
+    //Internally-used matrices
+    Eigen::SparseMatrix<double> DShape, DClose, DSmooth, Q, A, At, W;
+          
+    min_quad_with_fixed_data<double> solver_data;
+          
+    ShapeupData():
+    maxIterations(50),
+    pTolerance(10e-6),
+    shapeCoeff(1.0),
+    closeCoeff(100.0),
+    smoothCoeff(0.0){}
+  };
+      
     
-    struct ShapeupData{
-        //input data
-        Eigen::MatrixXd P;
-        Eigen::VectorXi SC;
-        Eigen::MatrixXi S;
-        Eigen::VectorXi b;
-        //  maxIterations     referring to number of local-global pairs.
-        int maxIterations;
-     //  pTolerance        algorithm stops when max(|P_k-P_{k-1}|)<pTolerance.
-        double pTolerance;
-        double shapeCoeff, closeCoeff, smoothCoeff;
-        std::function<bool(const Eigen::MatrixXd&, const Eigen::VectorXi&, const Eigen::MatrixXi&,  Eigen::MatrixXd&)> local_projection;
-        
-        //Internally-used matrices
-        Eigen::SparseMatrix<double> DShape, DClose, DSmooth, Q, A, At, W;
-        
-         min_quad_with_fixed_data<double> solver_data;
-        
-        ShapeupData():
-        maxIterations(50),
-        pTolerance(10e-6),
-        shapeCoeff(1.0),
-        closeCoeff(100.0),
-        smoothCoeff(0.0),
-        local_projection(igl::shapeup_identity_projection)
-        {}
-    };
     
+  //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"
+  //  E   #E by 2             the "edges" of a mesh; used for the smoothness energy.
+  //  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.
     
-    //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"
-    //  E   #E by 2             the "edges" of a mesh; used for the smoothness energy.
-    //  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.
-    
-    // Output:
-    //  sudata struct ShapeupData     the data necessary to solve the system in shapeup_solve
+  // 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 bool shapeup_precomputation(const Eigen::PlainObjectBase<DerivedP>& P,
-                                           const Eigen::PlainObjectBase<DerivedSC>& SC,
-                                           const Eigen::PlainObjectBase<DerivedS>& S,
-                                           const Eigen::PlainObjectBase<DerivedS>& E,
-                                           const Eigen::PlainObjectBase<Derivedb>& b,
-                                           const Eigen::PlainObjectBase<Derivedw>& w,
-                                           const std::function<bool(const Eigen::PlainObjectBase<DerivedP>&, const Eigen::PlainObjectBase<DerivedSC>&, const Eigen::PlainObjectBase<DerivedS>&,  Eigen::PlainObjectBase<DerivedP>&)>& local_projection,
-                                           ShapeupData & sudata);
+  template <
+  typename DerivedP,
+  typename DerivedSC,
+  typename DerivedS,
+  typename Derivedw>
+  IGL_INLINE bool shapeup_precomputation(const Eigen::PlainObjectBase<DerivedP>& P,
+                                         const Eigen::PlainObjectBase<DerivedSC>& SC,
+                                         const Eigen::PlainObjectBase<DerivedS>& S,
+                                         const Eigen::PlainObjectBase<DerivedS>& E,
+                                         const Eigen::PlainObjectBase<DerivedSC>& b,
+                                         const Eigen::PlainObjectBase<Derivedw>& w,
+                                         ShapeupData & sudata);
     
     
     
-    //This function solve the shapeup project optimization. shapeup_precompute must be called before with the same sudata, or results are unpredictable
+  //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)
+  //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)
 
-    //Output:
-    //P                 the solution to the problem, corresponding to P0.
-    template <
-    typename Derivedbc,
-    typename DerivedP>
-    IGL_INLINE bool shapeup_solve(const Eigen::PlainObjectBase<Derivedbc>& bc,
-                                   const Eigen::PlainObjectBase<DerivedP>& P0,
-                                    const ShapeupData & sudata,
-                                    Eigen::PlainObjectBase<DerivedP>& P);
+  //Output:
+  //P                 the solution to the problem, corresponding to P0.
+  template <
+  typename DerivedP,
+  typename DerivedSC,
+  typename DerivedS>
+  IGL_INLINE bool shapeup_solve(const Eigen::PlainObjectBase<DerivedP>& bc,
+                                const std::function<bool(const Eigen::PlainObjectBase<DerivedP>&, const Eigen::PlainObjectBase<DerivedSC>&, const Eigen::PlainObjectBase<DerivedS>&,  Eigen::PlainObjectBase<DerivedP>&)>& local_projection,
+                                const Eigen::PlainObjectBase<DerivedP>& P0,
+                                const bool quiet,
+                                Eigen::PlainObjectBase<DerivedP>& P);
 }
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/shapeup_local_projections.cpp

@@ -18,61 +18,61 @@
 
 namespace igl
 {
-    
-    //This projection does nothing but render points into projP. Mostly used for "echoing" the global step
-    IGL_INLINE bool shapeup_identity_projection(const Eigen::PlainObjectBase<Eigen::MatrixXd>& P, const Eigen::PlainObjectBase<Eigen::VectorXi>& SC, const Eigen::PlainObjectBase<Eigen::MatrixXi>& S,  Eigen::PlainObjectBase<Eigen::MatrixXd>& projP){
-        projP.conservativeResize(SC.rows(), 3*SC.maxCoeff());
-        for (int i=0;i<S.rows();i++){
-            Eigen::RowVector3d avgCurrP=Eigen::RowVector3d::Zero();
-            for (int j=0;j<SC(i);j++)
-                avgCurrP+=P.row(S(i,j))/(double)(SC(i));
-            
-            for (int j=0;j<SC(i);j++)
-                projP.block(i,3*j,1,3)=P.row(S(i,j))-avgCurrP;
-        }
-        return true;
+  
+  //This projection does nothing but render points into projP. Mostly used for "echoing" the global step
+  IGL_INLINE bool shapeup_identity_projection(const Eigen::PlainObjectBase<Eigen::MatrixXd>& P, const Eigen::PlainObjectBase<Eigen::VectorXi>& SC, const Eigen::PlainObjectBase<Eigen::MatrixXi>& S,  Eigen::PlainObjectBase<Eigen::MatrixXd>& projP){
+    projP.conservativeResize(SC.rows(), 3*SC.maxCoeff());
+    for (int i=0;i<S.rows();i++){
+      Eigen::RowVector3d avgCurrP=Eigen::RowVector3d::Zero();
+      for (int j=0;j<SC(i);j++)
+        avgCurrP+=P.row(S(i,j))/(double)(SC(i));
+  
+      for (int j=0;j<SC(i);j++)
+        projP.block(i,3*j,1,3)=P.row(S(i,j))-avgCurrP;
     }
-    
-    
-    //the projection assumes that the sets are vertices of polygons in order
-    IGL_INLINE bool shapeup_regular_face_projection(const Eigen::PlainObjectBase<Eigen::MatrixXd>& P, const Eigen::PlainObjectBase<Eigen::VectorXi>& SC, const Eigen::PlainObjectBase<Eigen::MatrixXi>& S,  Eigen::PlainObjectBase<Eigen::MatrixXd>& projP){
-        projP.conservativeResize(SC.rows(), 3*SC.maxCoeff());
-        for (int currRow=0;currRow<SC.rows();currRow++){
-            //computing average
-            int N=SC(currRow);
-            const Eigen::RowVectorXi SRow=S.row(currRow);
-            Eigen::RowVector3d avgCurrP=Eigen::RowVector3d::Zero();
-            Eigen::MatrixXd targetPolygon(N, 3);
-            Eigen::MatrixXd sourcePolygon(N, 3);
-            for (int j=0;j<N;j++)
-                avgCurrP+=P.row(SRow(j))/(double)(N);
-            
-            for (int j=0;j<N;j++)
-                targetPolygon.row(j)=P.row(SRow(j))-avgCurrP;
-            
-            //creating perfectly regular source polygon
-            for (int j=0;j<N;j++)
-                sourcePolygon.row(j)<<cos(2*M_PI*(double)j/(double(N))), sin(2*M_PI*(double)j/(double(N))),0.0;
-            
-            //finding closest similarity transformation between source and target
-            Eigen::MatrixXd corrMat=sourcePolygon.transpose()*targetPolygon;
-            Eigen::JacobiSVD<Eigen::Matrix3d> svd(corrMat, Eigen::ComputeFullU | Eigen::ComputeFullV);
-            Eigen::MatrixXd R=svd.matrixU()*svd.matrixV().transpose();
-            //getting scale by edge length change average. TODO: by singular values
-            Eigen::VectorXd sourceEdgeLengths(N);
-            Eigen::VectorXd targetEdgeLengths(N);
-            for (int j=0;j<N;j++){
-                sourceEdgeLengths(j)=(sourcePolygon.row((j+1)%N)-sourcePolygon.row(j)).norm();
-                targetEdgeLengths(j)=(targetPolygon.row((j+1)%N)-targetPolygon.row(j)).norm();
-            }
-            double scale=(targetEdgeLengths.cwiseQuotient(sourceEdgeLengths)).mean();
-            
-            for (int j=0;j<N;j++)
-                projP.block(currRow,3*j,1,3)=sourcePolygon.row(j)*R*scale;
-        }
-        
-        return true;
+    return true;
+  }
+  
+  
+  //the projection assumes that the sets are vertices of polygons in order
+  IGL_INLINE bool shapeup_regular_face_projection(const Eigen::PlainObjectBase<Eigen::MatrixXd>& P, const Eigen::PlainObjectBase<Eigen::VectorXi>& SC, const Eigen::PlainObjectBase<Eigen::MatrixXi>& S,  Eigen::PlainObjectBase<Eigen::MatrixXd>& projP){
+    projP.conservativeResize(SC.rows(), 3*SC.maxCoeff());
+    for (int currRow=0;currRow<SC.rows();currRow++){
+    //computing average
+      int N=SC(currRow);
+      const Eigen::RowVectorXi SRow=S.row(currRow);
+      Eigen::RowVector3d avgCurrP=Eigen::RowVector3d::Zero();
+      Eigen::MatrixXd targetPolygon(N, 3);
+      Eigen::MatrixXd sourcePolygon(N, 3);
+      for (int j=0;j<N;j++)
+        avgCurrP+=P.row(SRow(j))/(double)(N);
+  
+      for (int j=0;j<N;j++)
+        targetPolygon.row(j)=P.row(SRow(j))-avgCurrP;
+  
+      //creating perfectly regular source polygon
+      for (int j=0;j<N;j++)
+        sourcePolygon.row(j)<<cos(2*M_PI*(double)j/(double(N))), sin(2*M_PI*(double)j/(double(N))),0.0;
+  
+      //finding closest similarity transformation between source and target
+      Eigen::MatrixXd corrMat=sourcePolygon.transpose()*targetPolygon;
+      Eigen::JacobiSVD<Eigen::Matrix3d> svd(corrMat, Eigen::ComputeFullU | Eigen::ComputeFullV);
+      Eigen::MatrixXd R=svd.matrixU()*svd.matrixV().transpose();
+      //getting scale by edge length change average. TODO: by singular values
+      Eigen::VectorXd sourceEdgeLengths(N);
+      Eigen::VectorXd targetEdgeLengths(N);
+      for (int j=0;j<N;j++){
+        sourceEdgeLengths(j)=(sourcePolygon.row((j+1)%N)-sourcePolygon.row(j)).norm();
+        targetEdgeLengths(j)=(targetPolygon.row((j+1)%N)-targetPolygon.row(j)).norm();
+      }
+      double scale=(targetEdgeLengths.cwiseQuotient(sourceEdgeLengths)).mean();
+  
+      for (int j=0;j<N;j++)
+        projP.block(currRow,3*j,1,3)=sourcePolygon.row(j)*R*scale;
     }
-    
+  
+    return true;
+  }
+
 }
 

include/igl/shapeup_local_projections.h

@@ -19,13 +19,13 @@
 
 namespace igl
 {
-    
-    //This projection does nothing but render points into projP. Mostly used for "echoing" the global step
-    IGL_INLINE bool shapeup_identity_projection(const Eigen::PlainObjectBase<Eigen::MatrixXd>& P, const Eigen::PlainObjectBase<Eigen::VectorXi>& SC, const Eigen::PlainObjectBase<Eigen::MatrixXi>& S,  Eigen::PlainObjectBase<Eigen::MatrixXd>& projP);
-    
-    //the projection assumes that the sets are vertices of polygons in order
-    IGL_INLINE bool shapeup_regular_face_projection(const Eigen::PlainObjectBase<Eigen::MatrixXd>& P, const Eigen::PlainObjectBase<Eigen::VectorXi>& SC, const Eigen::PlainObjectBase<Eigen::MatrixXi>& S,  Eigen::PlainObjectBase<Eigen::MatrixXd>& projP);
-    
+  
+  //This projection does nothing but render points into projP. Mostly used for "echoing" the global step
+  IGL_INLINE bool shapeup_identity_projection(const Eigen::PlainObjectBase<Eigen::MatrixXd>& P, const Eigen::PlainObjectBase<Eigen::VectorXi>& SC, const Eigen::PlainObjectBase<Eigen::MatrixXi>& S,  Eigen::PlainObjectBase<Eigen::MatrixXd>& projP);
+  
+  //the projection assumes that the sets are vertices of polygons in order
+  IGL_INLINE bool shapeup_regular_face_projection(const Eigen::PlainObjectBase<Eigen::MatrixXd>& P, const Eigen::PlainObjectBase<Eigen::VectorXi>& SC, const Eigen::PlainObjectBase<Eigen::MatrixXi>& S,  Eigen::PlainObjectBase<Eigen::MatrixXd>& projP);
+  
 }
 
 #ifndef IGL_STATIC_LIBRARY

tutorial/801_ShapeUp/main.cpp

@@ -7,7 +7,8 @@
 #include <igl/readDMAT.h>
 #include <igl/readOFF.h>
 #include <igl/slice.h>
-#include <igl/viewer/Viewer.h>
+#include <igl/viewer/Viewer.h>
+#include <igl/PI.h>
 #include <vector>
 #include <cstdlib>
 
@@ -29,7 +30,21 @@ igl::ShapeupData su_data;
 
 
 // Scale for visualizing the fields
-double global_scale; //TODO: not used
+double global_scale; //TODO: not used
+
+void quadAngleRegularity(const Eigen::MatrixXd& V, const Eigen::MatrixXi& Q, Eigen::VectorXd& angleRegularity)
+{
+  angleRegularity.conservativeResize(Q.rows());
+  angleRegularity.setZero();
+  for (int i=0;i<Q.rows();i++){
+    for (int j=0;j<4;j++){
+      Eigen::RowVectorXd v21=(V.row(Q(i,j))-V.row(Q(i,(j+1)%4))).normalized();
+      Eigen::RowVectorXd v23=(V.row(Q(i,(j+2)%4))-V.row(Q(i,(j+1)%4))).normalized();
+  
+      angleRegularity(i)+=(abs(acos(v21.dot(v23))-igl::PI/2.0)/(igl::PI/2.0))/4.0;
+    }
+  }
+}
 
 
 bool key_down(igl::viewer::Viewer& viewer, unsigned char key, int modifier)
@@ -37,18 +52,19 @@ bool key_down(igl::viewer::Viewer& viewer, unsigned char key, int modifier)
   using namespace std;
   using namespace Eigen;
 
-  // Plot the original quad mesh
+  // Plot the original quad mesh
+  
   if (key == '1')
-  {
-      cout<<"before setting mesh 1"<<endl;
+  {
+    viewer.data.clear();
     // Draw the triangulated quad mesh
     viewer.data.set_mesh(VQC, FQCtri);
 
-    // Assign a color to each quad that corresponds to its planarity
-    VectorXd planarity;
-    igl::quad_planarity( VQC, FQC, planarity);
+    // Assign a color to each quad that corresponds to the average deviation of each angle from pi/2
+    VectorXd angleRegularity(FQC.rows());
+    quadAngleRegularity( VQC, FQC, angleRegularity);
     MatrixXd Ct;
-    igl::jet(planarity, 0, 0.01, Ct);
+    igl::jet(angleRegularity, 0.0, 0.1, Ct);
     MatrixXd C(FQCtri.rows(),3);
     C << Ct, Ct;
     viewer.data.set_colors(C);
@@ -62,16 +78,16 @@ bool key_down(igl::viewer::Viewer& viewer, unsigned char key, int modifier)
 
   // Plot the planarized quad mesh
   if (key == '2')
-  {
+  {
+    viewer.data.clear();
     // Draw the triangulated quad mesh
-      cout<<"before setting mesh 2"<<endl;
     viewer.data.set_mesh(VQCregular, FQCtri);
 
     // Assign a color to each quad that corresponds to its planarity
-    VectorXd planarity;
-    igl::quad_planarity( VQCregular, FQC, planarity);
+    VectorXd angleRegularity(FQC.rows());
+    quadAngleRegularity( VQCregular, FQC, angleRegularity);
     MatrixXd Ct;
-    igl::jet(planarity, 0, 0.01, Ct);
+    igl::jet(angleRegularity, 0, 0.1, Ct);
     MatrixXd C(FQCtri.rows(),3);
     C << Ct, Ct;
     viewer.data.set_colors(C);
@@ -105,24 +121,22 @@ int main(int argc, char *argv[])
 
   // Create a planar version with ShapeUp
   //igl::planarize_quad_mesh(VQC, FQC, 100, 0.005, VQCregular);
-    
+  
   E.resize(FQC.size(),2);
   E.col(0)<<FQC.col(0),FQC.col(1),FQC.col(2),FQC.col(3);
   E.col(1)<<FQC.col(1),FQC.col(2),FQC.col(3),FQC.col(0);
-    
-  VectorXi b;
-  VectorXd w(FQC.rows());
-  MatrixXd bc(0,3);
-    
+  
+  VectorXi b(1); b(0)=0;  //setting the first vertex to be the same.
+  
+  VectorXd w=VectorXd::Constant(FQC.rows(),1.0);
+  MatrixXd bc(1,3); bc<<VQC.row(0);
+  
   VectorXi array_of_fours=VectorXi::Constant(FQC.rows(),4);
-    cout<<"before pre-computation"<<endl;
-    std::function<bool(const Eigen::PlainObjectBase<MatrixXd>&, const Eigen::PlainObjectBase<VectorXi>&, const Eigen::PlainObjectBase<MatrixXi>&, Eigen::PlainObjectBase<MatrixXd>&)> localFunction(igl::shapeup_identity_projection);
-    
-    shapeup_precomputation(VQC, array_of_fours,FQC,E,b,w, localFunction,su_data);
-    cout<<"after pre-computation"<<endl;
-    shapeup_solve(bc,VQC,su_data,VQCregular);
-    cout<<"after computation"<<endl;
-    
+  std::function<bool(const Eigen::PlainObjectBase<MatrixXd>&, const Eigen::PlainObjectBase<VectorXi>&, const Eigen::PlainObjectBase<MatrixXi>&, Eigen::PlainObjectBase<MatrixXd>&)> localFunction(igl::shapeup_regular_face_projection);
+  
+  su_data.maxIterations=200;
+  shapeup_precomputation(VQC, array_of_fours,FQC,E,b,w,su_data);
+  shapeup_solve(bc,localFunction, VQC,su_data, false,VQCregular);
 
   // Convert the planarized mesh to triangles
   igl::slice( VQCregular, FQC.col(0).eval(), 1, PQC0regular);