working code, incl. rendering localisation volume clipping and skull masking

LocalizationService.cpp

@@ -15,6 +15,7 @@
 #include <vtkSphereSource.h>
 #include <vtkThreshold.h>
 #include <vtkWindowedSincPolyDataFilter.h>
+#include <vtkCleanPolyData.h>
 
 #include <cmath>
 #include <iostream>
@@ -146,62 +147,72 @@ MarkerList LocalizationService::localize(
     normals->SetInputConnection(smoother->GetOutputPort());
     normals->Update();
 
-    // Connectivity / region labeling
-    vtkNew<vtkPolyDataConnectivityFilter> connectivity;
-    connectivity->SetInputConnection(normals->GetOutputPort());
-    connectivity->SetExtractionModeToAllRegions();
-    connectivity->ColorRegionsOn();
-    connectivity->Update();
+    // 1) Label / count regions
+    vtkNew<vtkPolyDataConnectivityFilter> connAll;
+    connAll->SetInputConnection(normals->GetOutputPort());
+    connAll->SetExtractionModeToAllRegions();
+    connAll->ColorRegionsOn();
+    connAll->Update();
 
-    vtkPolyData* labeled = connectivity->GetOutput();
-    if (!labeled)
-        return result;
-
-    // VTK has changed over time whether "RegionId" is stored on cells or points.
-    // For stable per-region extraction we want it on *cells*.
-    vtkPolyData* labeledForThreshold = labeled;
-    vtkNew<vtkPointDataToCellData> p2c;
-    if (!labeled->GetCellData()->GetArray("RegionId") && labeled->GetPointData()->GetArray("RegionId"))
-    {
-        p2c->SetInputData(labeled);
-        p2c->PassPointDataOn();
-        p2c->Update();
-        labeledForThreshold = vtkPolyData::SafeDownCast(p2c->GetOutput());
-    }
-
-    const int nRegions = connectivity->GetNumberOfExtractedRegions();
+    const int nRegions = connAll->GetNumberOfExtractedRegions();
+    std::cout << "nRegions " << nRegions << std::endl;
     if (nRegions <= 0)
         return result;
 
-    // Threshold each region out of labeled polydata. This is robust across recent VTK versions.
+    // 2) Extract each region
+    vtkNew<vtkPolyDataConnectivityFilter> connOne;
+    connOne->SetInputConnection(normals->GetOutputPort());
+    connOne->SetExtractionModeToSpecifiedRegions();
+    connOne->ColorRegionsOff();
+    connOne->Update();
+
     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));
+            onProgress(100.0 * double(regionId + 1) / double(nRegions));
 
-        vtkNew<vtkThreshold> th;
-        th->SetInputData(labeledForThreshold);
-        th->SetInputArrayToProcess(0, 0, 0, vtkDataObject::FIELD_ASSOCIATION_CELLS, "RegionId");
-        th->SetLowerThreshold(regionId);
-        th->SetUpperThreshold(regionId);
-        th->SetThresholdFunction(vtkThreshold::THRESHOLD_BETWEEN);
+        connOne->InitializeSpecifiedRegionList();
+        connOne->AddSpecifiedRegion(regionId);
+        connOne->Modified();
+        connOne->Update();
 
-        vtkNew<vtkGeometryFilter> geom;
-        geom->SetInputConnection(th->GetOutputPort());
-        geom->Update();
-
-        vtkPolyData* out = geom->GetOutput();
-        if (!out || out->GetNumberOfCells() == 0)
+        vtkPolyData* raw = connOne->GetOutput();
+        if (!raw || raw->GetNumberOfCells() == 0)
+        {
+            std::cout << "regionId " << regionId << " => cells 0\n";
             continue;
+        }
+
+        // ✅ IMPORTANT: compact points so bounds/length reflect ONLY this region
+        vtkNew<vtkCleanPolyData> clean;
+        clean->SetInputData(raw);
+        clean->PointMergingOff();     // keep geometry identical; just drop unused points
+        clean->Update();
+
+        vtkPolyData* out = clean->GetOutput();
+        if (!out || out->GetNumberOfCells() == 0 || out->GetNumberOfPoints() == 0)
+        {
+            std::cout << "regionId " << regionId << " => empty after clean\n";
+            continue;
+        }
 
         out->ComputeBounds();
         double b[6];
         out->GetBounds(b);
         const double len = out->GetLength();
 
+        std::cout << "regionId " << regionId
+                  << " => cells " << out->GetNumberOfCells()
+                  << " points " << out->GetNumberOfPoints()
+                  << " len " << len
+                  << " bounds "
+                  << b[0] << " " << b[1] << " "
+                  << b[2] << " " << b[3] << " "
+                  << b[4] << " " << b[5] << "\n";
+
         bool ok = true;
 
         // 0: bounding-box diagonal length