softconstraints and stuff

NSolver/SLIMProblemInterface.cc

@@ -65,6 +65,18 @@ void SLIMProblemInterface::getFace(size_t i, int &_a, int &_b, int &_c)
     _c = this->i_faces[i](2);
 }
 
+void SLIMProblemInterface::fixVertex(int i)
+{
+    this->softconstraintIndices.push_back(i);
+}
+
+void SLIMProblemInterface::addRelation(int a, int b, Eigen::Matrix2d transformation, Eigen::Vector2d translation)
+{
+    this->relationIndices.push_back(Eigen::Vector2i(a, b));
+    this->transformation.push_back(transformation);
+    this->translation.push_back(translation);
+}
+
 void SLIMProblemInterface::storeResult(Eigen::MatrixXd &res)
 {
     this->res = res;

NSolver/SLIMProblemInterface.hh

@@ -68,7 +68,8 @@ public:
   void getParamVertex(size_t i, double& _x, double& _y);
   void getFace(size_t i, int& _a, int& _b, int& _c);
 
-
+  void fixVertex(int i);
+  void addRelation(int a, int b, Eigen::Matrix2d transformation, Eigen::Vector2d translation);
 
   virtual double computeEnergy(Eigen::MatrixXd& _j, Eigen::VectorXd& _a) = 0; //Evaluates the energy. _j #Faces x 4
   virtual double computeEnergySV(Eigen::MatrixXd& _sv, Eigen::VectorXd& _a) = 0;
@@ -122,6 +123,14 @@ public:
   Eigen::MatrixXd sv; //Singular Values of the jacobians, #F x 2
 
   Eigen::VectorXd area;
+
+  //Relation between vertices
+  std::vector<Eigen::Matrix2d> transformation;
+  std::vector<Eigen::Vector2d> translation;
+  std::vector<Eigen::Vector2i> relationIndices;
+
+  //Fixation/Softconstraints
+  std::vector<int> softconstraintIndices;
 };
 
 

NSolver/SLIMSolver.cc

@@ -18,13 +18,17 @@ int SLIMSolver::solve(SLIMProblemInterface &problem, int iterations)
 {
     std::cout << "Init start" << std::endl;
     problem.init(); //Some precomputations if this is the first iteration
+
     std::cout << "Init finished" << std::endl;
     Eigen::MatrixXd j;
     Eigen::SparseMatrix<double> ata;
     Eigen::VectorXd atb;
     Eigen::VectorXd p;
     Eigen::MatrixXd dir;
+    if (!checkFlipFree(problem, j))
+        return -1;
 
+    calculateJacobians(problem, problem.par, j);
     for (int i = 0; i < iterations; i++)
     {
         std::cout << "solve():calculateJacobians()" << std::endl;
@@ -166,7 +170,7 @@ void SLIMSolver::constructLinearSystem(SLIMProblemInterface &p, Eigen::SparseMat
 
     //Building of through triplets, since i have no idea how else without mega mat. multiplication
     std::vector<Eigen::Triplet<double>> a_trip;
-    a_trip.reserve(p.dx.size()*4 + p.dy.size()*4); //The D Matrices are only weighted, therefore no chnage in size.
+    //a_trip.reserve(p.dx.size()*4 + p.dy.size()*4); //The D Matrices are only weighted, therefore no chnage in size.
 
     for (int k = 0; k < p.dx.outerSize(); k++)
     {
@@ -241,7 +245,7 @@ void SLIMSolver::constructLinearSystem(SLIMProblemInterface &p, Eigen::SparseMat
 void SLIMSolver::addProximalTerm(SLIMProblemInterface &p, Eigen::SparseMatrix<double> &ata, Eigen::VectorXd &atb)
 {
     int v = p.n_vertices();
-    double l = 10e-4;
+    double l = 10e-8;
 
     Eigen::SparseMatrix<double> id;
     id.resize(2*v,2*v);
@@ -252,6 +256,76 @@ void SLIMSolver::addProximalTerm(SLIMProblemInterface &p, Eigen::SparseMatrix<do
     atb.block(v,0,v,1) += p.par.col(1) * l;
 }
 
+void SLIMSolver::addSoftConstraint(SLIMProblemInterface &p, Eigen::SparseMatrix<double> &ata, Eigen::VectorXd &atb)
+{
+    double penalty = 10e4;
+
+    int v = p.n_vertices();
+    int n = p.softconstraintIndices.size();
+
+    Eigen::SparseMatrix<double> sc;
+    sc.resize(2*v, 2*v);
+
+    std::vector<Eigen::Triplet<double>> trip_sc;
+
+    for (int i = 0; i < n; i++)
+    {
+        trip_sc.push_back(Eigen::Triplet<double>(i+0, i+0, penalty));
+        trip_sc.push_back(Eigen::Triplet<double>(i+v, i+v, penalty));
+
+        atb(i+0) += p.par(i, 0) * penalty;
+        atb(i+v) += p.par(i, 1) * penalty;
+    }
+    sc.setFromTriplets(trip_sc.begin(), trip_sc.end());
+    ata += sc;
+}
+
+void SLIMSolver::addRelationConstraint(SLIMProblemInterface &p, Eigen::SparseMatrix<double> &ata, Eigen::VectorXd &atb)
+{
+
+    double penalty = 10e4;
+
+    int v = p.n_vertices();
+    int n = p.relationIndices.size();
+
+    Eigen::SparseMatrix<double> sc;
+    sc.resize(2*n, 2*v);
+
+    Eigen::VectorXd t;
+    t.resize(2*n);
+
+    std::vector<Eigen::Triplet<double>> trip_sc;
+
+    for (int i = 0; i < n; i++)
+    {
+        int a = p.relationIndices[i](0);
+        int b = p.relationIndices[i](1);
+
+        double m00 = p.transformation[i](0, 0);
+        double m01 = p.transformation[i](0, 1);
+        double m10 = p.transformation[i](1, 0);
+        double m11 = p.transformation[i](1, 1);
+
+        //x_a*m00 + y_a*m01 - x_b = t_x
+        trip_sc.push_back(Eigen::Triplet<double>(i, a+0, m00)); //x_a*m00
+        trip_sc.push_back(Eigen::Triplet<double>(i, a+v, m01)); //y_a*m01
+        trip_sc.push_back(Eigen::Triplet<double>(i, b+0, -1));  //x_b* -1
+        t(i) = p.translation[i](0);
+
+        //x_a*m10 + y_a*m11 - y_b = t_y
+        trip_sc.push_back(Eigen::Triplet<double>(i+n, a+0, m10)); //x_a*m10
+        trip_sc.push_back(Eigen::Triplet<double>(i+n, a+v, m11)); //y_a*m11
+        trip_sc.push_back(Eigen::Triplet<double>(i+n, b+v, -1));  //y_b* -1
+        t(i+n) = p.translation[i](1);
+
+    }
+    sc.setFromTriplets(trip_sc.begin(), trip_sc.end());
+    atb += sc.transpose() * t;
+    ata += sc.transpose() * sc;
+}
+
+
+
 void SLIMSolver::solveLinearSystem(Eigen::SparseMatrix<double> &ata, Eigen::VectorXd &atb, Eigen::VectorXd &res)
 {
     Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> solver;
@@ -281,8 +355,9 @@ void SLIMSolver::lineSearch(SLIMProblemInterface &p, Eigen::MatrixXd &dir, doubl
         newEnergy = p.computeEnergy(j, p.area);
         std::printf("lineSearch(), Iteration %d, Energy %f \n", i, newEnergy);
         alphaMax *= 0.5;
-    } while (newEnergy > curEnergy && ++i < 100);
-    p.par = newParam;
+    } while (newEnergy > curEnergy && ++i < 200);
+    if (newEnergy < curEnergy)
+        p.par = newParam;
 }
 
 void SLIMSolver::calculateSearchDirection(SLIMProblemInterface &p, Eigen::VectorXd &proxy, Eigen::MatrixXd &dir)
@@ -295,8 +370,23 @@ void SLIMSolver::calculateSearchDirection(SLIMProblemInterface &p, Eigen::Vector
     dir.col(1) = proxy.block(v,0,v,1);
 
     dir = dir - p.par;
+
 }
 
+bool SLIMSolver::checkFlipFree(SLIMProblemInterface &p, Eigen::MatrixXd &j)
+{
+    calculateJacobians(p, p.par, j);
+    for (int i = 0; i < j.rows(); i++)
+    {
+        double det = (j(i, 0)*j(i,3))-(j(i,1)*j(i,2));
+        if (det <= 0)
+        {
+            std::cerr << "Initial Param not flip-free: det:" << det << std::endl;
+            return false;
+        }
+    }
+    return true;
+}
 double SLIMSolver::findAlphaMax(SLIMProblemInterface &p, Eigen::VectorXd& dir)
 {
 #warning not implemented

NSolver/SLIMSolver.hh

@@ -68,12 +68,17 @@ private:
   void calculateJacobians(COMISO::SLIMProblemInterface& p, Eigen::MatrixXd& param, Eigen::MatrixXd& j);
   void calculateWeightAndLambda(COMISO::SLIMProblemInterface& p, Eigen::MatrixXd& j, Eigen::MatrixXd &lambda);
   void constructLinearSystem(COMISO::SLIMProblemInterface& p, Eigen::SparseMatrix<double>& ata, Eigen::VectorXd& atb);
+
   void addProximalTerm(COMISO::SLIMProblemInterface& p, Eigen::SparseMatrix<double>& ata, Eigen::VectorXd& atb);
+  void addSoftConstraint(COMISO::SLIMProblemInterface& p, Eigen::SparseMatrix<double>& ata, Eigen::VectorXd& atb);
+  void addRelationConstraint(COMISO::SLIMProblemInterface& p, Eigen::SparseMatrix<double>& ata, Eigen::VectorXd& atb);
+
   void solveLinearSystem(Eigen::SparseMatrix<double>& ata, Eigen::VectorXd& atb, Eigen::VectorXd& res);
   double findAlphaMax(COMISO::SLIMProblemInterface& p, Eigen::VectorXd& dir);
   void lineSearch(COMISO::SLIMProblemInterface&p, Eigen::MatrixXd& dir, double alphaMax);
   void calculateSearchDirection(COMISO::SLIMProblemInterface& p, Eigen::VectorXd& proxy, Eigen::MatrixXd& dir);
 
+  bool checkFlipFree(COMISO::SLIMProblemInterface& p, Eigen::MatrixXd& j);
 private:
 
   Eigen::Matrix3Xd vertices;