Everything in shapeup is written, but nothing is compiled.

Former-commit-id: 301d420633dd40d49d80036a06b67c35d5875c74

include/igl/shapeup.cpp

@@ -9,6 +9,7 @@
 #define IGL_SHAPEUP_H
 
 #include "shapeup.h"
+#include "igl/min_quad_with_fixed.h"
 #include <igl/igl_inline.h>
 #include <igl/setdiff.h>
 #include <igl/cat.h>
@@ -16,120 +17,134 @@
 #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)
+    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<DerivedS>& E,
+                                           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,
+                                           struct ShapeupData& sudata)
     {
+        using namespace std;
         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
+        sudata.P=P;
+        sudata.SC=SC;
+        sudata.S=S;
+        sudata.E=E;
+        sudata.b=b;
+        sudata.local_projection=local_projection;
         
-        std::vector<Triplet<double> > QTriplets;
+        sudata.DShape.conservativeResize(SC.sum(), P.rows());  //Shape matrix (integration);
+        sudata.DClose.conservativeResize(h.rows(), P.rows());  //Closeness matrix for positional constraints
+        sudata.DSmooth.conservativeResize(EV.rows(), P.rows());  //smoothness matrix
+        
+        //Building shape matrix
+        std::vector<Triplet<double> > DShapeTriplets;
         int currRow=0;
         for (int i=0;i<S.rows();i++){
-            double avgCoeff=1.0/(double)SD(i);
+            double avgCoeff=1.0/(double)SC(i);
             
-            for (int j=0;j<SD(i);j++){
-                for (int k=0;k<SD(i);k++){
+            for (int j=0;j<SC(i);j++){
+                for (int k=0;k<SC(i);k++){
                     if (j==k)
-                        QTriplets.push_back(Triplet<double>(currRow+j, S(i,k), (1.0-avgCoeff)));
+                        DShapeTriplets.push_back(Triplet<double>(currRow+j, S(i,k), (1.0-avgCoeff)));
                     else
-                        QTriplets.push_back(Triplet<double>(currRow+j, S(i,k), (-avgCoeff)));
+                        DShapeTriplets.push_back(Triplet<double>(currRow+j, S(i,k), (-avgCoeff)));
                 }
             }
-            currRow+=SD(i);
+            currRow+=SC(i);
         }
         
-        sudata.Q.setFromTriplets(QTriplets.begin(), QTriplets.end());
+        sudata.DShape.setFromTriplets(DShapeTriplets.begin(), DShapeTriplets.end());
         
-        std::vector<Triplet<double> > CTriplets;
-        for (int i=0;i<h.size();i++)
-            CTriplets.push_back(Triplet<double>(i,h(i), 1.0));
+        //Building closeness matrix
+        std::vector<Triplet<double> > DCloseTriplets;
+        for (int i=0;i<b.size();i++)
+            DCloseTriplets.push_back(Triplet<double>(i,h(i), 1.0));
         
-        sudata.C.setFromTriplets(CTriplets.begin(), CTriplets.end());
+        sudata.DClose.setFromTriplets(DCloseTriplets.begin(), DCloseTriplets.end());
         
-        igl::cat(1, sudata.Q, sudata.C, sudata.A);
-        sudata.At=sudata.A.transpose();  //to save up this expensive computation.
+        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.
         
         //weight matrix
         vector<Triplet<double> > WTriplets;
-        //std::cout<<"w: "<<w<<std::endl;
+        
+        //one weight per set in S.
         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++)
+        
+        for (int i=0;i<b.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());
+        
+        for (int i=0;i<EV.rows();i++)
+            WTriplets.push_back(Triplet<double>(SD.sum()+b.size()+i, SD.sum()+b.size()+i, smoothCoeff));
+        
+        
+        sudata.W.conserativeResize(SD.sum()+b.size()+EV.rows(), SD.sum()+b.size()+EV.rows());
         sudata.W.setFromTriplets(WTriplets.begin(), WTriplets.end());
         
-        sudata.E=sudata.At*sudata.W*sudata.A;
-        sudata.solver.compute(sudata.E);
+        sudata.At=sudata.A.transpose();
+        sudata.Q=sudata.At*sudata.W*sudata.A;
+        
+        return min_quad_with_fixed_precompute(sudata.Q,VectorXi(),SparseMatrix<double>(),true,solver_data);
     }
     
     
     
-    IGL_INLINE void shapeup_compute(void (*projection)(int , const hedra::ShapeupData&, const Eigen::MatrixXd& , Eigen::MatrixXd&),
-                                    const Eigen::MatrixXd& vh,
+    template <
+    typename Derivedbc,
+    typename DerivedP>
+    IGL_INLINE void shapeup_solve(const Eigen::PlainObjectBase<Derivedbc>& bc,
+                                    const Eigen::PlainObjectBase<DerivedP>& P0,
                                     const struct ShapeupData& sudata,
-                                    Eigen::MatrixXd& currV,
-                                    const int maxIterations=50,
-                                    const double vTolerance=10e-6)
+                                    Eigen::PlainObjectBase<DerivedP>& P)
     {
         using namespace Eigen;
-        MatrixXd prevV=currV;
-        MatrixXd PV;
+        using namespace std;
+        MatrixXd currP;
+        MatrixXd prevP=P0;
+        MatrixXd projP;
         MatrixXd b(sudata.A.rows(),3);
-        b.block(sudata.Q.rows(), 0, sudata.h.rows(),3)=vh;  //this stays constant throughout the iterations
+        b.block(sudata.Q.rows(), 0, sudata.b.rows(),3)=bc;  //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());
+        projP.conservativeResize(sudata.SD.rows(), 3*sudata.SC.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);
+            
+            for (int j=0;j<sudata.SC.rows();j++)
+                sudata.local_projection(currV, SC,S,projP);
             //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);
+                for (int j=0;j<sudata.SC(i);j++){
+                    b.row(currRow++)=projP.block(i, 3*j, 1,3);
+                }
             }
-            //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;
+   
+            //the global solve is independent per dimension
+            Eigen::PlainObjectBase<DerivedP> rhs=-sudata.At*sudata.W*b;
+            min_quad_with_fixed_solve(sudata.solver_data, rhs,Eigen::PlainObjectBase<DerivedP>(),Eigen::PlainObjectBase<DerivedP>(), currP);
+
+            currV=sudata.solver.solve();
+            double currChange=(currP-prevP).lpNorm<Infinity>();
+            prevP=currP;
             if (currChange<vTolerance)
                 break;
             
@@ -137,8 +152,7 @@ namespace igl
     }
 }
 
-#ifndef IGL_STATIC_LIBRARY
-#include "shapeup.cpp"
-#endif
-
-#endif
+/*#ifdef IGL_STATIC_LIBRARY
+template bool igl::shapeup_solve<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, igl::ARAPData&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template bool igl::arap_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::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, igl::ARAPData&);
+#endif*/

include/igl/shapeup.h

@@ -8,6 +8,7 @@
 #ifndef IGL_SHAPEUP_H
 #define IGL_SHAPEUP_H
 
+#include "shapeup_local_projections.h"
 #include <igl/igl_inline.h>
 #include <igl/setdiff.h>
 #include <igl/cat.h>
@@ -28,16 +29,29 @@ namespace igl
 
         //input data
         Eigen::MatrixXd P;
-        Eigen::VectorXi SD;
+        Eigen::VectorXi SC;
         Eigen::MatrixXi S;
         Eigen::VectorXi b;
-        double shapeCoeff, closeCoeff;
-        std::function<bool(const MatrixXd&, const VectorXi&, const MatrixXi&,  Eigen::MatrixXd&)> local_projection;
+        //  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> A, Q, C, E, At, W;
+        Eigen::SparseMatrix<double> DShape, DClose, DSmooth, Q, A, At;
         
          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),
+        {}
     };
     
     
@@ -48,13 +62,11 @@ namespace igl
     //  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.
-    //  shapecoeff,
-    //  closeCoeff,
-    //  smoothCoeff             energy coefficients as mentioned in the paper
     
     // Output:
     //  sudata struct ShapeupData     the data necessary to solve the system in shapeup_solve
@@ -69,12 +81,10 @@ namespace igl
                                            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<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
@@ -83,22 +93,16 @@ namespace igl
     //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,
+    IGL_INLINE void shapeup_solve(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)
-    {
-      
-    }
+                                    Eigen::PlainObjectBase<DerivedP>& P);
 }
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/shapeup_local_projections.h

@@ -0,0 +1,84 @@
+// 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_LOCAL_PROJECTIONS_H
+#define IGL_SHAPEUP_LOCAL_PROJECTIONS_H
+
+#include <igl/igl_inline.h>
+#include <igl/setdiff.h>
+#include <igl/cat.h>
+#include <Eigen/Core>
+#include <vector>
+
+
+//This file implements several basic lcaol projection functions for the shapeup algorithm in shapeup.h
+
+namespace igl
+{
+    
+    //This projection does nothing but render points into projP. Mostly used for "echoing" the global step
+    bool shapeup_identity_projection(const Eigen::MatrixXd& P, const Eigen::VectorXi& SC, const Eigen::MatrixXi& S,  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<sudata.SC(i);j++)
+                avgCurrP+=currP.row(sudata.S(i,j))/(double)(sudata.SC(i));
+            
+            for (int j=0;j<sudata.SC(i);j++)
+                projP.block(i,3*j,1,3)=currP.row(sudata.S(i,j))-avgCurrP;
+        }
+    }
+    
+    
+    //the projection assumes that the sets are vertices of polygons in order
+    IGL_INLINE void shapeup_regular_face_projection(const Eigen::MatrixXd& P, const Eigen::VectorXi& SC, const Eigen::MatrixXi& S,  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 targetPos(N, 3);
+            Eigen::MatrixXd sourcePos(N, 3);
+            for (int j=0;j<N;j++)
+                avgCurrP+=currP.row(SRow(j))/(double)(N);
+            
+            for (int j=0;j<N;j++)
+                targetPolygon.row(j)=currP.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;
+        }
+        
+        
+    }
+    
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#include "shapeup_local_projections.cpp"
+#endif
+
+#endif