lscm.cpp 2.1 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283
  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. // Bug in unsupported/Eigen/SparseExtra needs iostream first
  10. #include <iostream>
  11. #include <unsupported/Eigen/SparseExtra>
  12. #include <igl/areamatrix.h>
  13. #include <igl/cotmatrix.h>
  14. #include <igl/kronecker_product.h>
  15. #include <igl/boundary_vertices_sorted.h>
  16. #include <igl/min_quad_with_fixed.h>
  17. IGL_INLINE Eigen::MatrixXd igl::lscm(
  18. const Eigen::MatrixXd& V,
  19. const Eigen::MatrixXi& F,
  20. const Eigen::VectorXi& b,
  21. const Eigen::MatrixXd& bc)
  22. {
  23. using namespace Eigen;
  24. // Assemble the area matrix (note that A is #Vx2 by #Vx2)
  25. SparseMatrix<double> A;
  26. igl::areamatrix(V,F,A);
  27. // Assemble the cotan laplacian matrix
  28. SparseMatrix<double> L;
  29. igl::cotmatrix(V,F,L);
  30. SparseMatrix<double> I2(2,2);
  31. I2.insert(0,0) = 1;
  32. I2.insert(1,1) = 1;
  33. SparseMatrix<double> L_flat = igl::kronecker_product(I2,L);
  34. VectorXi b_flat = b;
  35. MatrixXd bc_flat = bc;
  36. if (b.size() < 2)
  37. {
  38. // if no boundary conditions are provided, fix two boundary points
  39. Eigen::VectorXi bnd;
  40. igl::boundary_vertices_sorted(V,F,bnd);
  41. int v1 = bnd(0);
  42. int v2 = bnd(round(bnd.size()/2));
  43. b_flat.resize(4);
  44. b_flat << v1, v1+1, v2, v2+1;
  45. bc_flat.resize(4,1);
  46. bc_flat << 0, 1, 0, 0;
  47. }
  48. // Minimize the LSCM energy
  49. SparseMatrix<double> Q = -0.5*L_flat + A;
  50. const VectorXd B_flat = VectorXd::Zero(V.rows()*2);
  51. igl::min_quad_with_fixed_data<double> data;
  52. igl::min_quad_with_fixed_precompute(Q,b_flat,SparseMatrix<double>(),true,data);
  53. MatrixXd W_flat;
  54. if(!min_quad_with_fixed_solve(data,B_flat,bc_flat,VectorXd(),W_flat))
  55. assert(false);
  56. MatrixXd V_uv(V.rows(),2);
  57. for (unsigned i=0;i<V.rows();++i)
  58. {
  59. V_uv(i,0) = W_flat(i);
  60. V_uv(i,1) = W_flat(i+V.rows());
  61. }
  62. return V_uv;
  63. }
  64. #ifndef IGL_HEADER_ONLY
  65. // Explicit template specialization
  66. // generated by autoexplicit.sh
  67. #endif