[src] Adding buildTreeAutoRoot method in SkeletonGraph

CMakeLists.txt

@@ -18,6 +18,7 @@ set(HEADER_FILES
     ${PLUGIN_SKELETONIZATION_SRC_DIR}/SkeletonGraph/SkeletonGraph.h
     ${PLUGIN_SKELETONIZATION_SRC_DIR}/SkeletonGraph/SkeletonReader.h
     ${PLUGIN_SKELETONIZATION_SRC_DIR}/SkeletonGraph/SkeletonReader.inl
+    ${PLUGIN_SKELETONIZATION_SRC_DIR}/CGALMeshUtils.h
 )
 set(SOURCE_FILES
     ${PLUGIN_SKELETONIZATION_SRC_DIR}/init.cpp

src/MeshSkeletonizationPlugin/CGALMeshUtils.h

@@ -0,0 +1,89 @@
+#pragma once
+// Reuses the Kernel/Polyhedron/Point/HalfedgeDS typedefs already declared (at global scope) in MeshSkeletonization.h.
+#include <MeshSkeletonizationPlugin/MeshSkeletonization.h>
+#include <CGAL/AABB_tree.h>
+#include <CGAL/AABB_traits.h>
+#include <CGAL/AABB_face_graph_triangle_primitive.h>
+#include <CGAL/Side_of_triangle_mesh.h>
+#include <CGAL/Polyhedron_incremental_builder_3.h>
+
+#include <memory>
+
+namespace meshskeletonizationplugin
+{
+namespace cgalutils
+{
+    /// Builds a CGAL Polyhedron_3 from a flat vertex/triangle mesh, mirroring
+    /// MeshSkeletonization::geometryToPolyhedronOp but as a free/shared helper
+    /// so multiple components can build closed meshes without duplicating it.
+    template <class VecCoord, class SeqTriangles>
+    class MeshToPolyhedronOp : public CGAL::Modifier_base<HalfedgeDS>
+    {
+    public:
+        MeshToPolyhedronOp(const VecCoord& vertices, const SeqTriangles& triangles)
+            : m_vertices(vertices), m_triangles(triangles)
+        {
+        }
+
+        void operator()(HalfedgeDS& hds) override
+        {
+            CGAL::Polyhedron_incremental_builder_3<HalfedgeDS> builder(hds, true);
+            builder.begin_surface(m_vertices.size(), m_triangles.size());
+
+            for (const auto& v : m_vertices)
+                builder.add_vertex(Point(v[0], v[1], v[2]));
+
+            for (const auto& tri : m_triangles)
+            {
+                builder.begin_facet();
+                for (int j = 0; j < 3; ++j)
+                    builder.add_vertex_to_facet(tri[j]);
+                builder.end_facet();
+            }
+
+            if (builder.check_unconnected_vertices())
+                builder.remove_unconnected_vertices();
+
+            builder.end_surface();
+        }
+
+    private:
+        const VecCoord& m_vertices;
+        const SeqTriangles& m_triangles;
+    };
+
+    using AABBTraits = CGAL::AABB_traits<Kernel, CGAL::AABB_face_graph_triangle_primitive<Polyhedron>>;
+    using AABBTree = CGAL::AABB_tree<AABBTraits>;
+    using PointInsideTest = CGAL::Side_of_triangle_mesh<Polyhedron, Kernel>;
+
+    /// A closed mesh (segment, tumor, ...) ready for point-in-mesh and
+    /// closest-point queries: build a Polyhedron_3 from raw vertices/
+    /// triangles, then an AABB tree + inside/outside test backed by it.
+    struct ClosedMeshQuery
+    {
+        Polyhedron polyhedron;
+        std::unique_ptr<AABBTree> tree;
+        std::unique_ptr<PointInsideTest> insideTest;
+
+        void build()
+        {
+            tree = std::make_unique<AABBTree>(CGAL::faces(polyhedron).first, CGAL::faces(polyhedron).second, polyhedron);
+            tree->accelerate_distance_queries();
+            insideTest = std::make_unique<PointInsideTest>(*tree);
+        }
+
+        /// Convenience: build the polyhedron from flat vertices/triangles and
+        /// immediately build the tree/inside-test on top of it. No-ops (and
+        /// leaves tree/insideTest null) if the resulting polyhedron is empty.
+        template <class VecCoord, class SeqTriangles>
+        void buildFrom(const VecCoord& vertices, const SeqTriangles& triangles)
+        {
+            MeshToPolyhedronOp<VecCoord, SeqTriangles> op(vertices, triangles);
+            polyhedron.delegate(op);
+            if (!polyhedron.is_empty())
+                build();
+        }
+    };
+
+} // namespace cgalutils
+} // namespace meshskeletonizationplugin

src/MeshSkeletonizationPlugin/SkeletonGraph/SkeletonGraph.cpp

@@ -1,6 +1,5 @@
 #include <MeshSkeletonizationPlugin/SkeletonGraph/SkeletonGraph.h>
 
-
 #include <algorithm>
 #include <cmath>
 #include <fstream>
@@ -66,8 +65,7 @@ bool SkeletonGraph::loadFromFile(const std::string& filename, double mergeTolera
 
     while (std::getline(in, line))
     {
-        // Trim trailing whitespace/CR so blank-line detection works on
-        // files written or edited on Windows too.
+        // Trim trailing whitespace/CR so blank-line detection works on files written or edited on Windows too.
         while (!line.empty() && (line.back() == '\r' || line.back() == ' ' || line.back() == '\t'))
             line.pop_back();
 
@@ -200,6 +198,57 @@ int SkeletonGraph::closestNodeId(const std::array<double, 3>& p) const
     return best;
 }
 
+void SkeletonGraph::buildTreeAutoRoot()
+{
+    if (m_nodes.empty())
+        return;
+
+    std::vector<bool> visited(m_nodes.size(), false);
+    int bestRoot = -1;
+    std::size_t bestSize = 0;
+
+    for (const SkeletonNode& start : m_nodes)
+    {
+        if (visited[start.id()])
+            continue;
+
+        // BFS this component, just to measure it and grab a representative node.
+        std::vector<int> component;
+        std::queue<int> q;
+        visited[start.id()] = true;
+        q.push(start.id());
+
+        while (!q.empty())
+        {
+            int u = q.front();
+            q.pop();
+            component.push_back(u);
+
+            auto it = m_adjacency.find(u);
+            if (it == m_adjacency.end())
+                continue;
+
+            for (int v : it->second)
+            {
+                if (!visited[v])
+                {
+                    visited[v] = true;
+                    q.push(v);
+                }
+            }
+        }
+
+        if (component.size() > bestSize)
+        {
+            bestSize = component.size();
+            bestRoot = component.front();
+        }
+    }
+
+    if (bestRoot >= 0)
+        buildTree(bestRoot);
+}
+
 void SkeletonGraph::buildTree(const std::array<double, 3>& entryPoint)
 {
     int root = closestNodeId(entryPoint);

src/MeshSkeletonizationPlugin/SkeletonGraph/SkeletonGraph.h

@@ -38,6 +38,10 @@ class SOFA_MESHSKELETONIZATIONPLUGIN_API SkeletonGraph
     void buildTree(const std::array<double, 3>& entryPoint);
     void buildTree(int rootId);
 
+    /// Picks a root automatically, with no entry point needed: the node belonging to the LARGEST connected component of the raw connectivity graph.  
+    /// This guarantees the tree is rooted in the dominant structure instead of risking a tiny isolated fragment.
+    void buildTreeAutoRoot();
+
     /// Best-effort correspondence between each skeleton node and the closest
     /// vertex of the input (vessel) mesh; fills meshVertexId/distanceToMesh.
     void computeMeshCorrespondence(const std::vector<std::array<double, 3>>& meshVertices);

src/MeshSkeletonizationPlugin/SkeletonGraph/SkeletonReader.inl

@@ -10,7 +10,7 @@ template <class DataTypes>
 SkeletonReader<DataTypes>::SkeletonReader()
     : d_inSkeletonFilename(initData(&d_inSkeletonFilename, "filename", "Skeleton polyline file to read (e.g. skeleton.txt)"))
     , d_inVertices(initData(&d_inVertices, "inputVertices", "Optional input mesh vertices, to link skeleton nodes to the mesh"))
-    , d_inEntryPoint(initData(&d_inEntryPoint, Vec3(0, 0, 0), "entryPoint", "Approx. entry point; closest node becomes the tree root"))
+    , d_inEntryPoint(initData(&d_inEntryPoint, Vec3(0, 0, 0), "entryPoint", "Approx. entry point; closest node becomes the tree root. If left unset, a root is instead picked automatically from the largest connected component of the skeleton."))
     , d_outVTKFilename(initData(&d_outVTKFilename, "outputVTK", "File path to export the rooted tree (.vtk)"))
     , d_outReportFilename(initData(&d_outReportFilename, "outputReport", "File path to export a per-node CSV report (id, x, y, z, parentId, childrenIds, pathFromRoot)"))
     , d_outNodeCount(initData(&d_outNodeCount, 0, "nodeCount", "Number of skeleton nodes read"))
@@ -61,8 +61,18 @@ void SkeletonReader<DataTypes>::doUpdate()
     msg_info() << "Skeleton loaded: " << m_graph.nodes().size() << " node(s).";
     d_outNodeCount.setValue(static_cast<int>(m_graph.nodes().size()));
 
-    const Vec3& entry = d_inEntryPoint.getValue();
-    m_graph.buildTree({ double(entry[0]), double(entry[1]), double(entry[2]) });
+    if (d_inEntryPoint.isSet())
+    {
+        const Vec3& entry = d_inEntryPoint.getValue();
+        m_graph.buildTree({ double(entry[0]), double(entry[1]), double(entry[2]) });
+        msg_info() << "Rooted using explicit entryPoint " << entry << ".";
+    }
+    else
+    {
+        // No explicit entry point given . Pick a root from the largest connected component instead.
+        m_graph.buildTreeAutoRoot();
+        msg_info() << "No entryPoint set; auto-rooted from the largest connected component.";
+    }
 
     if (!m_graph.hasRoot())
     {