diff --git a/CMakeLists.txt b/CMakeLists.txt index 94d5128..3914fc8 100644 --- a/CMakeLists.txt +++ b/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 diff --git a/TestScene.scn b/TestScene.scn deleted file mode 100644 index d077089..0000000 --- a/TestScene.scn +++ /dev/null @@ -1,80 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - \ No newline at end of file diff --git a/src/MeshSkeletonizationPlugin/CGALMeshUtils.h b/src/MeshSkeletonizationPlugin/CGALMeshUtils.h new file mode 100644 index 0000000..9ac9857 --- /dev/null +++ b/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 +#include +#include +#include +#include +#include + +#include + +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 MeshToPolyhedronOp : public CGAL::Modifier_base + { + public: + MeshToPolyhedronOp(const VecCoord& vertices, const SeqTriangles& triangles) + : m_vertices(vertices), m_triangles(triangles) + { + } + + void operator()(HalfedgeDS& hds) override + { + CGAL::Polyhedron_incremental_builder_3 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>; + using AABBTree = CGAL::AABB_tree; + using PointInsideTest = CGAL::Side_of_triangle_mesh; + + /// 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 tree; + std::unique_ptr insideTest; + + void build() + { + tree = std::make_unique(CGAL::faces(polyhedron).first, CGAL::faces(polyhedron).second, polyhedron); + tree->accelerate_distance_queries(); + insideTest = std::make_unique(*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 + void buildFrom(const VecCoord& vertices, const SeqTriangles& triangles) + { + MeshToPolyhedronOp op(vertices, triangles); + polyhedron.delegate(op); + if (!polyhedron.is_empty()) + build(); + } + }; + +} // namespace cgalutils +} // namespace meshskeletonizationplugin diff --git a/src/MeshSkeletonizationPlugin/SkeletonGraph/SkeletonGraph.cpp b/src/MeshSkeletonizationPlugin/SkeletonGraph/SkeletonGraph.cpp index 63a5d8c..256cb1e 100644 --- a/src/MeshSkeletonizationPlugin/SkeletonGraph/SkeletonGraph.cpp +++ b/src/MeshSkeletonizationPlugin/SkeletonGraph/SkeletonGraph.cpp @@ -1,6 +1,5 @@ #include - #include #include #include @@ -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& p) const return best; } +void SkeletonGraph::buildTreeAutoRoot() +{ + if (m_nodes.empty()) + return; + + std::vector 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 component; + std::queue 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& entryPoint) { int root = closestNodeId(entryPoint); diff --git a/src/MeshSkeletonizationPlugin/SkeletonGraph/SkeletonGraph.h b/src/MeshSkeletonizationPlugin/SkeletonGraph/SkeletonGraph.h index a61aaf6..f621e0c 100644 --- a/src/MeshSkeletonizationPlugin/SkeletonGraph/SkeletonGraph.h +++ b/src/MeshSkeletonizationPlugin/SkeletonGraph/SkeletonGraph.h @@ -38,6 +38,10 @@ class SOFA_MESHSKELETONIZATIONPLUGIN_API SkeletonGraph void buildTree(const std::array& 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>& meshVertices); diff --git a/src/MeshSkeletonizationPlugin/SkeletonGraph/SkeletonReader.inl b/src/MeshSkeletonizationPlugin/SkeletonGraph/SkeletonReader.inl index 1745b8f..eeb5202 100644 --- a/src/MeshSkeletonizationPlugin/SkeletonGraph/SkeletonReader.inl +++ b/src/MeshSkeletonizationPlugin/SkeletonGraph/SkeletonReader.inl @@ -10,7 +10,7 @@ template SkeletonReader::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::doUpdate() msg_info() << "Skeleton loaded: " << m_graph.nodes().size() << " node(s)."; d_outNodeCount.setValue(static_cast(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()) {