Merge pull request #2126 from SCIInstitute/2081-mesh_based_attributes

Mesh field based correspondence attributes

Author: akenmorris

Applications/shapeworks/Commands.cpp

@@ -126,7 +126,13 @@ bool OptimizeCommand::execute(const optparse::Values& options, SharedCommandData
     try {
       // load spreadsheet project
       ProjectHandle project = std::make_shared<Project>();
-      project->load(projectFile);
+      try {
+        project->load(projectFile);
+      } catch (std::exception& e) {
+        SW_ERROR("Project failed to load: {}", e.what());
+        return false;
+      }
+
 
       const auto oldBasePath = boost::filesystem::current_path();
       auto base = StringUtils::getPath(projectFile);

Libs/Mesh/Mesh.cpp

@@ -18,6 +18,7 @@
 
 // vtk
 #include <vtkAppendPolyData.h>
+#include <vtkGradientFilter.h>
 #include <vtkButterflySubdivisionFilter.h>
 #include <vtkCenterOfMass.h>
 #include <vtkCleanPolyData.h>
@@ -863,6 +864,150 @@ Field Mesh::curvature(const CurvatureType type) const {
   return curv;
 }
 
+
+void computeGradient(vtkDataSet* inputDataSet, const char* scalarFieldName, const char* gradientFieldName) {
+  vtkSmartPointer<vtkGradientFilter> gradientFilter = vtkSmartPointer<vtkGradientFilter>::New();
+  gradientFilter->SetInputData(inputDataSet);
+  gradientFilter->SetInputScalars(vtkDataSet::FIELD_ASSOCIATION_POINTS, scalarFieldName);
+  gradientFilter->SetResultArrayName(gradientFieldName);
+  gradientFilter->Update();
+  /*
+
+  vtkSmartPointer<vtkDoubleArray> gradientData = vtkSmartPointer<vtkDoubleArray>::New();
+  gradientData->SetNumberOfComponents(3);
+  gradientData->SetNumberOfTuples(inputDataSet->GetNumberOfPoints());
+  gradientData->SetName(gradientFieldName);
+*/
+  vtkDoubleArray* gradPointArray = vtkArrayDownCast<vtkDoubleArray>(
+      vtkDataSet::SafeDownCast(gradientFilter->GetOutput())->GetPointData()->GetArray(gradientFieldName));
+
+
+/*
+  for (vtkIdType pointId = 0; pointId < inputDataSet->GetNumberOfPoints(); ++pointId) {
+    double gradient[3];
+    gradientFilter->GetOutput()->GetPointData()->GetTuple(pointId, gradient);
+    gradientData->SetTuple(pointId, gradient);
+  }
+*/
+
+  //  inputDataSet->GetPointData()->AddArray(gradientData);
+  inputDataSet->GetPointData()->AddArray(gradPointArray);
+}
+
+
+void Mesh::computeFieldGradient(const std::string& field) const {
+  computeGradient(poly_data_, field.c_str(), (std::string("gradient_") + field).c_str());
+  return;
+  auto arr = poly_data_->GetPointData()->GetArray(field.c_str());
+
+  // for each vertex, compute the gradient of the field and store x,y,z
+  // components in a vector
+  auto gradient = vtkSmartPointer<vtkDoubleArray>::New();
+  gradient->SetNumberOfComponents(3);
+  gradient->SetNumberOfTuples(numPoints());
+  gradient->SetName((std::string("gradient_") + field).c_str());
+
+  for (int i = 0; i < numPoints(); i++) {
+    // get the field value at this vertex
+    double value = arr->GetTuple1(i);
+
+    // collect all neighboring vertices using vtk
+    auto cell = vtkSmartPointer<vtkGenericCell>::New();
+    auto neighbors = vtkSmartPointer<vtkIdList>::New();
+    poly_data_->GetPointCells(i, neighbors);
+    int num_neighbors = neighbors->GetNumberOfIds();
+
+    // for each neighbor vertex
+    double x = 0.0;
+    double y = 0.0;
+    double z = 0.0;
+
+    for (int j = 0; j < num_neighbors; j++) {
+      // get the neighbor vertex id
+      int neighbor_id = neighbors->GetId(j);
+
+      // get the neighbor vertex
+      auto neighbor = poly_data_->GetPoint(neighbor_id);
+
+      // get the neighbor vertex value
+      double neighbor_value = arr->GetTuple1(neighbor_id);
+
+      // compute the gradient
+      x += (neighbor[0] - getPoint(i)[0]) * (neighbor_value - value);
+      y += (neighbor[1] - getPoint(i)[1]) * (neighbor_value - value);
+      z += (neighbor[2] - getPoint(i)[2]) * (neighbor_value - value);
+    }
+
+    gradient->SetTuple3(i, x, y, z);
+  }
+
+  poly_data_->GetPointData()->AddArray(gradient);
+}
+
+Eigen::Vector3d Mesh::computeFieldGradientAtPoint(const std::string& field, const Point3& query) const {
+  this->updateCellLocator();
+
+  // compute gradient if not already computed
+  if (poly_data_->GetPointData()->GetArray((std::string("gradient_") + field).c_str()) == nullptr) {
+    computeFieldGradient(field);
+  }
+
+  double closestPoint[3];
+  vtkIdType cellId;
+  int subId;
+  double dist;
+  cellLocator->FindClosestPoint(query.data(), closestPoint, cellId, subId, dist);
+
+  auto cell = poly_data_->GetCell(cellId);
+
+  size_t v1 = cell->GetPointId(0);
+  size_t v2 = cell->GetPointId(1);
+  size_t v3 = cell->GetPointId(2);
+
+  auto gradient = poly_data_->GetPointData()->GetArray((std::string("gradient_") + field).c_str());
+
+  Eigen::Vector3d grad1(gradient->GetTuple3(v1)[0], gradient->GetTuple3(v1)[1], gradient->GetTuple3(v1)[2]);
+  Eigen::Vector3d grad2(gradient->GetTuple3(v2)[0], gradient->GetTuple3(v2)[1], gradient->GetTuple3(v2)[2]);
+  Eigen::Vector3d grad3(gradient->GetTuple3(v3)[0], gradient->GetTuple3(v3)[1], gradient->GetTuple3(v3)[2]);
+
+  // Compute barycentric distances
+  Eigen::Vector3d cp(closestPoint[0], closestPoint[1], closestPoint[2]);
+  Eigen::Vector3d bary = computeBarycentricCoordinates(cp, cellId);
+
+  bary = bary / bary.sum();
+
+  Eigen::Vector3d result;
+  result = bary[0] * grad1 + bary[1] * grad2 + bary[2] * grad3;
+  return result;
+}
+
+double Mesh::interpolateFieldAtPoint(const std::string& field, const Point3& query) const {
+  this->updateCellLocator();
+
+  double closestPoint[3];
+  vtkIdType cellId;
+  int subId;
+  double dist;
+  this->cellLocator->FindClosestPoint(query.data(), closestPoint, cellId, subId, dist);
+
+  auto cell = this->poly_data_->GetCell(cellId);
+
+  size_t v1 = cell->GetPointId(0);
+  size_t v2 = cell->GetPointId(1);
+  size_t v3 = cell->GetPointId(2);
+
+  // Compute barycentric distances
+  Eigen::Vector3d cp(closestPoint[0], closestPoint[1], closestPoint[2]);
+  Eigen::Vector3d bary = computeBarycentricCoordinates(cp, cellId);
+
+  bary = bary / bary.sum();
+
+  Eigen::Vector3d values(this->getFieldValue(field, v1), this->getFieldValue(field, v2),
+                         this->getFieldValue(field, v3));
+
+  return (bary * values.transpose()).mean();
+}
+
 Mesh& Mesh::applySubdivisionFilter(const SubdivisionType type, int subdivision) {
   if (type == Mesh::SubdivisionType::Loop) {
     auto filter = vtkSmartPointer<vtkLoopSubdivisionFilter>::New();
@@ -1466,30 +1611,8 @@ Eigen::Vector3d Mesh::computeBarycentricCoordinates(const Eigen::Vector3d& pt, i
 }
 
 double Mesh::getFFCValue(Eigen::Vector3d query) const {
-  this->updateCellLocator();
-
-  double closestPoint[3];
-  vtkIdType cellId;
-  int subId;
-  double dist;
-  this->cellLocator->FindClosestPoint(query.data(), closestPoint, cellId, subId, dist);
-
-  auto cell = this->poly_data_->GetCell(cellId);
-
-  size_t v1 = cell->GetPointId(0);
-  size_t v2 = cell->GetPointId(1);
-  size_t v3 = cell->GetPointId(2);
-
-  // Compute barycentric distances
-  Eigen::Vector3d cp(closestPoint[0], closestPoint[1], closestPoint[2]);
-  Eigen::Vector3d bary = computeBarycentricCoordinates(cp, cellId);
-
-  bary = bary / bary.sum();
-
-  Eigen::Vector3d values(this->getFieldValue("value", v1), this->getFieldValue("value", v2),
-                         this->getFieldValue("value", v3));
-
-  return (bary * values.transpose()).mean();
+  Point3 point(query.data());
+  return interpolateFieldAtPoint("value", point);
 }
 
 Eigen::Vector3d Mesh::getFFCGradient(Eigen::Vector3d query) const {

Libs/Mesh/Mesh.h

@@ -162,6 +162,15 @@ class Mesh {
   /// computes curvature using principal (default) or gaussian or mean algorithms
   Field curvature(const CurvatureType type = Principal) const;
 
+  /// compute the gradient of a scalar field for all vertices
+  void computeFieldGradient(const std::string& field) const;
+
+  /// compute the gradient of a scalar field at a point
+  Eigen::Vector3d computeFieldGradientAtPoint(const std::string& field, const Point3& query) const;
+
+  /// interpolate a scalar field at a given point
+  double interpolateFieldAtPoint(const std::string& field, const Point3& query) const;
+
   /// applies subdivision filter (butterfly (default) or loop)
   Mesh& applySubdivisionFilter(const SubdivisionType type = Butterfly, int subdivision = 1);
 

Libs/Optimize/Domain/ImplicitSurfaceDomain.h

@@ -1,9 +1,6 @@
 #pragma once
 
 #include "ImageDomainWithCurvature.h"
-#include "Libs/Mesh/meshFIM.h"
-#include "TriMesh.h"
-#include "TriMesh_algo.h"
 #include "vnl/vnl_cross.h"
 #include "vnl/vnl_inverse.h"
 #include "vnl/vnl_math.h"
@@ -105,34 +102,6 @@ class ImplicitSurfaceDomain : public ImageDomainWithCurvature<T> {
     return newpoint;
   }
 
-  void SetMesh(TriMesh* mesh) {
-    m_mesh = new meshFIM();
-    m_mesh->SetMesh(mesh);
-  };
-
-  void SetFeaMesh(const char* feaFile) { m_mesh->ReadFeatureFromFile(feaFile); };
-  void SetFeaGrad(const char* feaGradFile) { m_mesh->ReadFeatureGradientFromFile(feaGradFile); };
-  void SetFids(const char* fidsFile) {
-    m_mesh->ReadFaceIndexMap(fidsFile);
-    const typename ImageType::PointType orgn = this->GetOrigin();
-    m_mesh->imageOrigin[0] = orgn[0];
-    m_mesh->imageOrigin[1] = orgn[1];
-    m_mesh->imageOrigin[2] = orgn[2];
-    typename ImageType::RegionType::SizeType sz = this->GetSize();
-    m_mesh->imageSize[0] = sz[0];
-    m_mesh->imageSize[1] = sz[1];
-    m_mesh->imageSize[2] = sz[2];
-    typename ImageType::SpacingType sp = this->GetSpacing();
-    m_mesh->imageSpacing[0] = sp[0];
-    m_mesh->imageSpacing[1] = sp[1];
-    m_mesh->imageSpacing[2] = sp[2];
-    typename ImageType::RegionType::IndexType idx = this->GetIndex();
-    m_mesh->imageIndex[0] = idx[0];
-    m_mesh->imageIndex[1] = idx[1];
-    m_mesh->imageIndex[2] = idx[2];
-  };
-  meshFIM* GetMesh() { return m_mesh; }
-  meshFIM* GetMesh() const { return m_mesh; }
 
   /** Get any valid point on the domain. This is used to place the first particle. */
   PointType GetZeroCrossingPoint() const override {
@@ -142,7 +111,7 @@ class ImplicitSurfaceDomain : public ImageDomainWithCurvature<T> {
     return p;
   }
 
-  ImplicitSurfaceDomain() : m_Tolerance(1.0e-4) { m_mesh = NULL; }
+  ImplicitSurfaceDomain() : m_Tolerance(1.0e-4) { }
   void PrintSelf(std::ostream& os, itk::Indent indent) const {
     Superclass::PrintSelf(os, indent);
     os << indent << "m_Tolerance = " << m_Tolerance << std::endl;
@@ -152,7 +121,7 @@ class ImplicitSurfaceDomain : public ImageDomainWithCurvature<T> {
  private:
   T m_Tolerance;
 
-  meshFIM* m_mesh;
+
 };
 
 }  // end namespace shapeworks

Libs/Optimize/Domain/MeshDomain.cpp

@@ -12,10 +12,10 @@
 namespace shapeworks {
 
 bool MeshDomain::ApplyConstraints(PointType &p, int idx, bool dbg) const {
-  if (!meshWrapper) {
+  if (!mesh_wrapper_) {
     return true;
   }
-  p = meshWrapper->SnapToMesh(p, idx);
+  p = mesh_wrapper_->SnapToMesh(p, idx);
   return true;
 }
 
@@ -27,7 +27,7 @@ bool MeshDomain::ApplyVectorConstraints(vnl_vector_fixed<double, DIMENSION> &gra
 
 vnl_vector_fixed<double, DIMENSION> MeshDomain::ProjectVectorToSurfaceTangent(
     vnl_vector_fixed<double, DIMENSION> &gradE, const PointType &pos, int idx) const {
-  return meshWrapper->ProjectVectorToSurfaceTangent(pos, idx, gradE);
+  return mesh_wrapper_->ProjectVectorToSurfaceTangent(pos, idx, gradE);
 }
 
 MeshDomain::PointType MeshDomain::UpdateParticlePosition(const PointType &point, int idx,
@@ -36,20 +36,20 @@ MeshDomain::PointType MeshDomain::UpdateParticlePosition(const PointType &point,
   for (unsigned int i = 0; i < DIMENSION; i++) {
     negativeUpdate[i] = -update[i];
   }
-  PointType newPoint = meshWrapper->GeodesicWalk(point, idx, negativeUpdate);
+  PointType newPoint = mesh_wrapper_->GeodesicWalk(point, idx, negativeUpdate);
   return newPoint;
 }
 
 double MeshDomain::GetMaxDiameter() const {
   // todo should this not be the length of the bounding box diagonal?
-  PointType boundingBoxSize = meshWrapper->GetMeshUpperBound() - meshWrapper->GetMeshLowerBound();
+  PointType boundingBoxSize = mesh_wrapper_->GetMeshUpperBound() - mesh_wrapper_->GetMeshLowerBound();
   double max = 0;
   for (int d = 0; d < 3; d++) {
     max = max > boundingBoxSize[d] ? max : boundingBoxSize[d];
   }
   return max;
 }
 
-void MeshDomain::InvalidateParticlePosition(int idx) const { this->meshWrapper->InvalidateParticle(idx); }
+void MeshDomain::InvalidateParticlePosition(int idx) const { this->mesh_wrapper_->InvalidateParticle(idx); }
 
 }  // namespace shapeworks