Merge pull request #2143 from SCIInstitute/cortical_thickness

Cortical thickness updates
SCIInstitute · Sep 27, 2023 · 160f306 · 160f306
2 parents 18da3d7 + e77b194
commit 160f306
Show file tree

Hide file tree

Showing 7 changed files with 370 additions and 180 deletions.
diff --git a/Applications/shapeworks/MeshCommands.cpp b/Applications/shapeworks/MeshCommands.cpp
@@ -1541,6 +1541,11 @@ void ComputeThickness::buildParser() {
       .type("double")
       .set_default(100000.0)
       .help("Maximum distance to determine thickness");
+  parser.add_option("--median_radius")
+      .action("store")
+      .type("double")
+      .set_default(5.0)
+      .help("Median radius for smoothing, multiplier of average edge length");
   parser.add_option("--distance_mesh")
       .action("store")
       .type("string")
@@ -1565,15 +1570,17 @@ bool ComputeThickness::execute(const optparse::Values &options, SharedCommandDat
 
   double max_dist = static_cast<double>(options.get("max_dist"));
 
+  double median_radius = static_cast<double>(options.get("median_radius"));
+
   std::string dt_filename = static_cast<std::string>(options.get("distance_transform"));
 
   std::string distance_mesh_filename = static_cast<std::string>(options.get("distance_mesh"));
 
   if (dt_filename == "") {
-    sharedData.mesh->computeThickness(img, nullptr, max_dist, distance_mesh_filename);
+    sharedData.mesh->computeThickness(img, nullptr, max_dist, median_radius, distance_mesh_filename);
   } else {
     Image dt(dt_filename);
-    sharedData.mesh->computeThickness(img, &dt, max_dist, distance_mesh_filename);
+    sharedData.mesh->computeThickness(img, &dt, max_dist, median_radius, distance_mesh_filename);
   }
 
   return sharedData.validMesh();

diff --git a/Libs/Mesh/Mesh.cpp b/Libs/Mesh/Mesh.cpp
@@ -1128,8 +1128,8 @@ Image Mesh::toDistanceTransform(PhysicalRegion region, const Point3 spacing, con
   return img;
 }
 
-Mesh& Mesh::computeThickness(Image& image, Image* dt, double max_dist, std::string distance_mesh) {
-  mesh::compute_thickness(*this, image, dt, max_dist, distance_mesh);
+Mesh& Mesh::computeThickness(Image& image, Image* dt, double max_dist, double median_radius, std::string distance_mesh) {
+  mesh::compute_thickness(*this, image, dt, max_dist, median_radius, distance_mesh);
   return *this;
 }
 

diff --git a/Libs/Mesh/Mesh.h b/Libs/Mesh/Mesh.h
@@ -182,7 +182,8 @@ class Mesh {
                             const Dims padding = Dims({1, 1, 1})) const;
 
   /// assign cortical thickness values from mesh points
-  Mesh& computeThickness(Image& image, Image* dt = nullptr, double max_dist = 10000, std::string distance_mesh = "");
+  Mesh& computeThickness(Image& image, Image* dt = nullptr, double max_dist = 10000, double median_radius = 5.0,
+                         std::string distance_mesh = "");
 
   /// compute geodesic distances to landmarks and assign as fields
   Mesh& computeLandmarkGeodesics(const std::vector<Point3>& landmarks);
@@ -275,6 +276,12 @@ class Mesh {
   //! Clips the mesh according to a field value
   vtkSmartPointer<vtkPolyData> clipByField(const std::string& name, double value);
 
+  //! Returns the cell locator
+  vtkSmartPointer<vtkStaticCellLocator> getCellLocator() const {
+    updateCellLocator();
+    return cellLocator;
+  }
+
  private:
   friend struct SharedCommandData;
   Mesh()