Fix T53581: remesh modifier artifacts in sharp mode

Replace relative threshold for pseudo inverse in sharp remeshing modifier with
0.1 as proposed in the original paper.

Also change pseudo-inverse implementation that works with dynamic heap-allocated
matrix to static 3x3 version, for performance

Differential Revision: https://developer.blender.org/D5078

1 files changed, 6 insertions, 12 deletions

diff --git a/intern/dualcon/intern/octree.cpp b/intern/dualcon/intern/octree.cpp
index a6fac6fed33..70b3b8bb457 100644
--- a/intern/dualcon/intern/octree.cpp
+++ b/intern/dualcon/intern/octree.cpp
@@ -2180,19 +2180,13 @@ void Octree::countIntersection(Node *node, int height, int &nedge, int &ncell, i
 }
 
 /* from http://eigen.tuxfamily.org/bz/show_bug.cgi?id=257 */
-template<typename _Matrix_Type_>
-void pseudoInverse(const _Matrix_Type_ &a,
-                   _Matrix_Type_ &result,
-                   double epsilon = std::numeric_limits<typename _Matrix_Type_::Scalar>::epsilon())
+static void pseudoInverse(const Eigen::Matrix3f &a, Eigen::Matrix3f &result, float tolerance)
 {
-  Eigen::JacobiSVD<_Matrix_Type_> svd = a.jacobiSvd(Eigen::ComputeFullU | Eigen::ComputeFullV);
-
-  typename _Matrix_Type_::Scalar tolerance = epsilon * std::max(a.cols(), a.rows()) *
-                                             svd.singularValues().array().abs().maxCoeff();
+  Eigen::JacobiSVD<Eigen::Matrix3f> svd = a.jacobiSvd(Eigen::ComputeFullU | Eigen::ComputeFullV);
 
   result = svd.matrixV() *
-           _Matrix_Type_((svd.singularValues().array().abs() > tolerance)
-                             .select(svd.singularValues().array().inverse(), 0))
+           Eigen::Vector3f((svd.singularValues().array().abs() > tolerance)
+                               .select(svd.singularValues().array().inverse(), 0))
                .asDiagonal() *
            svd.matrixU().adjoint();
 }
@@ -2203,9 +2197,9 @@ static void solve_least_squares(const float halfA[],
                                 float rvalue[])
 {
   /* calculate pseudo-inverse */
-  Eigen::MatrixXf A(3, 3), pinv(3, 3);
+  Eigen::Matrix3f A, pinv;
   A << halfA[0], halfA[1], halfA[2], halfA[1], halfA[3], halfA[4], halfA[2], halfA[4], halfA[5];
-  pseudoInverse(A, pinv);
+  pseudoInverse(A, pinv, 0.1f);
 
   Eigen::Vector3f b2(b), mp(midpoint), result;
   b2 = b2 + A * -mp;