#2189 internal image summary (#2193)

* toDistanceTransform needs to invalidate the cell locator in case it's not properly reset.

* Add Mesh isPointInside (currently unused)

* Add helpers to return minimum image spacing.

* Add internal image summary information fo mesh compute thickness functionality (#2189)

* Catch Image exception

* Check DT by spacing

* Cleanup, add smoothing

* Update baseline

* Fix bug causing undeterministic behavior

* Check inside image

* Update baseline

Author: akenmorris

Libs/Analyze/Shape.cpp

@@ -568,9 +568,13 @@ std::shared_ptr<Image> Shape::get_image_volume(std::string image_volume_name) {
     auto filename = filenames[image_volume_name];
 
     if (image_volume_filename_ != filename) {
-      std::shared_ptr<Image> image = std::make_shared<Image>(filename);
-      image_volume_ = image;
-      image_volume_filename_ = filename;
+      try {
+        std::shared_ptr<Image> image = std::make_shared<Image>(filename);
+        image_volume_ = image;
+        image_volume_filename_ = filename;
+      } catch (std::exception &ex) {
+        SW_ERROR("Unable to open file: {}", filename);
+      }
     }
 
     return image_volume_;

Libs/Image/Image.cpp

@@ -859,6 +859,11 @@ Image& Image::isolate() {
   return *this;
 }
 
+double Image::get_minimum_spacing() const {
+  auto spacing = this->spacing();
+  return std::min(spacing[0], std::min(spacing[1], spacing[2]));
+}
+
 Point3 Image::centerOfMass(PixelType minVal, PixelType maxVal) const {
   itk::ImageRegionIteratorWithIndex<ImageType> imageIt(this->itk_image_, itk_image_->GetLargestPossibleRegion());
   int numPixels = 0;

Libs/Image/Image.h

@@ -216,6 +216,9 @@ class Image {
   /// physical spacing of the image
   Vector spacing() const { return itk_image_->GetSpacing(); }
 
+  /// minimum physical spacing of the image
+  double get_minimum_spacing() const;
+
   /// physical coordinates of image origin
   Point3 origin() const { return itk_image_->GetOrigin(); }
 

Libs/Mesh/Mesh.cpp

@@ -763,6 +763,23 @@ int Mesh::closestPointId(const Point3 point) const {
   return closestPointId;
 }
 
+bool Mesh::isPointInside(const Point3 point) const {
+  // create point set
+  auto points = vtkSmartPointer<vtkPoints>::New();
+  points->InsertNextPoint(point.GetDataPointer());
+  auto polydata = vtkSmartPointer<vtkPolyData>::New();
+  polydata->SetPoints(points);
+
+  auto select = vtkSmartPointer<vtkSelectEnclosedPoints>::New();
+
+  select->SetInputData(polydata);
+  select->SetSurfaceData(this->poly_data_);
+  select->SetTolerance(0.0001);
+  select->Update();
+
+  return select->IsInside(0);
+}
+
 double Mesh::geodesicDistance(int source, int target) const {
   if (source < 0 || target < 0 || numPoints() < source || numPoints() < target) {
     throw std::invalid_argument("requested point ids outside range of points available in mesh");
@@ -1066,6 +1083,7 @@ Image Mesh::toImage(PhysicalRegion region, Point3 spacing) const {
 }
 
 Image Mesh::toDistanceTransform(PhysicalRegion region, const Point3 spacing, const Dims padding) const {
+  invalidateLocators();
   this->updateCellLocator();
 
   // if no region, use mesh bounding box
@@ -1106,14 +1124,20 @@ Image Mesh::toDistanceTransform(PhysicalRegion region, const Point3 spacing, con
   enclosed->SetSurfaceData(this->poly_data_);
   enclosed->Update();
 
+  std::mutex cell_mutex;
+
   tbb::parallel_for(tbb::blocked_range<size_t>(0, indices.size()), [&](const tbb::blocked_range<size_t>& r) {
     for (size_t i = r.begin(); i != r.end(); ++i) {
       Image::ImageType::PointType p;
       itkimg->TransformIndexToPhysicalPoint(indices[i], p);
 
       double distance = 0.0;
       vtkIdType face_id = 0;
-      closestPoint(p, distance, face_id);
+
+      {
+        // std::lock_guard<std::mutex> lock(cell_mutex);
+        closestPoint(p, distance, face_id);
+      }
 
       bool outside = !enclosed->IsInside(i);
 
@@ -1125,7 +1149,8 @@ Image Mesh::toDistanceTransform(PhysicalRegion region, const Point3 spacing, con
   return img;
 }
 
-Mesh& Mesh::computeThickness(Image& image, Image* dt, double max_dist, double median_radius, std::string 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;
 }

Libs/Mesh/Mesh.h

@@ -154,6 +154,9 @@ class Mesh {
   /// returns closest point id in this mesh to the given point in space
   int closestPointId(const Point3 point) const;
 
+  /// returns if the given point is inside the mesh
+  bool isPointInside(const Point3 point) const;
+
   /// computes geodesic distance between two vertices (specified by their indices) on mesh
   double geodesicDistance(int source, int target) const;
 

Libs/Mesh/MeshComputeThickness.cpp

@@ -4,6 +4,7 @@
 #include <itkVectorLinearInterpolateImageFunction.h>
 #include <tbb/parallel_for.h>
 #include <vtkPointData.h>
+#include <vtkSelectEnclosedPoints.h>
 #include <vtkStaticCellLocator.h>
 #include <vtkStaticPointLocator.h>
 
@@ -126,7 +127,9 @@ static std::vector<double> median_smooth_signal_intensities(std::vector<double>
   for (int i = 0; i < intensities.size(); i++) {
     std::vector<double> local_intensities;
     for (int j = -2; j <= 2; j++) {
-      local_intensities.push_back(intensities[i + j]);
+      if (i + j >= 0 && i + j < intensities.size()) {
+        local_intensities.push_back(intensities[i + j]);
+      }
     }
     // compute median
     std::sort(local_intensities.begin(), local_intensities.end());
@@ -174,6 +177,58 @@ static double get_distance_to_opposite_side(Mesh& mesh, int point_id) {
   return distance;
 }
 
+/*
+
+//---------------------------------------------------------------------------
+static void get_distances_to_inner_mesh(Mesh &outer_mesh, int point_id, Mesh &inner_mesh, double &p1, double &p2) {
+// using the surface normal from the points on the outer mesh, find the two intersections of the inner mesh
+
+    vtkSmartPointer<vtkPolyData> poly_data = outer_mesh.getVTKMesh();
+
+    // Get the surface normal from the given point
+    double normal[3];
+    poly_data->GetPointData()->GetNormals()->GetTuple(point_id, normal);
+
+    // Create a ray in the opposite direction of the normal
+    double ray_start[3];
+    double ray_end[3];
+    poly_data->GetPoint(point_id, ray_start);
+    for (int i = 0; i < 3; ++i) {
+        ray_start[i] = ray_start[i] - normal[i] * 0.001;
+        ray_end[i] = ray_start[i] - normal[i] * 100.0;
+    }
+
+    int result = 0;
+    // Find the intersection point with the mesh
+
+
+
+
+}
+
+*/
+
+/*
+//---------------------------------------------------------------------------
+bool is_inside(Mesh &mesh, Point3 start, Point3 current) {
+
+    // use the sign of t to determine if the point is inside or outside the mesh
+
+    // lock mutex (IntersectWithLine is not thread safe)
+    std::lock_guard<std::mutex> lock(cell_mutex);
+    auto cellLocator = mesh.getCellLocator();
+    double t;
+    double intersectionPoint[3];
+    double pcoords[3];
+    int subId;
+    int result = cellLocator->IntersectWithLine(start.GetDataPointer(), current.GetDataPointer(), 0.0, t,
+intersectionPoint, pcoords, subId);
+
+    return result != 0/;
+
+}
+*/
+
 //---------------------------------------------------------------------------
 double compute_thickness_from_signal(const std::vector<double>& intensities_input, double step_size) {
   auto smoothed = median_smooth_signal_intensities(intensities_input);
@@ -358,14 +413,7 @@ void compute_thickness(Mesh& mesh, Image& image, Image* dt, double max_dist, dou
     }
   };
 
-  // get minimum spacing
-  double spacing = image.spacing()[0];
-  for (int i = 1; i < 3; i++) {
-    if (image.spacing()[i] < spacing) {
-      spacing = image.spacing()[i];
-    }
-  }
-
+  double spacing = image.get_minimum_spacing();
   double step_size = spacing / 10.0;
 
   auto step = [&](Point3 point, VectorPixelType gradient, double multiplier) -> Point3 {
@@ -518,6 +566,10 @@ void compute_thickness(Mesh& mesh, Image& image, Image* dt, double max_dist, dou
 
     d_mesh.write(distance_mesh);
   }
+
+  // compute inner mesh
+  Mesh inner_mesh = compute_inner_mesh(mesh, "thickness");
+  summarize_internal_intensities(mesh, inner_mesh, image);
 }
 
 //---------------------------------------------------------------------------
@@ -564,4 +616,106 @@ Mesh compute_inner_mesh(const Mesh& mesh, std::string array_name) {
   return new_mesh;
 }
 
+//---------------------------------------------------------------------------
+void summarize_internal_intensities(Mesh& outer_mesh, Mesh& inner_mesh, Image& image) {
+  outer_mesh.computeNormals();
+  inner_mesh.computeNormals();
+
+  vtkSmartPointer<vtkPolyData> outerMeshWithNormals = outer_mesh.getVTKMesh();
+  vtkSmartPointer<vtkPolyData> innerMesh = inner_mesh.getVTKMesh();
+
+  // Get the normal array
+  vtkDataArray* normalArray = outerMeshWithNormals->GetPointData()->GetNormals();
+
+  double spacing = image.get_minimum_spacing();
+  double step_size = spacing / 10.0;
+
+  const unsigned long num_points = outerMeshWithNormals->GetNumberOfPoints();
+
+  // create arrays for min, max, mean
+  auto values_min = vtkSmartPointer<vtkDoubleArray>::New();
+  values_min->SetNumberOfComponents(1);
+  values_min->SetNumberOfTuples(num_points);
+  values_min->SetName("intensity_min");
+
+  auto values_max = vtkSmartPointer<vtkDoubleArray>::New();
+  values_max->SetNumberOfComponents(1);
+  values_max->SetNumberOfTuples(num_points);
+  values_max->SetName("intensity_max");
+
+  auto values_mean = vtkSmartPointer<vtkDoubleArray>::New();
+  values_mean->SetNumberOfComponents(1);
+  values_mean->SetNumberOfTuples(num_points);
+  values_mean->SetName("intensity_mean");
+
+  auto inner_dt = inner_mesh.toDistanceTransform(PhysicalRegion(), image.spacing(), {1, 1, 1});
+  auto outer_dt = outer_mesh.toDistanceTransform(PhysicalRegion(), image.spacing(), {1, 1, 1});
+
+  // Iterate over each point in the outer mesh
+  for (vtkIdType pointId = 0; pointId < outerMeshWithNormals->GetNumberOfPoints(); ++pointId) {
+    // parallel loop over all points using TBB
+    ////tbb::parallel_for(tbb::blocked_range<size_t>{0, num_points}, [&](const tbb::blocked_range<size_t>& r) {
+    ////for (size_t i = r.begin(); i < r.end(); ++i) {
+    // int pointId = i;
+
+    // Get the surface normal at this point
+    double normal[3];
+    normalArray->GetTuple(pointId, normal);
+
+    // Get the coordinates of the current point
+    double point[3];
+    outerMeshWithNormals->GetPoint(pointId, point);
+
+    // Move to the next voxel along the surface normal by the step_size
+    for (int i = 0; i < 3; ++i) {
+      point[i] -= normal[i] * step_size * 10;
+    }
+
+    double max_value = 0;
+    double min_value = std::numeric_limits<double>::max();
+    double sum = 0;
+    double count = 0;
+
+    while (outer_dt.evaluate(point) > 0) {
+        if (inner_dt.isInside(point) && inner_dt.evaluate(point) > spacing && image.isInside(point)) {
+        double value = image.evaluate(point);
+        max_value = std::max<double>(max_value, value);
+        min_value = std::min<double>(min_value, value);
+        sum += value;
+        count++;
+      }
+
+      // Move to the next voxel along the surface normal by the step_size
+      for (int i = 0; i < 3; ++i) {
+        point[i] -= normal[i] * step_size;
+      }
+    }
+
+    if (count == 0) {
+      min_value = max_value = sum = 0;  // no points found
+      count = 1;
+    }
+
+    // Set the values in the output arrays
+    values_min->SetValue(pointId, min_value);
+    values_max->SetValue(pointId, max_value);
+    values_mean->SetValue(pointId, sum / count);
+
+    // std::cout << "point " << pointId << " min " << min_value << " max " << max_value << " mean " << sum / count
+    //          << std::endl;
+  }
+  //});
+
+  outer_mesh.setField("intensity_min", values_min, Mesh::Point);
+  outer_mesh.setField("intensity_max", values_max, Mesh::Point);
+  outer_mesh.setField("intensity_mean", values_mean, Mesh::Point);
+
+  double average_edge_length = compute_average_edge_length(outer_mesh.getVTKMesh());
+
+  double median_radius = 5.0;
+  median_smooth(outer_mesh.getVTKMesh(), "intensity_min", average_edge_length * median_radius);
+  median_smooth(outer_mesh.getVTKMesh(), "intensity_max", average_edge_length * median_radius);
+  median_smooth(outer_mesh.getVTKMesh(), "intensity_mean", average_edge_length * median_radius);
+}
+
 }  // namespace shapeworks::mesh

Libs/Mesh/MeshComputeThickness.h

@@ -9,6 +9,10 @@ namespace shapeworks::mesh {
 void compute_thickness(Mesh &mesh, Image &image, Image *dt, double max_dist, double median_radius,
                        std::string distance_mesh);
 
+//! Compute an internal mesh based on the thickness values
 Mesh compute_inner_mesh(const Mesh &mesh, std::string array_name);
 
+//! Summarize the internal intensities values of the area inside the inner mesh on the outer mesh
+void summarize_internal_intensities(Mesh &outer_mesh, Mesh &inner_mesh, Image &image);
+
 }  // namespace shapeworks::mesh

Testing/data/thickness/thickness.vtk

@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:a38d49a33baf9186fc640e216a316f85e3dd180ff92d00ff6bb4a17c3683afca
-size 687238
+oid sha256:c98a3f196ecd875435294b2323d96d7624e6da879fe48db2a211db57430564b8
+size 743076