From 2323fec1b0a1f3f3c73f021c7e99e2e618994aea Mon Sep 17 00:00:00 2001 From: Alan Morris Date: Thu, 31 Aug 2023 17:31:50 -0600 Subject: [PATCH 1/7] Restrict distance geometrically. The thickness shouldn't be more than half the distance to the other side. --- Libs/Mesh/Mesh.h | 6 +++++ Libs/Mesh/MeshComputeThickness.cpp | 43 +++++++++++++++++++++++++++--- 2 files changed, 46 insertions(+), 3 deletions(-) diff --git a/Libs/Mesh/Mesh.h b/Libs/Mesh/Mesh.h index a6eaad42cf..ed5aa6950b 100644 --- a/Libs/Mesh/Mesh.h +++ b/Libs/Mesh/Mesh.h @@ -275,6 +275,12 @@ class Mesh { //! Clips the mesh according to a field value vtkSmartPointer clipByField(const std::string& name, double value); + //! Returns the cell locator + vtkSmartPointer getCellLocator() const { + updateCellLocator(); + return cellLocator; + } + private: friend struct SharedCommandData; Mesh() diff --git a/Libs/Mesh/MeshComputeThickness.cpp b/Libs/Mesh/MeshComputeThickness.cpp index 31a911d5e4..e077c9d167 100644 --- a/Libs/Mesh/MeshComputeThickness.cpp +++ b/Libs/Mesh/MeshComputeThickness.cpp @@ -3,6 +3,7 @@ #include #include #include +#include #include "Logging.h" namespace shapeworks::mesh { @@ -13,7 +14,7 @@ using GradientFilterType = itk::GradientImageFilter; using GradientInterpolatorType = itk::VectorLinearInterpolateImageFunction; //--------------------------------------------------------------------------- -std::vector smoothIntensities(std::vector intensities) { +std::vector smooth_intensities(std::vector intensities) { // smooth using average of neighbors std::vector smoothed_intensities; @@ -32,6 +33,41 @@ std::vector smoothIntensities(std::vector intensities) { return smoothed_intensities; } +//--------------------------------------------------------------------------- +static double get_distance_to_opposite_side(Mesh& mesh, int point_id) { + vtkSmartPointer poly_data = 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; + } + + auto cellLocator = mesh.getCellLocator(); + + // Find the intersection point with the mesh + double intersectionPoint[3]; + double t; + int subId; + double pcoords[3]; + int result = cellLocator->IntersectWithLine(ray_start, ray_end, 0.0, t, intersectionPoint, pcoords, subId); + + if (result == 0) { + // return large number if no intersection found + return std::numeric_limits::max(); + } + // Calculate the distance between the input point and the intersection point + double distance = std::sqrt(vtkMath::Distance2BetweenPoints(ray_start, intersectionPoint)); + return distance; +} + //--------------------------------------------------------------------------- void compute_thickness(Mesh& mesh, Image& image, Image* dt, double max_dist, std::string distance_mesh) { SW_DEBUG("Computing thickness with max_dist {}", max_dist); @@ -162,9 +198,9 @@ void compute_thickness(Mesh& mesh, Image& image, Image* dt, double max_dist, std point[2] += gradient[2]; } - auto smoothed = smoothIntensities(intensities); + auto smoothed = smooth_intensities(intensities); for (int k = 0; k < 10; k++) { - smoothed = smoothIntensities(smoothed); + smoothed = smooth_intensities(smoothed); } intensities = smoothed; @@ -291,6 +327,7 @@ void compute_thickness(Mesh& mesh, Image& image, Image* dt, double max_dist, std #endif distance = std::min(distance, max_dist); + distance = std::min(distance, get_distance_to_opposite_side(mesh, i) / 2.0); values->InsertValue(i, distance); } From 78ae509a41632b57543ce28268a0d405b93149e3 Mon Sep 17 00:00:00 2001 From: Alan Morris Date: Thu, 31 Aug 2023 21:37:31 -0600 Subject: [PATCH 2/7] Add geometric distance restriction, add median filtering of scalar values. --- Applications/shapeworks/MeshCommands.cpp | 11 ++- Libs/Mesh/Mesh.cpp | 4 +- Libs/Mesh/Mesh.h | 3 +- Libs/Mesh/MeshComputeThickness.cpp | 93 +++++++++++++++++++++++- Libs/Mesh/MeshComputeThickness.h | 3 +- Libs/Python/ShapeworksPython.cpp | 2 +- 6 files changed, 106 insertions(+), 10 deletions(-) diff --git a/Applications/shapeworks/MeshCommands.cpp b/Applications/shapeworks/MeshCommands.cpp index fccb346ab9..b9a8ace8ed 100644 --- 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(options.get("max_dist")); + double median_radius = static_cast(options.get("median_radius")); + std::string dt_filename = static_cast(options.get("distance_transform")); std::string distance_mesh_filename = static_cast(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 index 74b20cf6c9..f44fdc9652 100644 --- 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 index ed5aa6950b..29add84380 100644 --- 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& landmarks); diff --git a/Libs/Mesh/MeshComputeThickness.cpp b/Libs/Mesh/MeshComputeThickness.cpp index e077c9d167..b1e714a6a6 100644 --- a/Libs/Mesh/MeshComputeThickness.cpp +++ b/Libs/Mesh/MeshComputeThickness.cpp @@ -4,6 +4,7 @@ #include #include #include +#include #include "Logging.h" namespace shapeworks::mesh { @@ -14,7 +15,87 @@ using GradientFilterType = itk::GradientImageFilter; using GradientInterpolatorType = itk::VectorLinearInterpolateImageFunction; //--------------------------------------------------------------------------- -std::vector smooth_intensities(std::vector intensities) { +static double compute_average_edge_length(vtkSmartPointer poly_data) { + vtkPoints* points = poly_data->GetPoints(); + vtkCellArray* cells = poly_data->GetPolys(); + + vtkIdType num_edges = 0; + double total_edge_length = 0.0; + + vtkIdType cell_size; + vtkIdType const* cell_points; + for (cells->InitTraversal(); cells->GetNextCell(cell_size, cell_points);) { + if (cell_size == 3 || cell_size == 4) { // Assuming triangles or quads + for (vtkIdType i = 0; i < cell_size; ++i) { + vtkIdType start_point_id = cell_points[i]; + vtkIdType end_point_id = cell_points[(i + 1) % cell_size]; + + double start_point[3]; + double end_point[3]; + points->GetPoint(start_point_id, start_point); + points->GetPoint(end_point_id, end_point); + + double edge_length = vtkMath::Distance2BetweenPoints(start_point, end_point); + total_edge_length += std::sqrt(edge_length); + + num_edges++; + } + } + } + + if (num_edges > 0) { + return total_edge_length / num_edges; + } else { + return 0.0; // Return 0 if no edges found + } +} + +//--------------------------------------------------------------------------- +static void median_smooth(vtkSmartPointer poly_data, const char* scalar_name, double radius) { + vtkSmartPointer smoothed_scalars = vtkSmartPointer::New(); + smoothed_scalars->SetName(scalar_name); + + auto locator = vtkSmartPointer::New(); + locator->SetDataSet(poly_data); + locator->BuildLocator(); + + vtkSmartPointer point_ids = vtkSmartPointer::New(); + vtkSmartPointer neighbor_point_ids = vtkSmartPointer::New(); + + for (vtkIdType point_id = 0; point_id < poly_data->GetNumberOfPoints(); ++point_id) { + double queryPoint[3]; + poly_data->GetPoint(point_id, queryPoint); + + locator->FindPointsWithinRadius(radius, queryPoint, point_ids); + + int num_neighbors = point_ids->GetNumberOfIds(); + if (num_neighbors == 0) { + smoothed_scalars->InsertNextValue(poly_data->GetPointData()->GetScalars(scalar_name)->GetComponent(point_id, 0)); + continue; + } + + std::vector neighbor_scalars; + neighbor_scalars.reserve(num_neighbors); + + for (int i = 0; i < num_neighbors; ++i) { + vtkIdType neighbor_id = point_ids->GetId(i); + double neighbor_scalar = poly_data->GetPointData()->GetScalars(scalar_name)->GetComponent(neighbor_id, 0); + neighbor_scalars.push_back(neighbor_scalar); + } + + std::sort(neighbor_scalars.begin(), neighbor_scalars.end()); + + int median_index = num_neighbors / 2; + double smoothed_value = neighbor_scalars[median_index]; + + smoothed_scalars->InsertNextValue(smoothed_value); + } + + poly_data->GetPointData()->AddArray(smoothed_scalars); +} + +//--------------------------------------------------------------------------- +static std::vector smooth_intensities(std::vector intensities) { // smooth using average of neighbors std::vector smoothed_intensities; @@ -69,7 +150,8 @@ static double get_distance_to_opposite_side(Mesh& mesh, int point_id) { } //--------------------------------------------------------------------------- -void compute_thickness(Mesh& mesh, Image& image, Image* dt, double max_dist, std::string distance_mesh) { +void compute_thickness(Mesh& mesh, Image& image, Image* dt, double max_dist, double median_radius, + std::string distance_mesh) { SW_DEBUG("Computing thickness with max_dist {}", max_dist); bool use_dt = dt != nullptr; @@ -333,12 +415,17 @@ void compute_thickness(Mesh& mesh, Image& image, Image* dt, double max_dist, std } mesh.setField("thickness", values, Mesh::Point); + double average_edge_length = compute_average_edge_length(poly_data); + SW_DEBUG("average edge length: {}", average_edge_length); + + median_smooth(mesh.getVTKMesh(), "thickness", average_edge_length * median_radius); + values = vtkDoubleArray::SafeDownCast(mesh.getVTKMesh()->GetPointData()->GetArray("thickness")); + if (distance_mesh != "") { // create a copy Mesh d_mesh = mesh; // for each vertex, move the particle the distance scalar - vtkSmartPointer poly_data = d_mesh.getVTKMesh(); for (int i = 0; i < mesh.numPoints(); i++) { Point3 point; diff --git a/Libs/Mesh/MeshComputeThickness.h b/Libs/Mesh/MeshComputeThickness.h index c0f720b277..7b41255b99 100644 --- a/Libs/Mesh/MeshComputeThickness.h +++ b/Libs/Mesh/MeshComputeThickness.h @@ -6,6 +6,7 @@ namespace shapeworks::mesh { //! Compute the cortical thickness of a mesh and image (e.g. CT) -void compute_thickness(Mesh &mesh, Image &image, Image *dt, double max_dist, std::string distance_mesh); +void compute_thickness(Mesh &mesh, Image &image, Image *dt, double max_dist, double median_radius, + std::string distance_mesh); } // namespace shapeworks::mesh diff --git a/Libs/Python/ShapeworksPython.cpp b/Libs/Python/ShapeworksPython.cpp index 83d7a51a5a..1c68c6f886 100644 --- a/Libs/Python/ShapeworksPython.cpp +++ b/Libs/Python/ShapeworksPython.cpp @@ -1246,7 +1246,7 @@ PYBIND11_MODULE(shapeworks_py, m) .def("computeThickness", &Mesh::computeThickness, "Computes cortical thickness", - "ct"_a, "dt"_a = nullptr, "maxDist"_a=10000, "distanceMesh"_a = "") + "ct"_a, "dt"_a = nullptr, "maxDist"_a=10000, "medianRadius"_a=5.0, "distanceMesh"_a = "") ; From acc0469576f3201a62d2ac245da95f780341b7d2 Mon Sep 17 00:00:00 2001 From: Alan Morris Date: Thu, 7 Sep 2023 12:55:14 -0600 Subject: [PATCH 3/7] Add median smoothing of 1D signal profiles --- Libs/Mesh/MeshComputeThickness.cpp | 24 ++++++++++++++++++++++-- 1 file changed, 22 insertions(+), 2 deletions(-) diff --git a/Libs/Mesh/MeshComputeThickness.cpp b/Libs/Mesh/MeshComputeThickness.cpp index b1e714a6a6..9816aed86a 100644 --- a/Libs/Mesh/MeshComputeThickness.cpp +++ b/Libs/Mesh/MeshComputeThickness.cpp @@ -114,6 +114,26 @@ static std::vector smooth_intensities(std::vector intensities) { return smoothed_intensities; } +//--------------------------------------------------------------------------- +static std::vector median_smooth_signal_intensities(std::vector intensities) { + // smooth using average of neighbors + std::vector smoothed_intensities; + + for (int i = 0; i < intensities.size(); i++) { + double sum = 0; + int count = 0; + std::vector local_intensities; + for (int j = -2; j <= 2; j++) { + local_intensities.push_back(intensities[i + j]); + } + // compute median + std::sort(local_intensities.begin(), local_intensities.end()); + smoothed_intensities.push_back(local_intensities[2]); + } + + return smoothed_intensities; +} + //--------------------------------------------------------------------------- static double get_distance_to_opposite_side(Mesh& mesh, int point_id) { vtkSmartPointer poly_data = mesh.getVTKMesh(); @@ -280,9 +300,9 @@ void compute_thickness(Mesh& mesh, Image& image, Image* dt, double max_dist, dou point[2] += gradient[2]; } - auto smoothed = smooth_intensities(intensities); + auto smoothed = median_smooth_signal_intensities(intensities); for (int k = 0; k < 10; k++) { - smoothed = smooth_intensities(smoothed); + smoothed = median_smooth_signal_intensities(smoothed); } intensities = smoothed; From 764daf8dcfac55acc9f8a4fd660c1cf7f6b3ca8a Mon Sep 17 00:00:00 2001 From: Alan Morris Date: Fri, 8 Sep 2023 10:15:29 -0600 Subject: [PATCH 4/7] Calibrate outside values using average across surface. This fixes the problem of nearby bones. --- Libs/Mesh/MeshComputeThickness.cpp | 89 +++++++++++++++++++----------- 1 file changed, 58 insertions(+), 31 deletions(-) diff --git a/Libs/Mesh/MeshComputeThickness.cpp b/Libs/Mesh/MeshComputeThickness.cpp index 9816aed86a..5208e4da16 100644 --- a/Libs/Mesh/MeshComputeThickness.cpp +++ b/Libs/Mesh/MeshComputeThickness.cpp @@ -229,38 +229,78 @@ void compute_thickness(Mesh& mesh, Image& image, Image* dt, double max_dist, dou double step_size = spacing / 10.0; + auto step = [&](Point3 point, VectorPixelType gradient, double multiplier) -> Point3 { + // normalize the gradient + float norm = std::sqrt(gradient[0] * gradient[0] + gradient[1] * gradient[1] + gradient[2] * gradient[2]); + Point3 new_point; + for (int i = 0; i < 3; i++) { + gradient[i] /= norm; + gradient[i] *= step_size; + new_point[i] = point[i] + gradient[i] * multiplier; + } + return new_point; + }; + + const double distance_outside = 4.0; + const double distance_inside = 12.0; + + // first establish average 'outside' values + + std::vector average_intensities; + std::vector all_intensities; for (int i = 0; i < mesh.numPoints(); i++) { Point3 point; poly_data->GetPoint(i, point.GetDataPointer()); - std::vector intensities; - double distance_outside = 4.0; - double distance_inside = 12.0; - for (int j = 0; j <= distance_outside / step_size; j++) { double intensity = 0; // check if point is inside image if (check_inside(point)) { intensity = image.evaluate(point); + all_intensities.push_back(intensity); + } + if (j < average_intensities.size()) { + average_intensities[j] += intensity; + } else { + average_intensities.push_back(intensity); } - intensities.push_back(intensity); VectorPixelType gradient = get_gradient(i, point); + point = step(point, gradient, -1.0); + } + } - // normalize the gradient - float norm = std::sqrt(gradient[0] * gradient[0] + gradient[1] * gradient[1] + gradient[2] * gradient[2]); - gradient[0] /= norm; - gradient[1] /= norm; - gradient[2] /= norm; + // divide by mean to compute average + for (int i = 0; i < average_intensities.size(); i++) { + average_intensities[i] /= mesh.numPoints(); + } - gradient[0] *= step_size; - gradient[1] *= step_size; - gradient[2] *= step_size; + // compute median of all intensities + std::sort(all_intensities.begin(), all_intensities.end()); + double median_intensity = all_intensities[all_intensities.size() / 2]; - // take step - point[0] -= gradient[0]; - point[1] -= gradient[1]; - point[2] -= gradient[2]; + // parallel loop over all points using TBB + + for (int i = 0; i < mesh.numPoints(); i++) { + Point3 point; + poly_data->GetPoint(i, point.GetDataPointer()); + + std::vector intensities; + + for (int j = 0; j <= distance_outside / step_size; j++) { + double intensity = 0; + // check if point is inside image + if (check_inside(point)) { + intensity = image.evaluate(point); + } + intensities.push_back(average_intensities[j]); + // intensities.push_back(median_intensity); + // intensities.push_back(intensity); + //SW_LOG("median: {}, average: {}", median_intensity, average_intensities[j]); + //std::cerr << "median: " << median_intensity << ", average: " << average_intensities[j] << "\n"; + + VectorPixelType gradient = get_gradient(i, point); + point = step(point, gradient, -1.0); } // reverse the intensities vector @@ -284,20 +324,7 @@ void compute_thickness(Mesh& mesh, Image& image, Image* dt, double max_dist, dou // evaluate gradient VectorPixelType gradient = get_gradient(i, point); - // normalize the gradient - float norm = std::sqrt(gradient[0] * gradient[0] + gradient[1] * gradient[1] + gradient[2] * gradient[2]); - gradient[0] /= norm; - gradient[1] /= norm; - gradient[2] /= norm; - - gradient[0] *= step_size; - gradient[1] *= step_size; - gradient[2] *= step_size; - - // take step - point[0] += gradient[0]; - point[1] += gradient[1]; - point[2] += gradient[2]; + point = step(point, gradient, 1.0); } auto smoothed = median_smooth_signal_intensities(intensities); From 93464a5648aa4fe03808580c8053afe9dcce6362 Mon Sep 17 00:00:00 2001 From: Alan Morris Date: Fri, 8 Sep 2023 10:38:08 -0600 Subject: [PATCH 5/7] Parallelize mesh compute thickness --- Libs/Mesh/MeshComputeThickness.cpp | 338 ++++++++++++++--------------- 1 file changed, 168 insertions(+), 170 deletions(-) diff --git a/Libs/Mesh/MeshComputeThickness.cpp b/Libs/Mesh/MeshComputeThickness.cpp index 5208e4da16..e4652c4f4b 100644 --- a/Libs/Mesh/MeshComputeThickness.cpp +++ b/Libs/Mesh/MeshComputeThickness.cpp @@ -2,6 +2,7 @@ #include #include +#include #include #include #include @@ -120,8 +121,6 @@ static std::vector median_smooth_signal_intensities(std::vector std::vector smoothed_intensities; for (int i = 0; i < intensities.size(); i++) { - double sum = 0; - int count = 0; std::vector local_intensities; for (int j = -2; j <= 2; j++) { local_intensities.push_back(intensities[i + j]); @@ -169,6 +168,138 @@ static double get_distance_to_opposite_side(Mesh& mesh, int point_id) { return distance; } +double compute_thickness_from_signal(const std::vector& intensities_input, double step_size) { + auto smoothed = median_smooth_signal_intensities(intensities_input); + for (int k = 0; k < 10; k++) { + smoothed = median_smooth_signal_intensities(smoothed); + } + auto intensities = smoothed; + + // find the highest intensity and its index + double max_intensity = 0; + int point_a = 0; + for (int j = 0; j < intensities.size(); j++) { + if (intensities[j] > max_intensity) { + max_intensity = intensities[j]; + point_a = j; + } + } + + // compute derivative of intensities + std::vector derivatives; + for (int j = 0; j < intensities.size() - 1; j++) { + derivatives.push_back((intensities[j + 1] - intensities[j]) / step_size); + } + + // compute mean of DHU (derivative houndsfield unit) of the first 2mm (point x) + int count_first_2mm = (1 / step_size) * 2.0; + int point_x = count_first_2mm; + int point_o = (1 / step_size) * 3.0; + + double mean_dhu = 0; + for (int j = 0; j < count_first_2mm; j++) { + mean_dhu += derivatives[j]; + } + mean_dhu /= count_first_2mm; + + // compute standard deviation of DHU of the first 2mm + double std_dhu = 0; + for (int j = 0; j < count_first_2mm; j++) { + std_dhu += (derivatives[j] - mean_dhu) * (derivatives[j] - mean_dhu); + } + + std_dhu = std::sqrt(std_dhu / count_first_2mm); + + // equation 2 + double dhu_threshold = mean_dhu + 4.27 * std_dhu; + + // find the first index where the derivative is greater than the threshold (point b) + int point_b = -1; + for (int j = point_o; j < derivatives.size(); j++) { + if (derivatives[j] > dhu_threshold) { + point_b = j; + break; + } + } + + bool not_found = false; + if (point_b < 0) { + not_found = true; + point_b = 0; + } + double hu_a = intensities[point_a]; + double hu_c = intensities[point_b]; + + // constant from paper + const double k_cor = 0.605; + + double hu_threshold = (hu_a - hu_c) * k_cor + hu_c; + + // find the first and last point above hu_threshold + int point_d = -1; + for (int j = point_b; j < intensities.size(); j++) { + if (intensities[j] > hu_threshold) { + point_d = j; + break; + } + } + int point_e = -1; + if (point_d != -1) { + // find the next point below hu_threshold + for (int j = point_d; j < intensities.size(); j++) { + if (intensities[j] < hu_threshold) { + point_e = j; + break; + } + } + } + + // compute the distance between the first and last point above the threshold + double distance = (point_e - point_d) * step_size; + + if (point_d == -1 || point_e == -1 || not_found) { + distance = 0; + } + +#if 0 // debug + if (i < 10) { + // write out the intensities to a file named with the point id + std::ofstream out("intensities_" + std::to_string(i) + ".txt"); + for (auto intensity : intensities) { + out << intensity << std::endl; + } + std::ofstream out2("derivatives_" + std::to_string(i) + ".txt"); + for (auto derivative : derivatives) { + out2 << derivative << std::endl; + } + + std::ofstream out4("smoothed_" + std::to_string(i) + ".txt"); + for (auto intensity : smoothed) { + out4 << intensity << std::endl; + } + + // write out json file with the parameters + std::ofstream out3("parameters_" + std::to_string(i) + ".json"); + out3 << "{\n"; + out3 << " \"hu_a\": " << hu_a << ",\n"; + out3 << " \"hu_c\": " << hu_c << ",\n"; + out3 << " \"k_cor\": " << k_cor << ",\n"; + out3 << " \"hu_threshold\": " << hu_threshold << ",\n"; + out3 << " \"distance\": " << distance << ",\n"; + out3 << " \"point_a\": " << point_a << ",\n"; + out3 << " \"point_b\": " << point_b << ",\n"; + out3 << " \"point_d\": " << point_d << ",\n"; + out3 << " \"point_e\": " << point_e << ",\n"; + out3 << " \"mean_dhu\": " << mean_dhu << ",\n"; + out3 << " \"std_dhu\": " << std_dhu << ",\n"; + out3 << " \"dhu_threshold\": " << dhu_threshold << "\n"; + out3 << "}\n"; + } +#endif + + return distance; +} + //--------------------------------------------------------------------------- void compute_thickness(Mesh& mesh, Image& image, Image* dt, double max_dist, double median_radius, std::string distance_mesh) { @@ -177,6 +308,7 @@ void compute_thickness(Mesh& mesh, Image& image, Image* dt, double max_dist, dou bool use_dt = dt != nullptr; mesh.computeNormals(); + mesh.getCellLocator(); auto poly_data = mesh.getVTKMesh(); @@ -247,7 +379,6 @@ void compute_thickness(Mesh& mesh, Image& image, Image* dt, double max_dist, dou // first establish average 'outside' values std::vector average_intensities; - std::vector all_intensities; for (int i = 0; i < mesh.numPoints(); i++) { Point3 point; poly_data->GetPoint(i, point.GetDataPointer()); @@ -257,7 +388,6 @@ void compute_thickness(Mesh& mesh, Image& image, Image* dt, double max_dist, dou // check if point is inside image if (check_inside(point)) { intensity = image.evaluate(point); - all_intensities.push_back(intensity); } if (j < average_intensities.size()) { average_intensities[j] += intensity; @@ -275,191 +405,59 @@ void compute_thickness(Mesh& mesh, Image& image, Image* dt, double max_dist, dou average_intensities[i] /= mesh.numPoints(); } - // compute median of all intensities - std::sort(all_intensities.begin(), all_intensities.end()); - double median_intensity = all_intensities[all_intensities.size() / 2]; + // mutex + std::mutex mutex; + unsigned long num_points = mesh.numPoints(); // parallel loop over all points using TBB + tbb::parallel_for(tbb::blocked_range{0, num_points}, [&](const tbb::blocked_range& r) { + for (size_t i = r.begin(); i < r.end(); ++i) { + // for (int i = 0; i < mesh.numPoints(); i++) { + Point3 point; + poly_data->GetPoint(i, point.GetDataPointer()); - for (int i = 0; i < mesh.numPoints(); i++) { - Point3 point; - poly_data->GetPoint(i, point.GetDataPointer()); - - std::vector intensities; - - for (int j = 0; j <= distance_outside / step_size; j++) { - double intensity = 0; - // check if point is inside image - if (check_inside(point)) { - intensity = image.evaluate(point); - } - intensities.push_back(average_intensities[j]); - // intensities.push_back(median_intensity); - // intensities.push_back(intensity); - //SW_LOG("median: {}, average: {}", median_intensity, average_intensities[j]); - //std::cerr << "median: " << median_intensity << ", average: " << average_intensities[j] << "\n"; - - VectorPixelType gradient = get_gradient(i, point); - point = step(point, gradient, -1.0); - } - - // reverse the intensities vector - std::reverse(intensities.begin(), intensities.end()); - - // drop the last intensity, since we're already there - intensities.pop_back(); - - // reset point position back to the surface - poly_data->GetPoint(i, point.GetDataPointer()); - - for (int j = 0; j < distance_inside / step_size; j++) { - double intensity = 0; - // check if point is inside image - if (check_inside(point)) { - intensity = image.evaluate(point); - } - - intensities.push_back(intensity); - - // evaluate gradient - VectorPixelType gradient = get_gradient(i, point); - - point = step(point, gradient, 1.0); - } + std::vector intensities; - auto smoothed = median_smooth_signal_intensities(intensities); - for (int k = 0; k < 10; k++) { - smoothed = median_smooth_signal_intensities(smoothed); - } - intensities = smoothed; - - // find the highest intensity and its index - double max_intensity = 0; - int point_a = 0; - for (int j = 0; j < intensities.size(); j++) { - if (intensities[j] > max_intensity) { - max_intensity = intensities[j]; - point_a = j; + for (int j = 0; j <= distance_outside / step_size; j++) { + intensities.push_back(average_intensities[j]); } - } - - // compute derivative of intensities - std::vector derivatives; - for (int j = 0; j < intensities.size() - 1; j++) { - derivatives.push_back((intensities[j + 1] - intensities[j]) / step_size); - } - - // compute mean of DHU (derivative houndsfield unit) of the first 2mm (point x) - int count_first_2mm = (1 / step_size) * 2.0; - int point_x = count_first_2mm; - int point_o = (1 / step_size) * 3.0; - - double mean_dhu = 0; - for (int j = 0; j < count_first_2mm; j++) { - mean_dhu += derivatives[j]; - } - mean_dhu /= count_first_2mm; - - // compute standard deviation of DHU of the first 2mm - double std_dhu = 0; - for (int j = 0; j < count_first_2mm; j++) { - std_dhu += (derivatives[j] - mean_dhu) * (derivatives[j] - mean_dhu); - } - std_dhu = std::sqrt(std_dhu / count_first_2mm); + // reverse the intensities vector + std::reverse(intensities.begin(), intensities.end()); - // equation 2 - double dhu_threshold = mean_dhu + 4.27 * std_dhu; + // drop the last intensity, since we're already there + intensities.pop_back(); - // find the first index where the derivative is greater than the threshold (point b) - int point_b = -1; - for (int j = point_o; j < derivatives.size(); j++) { - if (derivatives[j] > dhu_threshold) { - point_b = j; - break; - } - } + // reset point position back to the surface + poly_data->GetPoint(i, point.GetDataPointer()); - bool not_found = false; - if (point_b < 0) { - not_found = true; - point_b = 0; - } - double hu_a = intensities[point_a]; - double hu_c = intensities[point_b]; + for (int j = 0; j < distance_inside / step_size; j++) { + double intensity = 0; + // check if point is inside image + if (check_inside(point)) { + intensity = image.evaluate(point); + } - // constant from paper - const double k_cor = 0.605; + intensities.push_back(intensity); - double hu_threshold = (hu_a - hu_c) * k_cor + hu_c; + // evaluate gradient + VectorPixelType gradient = get_gradient(i, point); - // find the first and last point above hu_threshold - int point_d = -1; - for (int j = point_b; j < intensities.size(); j++) { - if (intensities[j] > hu_threshold) { - point_d = j; - break; - } - } - int point_e = -1; - if (point_d != -1) { - // find the next point below hu_threshold - for (int j = point_d; j < intensities.size(); j++) { - if (intensities[j] < hu_threshold) { - point_e = j; - break; - } + point = step(point, gradient, 1.0); } - } - - // compute the distance between the first and last point above the threshold - double distance = (point_e - point_d) * step_size; - if (point_d == -1 || point_e == -1 || not_found) { - distance = 0; - } + auto distance = compute_thickness_from_signal(intensities, step_size); -#if 0 // debug - if (i < 10) { - // write out the intensities to a file named with the point id - std::ofstream out("intensities_" + std::to_string(i) + ".txt"); - for (auto intensity : intensities) { - out << intensity << std::endl; - } - std::ofstream out2("derivatives_" + std::to_string(i) + ".txt"); - for (auto derivative : derivatives) { - out2 << derivative << std::endl; - } + distance = std::min(distance, max_dist); + distance = std::min(distance, get_distance_to_opposite_side(mesh, i) / 2.0); - std::ofstream out4("smoothed_" + std::to_string(i) + ".txt"); - for (auto intensity : smoothed) { - out4 << intensity << std::endl; + { + std::lock_guard lock(mutex); + values->InsertValue(i, distance); } - - // write out json file with the parameters - std::ofstream out3("parameters_" + std::to_string(i) + ".json"); - out3 << "{\n"; - out3 << " \"hu_a\": " << hu_a << ",\n"; - out3 << " \"hu_c\": " << hu_c << ",\n"; - out3 << " \"k_cor\": " << k_cor << ",\n"; - out3 << " \"hu_threshold\": " << hu_threshold << ",\n"; - out3 << " \"distance\": " << distance << ",\n"; - out3 << " \"point_a\": " << point_a << ",\n"; - out3 << " \"point_b\": " << point_b << ",\n"; - out3 << " \"point_d\": " << point_d << ",\n"; - out3 << " \"point_e\": " << point_e << ",\n"; - out3 << " \"mean_dhu\": " << mean_dhu << ",\n"; - out3 << " \"std_dhu\": " << std_dhu << ",\n"; - out3 << " \"dhu_threshold\": " << dhu_threshold << "\n"; - out3 << "}\n"; } -#endif + }); - distance = std::min(distance, max_dist); - distance = std::min(distance, get_distance_to_opposite_side(mesh, i) / 2.0); - - values->InsertValue(i, distance); - } mesh.setField("thickness", values, Mesh::Point); double average_edge_length = compute_average_edge_length(poly_data); From c49204c690f2942d8916711ea7cbe3013d2a9c42 Mon Sep 17 00:00:00 2001 From: Alan Morris Date: Mon, 25 Sep 2023 13:52:11 -0600 Subject: [PATCH 6/7] Update baseline --- Testing/data/thickness/thickness.vtk | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Testing/data/thickness/thickness.vtk b/Testing/data/thickness/thickness.vtk index 55b5b63c14..6ca5e15a37 100644 --- a/Testing/data/thickness/thickness.vtk +++ b/Testing/data/thickness/thickness.vtk @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:3e9a4c205e2d8099c3ad0c5dc0bff2e21878abd1eb3e24d74c950f2560539bda -size 687380 +oid sha256:a38d49a33baf9186fc640e216a316f85e3dd180ff92d00ff6bb4a17c3683afca +size 687238 From e77b19447912cc8ee837a39c6abdde35db671c3f Mon Sep 17 00:00:00 2001 From: Alan Morris Date: Tue, 26 Sep 2023 14:24:10 -0600 Subject: [PATCH 7/7] Lock guard around non-threadsafe cell locator IntersectWithLine --- Libs/Mesh/MeshComputeThickness.cpp | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/Libs/Mesh/MeshComputeThickness.cpp b/Libs/Mesh/MeshComputeThickness.cpp index e4652c4f4b..66dc605f6d 100644 --- a/Libs/Mesh/MeshComputeThickness.cpp +++ b/Libs/Mesh/MeshComputeThickness.cpp @@ -15,6 +15,9 @@ using GradientImageType = itk::Image; using GradientFilterType = itk::GradientImageFilter; using GradientInterpolatorType = itk::VectorLinearInterpolateImageFunction; +// mutex for cell locator (not thread safe) +static std::mutex cell_mutex; + //--------------------------------------------------------------------------- static double compute_average_edge_length(vtkSmartPointer poly_data) { vtkPoints* points = poly_data->GetPoints(); @@ -150,15 +153,18 @@ static double get_distance_to_opposite_side(Mesh& mesh, int point_id) { ray_end[i] = ray_start[i] - normal[i] * 100.0; } - auto cellLocator = mesh.getCellLocator(); - + int result = 0; // Find the intersection point with the mesh double intersectionPoint[3]; double t; int subId; double pcoords[3]; - int result = cellLocator->IntersectWithLine(ray_start, ray_end, 0.0, t, intersectionPoint, pcoords, subId); - + { + // lock mutex (IntersectWithLine is not thread safe) + std::lock_guard lock(cell_mutex); + auto cellLocator = mesh.getCellLocator(); + result = cellLocator->IntersectWithLine(ray_start, ray_end, 0.0, t, intersectionPoint, pcoords, subId); + } if (result == 0) { // return large number if no intersection found return std::numeric_limits::max(); @@ -409,10 +415,10 @@ void compute_thickness(Mesh& mesh, Image& image, Image* dt, double max_dist, dou std::mutex mutex; unsigned long num_points = mesh.numPoints(); + // parallel loop over all points using TBB tbb::parallel_for(tbb::blocked_range{0, num_points}, [&](const tbb::blocked_range& r) { for (size_t i = r.begin(); i < r.end(); ++i) { - // for (int i = 0; i < mesh.numPoints(); i++) { Point3 point; poly_data->GetPoint(i, point.GetDataPointer());