2349 - Improve mesh domain memory usage (#2353)

Fix #2349 by switching to LibIGL style storage of polydata and associated functions.

It also generates gradient normals on demand, so if they are not needed, they are not computed, saving further memory and computation time.

Author: akenmorris

ExternalLibs/robin_hood/robin_hood.h

@@ -192,13 +192,7 @@ static Counts& counts() {
 #    define ROBIN_HOOD_PRIVATE_DEFINITION_HAS_NATIVE_WCHART() 1
 #endif
 
-// workaround missing "is_trivially_copyable" in g++ < 5.0
-// See https://stackoverflow.com/a/31798726/48181
-#if defined(__GNUC__) && __GNUC__ < 5
-#    define ROBIN_HOOD_IS_TRIVIALLY_COPYABLE(...) __has_trivial_copy(__VA_ARGS__)
-#else
-#    define ROBIN_HOOD_IS_TRIVIALLY_COPYABLE(...) std::is_trivially_copyable<__VA_ARGS__>::value
-#endif
+#define ROBIN_HOOD_IS_TRIVIALLY_COPYABLE(...) std::is_trivially_copyable<__VA_ARGS__>::value
 
 // helpers for C++ versions, see https://gcc.gnu.org/onlinedocs/cpp/Standard-Predefined-Macros.html
 #define ROBIN_HOOD_PRIVATE_DEFINITION_CXX() __cplusplus

Libs/Analyze/Analyze.cpp

@@ -75,7 +75,7 @@ static void write_good_bad_angles(json& json_object, ProjectHandle project, Anal
   std::vector<json> good_bad_angles;
 
   for (int d = 0; d < project->get_number_of_domains_per_subject(); d++) {
-    std::vector<std::shared_ptr<MeshWrapper>> meshes;
+    std::vector<std::shared_ptr<Surface>> meshes;
     for (auto& shape : analyze.get_shapes()) {
       meshes.push_back(shape->get_groomed_mesh_wrappers()[d]);
     }

Libs/Analyze/Shape.cpp

@@ -16,7 +16,7 @@
 #include <vtkKdTreePointLocator.h>
 #include <vtkPointData.h>
 
-#include "Libs/Optimize/Domain/MeshWrapper.h"
+#include "Libs/Optimize/Domain/Surface.h"
 
 using ReaderType = itk::ImageFileReader<ImageType>;
 
@@ -931,14 +931,14 @@ bool Shape::has_planes() {
 }
 
 //---------------------------------------------------------------------------
-std::vector<std::shared_ptr<MeshWrapper>> Shape::get_groomed_mesh_wrappers() {
+std::vector<std::shared_ptr<Surface>> Shape::get_groomed_mesh_wrappers() {
   if (!groomed_mesh_wrappers_.empty()) {
     return groomed_mesh_wrappers_;
   }
 
   auto group = get_groomed_meshes(true /* wait */);
   for (auto& mesh : group.meshes()) {
-    auto wrapper = std::make_shared<MeshWrapper>(mesh->get_poly_data());
+    auto wrapper = std::make_shared<Surface>(mesh->get_poly_data());
     groomed_mesh_wrappers_.push_back(wrapper);
   }
   return groomed_mesh_wrappers_;

Libs/Analyze/Shape.h

@@ -22,7 +22,7 @@ namespace shapeworks {
 class Shape;
 using ShapeHandle = std::shared_ptr<Shape>;
 using ShapeList = std::vector<ShapeHandle>;
-class MeshWrapper;
+class Surface;
 
 //! Representation of a single shape/patient/subject.
 class Shape {
@@ -187,7 +187,7 @@ class Shape {
 
   bool has_planes();
 
-  std::vector<std::shared_ptr<MeshWrapper>> get_groomed_mesh_wrappers();
+  std::vector<std::shared_ptr<Surface>> get_groomed_mesh_wrappers();
 
   void recompute_original_surface();
 
@@ -208,7 +208,7 @@ class Shape {
   MeshGroup original_meshes_;
   MeshGroup groomed_meshes_;
   MeshGroup reconstructed_meshes_;
-  std::vector<std::shared_ptr<MeshWrapper>> groomed_mesh_wrappers_;
+  std::vector<std::shared_ptr<Surface>> groomed_mesh_wrappers_;
 
   std::string override_feature_;
 

Libs/Mesh/Mesh.cpp

@@ -60,7 +60,7 @@
 
 #include "FEFixMesh.h"
 #include "Image.h"
-#include "Libs/Optimize/Domain/MeshWrapper.h"
+#include "Libs/Optimize/Domain/Surface.h"
 #include "Logging.h"
 #include "MeshComputeThickness.h"
 #include "MeshUtils.h"
@@ -785,8 +785,8 @@ double Mesh::geodesicDistance(int source, int target) const {
     throw std::invalid_argument("requested point ids outside range of points available in mesh");
   }
 
-  MeshWrapper wrap(this->poly_data_, true);
-  return wrap.ComputeDistance(getPoint(source), -1, getPoint(target), -1);
+  Surface surface(this->poly_data_, true);
+  return surface.compute_distance(getPoint(source), -1, getPoint(target), -1);
 }
 
 Field Mesh::geodesicDistance(const Point3 landmark) const {
@@ -795,10 +795,10 @@ Field Mesh::geodesicDistance(const Point3 landmark) const {
   distance->SetNumberOfTuples(numPoints());
   distance->SetName("GeodesicDistanceToLandmark");
 
-  MeshWrapper wrap(this->poly_data_, true);
+  Surface surface(this->poly_data_, true);
 
   for (int i = 0; i < numPoints(); i++) {
-    distance->SetValue(i, wrap.ComputeDistance(landmark, -1, getPoint(i), -1));
+    distance->SetValue(i, surface.compute_distance(landmark, -1, getPoint(i), -1));
   }
 
   return distance;
@@ -1767,7 +1767,7 @@ Eigen::Vector3d Mesh::getFFCGradient(Eigen::Vector3d query) const {
   return grad;
 }
 
-// WARNING: Copied directly from Meshwrapper. TODO: When refactoring, take this into account.
+// WARNING: Copied directly from Surface. TODO: When refactoring, take this into account.
 vtkSmartPointer<vtkPoints> Mesh::getIGLMesh(Eigen::MatrixXd& V, Eigen::MatrixXi& F) const {
   const int n_verts = this->poly_data_->GetNumberOfPoints();
   const int n_faces = this->poly_data_->GetNumberOfCells();

Libs/Mesh/Mesh.h

@@ -277,8 +277,8 @@ class Mesh {
 
   //! Formats mesh into an IGL format
   vtkSmartPointer<vtkPoints> getIGLMesh(Eigen::MatrixXd& V, Eigen::MatrixXi& F)
-      const;  // Copied directly from MeshWrapper. this->poly_data_ becomes this->mesh. // WARNING: Copied directly
-              // from Meshwrapper. TODO: When refactoring, take this into account.
+      const;  // Copied directly from Surface. this->poly_data_ becomes this->mesh. // WARNING: Copied directly
+              // from Surface. TODO: When refactoring, take this into account.
 
   //! Clips the mesh according to a field value
   vtkSmartPointer<vtkPolyData> clipByField(const std::string& name, double value);
@@ -293,7 +293,7 @@ class Mesh {
 
   /// Computes baricentric coordinates given a query point and a face number
   Eigen::Vector3d computeBarycentricCoordinates(const Eigen::Vector3d& pt, int face)
-      const;  // // WARNING: Copied directly from Meshwrapper. TODO: When refactoring, take this into account.
+      const;  // // WARNING: Copied directly from Surface. TODO: When refactoring, take this into account.
 
   //! Interpolates scalar values at points (e.g. correspondence points) to the mesh, assign/create a field with the
   //! given name