shapeup.cpp 7.7 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164
  1. // This file is part of libigl, a simple c++ geometry processing library.
  2. //
  3. // Copyright (C) 2017 Amir Vaxman <avaxman@gmail.com>
  4. //
  5. // This Source Code Form is subject to the terms of the Mozilla Public License
  6. // v. 2.0. If a copy of the MPL was not distributed with this file, You can
  7. // obtain one at http://mozilla.org/MPL/2.0/.
  8. #include <igl/shapeup.h>
  9. #include <igl/min_quad_with_fixed.h>
  10. #include <igl/igl_inline.h>
  11. #include <igl/setdiff.h>
  12. #include <igl/cat.h>
  13. #include <Eigen/Core>
  14. #include <vector>
  15. namespace igl
  16. {
  17. template <
  18. typename DerivedP,
  19. typename DerivedSC,
  20. typename DerivedS,
  21. typename Derivedb,
  22. typename Derivedw>
  23. IGL_INLINE bool shapeup_precomputation(const Eigen::PlainObjectBase<DerivedP>& P,
  24. const Eigen::PlainObjectBase<DerivedSC>& SC,
  25. const Eigen::PlainObjectBase<DerivedS>& S,
  26. const Eigen::PlainObjectBase<DerivedS>& E,
  27. const Eigen::PlainObjectBase<Derivedb>& b,
  28. const Eigen::PlainObjectBase<Derivedw>& w,
  29. const std::function<bool(const Eigen::PlainObjectBase<DerivedP>&, const Eigen::PlainObjectBase<DerivedSC>&, const Eigen::PlainObjectBase<DerivedS>&, Eigen::PlainObjectBase<DerivedP>&)>& local_projection,
  30. ShapeupData & sudata)
  31. {
  32. using namespace std;
  33. using namespace Eigen;
  34. sudata.P=P;
  35. sudata.SC=SC;
  36. sudata.S=S;
  37. sudata.b=b;
  38. //sudata.local_projection=local_projection;
  39. sudata.DShape.conservativeResize(SC.sum(), P.rows()); //Shape matrix (integration);
  40. sudata.DClose.conservativeResize(b.rows(), P.rows()); //Closeness matrix for positional constraints
  41. sudata.DSmooth.conservativeResize(E.rows(), P.rows()); //smoothness matrix
  42. //Building shape matrix
  43. std::vector<Triplet<double> > DShapeTriplets;
  44. int currRow=0;
  45. for (int i=0;i<S.rows();i++){
  46. double avgCoeff=1.0/(double)SC(i);
  47. for (int j=0;j<SC(i);j++){
  48. for (int k=0;k<SC(i);k++){
  49. if (j==k)
  50. DShapeTriplets.push_back(Triplet<double>(currRow+j, S(i,k), (1.0-avgCoeff)));
  51. else
  52. DShapeTriplets.push_back(Triplet<double>(currRow+j, S(i,k), (-avgCoeff)));
  53. }
  54. }
  55. currRow+=SC(i);
  56. }
  57. sudata.DShape.setFromTriplets(DShapeTriplets.begin(), DShapeTriplets.end());
  58. //Building closeness matrix
  59. std::vector<Triplet<double> > DCloseTriplets;
  60. for (int i=0;i<b.size();i++)
  61. DCloseTriplets.push_back(Triplet<double>(i,b(i), 1.0));
  62. sudata.DClose.setFromTriplets(DCloseTriplets.begin(), DCloseTriplets.end());
  63. //Building smoothness matrix
  64. std::vector<Triplet<double> > DSmoothTriplets;
  65. for (int i = 0; i < E.rows(); i++) {
  66. DSmoothTriplets.push_back(Triplet<double>(i, E(i, 0), -1));
  67. DSmoothTriplets.push_back(Triplet<double>(i, E(i, 1), 1));
  68. }
  69. igl::cat(1, sudata.DShape, sudata.DClose, sudata.A);
  70. //is this allowed? repeating A.
  71. igl::cat(1, sudata.A, sudata.DSmooth, sudata.A);
  72. //sudata.At=sudata.A.transpose(); //to save up this expensive computation.
  73. //weight matrix
  74. vector<Triplet<double> > WTriplets;
  75. //one weight per set in S.
  76. currRow=0;
  77. for (int i=0;i<SC.rows();i++){
  78. for (int j=0;j<SC(i);j++)
  79. WTriplets.push_back(Triplet<double>(currRow+j,currRow+j,sudata.shapeCoeff*w(i)));
  80. currRow+=SC(i);
  81. }
  82. for (int i=0;i<b.size();i++)
  83. WTriplets.push_back(Triplet<double>(SC.sum()+i, SC.sum()+i, sudata.closeCoeff));
  84. for (int i=0;i<E.rows();i++)
  85. WTriplets.push_back(Triplet<double>(SC.sum()+b.size()+i, SC.sum()+b.size()+i, sudata.smoothCoeff));
  86. sudata.W.conservativeResize(SC.sum()+b.size()+E.rows(), SC.sum()+b.size()+E.rows());
  87. sudata.W.setFromTriplets(WTriplets.begin(), WTriplets.end());
  88. sudata.At=sudata.A.transpose(); //for efficieny, as we use the transpose a lot in the least squares
  89. sudata.Q=sudata.At*sudata.W*sudata.A;
  90. return min_quad_with_fixed_precompute(sudata.Q,VectorXi(),SparseMatrix<double>(),true,sudata.solver_data);
  91. }
  92. template <
  93. typename Derivedbc,
  94. typename DerivedP>
  95. IGL_INLINE bool shapeup_solve(const Eigen::PlainObjectBase<Derivedbc>& bc,
  96. const Eigen::PlainObjectBase<DerivedP>& P0,
  97. const ShapeupData & sudata,
  98. Eigen::PlainObjectBase<DerivedP>& P)
  99. {
  100. using namespace Eigen;
  101. using namespace std;
  102. MatrixXd currP=P0;
  103. MatrixXd prevP=P0;
  104. MatrixXd projP;
  105. MatrixXd rhs(sudata.A.rows(), 3); rhs.setZero();
  106. rhs.block(sudata.DShape.rows(), 0, sudata.b.rows(),3)=bc; //this stays constant throughout the iterations
  107. projP.conservativeResize(sudata.SC.rows(), 3*sudata.SC.maxCoeff());
  108. for (int i=0;i<sudata.maxIterations;i++){
  109. for (int j=0;j<sudata.SC.rows();j++)
  110. sudata.local_projection(currP, sudata.SC,sudata.S,projP);
  111. //constructing the projection part of the (DShape rows of the) right hand side
  112. int currRow=0;
  113. for (int i=0;i<sudata.S.rows();i++){
  114. for (int j=0;j<sudata.SC(i);j++){
  115. rhs.row(currRow++)=projP.block(i, 3*j, 1,3);
  116. }
  117. }
  118. //the global solve is independent per dimension
  119. Eigen::PlainObjectBase<DerivedP> lsrhs=-sudata.At*sudata.W*rhs;
  120. min_quad_with_fixed_solve(sudata.solver_data, rhs,Eigen::PlainObjectBase<DerivedP>(),Eigen::PlainObjectBase<DerivedP>(), currP);
  121. double currChange=(currP-prevP).lpNorm<Infinity>();
  122. cout << "currChange: " << currChange << endl;
  123. prevP=currP;
  124. if (currChange<sudata.pTolerance)
  125. break;
  126. }
  127. return true;
  128. }
  129. }
  130. #ifdef IGL_STATIC_LIBRARY
  131. 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&);
  132. 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);
  133. #endif