lscm.cpp 1.9 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869
  1. // This file is part of libigl, a simple c++ geometry processing library.
  2. //
  3. // Copyright (C) 2014 Daniele Panozzo <daniele.panozzo@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 "lscm.h"
  9. #include <igl/vector_area_matrix.h>
  10. #include <igl/cotmatrix.h>
  11. //#include <igl/kronecker_product.h>
  12. #include <igl/repdiag.h>
  13. #include <igl/min_quad_with_fixed.h>
  14. #include <igl/matlab_format.h>
  15. #include <iostream>
  16. IGL_INLINE void igl::lscm(
  17. const Eigen::MatrixXd& V,
  18. const Eigen::MatrixXi& F,
  19. const Eigen::VectorXi& b,
  20. const Eigen::MatrixXd& bc,
  21. Eigen::MatrixXd & V_uv)
  22. {
  23. using namespace Eigen;
  24. using namespace std;
  25. // Assemble the area matrix (note that A is #Vx2 by #Vx2)
  26. SparseMatrix<double> A;
  27. igl::vector_area_matrix(F,A);
  28. // Assemble the cotan laplacian matrix
  29. SparseMatrix<double> L;
  30. igl::cotmatrix(V,F,L);
  31. SparseMatrix<double> L_flat;
  32. repdiag(L,2,L_flat);
  33. VectorXi b_flat(b.size()*bc.cols(),1);
  34. VectorXd bc_flat(bc.size(),1);
  35. for(int c = 0;c<bc.cols();c++)
  36. {
  37. b_flat.block(c*b.size(),0,b.rows(),1) = c*V.rows() + b.array();
  38. bc_flat.block(c*bc.rows(),0,bc.rows(),1) = bc.col(c);
  39. }
  40. // Minimize the LSCM energy
  41. SparseMatrix<double> Q = -L_flat + 2.*A;
  42. const VectorXd B_flat = VectorXd::Zero(V.rows()*2);
  43. igl::min_quad_with_fixed_data<double> data;
  44. igl::min_quad_with_fixed_precompute(Q,b_flat,SparseMatrix<double>(),true,data);
  45. MatrixXd W_flat;
  46. if(!min_quad_with_fixed_solve(data,B_flat,bc_flat,VectorXd(),W_flat))
  47. assert(false);
  48. assert(W_flat.rows() == V.rows()*2);
  49. V_uv.resize(V.rows(),2);
  50. for (unsigned i=0;i<V_uv.cols();++i)
  51. {
  52. V_uv.col(i) = W_flat.block(V_uv.rows()*i,0,V_uv.rows(),1);
  53. }
  54. }
  55. #ifdef IGL_STATIC_LIBRARY
  56. // Explicit template specialization
  57. #endif