Fix issue with anisotropic surface reconstruction (#132)

This revision fixes Issue #130 which was caused by the singular value clamping (taking max instead of absmax). It also includes VSCode setting fix which was outdated.

.vscode/settings.json

@@ -1,14 +1,8 @@
 // Place your settings in this file to overwrite the default settings
 {
 	"editor.tabSize": 4,
-
 	"editor.insertSpaces": true,
-
 	"files.trimTrailingWhitespace": true,
-
-	"files.insertFinalNewline": false,
-
-    "C_Cpp.clang_format_style": "file",
-
-    "C_Cpp.clang_format_formatOnSave": true
+	"files.insertFinalNewline": true,
+	"C_Cpp.clang_format_style": "file"
 }

scripts/travis_build.sh

@@ -7,7 +7,13 @@ cd build
 cmake ..
 make
 bin/unit_tests
-cd ..
-pip install --user -r requirements.txt
-pip install --user .
-python src/tests/python_tests/main.py
+
+unamestr=`uname`
+if [[ "$unamestr" == 'Darwin' ]]; then
+    echo "Disabling pip test for macOS"
+else
+    cd ..
+    pip install --user -r requirements.txt
+    pip install --user .
+    python src/tests/python_tests/main.py
+fi

src/jet/anisotropic_points_to_implicit2.cpp

@@ -82,6 +82,8 @@ void AnisotropicPointsToImplicit2::convert(
         PointKdTreeSearcher2::builder().makeShared();
     meanNeighborSearcher->build(points);
 
+    JET_INFO << "Built neighbor searcher.";
+
     SphSystemData2 meanParticles;
     meanParticles.addParticles(points);
     meanParticles.setNeighborSearcher(meanNeighborSearcher);
@@ -136,21 +138,26 @@ void AnisotropicPointsToImplicit2::convert(
             Matrix2x2D w;
             svd(cov, u, v, w);
 
+            // Take off the sign
+            v.x = std::fabs(v.x);
+            v.y = std::fabs(v.y);
+
             // Constrain Sigma
-            const double maxSingularVal = v.absmax();
+            const double maxSingularVal = v.max();
             const double kr = 4.0;
-            v[0] = std::max(v[0], maxSingularVal / kr);
-            v[1] = std::max(v[1], maxSingularVal / kr);
+            v.x = std::max(v.x, maxSingularVal / kr);
+            v.y = std::max(v.y, maxSingularVal / kr);
             const auto invSigma = Matrix2x2D::makeScaleMatrix(1.0 / v);
 
             // Compute G
-            const double relA = v[0] * v[1];  // area preservation
-            const Matrix2x2D g =
-                invH * std::sqrt(relA) * (w * invSigma * u.transposed());
+            const double scale = std::sqrt(v.x * v.y);  // area preservation
+            const Matrix2x2D g = invH * scale * (w * invSigma * u.transposed());
             gs[i] = g;
         }
     });
 
+    JET_INFO << "Computed G and means.";
+
     // SPH estimator
     meanParticles.setKernelRadius(h);
     meanParticles.updateDensities();
@@ -173,10 +180,16 @@ void AnisotropicPointsToImplicit2::convert(
         return _cutOffDensity - sum;
     });
 
+    JET_INFO << "Computed SDF.";
+
     if (_isOutputSdf) {
         FmmLevelSetSolver2 solver;
         solver.reinitialize(*temp, kMaxD, output);
+
+        JET_INFO << "Completed einitialization.";
     } else {
         temp->swap(output);
     }
+
+    JET_INFO << "Done converting points to implicit surface.";
 }

src/jet/anisotropic_points_to_implicit3.cpp

@@ -85,6 +85,8 @@ void AnisotropicPointsToImplicit3::convert(
         PointKdTreeSearcher3::builder().makeShared();
     meanNeighborSearcher->build(points);
 
+    JET_INFO << "Built neighbor searcher.";
+
     SphSystemData3 meanParticles;
     meanParticles.addParticles(points);
     meanParticles.setNeighborSearcher(meanNeighborSearcher);
@@ -139,22 +141,30 @@ void AnisotropicPointsToImplicit3::convert(
             Matrix3x3D w;
             svd(cov, u, v, w);
 
+            // Take off the sign
+            v.x = std::fabs(v.x);
+            v.y = std::fabs(v.y);
+            v.z = std::fabs(v.z);
+
             // Constrain Sigma
-            const double maxSingularVal = v.absmax();
+            const double maxSingularVal = v.max();
             const double kr = 4.0;
-            v[0] = std::max(v[0], maxSingularVal / kr);
-            v[1] = std::max(v[1], maxSingularVal / kr);
-            v[2] = std::max(v[2], maxSingularVal / kr);
+            v.x = std::max(v.x, maxSingularVal / kr);
+            v.y = std::max(v.y, maxSingularVal / kr);
+            v.z = std::max(v.z, maxSingularVal / kr);
+
             const auto invSigma = Matrix3x3D::makeScaleMatrix(1.0 / v);
 
             // Compute G
-            const double relV = v[0] * v[1] * v[2];  // area preservation
-            const Matrix3x3D g = invH * std::pow(relV, 1.0 / 3.0) *
-                                 (w * invSigma * u.transposed());
+            const double scale =
+                std::pow(v.x * v.y * v.z, 1.0 / 3.0);  // volume preservation
+            const Matrix3x3D g = invH * scale * (w * invSigma * u.transposed());
             gs[i] = g;
         }
     });
 
+    JET_INFO << "Computed G and means.";
+
     // SPH estimator
     meanParticles.setKernelRadius(h);
     meanParticles.updateDensities();
@@ -177,9 +187,13 @@ void AnisotropicPointsToImplicit3::convert(
         return _cutOffDensity - sum;
     });
 
+    JET_INFO << "Computed SDF.";
+
     if (_isOutputSdf) {
         FmmLevelSetSolver3 solver;
         solver.reinitialize(*temp, kMaxD, output);
+
+        JET_INFO << "Completed einitialization.";
     } else {
         temp->swap(output);
     }