trying a code refactoring

LocalizationService.cpp

@@ -0,0 +1,238 @@
+#include "LocalizationService.h"
+
+#include <vtkCellData.h>
+#include <vtkExtractEdges.h>
+#include <vtkGeometryFilter.h>
+#include <vtkImageData.h>
+#include <vtkMarchingCubes.h>
+#include <vtkNew.h>
+#include <vtkPoints.h>
+#include <vtkPolyData.h>
+#include <vtkPolyDataConnectivityFilter.h>
+#include <vtkPolyDataNormals.h>
+#include <vtkSphereSource.h>
+#include <vtkThreshold.h>
+#include <vtkWindowedSincPolyDataFilter.h>
+
+#include <cmath>
+#include <iostream>
+
+namespace
+{
+struct Surf
+{
+    double length{0.0};
+    double hausdorff{0.0};
+    double bounds[6]{0, 0, 0, 0, 0, 0};
+    double center[3]{0, 0, 0};
+    double confronto_lati{0.0};
+    int regionId{-1};
+    Point3d centroid;
+};
+
+static Point3d centroidFromPoints(const vtkPoints* pts)
+{
+    Point3d c;
+    if (!pts) return c;
+    const vtkIdType n = pts->GetNumberOfPoints();
+    if (n <= 0) return c;
+
+    double sx = 0.0, sy = 0.0, sz = 0.0;
+    double p[3];
+    for (vtkIdType i = 0; i < n; ++i)
+    {
+        pts->GetPoint(i, p);
+        sx += p[0];
+        sy += p[1];
+        sz += p[2];
+    }
+    c.x = sx / static_cast<double>(n);
+    c.y = sy / static_cast<double>(n);
+    c.z = sz / static_cast<double>(n);
+    return c;
+}
+
+static double sideDifferenceMetric(const double b[6])
+{
+    // Legacy "lati" metric: compare extents.
+    const double dx = std::abs(b[1] - b[0]);
+    const double dy = std::abs(b[3] - b[2]);
+    const double dz = std::abs(b[5] - b[4]);
+
+    // difference between the two largest extents (simple, stable)
+    const double a = std::max({dx, dy, dz});
+    const double c = std::min({dx, dy, dz});
+    const double mid = dx + dy + dz - a - c;
+    return std::max({std::abs(a - mid), std::abs(mid - c), std::abs(a - c)});
+}
+
+static double hausdorffToSphere(vtkPolyData* surface, const double center[3])
+{
+    if (!surface || surface->GetNumberOfPoints() == 0) return 0.0;
+
+    vtkNew<vtkSphereSource> source_sphere;
+    source_sphere->SetCenter(center[0], center[1], center[2]);
+    source_sphere->SetRadius(5.0);
+    source_sphere->SetThetaResolution(24);
+    source_sphere->SetPhiResolution(24);
+    source_sphere->Update();
+
+    vtkPolyData* sphere = source_sphere->GetOutput();
+    if (!sphere || sphere->GetNumberOfPoints() == 0) return 0.0;
+
+    const vtkIdType nSource = surface->GetNumberOfPoints();
+    const vtkIdType nTarget = sphere->GetNumberOfPoints();
+
+    double maxMin = 0.0;
+    double p1[3], p2[3];
+
+    for (vtkIdType i = 0; i < nSource; ++i)
+    {
+        surface->GetPoint(i, p1);
+        double smallest = std::numeric_limits<double>::infinity();
+
+        for (vtkIdType j = 0; j < nTarget; ++j)
+        {
+            sphere->GetPoint(j, p2);
+            const double dx = p1[0] - p2[0];
+            const double dy = p1[1] - p2[1];
+            const double dz = p1[2] - p2[2];
+            const double dist2 = dx * dx + dy * dy + dz * dz;
+            if (dist2 < smallest) smallest = dist2;
+        }
+
+        if (smallest > maxMin) maxMin = smallest;
+    }
+
+    return maxMin;
+}
+
+} // namespace
+
+MarkerList LocalizationService::localize(
+    vtkImageData* volume,
+    const LocalizationParams& params,
+    const ProgressFn& onProgress,
+    const AbortFn& shouldAbort)
+{
+    MarkerList result;
+    if (!volume)
+        return result;
+
+    // Marching cubes
+    vtkNew<vtkMarchingCubes> marchingcubes;
+    marchingcubes->SetInputData(volume);
+    marchingcubes->SetValue(0, params.marchThreshold);
+    marchingcubes->SetComputeNormals(0);
+    marchingcubes->Update();
+
+    if (marchingcubes->GetNumberOfContours() <= 0)
+        return result;
+
+    // Smoothing
+    vtkNew<vtkWindowedSincPolyDataFilter> smoother;
+    smoother->SetInputConnection(marchingcubes->GetOutputPort());
+    smoother->SetNumberOfIterations(20);
+    smoother->NormalizeCoordinatesOn();
+    smoother->SetPassBand(0.01);
+    smoother->BoundarySmoothingOn();
+    smoother->Update();
+
+    // Normals
+    vtkNew<vtkPolyDataNormals> normals;
+    normals->SetInputConnection(smoother->GetOutputPort());
+    normals->Update();
+
+    // Connectivity / region labeling
+    vtkNew<vtkPolyDataConnectivityFilter> connectivity;
+    connectivity->SetInputConnection(normals->GetOutputPort());
+    connectivity->SetExtractionModeToAllRegions();
+    connectivity->ColorRegionsOn();
+    connectivity->Update();
+
+    vtkPolyData* labeled = connectivity->GetOutput();
+    if (!labeled)
+        return result;
+
+    const int nRegions = connectivity->GetNumberOfExtractedRegions();
+    if (nRegions <= 0)
+        return result;
+
+    // Threshold each region out of labeled polydata. This is robust across recent VTK versions.
+    for (int regionId = 0; regionId < nRegions; ++regionId)
+    {
+        if (shouldAbort && shouldAbort())
+            break;
+
+        if (onProgress)
+            onProgress(100.0 * static_cast<double>(regionId + 1) / static_cast<double>(nRegions));
+
+        vtkNew<vtkThreshold> th;
+        th->SetInputData(labeled);
+        th->SetInputArrayToProcess(0, 0, 0, vtkDataObject::FIELD_ASSOCIATION_CELLS, "RegionId");
+        th->SetLowerThreshold(regionId);
+        th->SetUpperThreshold(regionId);
+        th->SetThresholdFunction(vtkThreshold::THRESHOLD_BETWEEN);
+
+        vtkNew<vtkGeometryFilter> geom;
+        geom->SetInputConnection(th->GetOutputPort());
+        geom->Update();
+
+        vtkPolyData* out = geom->GetOutput();
+        if (!out || out->GetNumberOfCells() == 0)
+            continue;
+
+        out->ComputeBounds();
+        double b[6];
+        out->GetBounds(b);
+        const double len = out->GetLength();
+
+        bool ok = true;
+
+        // 0: bounding-box diagonal length
+        if (params.selectedFilters.size() > 0 && params.selectedFilters[0])
+        {
+            if (len < params.thrDown_D || len > params.thrUp_D)
+                ok = false;
+        }
+
+        // 1: Hausdorff
+        if (ok && params.selectedFilters.size() > 1 && params.selectedFilters[1])
+        {
+            double center[3];
+            out->GetCenter(center);
+            const double hd = hausdorffToSphere(out, center);
+            if (hd > params.thr_HAUSD)
+                ok = false;
+        }
+
+        // 2: side-difference
+        if (ok && params.selectedFilters.size() > 2 && params.selectedFilters[2])
+        {
+            const double metric = sideDifferenceMetric(b);
+            if (metric > params.thr_S)
+                ok = false;
+        }
+
+        if (!ok)
+            continue;
+
+        // Centroid from edges, to match legacy behavior.
+        vtkNew<vtkExtractEdges> edge_extractor;
+        edge_extractor->SetInputData(out);
+        edge_extractor->Update();
+
+        vtkPolyData* edges = edge_extractor->GetOutput();
+        if (!edges || edges->GetNumberOfPoints() == 0)
+            continue;
+
+        const Point3d c_mm = centroidFromPoints(edges->GetPoints());
+
+        Marker m;
+        m.regionId = regionId;
+        m.centroid = Point3d{c_mm.x / 1000.0, c_mm.y / 1000.0, c_mm.z / 1000.0};
+        result.push_back(m);
+    }
+
+    return result;
+}