|
@@ -0,0 +1,79 @@
|
|
|
+//
|
|
|
+// IGL Lib - Simple C++ mesh library
|
|
|
+//
|
|
|
+// Copyright 2011, Daniele Panozzo. All rights reserved.
|
|
|
+
|
|
|
+#ifndef COTMATRIX_H
|
|
|
+#define COTMATRIX_H
|
|
|
+
|
|
|
+#define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
|
|
|
+#include <Eigen/Dense>
|
|
|
+#include <Eigen/Sparse>
|
|
|
+
|
|
|
+
|
|
|
+namespace igl
|
|
|
+{
|
|
|
+ void computeCotWeights(Eigen::Vector3d& v1, Eigen::Vector3d& v2, Eigen::Vector3d& v3, double& cotAlpha, double& cotBeta, double& cotGamma)
|
|
|
+ {
|
|
|
+ Eigen::Vector3d v12 = v2-v1;
|
|
|
+ Eigen::Vector3d v13 = v3-v1;
|
|
|
+ Eigen::Vector3d v23 = v3-v2;
|
|
|
+
|
|
|
+ double halfArea = (v12.cross(v13)).norm();//squaredNorm();
|
|
|
+
|
|
|
+ //improve numerical stability
|
|
|
+ const double cotTolerance = 1e-10;
|
|
|
+ if(halfArea < cotTolerance)
|
|
|
+ {
|
|
|
+ std::cout << "Cot weights are close to singular!" << std::endl;
|
|
|
+ halfArea = cotTolerance;
|
|
|
+ }
|
|
|
+
|
|
|
+ cotAlpha = (v12.dot(v13)) / halfArea;
|
|
|
+ cotBeta = (v23.dot(-v12)) / halfArea;
|
|
|
+ cotGamma = (-v23.dot(-v13)) / halfArea;
|
|
|
+ }
|
|
|
+
|
|
|
+ void cotmatrix (Eigen::MatrixXd V, Eigen::MatrixXi F, Eigen::SparseMatrix<double>& L, bool cotan = true)
|
|
|
+ {
|
|
|
+ Eigen::DynamicSparseMatrix<double, Eigen::RowMajor> dyn_L (V.rows(), V.rows());
|
|
|
+
|
|
|
+ for (unsigned i = 0; i < F.rows(); i++)
|
|
|
+ {
|
|
|
+ int vi1 = F(i,0);
|
|
|
+ int vi2 = F(i,1);
|
|
|
+ int vi3 = F(i,2);
|
|
|
+
|
|
|
+ if (cotan)
|
|
|
+ {
|
|
|
+ Eigen::Vector3d v1 (V(vi1,0), V(vi1,1), V(vi1,2));
|
|
|
+ Eigen::Vector3d v2 (V(vi2,0), V(vi2,1), V(vi2,2));
|
|
|
+ Eigen::Vector3d v3 (V(vi3,0), V(vi3,1), V(vi3,2));
|
|
|
+ double cot_a, cot_b, cot_g;
|
|
|
+ computeCotWeights (v1, v2, v3, cot_a, cot_b, cot_g);
|
|
|
+
|
|
|
+ dyn_L.coeffRef (vi1, vi2) += cot_g;
|
|
|
+ dyn_L.coeffRef (vi2, vi3) += cot_a;
|
|
|
+ dyn_L.coeffRef (vi3, vi1) += cot_b;
|
|
|
+ }
|
|
|
+ else
|
|
|
+ {
|
|
|
+ dyn_L.coeffRef (vi1, vi2) += 1.0;
|
|
|
+ dyn_L.coeffRef (vi2, vi3) += 1.0;
|
|
|
+ dyn_L.coeffRef (vi3, vi1) += 1.0;
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ for (int k=0; k < dyn_L.outerSize(); ++k)
|
|
|
+ {
|
|
|
+ double tmp = 0.0f;
|
|
|
+ for (Eigen::DynamicSparseMatrix<double, Eigen::RowMajor>::InnerIterator it (dyn_L, k); it; ++it)
|
|
|
+ tmp += it.value ();
|
|
|
+ dyn_L.coeffRef (k,k) = -tmp;
|
|
|
+ }
|
|
|
+
|
|
|
+ L = Eigen::SparseMatrix<double> (dyn_L);
|
|
|
+ }
|
|
|
+}
|
|
|
+
|
|
|
+#endif
|