diff --git a/.gitignore b/.gitignore index 73a7d6f..5d43e36 100644 --- a/.gitignore +++ b/.gitignore @@ -5,6 +5,5 @@ *.vtk *.obj *.sxn -*.scn *.csv diff --git a/CMakeLists.txt b/CMakeLists.txt index 3b0951f..1dc5547 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -18,14 +18,18 @@ 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}/SegmentMapping/SkeletonSegmentMapper.inl + ${PLUGIN_SKELETONIZATION_SRC_DIR}/SegmentMapping/SkeletonSegmentMapper.h + ${PLUGIN_SKELETONIZATION_SRC_DIR}/SkeletonResection/SkeletonResectionSimulator.h + ${PLUGIN_SKELETONIZATION_SRC_DIR}/SkeletonResection/SkeletonResectionSimulator.inl ) set(SOURCE_FILES ${PLUGIN_SKELETONIZATION_SRC_DIR}/init.cpp ${PLUGIN_SKELETONIZATION_SRC_DIR}/MeshSkeletonization.cpp ${PLUGIN_SKELETONIZATION_SRC_DIR}/SkeletonGraph/SkeletonGraph.cpp ${PLUGIN_SKELETONIZATION_SRC_DIR}/SkeletonGraph/SkeletonReader.cpp - - + ${PLUGIN_SKELETONIZATION_SRC_DIR}/SegmentMapping/SkeletonSegmentMapper.cpp + ${PLUGIN_SKELETONIZATION_SRC_DIR}/SkeletonResection/SkeletonResectionSimulator.cpp ) set(README_FILES README.md) diff --git a/TestScene.scn b/TestScene.scn new file mode 100644 index 0000000..d077089 --- /dev/null +++ b/TestScene.scn @@ -0,0 +1,80 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/scenes/BeamFEMSkeleton.scn b/scenes/BeamFEMSkeleton.scn index 9c86bd9..b827dab 100644 --- a/scenes/BeamFEMSkeleton.scn +++ b/scenes/BeamFEMSkeleton.scn @@ -12,12 +12,7 @@ - + diff --git a/src/MeshSkeletonizationPlugin/SegmentMapping/SkeletonSegmentMapper.cpp b/src/MeshSkeletonizationPlugin/SegmentMapping/SkeletonSegmentMapper.cpp new file mode 100644 index 0000000..dbe59d5 --- /dev/null +++ b/src/MeshSkeletonizationPlugin/SegmentMapping/SkeletonSegmentMapper.cpp @@ -0,0 +1,22 @@ +#define SKELETONSEGMENTMAPPER_CPP +#include + +#include +#include + +namespace meshskeletonizationplugin +{ +using namespace sofa::defaulttype; + +void registerSkeletonSegmentMapper(sofa::core::ObjectFactory* factory) +{ + factory->registerObjects(sofa::core::ObjectRegistrationData( + "Maps a previously-read skeleton (SkeletonReader) onto liver segments (e.g. " + "Couinaud I-VIII), using a per-mesh-vertex segment label and smoothing " + "per-branch by majority vote") + .add< SkeletonSegmentMapper >()); +} + +template class SOFA_MESHSKELETONIZATIONPLUGIN_API SkeletonSegmentMapper; + +} // namespace meshskeletonizationplugin diff --git a/src/MeshSkeletonizationPlugin/SegmentMapping/SkeletonSegmentMapper.h b/src/MeshSkeletonizationPlugin/SegmentMapping/SkeletonSegmentMapper.h new file mode 100644 index 0000000..c59ea3d --- /dev/null +++ b/src/MeshSkeletonizationPlugin/SegmentMapping/SkeletonSegmentMapper.h @@ -0,0 +1,79 @@ +#pragma once + +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +using namespace sofa; +using namespace sofa::defaulttype; + +namespace meshskeletonizationplugin +{ + +template +class SkeletonSegmentMapper : public sofa::core::DataEngine +{ +public: + SOFA_CLASS(SOFA_TEMPLATE(SkeletonSegmentMapper, DataTypes), sofa::core::DataEngine); + + typedef typename DataTypes::Coord Coord; + typedef typename DataTypes::VecCoord VecCoord; + typedef typename Coord::value_type Real; + typedef type::Vec<3, Real> Vec3; + + /// The SkeletonReader whose graph should be segmented. If left unset, the + /// first SkeletonReader found in the context is used. + sofa::core::objectmodel::SingleLink< + SkeletonSegmentMapper, + SkeletonReader, + sofa::core::objectmodel::BaseLink::FLAG_STOREPATH | sofa::core::objectmodel::BaseLink::FLAG_STRONGLINK> + l_skeletonReader; + + /// One closed mesh per liver segment (e.g. the Couinaud collision meshes), + /// in the same order as d_inSegmentNames. + sofa::core::objectmodel::MultiLink< + SkeletonSegmentMapper, + sofa::core::loader::MeshLoader, + sofa::core::objectmodel::BaseLink::FLAG_STOREPATH> + l_segmentMeshes; + + //Component parameters + // Inputs + sofa::core::objectmodel::Data> d_inSegmentNames; ///< Human-readable name per entry of l_segmentMeshes, e.g. "II", "IVa", "VIII" + sofa::core::objectmodel::DataFileName d_outSegmentReportFilename; ///< Optional CSV report: id, x, y, z, segment name + // Outputs + sofa::core::objectmodel::Data> d_outNodeSegments; ///< Segment index per skeleton node (into l_segmentMeshes/d_inSegmentNames, -1 = unknown), indexed like graph().nodes() + sofa::core::objectmodel::Data> d_outNodeSegmentNames; ///< Same as d_outNodeSegments, resolved to names ("unknown" if -1) + sofa::core::objectmodel::Data d_outSegmentCount; ///< Number of distinct segments found (excluding unknown) + + /// Direct access to the segmented graph, e.g. for another component to + /// query via getContext()->get>()->graph(). + const SkeletonGraph& graph() const { return m_graph; } + + void init() override; + void doUpdate() override; + void draw(const sofa::core::visual::VisualParams* vparams) override; + +private: + SkeletonSegmentMapper(); + virtual ~SkeletonSegmentMapper() = default; + + /// Local working copy of the linked reader's graph, augmented with + /// segment labels (SkeletonGraph is a plain value type, cheap to copy). + SkeletonGraph m_graph; +}; + +#if !defined(SKELETONSEGMENTMAPPER_CPP) +extern template class SOFA_MESHSKELETONIZATIONPLUGIN_API SkeletonSegmentMapper; +#endif + +} // namespace meshskeletonizationplugin diff --git a/src/MeshSkeletonizationPlugin/SegmentMapping/SkeletonSegmentMapper.inl b/src/MeshSkeletonizationPlugin/SegmentMapping/SkeletonSegmentMapper.inl new file mode 100644 index 0000000..9f1870b --- /dev/null +++ b/src/MeshSkeletonizationPlugin/SegmentMapping/SkeletonSegmentMapper.inl @@ -0,0 +1,273 @@ +#pragma once +#include + +// Reuses the Kernel/Polyhedron/Point typedefs already declared (at global +// scope) in MeshSkeletonization.h, so segment meshes are built the same way +// MeshSkeletonization builds its own polyhedron. +#include + +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include + +namespace meshskeletonizationplugin +{ + +namespace +{ + // Builds a CGAL Polyhedron_3 from a flat vertex/triangle mesh, mirroring + // MeshSkeletonization::geometryToPolyhedronOp but as a free function so it + // can be reused here without duplicating a private nested class. + 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; + + /// One segment mesh, ready for point-in-mesh + closest-point queries. + struct SegmentMeshQuery + { + 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); + } + }; +} // anonymous namespace + +template +SkeletonSegmentMapper::SkeletonSegmentMapper() + : l_skeletonReader(initLink("skeletonReader", + "SkeletonReader whose graph should be segmented; if empty, the first " + "SkeletonReader found in the scene is used")) + , l_segmentMeshes(initLink("segmentMeshes", + "One closed mesh loader per liver segment (e.g. Couinaud collision meshes), " + "in the same order as segmentNames")) + , d_inSegmentNames(initData(&d_inSegmentNames, "segmentNames", + "Human-readable name per entry of segmentMeshes, e.g. 'II', 'IVa', 'VIII'")) + , d_outSegmentReportFilename(initData(&d_outSegmentReportFilename, "outSegmentReportFilename", + "Optional path to export a CSV report: id, x, y, z, segment name")) + , d_outNodeSegments(initData(&d_outNodeSegments, "nodeSegments", + "Segment index per skeleton node (-1 if unknown), indexed like the graph's nodes")) + , d_outNodeSegmentNames(initData(&d_outNodeSegmentNames, "nodeSegmentNames", + "Same as nodeSegments, resolved to names ('unknown' if -1)")) + , d_outSegmentCount(initData(&d_outSegmentCount, 0, "segmentCount", + "Number of distinct segments found (excluding -1/unknown)")) +{ + addOutput(&d_outNodeSegments); + addOutput(&d_outNodeSegmentNames); + addOutput(&d_outSegmentCount); +} + +template +void SkeletonSegmentMapper::init() +{ + if (!l_skeletonReader) + { + l_skeletonReader.set(this->getContext()->template get>()); + } + if (!l_skeletonReader) + { + msg_error() << "No SkeletonReader found. Set 'skeletonReader' to a valid " + "component path, or add one earlier in the scene."; + } + + if (l_segmentMeshes.empty()) + { + msg_error() << "No segment meshes linked: set 'segmentMeshes' to a list of mesh " + "loader paths (one per liver segment)."; + } + else if (l_segmentMeshes.size() != d_inSegmentNames.getValue().size()) + { + msg_warning() << "segmentMeshes (" << l_segmentMeshes.size() << ") and segmentNames (" + << d_inSegmentNames.getValue().size() << ") have different sizes; " + "unnamed segments will be reported by index only."; + } + + setDirtyValue(); + update(); // compute eagerly, don't wait for something to read an output Data +} + +template +void SkeletonSegmentMapper::doUpdate() +{ + if (!l_skeletonReader || l_segmentMeshes.empty()) + return; + + // Work on our own copy of the reader's graph so we don't need a + // non-const accessor on SkeletonReader just for this. + m_graph = l_skeletonReader->graph(); + + // Build one point-in-mesh query per linked segment mesh. + std::vector segmentQueries(l_segmentMeshes.size()); + for (std::size_t i = 0; i < l_segmentMeshes.size(); ++i) + { + sofa::core::loader::MeshLoader* loader = l_segmentMeshes.get(i); + if (!loader) + continue; + + MeshToPolyhedronOp< + sofa::type::vector, + sofa::type::vector> + op(loader->d_positions.getValue(), loader->d_triangles.getValue()); + + segmentQueries[i].polyhedron.delegate(op); + if (!segmentQueries[i].polyhedron.is_empty()) + segmentQueries[i].build(); + } + + // Raw per-node label: which segment mesh contains this node, or (if none + // does) whichever segment mesh's surface is closest. + const auto& nodes = m_graph.nodes(); + std::vector raw(nodes.size(), -1); + + for (const SkeletonNode& node : nodes) + { + const auto& p = node.position(); + Point query(p[0], p[1], p[2]); + + int inside = -1; + for (std::size_t i = 0; i < segmentQueries.size() && inside == -1; ++i) + { + if (segmentQueries[i].insideTest && + (*segmentQueries[i].insideTest)(query) == CGAL::ON_BOUNDED_SIDE) + { + inside = static_cast(i); + } + } + + if (inside != -1) + { + raw[node.id()] = inside; + continue; + } + + // Fallback: closest segment surface. + double bestDist2 = std::numeric_limits::max(); + int bestSeg = -1; + for (std::size_t i = 0; i < segmentQueries.size(); ++i) + { + if (!segmentQueries[i].tree) + continue; + double d2 = CGAL::to_double(segmentQueries[i].tree->squared_distance(query)); + if (d2 < bestDist2) + { + bestDist2 = d2; + bestSeg = static_cast(i); + } + } + raw[node.id()] = bestSeg; + } + + m_graph.assignSegmentLabels(raw); + + const auto& names = d_inSegmentNames.getValue(); + + sofa::helper::WriteAccessor>> outSegments(d_outNodeSegments); + sofa::helper::WriteAccessor>> outSegmentNames(d_outNodeSegmentNames); + outSegments.resize(nodes.size()); + outSegmentNames.resize(nodes.size()); + + std::set distinct; + for (const SkeletonNode& n : nodes) + { + int seg = m_graph.segmentOf(n.id()); + outSegments[n.id()] = seg; + outSegmentNames[n.id()] = (seg >= 0 && seg < static_cast(names.size())) ? names[seg] : "unknown"; + if (seg != -1) + distinct.insert(seg); + } + + d_outSegmentCount.setValue(static_cast(distinct.size())); + + if (d_outSegmentReportFilename.isSet() && !d_outSegmentReportFilename.getFullPath().empty()) + { + m_graph.exportSegmentReportCSV(d_outSegmentReportFilename.getFullPath(), names); + } +} + +template +void SkeletonSegmentMapper::draw(const sofa::core::visual::VisualParams* vparams) +{ + if (!vparams->displayFlags().getShowBehaviorModels()) + return; + + static const std::vector palette = { + sofa::type::RGBAColor(0.90f, 0.10f, 0.10f, 1.0f), + sofa::type::RGBAColor(0.10f, 0.60f, 0.90f, 1.0f), + sofa::type::RGBAColor(0.20f, 0.80f, 0.20f, 1.0f), + sofa::type::RGBAColor(0.95f, 0.75f, 0.10f, 1.0f), + sofa::type::RGBAColor(0.60f, 0.20f, 0.80f, 1.0f), + sofa::type::RGBAColor(0.90f, 0.45f, 0.10f, 1.0f), + sofa::type::RGBAColor(0.10f, 0.80f, 0.70f, 1.0f), + sofa::type::RGBAColor(0.50f, 0.50f, 0.50f, 1.0f), + }; + static const sofa::type::RGBAColor unknownColor(0.7f, 0.7f, 0.7f, 1.0f); + + std::vector points; + std::vector colors; + points.reserve(m_graph.nodes().size()); + colors.reserve(m_graph.nodes().size()); + + for (const SkeletonNode& n : m_graph.nodes()) + { + const auto& p = n.position(); + points.emplace_back(p[0], p[1], p[2]); + + int seg = m_graph.segmentOf(n.id()); + colors.push_back(seg >= 0 ? palette[seg % palette.size()] : unknownColor); + } + + vparams->drawTool()->drawPoints(points, 6.0f, colors); +} + +} // namespace meshskeletonizationplugin diff --git a/src/MeshSkeletonizationPlugin/SkeletonGraph/SkeletonGraph.cpp b/src/MeshSkeletonizationPlugin/SkeletonGraph/SkeletonGraph.cpp index 2c650e4..f7e51be 100644 --- a/src/MeshSkeletonizationPlugin/SkeletonGraph/SkeletonGraph.cpp +++ b/src/MeshSkeletonizationPlugin/SkeletonGraph/SkeletonGraph.cpp @@ -16,6 +16,7 @@ void SkeletonGraph::clear() m_nodes.clear(); m_adjacency.clear(); m_loopEdges.clear(); + m_nodeSegment.clear(); m_rootId = -1; } @@ -91,6 +92,93 @@ bool SkeletonGraph::loadFromFile(const std::string& filename, double mergeTolera return true; } +bool SkeletonGraph::loadFromVTK(const std::string& filename) +{ + std::ifstream in(filename); + if (!in.is_open()) + return false; + + clear(); + + // Skip the 4-line ASCII VTK header: version, title, "ASCII", "DATASET POLYDATA". + std::string discard; + for (int i = 0; i < 4; ++i) + std::getline(in, discard); + + std::string tag; + int currentBlockCount = 0; // set by POINT_DATA/CELL_DATA, used to size SCALARS reads + + while (in >> tag) + { + if (tag == "POINTS") + { + int numPoints = 0; + std::string dataType; + in >> numPoints >> dataType; + + m_nodes.clear(); + m_nodes.reserve(numPoints); + for (int i = 0; i < numPoints; ++i) + { + std::array p{}; + in >> p[0] >> p[1] >> p[2]; + m_nodes.emplace_back(i, p[0], p[1], p[2]); + } + } + else if (tag == "LINES") + { + int numLines = 0, listSize = 0; + in >> numLines >> listSize; + for (int i = 0; i < numLines; ++i) + { + int count = 0, a = 0, b = 0; + in >> count >> a >> b; + if (count == 2 && a >= 0 && b >= 0 && + a < static_cast(m_nodes.size()) && b < static_cast(m_nodes.size())) + { + // exportToVTK() writes "parentId childId" per line. + m_nodes[a].addChildId(b); + m_nodes[b].addParentId(a); + + // Also record the raw undirected connectivity - this is + // what simulateResection()'s reachability BFS walks, so + // without it a cut would never propagate to children. + connect(a, b); + } + } + } + else if (tag == "POINT_DATA" || tag == "CELL_DATA") + { + in >> currentBlockCount; + } + else if (tag == "SCALARS") + { + std::string name, dataType, lookupTag, lookupName; + in >> name >> dataType >> lookupTag >> lookupName; // " \nLOOKUP_TABLE default" + + if (name == "type" && currentBlockCount == static_cast(m_nodes.size())) + { + for (int i = 0; i < currentBlockCount; ++i) + { + int t = 0; + in >> t; + if (t == 1) + m_rootId = i; + } + } + else + { + // Not needed to rebuild the tree (depth, is_loop, ...): consume and discard. + double dummy; + for (int i = 0; i < currentBlockCount; ++i) + in >> dummy; + } + } + } + + return hasRoot(); +} + int SkeletonGraph::closestNodeId(const std::array& p) const { int best = -1; @@ -409,4 +497,204 @@ void SkeletonGraph::exportReportCSV(const std::string& filename) const << "\"" << joinIds(loopParents) << "\"\n"; } } + +// --- Liver-segment mapping --------------------------------------------------- + +void SkeletonGraph::assignSegmentLabels(const std::vector& rawNodeLabels) +{ + const int n = static_cast(m_nodes.size()); + m_nodeSegment.assign(n, -1); + + if (!hasRoot() || rawNodeLabels.empty()) + return; + + const std::vector& raw = rawNodeLabels; + + // Walk the tree from the root, splitting it into maximal branches: a + // branch runs from the root (or right after a branch point) down to the + // next branch point or leaf. Every node in a branch gets that branch's + // majority raw label, so a few mislabeled nodes near a segment boundary + // don't fragment an otherwise-clear branch. + std::vector stack; + stack.push_back(m_rootId); + + while (!stack.empty()) + { + int start = stack.back(); + stack.pop_back(); + + std::vector chain; + int cur = start; + while (true) + { + chain.push_back(cur); + const SkeletonNode* curNode = node(cur); + const auto& children = curNode->childrenIds(); + + if (children.size() == 1) + { + cur = children.front(); + continue; + } + + // Reached a leaf (0 children) or a branch point (>1 children): + // end this chain here, and start a fresh chain at each child. + for (int c : children) + stack.push_back(c); + break; + } + + std::map counts; + for (int id : chain) + if (id < static_cast(raw.size()) && raw[id] != -1) + ++counts[raw[id]]; + + int majority = -1; + int best = 0; + for (const auto& kv : counts) + { + if (kv.second > best) + { + best = kv.second; + majority = kv.first; + } + } + + for (int id : chain) + m_nodeSegment[id] = majority; + } +} + +int SkeletonGraph::segmentOf(int nodeId) const +{ + if (nodeId < 0 || nodeId >= static_cast(m_nodeSegment.size())) + return -1; + return m_nodeSegment[nodeId]; +} + +std::vector SkeletonGraph::nodesInSegment(int segment) const +{ + std::vector result; + for (const SkeletonNode& n : m_nodes) + if (segmentOf(n.id()) == segment) + result.push_back(n.id()); + return result; +} + +void SkeletonGraph::exportSegmentReportCSV(const std::string& filename) const +{ + if (!hasRoot()) + return; + + std::ofstream out(filename, std::ofstream::out | std::ofstream::trunc); + out << "id,x,y,z,segment\n"; + out << std::fixed << std::setprecision(6); + + for (const SkeletonNode& n : m_nodes) + { + const auto& p = n.position(); + out << n.id() << "," + << p[0] << "," << p[1] << "," << p[2] << "," + << segmentOf(n.id()) << "\n"; + } +} + +void SkeletonGraph::exportSegmentReportCSV(const std::string& filename, const std::vector& segmentNames) const +{ + if (!hasRoot()) + return; + + std::set distinctSegments; + for (const SkeletonNode& n : m_nodes) + { + int seg = segmentOf(n.id()); + if (seg != -1) + distinctSegments.insert(seg); + } + + std::ofstream out(filename, std::ofstream::out | std::ofstream::trunc); + out << "# totalNodes=" << m_nodes.size() << ",distinctSegments=" << distinctSegments.size() << "\n"; + out << "id,x,y,z,segment\n"; + out << std::fixed << std::setprecision(6); + + for (const SkeletonNode& n : m_nodes) + { + const auto& p = n.position(); + int seg = segmentOf(n.id()); + + std::string name = "unknown"; + if (seg >= 0 && seg < static_cast(segmentNames.size())) + name = segmentNames[seg]; + + out << n.id() << "," + << p[0] << "," << p[1] << "," << p[2] << "," + << name << "\n"; + } +} + +std::vector SkeletonGraph::simulateResection(const std::vector& cutNodeIds) const +{ + std::vector affected; + if (!hasRoot()) + return affected; + + std::set cutSet(cutNodeIds.begin(), cutNodeIds.end()); + + // If the root itself is cut, nothing downstream can be perfused any more. + if (cutSet.count(m_rootId)) + { + affected.reserve(m_nodes.size()); + for (const SkeletonNode& n : m_nodes) + affected.push_back(n.id()); + return affected; + } + + // Reachability BFS over the raw connectivity graph (which already + // includes any collateral/loop edges): a cut node is simply never + // pushed/expanded, which is equivalent to deleting it and all of its + // edges from the graph. Anything still reachable from the root through + // some other path (e.g. an anastomosis bypassing the cut) stays + // perfused; everything else is affected. + std::vector reachable(m_nodes.size(), false); + std::queue q; + reachable[m_rootId] = true; + q.push(m_rootId); + + while (!q.empty()) + { + int u = q.front(); + q.pop(); + + auto it = m_adjacency.find(u); + if (it == m_adjacency.end()) + continue; + + for (int v : it->second) + { + if (cutSet.count(v) || reachable[v]) + continue; + reachable[v] = true; + q.push(v); + } + } + + for (const SkeletonNode& n : m_nodes) + if (cutSet.count(n.id()) || !reachable[n.id()]) + affected.push_back(n.id()); + + return affected; +} + +std::vector SkeletonGraph::affectedSegments(const std::vector& affectedNodeIds) const +{ + std::set segments; + for (int id : affectedNodeIds) + { + int seg = segmentOf(id); + if (seg != -1) + segments.insert(seg); + } + return std::vector(segments.begin(), segments.end()); +} + } // namespace meshskeletonizationplugin diff --git a/src/MeshSkeletonizationPlugin/SkeletonGraph/SkeletonGraph.h b/src/MeshSkeletonizationPlugin/SkeletonGraph/SkeletonGraph.h index e2d9e3f..200542f 100644 --- a/src/MeshSkeletonizationPlugin/SkeletonGraph/SkeletonGraph.h +++ b/src/MeshSkeletonizationPlugin/SkeletonGraph/SkeletonGraph.h @@ -24,6 +24,14 @@ class SkeletonGraph /// to populate parent/children ids on the nodes. bool loadFromFile(const std::string& filename, double mergeTolerance = 1e-6); + /// Loads a skeleton already exported as VTK POLYDATA by exportToVTK() (points + + /// "parent child" LINES + a "type" POINT_DATA scalar marking the root with value + /// 1). Unlike loadFromFile(), the tree structure (parent/children ids, root) is + /// read directly from the file instead of recomputed by buildTree(); other scalar + /// fields (depth, is_loop) are read past and ignored. Returns false if the file + /// could not be opened or no root ("type"==1) was found. + bool loadFromVTK(const std::string& filename); + /// Picks the node closest to entryPoint as root and (re)computes every /// node's parentIds/childrenIds via BFS over the raw connectivity. void buildTree(const std::array& entryPoint); @@ -41,7 +49,7 @@ class SkeletonGraph /// Ids of nodes connected by a loop/anastomosis (extra edge beyond the tree). const std::vector>& loopEdges() const { return m_loopEdges; } - //newest additions + //newzest additions const SkeletonNode* node(int nodeId) const; bool parentsOf(int nodeId, std::vector& parents) const; bool childrenOf(int nodeId, std::vector& children) const; @@ -51,6 +59,61 @@ class SkeletonGraph void exportToVTK(const std::string& filename) const; void clear(); + // --- Liver-segment mapping ------------------------------------------------- + + /// Assigns each node a liver-segment id (e.g. an index into a list of + /// Couinaud segments) from `rawNodeLabels`, indexed the same way as + /// nodes() (rawNodeLabels[nodeId], -1 = unknown/undetermined for that + /// node, e.g. because it fell outside every segment mesh). Typically + /// produced externally per node, e.g. via a point-in-mesh test against a + /// set of closed segment meshes (see SkeletonSegmentMapper). + /// + /// Raw per-node labels can be noisy near segment boundaries, so the tree + /// is split into maximal branches (runs of nodes between the root, + /// branch points, and leaves) and every node in a branch is set to that + /// branch's majority raw label. Requires the tree to have a root + /// (buildTree() or loadFromVTK()). Branches with no labeled node at all + /// get -1. + void assignSegmentLabels(const std::vector& rawNodeLabels); + + /// Segment id assigned to a node by the last assignSegmentLabels() call, + /// or -1 if it hasn't been run (or the node has no label). + int segmentOf(int nodeId) const; + + /// Ids of all nodes currently assigned to the given segment. + std::vector nodesInSegment(int segment) const; + + /// Writes a CSV report (id, x, y, z, segment) reflecting the last + /// assignSegmentLabels() call. + void exportSegmentReportCSV(const std::string& filename) const; + + /// Same as above, but resolves each segment id to a human-readable name + /// via `segmentNames` (segmentNames[segmentId]); ids without a matching + /// entry, and -1 (unknown), are written as "unknown". + void exportSegmentReportCSV(const std::string& filename, const std::vector& segmentNames) const; + + // --- Resection (devascularization) simulation ------------------------------- + + /// Simulates a resection (vessel transection) at each of `cutNodeIds`: + /// every cut node is fully removed from the connectivity graph (all of + /// its edges, both primary tree and any loop/anastomosis edges), then + /// reachability from the root is recomputed over what remains. This is + /// deliberately NOT the same as subtree(): a node below a cut can still + /// be reachable (and so still perfused) if a collateral/anastomosis + /// connects it back to the root through a path that avoids every cut + /// node - subtree() would incorrectly flag it as affected. + /// + /// Returns every node that can no longer reach the root (i.e. would lose + /// blood supply), plus the cut nodes themselves. If the root itself is + /// cut, every node is returned. Requires hasRoot(). + std::vector simulateResection(const std::vector& cutNodeIds) const; + + /// Convenience: the distinct liver segments touched by `affectedNodeIds` + /// (typically the result of simulateResection()) - i.e. the segments + /// that would lose blood supply. Requires assignSegmentLabels() to have + /// been called already; nodes with no segment (-1) are ignored. + std::vector affectedSegments(const std::vector& affectedNodeIds) const; + private: int findOrCreateNode(const std::array& p, double tol); int closestNodeId(const std::array& p) const; @@ -60,6 +123,10 @@ class SkeletonGraph std::map> m_adjacency; ///< raw undirected connectivity from the file std::vector> m_loopEdges; int m_rootId{ -1 }; + + /// Segment id per node (parallel to m_nodes, indexed by node id), filled by + /// assignSegmentLabels(). Empty until that has been called at least once. + std::vector m_nodeSegment; }; } // namespace meshskeletonizationplugin diff --git a/src/MeshSkeletonizationPlugin/SkeletonGraph/SkeletonNode.h b/src/MeshSkeletonizationPlugin/SkeletonGraph/SkeletonNode.h index 6deba3d..514e199 100644 --- a/src/MeshSkeletonizationPlugin/SkeletonGraph/SkeletonNode.h +++ b/src/MeshSkeletonizationPlugin/SkeletonGraph/SkeletonNode.h @@ -38,8 +38,6 @@ class SkeletonNode bool isLeaf() const { return m_childrenIds.empty(); } bool isBranchPoint() const { return m_childrenIds.size() > 1; } - // --- Optional relation to the input (vessel) mesh this skeleton came from --- - /// Index of the closest vertex in the input mesh, -1 if not computed. int meshVertexId() const { return m_meshVertexId; } void setMeshVertexId(int vId) { m_meshVertexId = vId; } diff --git a/src/MeshSkeletonizationPlugin/SkeletonGraph/SkeletonReader.inl b/src/MeshSkeletonizationPlugin/SkeletonGraph/SkeletonReader.inl index 30cf5ad..f6d360b 100644 --- a/src/MeshSkeletonizationPlugin/SkeletonGraph/SkeletonReader.inl +++ b/src/MeshSkeletonizationPlugin/SkeletonGraph/SkeletonReader.inl @@ -1,12 +1,11 @@ #pragma once #include -using namespace sofa::core::objectmodel; +using namespace sofa::core::objectmodel; namespace meshskeletonizationplugin { - template SkeletonReader::SkeletonReader() : d_inSkeletonFilename(initData(&d_inSkeletonFilename, "filename", "Skeleton polyline file to read (e.g. skeleton.txt)")) @@ -19,7 +18,6 @@ SkeletonReader::SkeletonReader() addInput(&d_inSkeletonFilename); addInput(&d_inVertices); addInput(&d_inEntryPoint); - addOutput(&d_outVTKFilename); addOutput(&d_outReportFilename); addOutput(&d_outNodeCount); @@ -42,7 +40,8 @@ void SkeletonReader::init() template void SkeletonReader::doUpdate() { - if (d_inSkeletonFilename.getFullPath().empty()){ + if (d_inSkeletonFilename.getFullPath().empty()) + { d_componentState.setValue(ComponentState::Invalid); return; } @@ -62,14 +61,15 @@ void SkeletonReader::doUpdate() if (!m_graph.hasRoot()) { + // The file parsed but produced no usable rooted tree, so every + // downstream output would be empty. msg_error() << "No root could be established from: " << d_inSkeletonFilename.getFullPath(); d_componentState.setValue(ComponentState::Invalid); return; } - if (m_graph.hasRoot()) - msg_info() << "Tree built, root id " << m_graph.rootId() - << ", " << m_graph.loopEdges().size() << " loop edge(s) detected."; + msg_info() << "Tree built, root id " << m_graph.rootId() + << ", " << m_graph.loopEdges().size() << " loop edge(s) detected."; if (!d_inVertices.getValue().empty()) { @@ -82,18 +82,22 @@ void SkeletonReader::doUpdate() if (d_outVTKFilename.isSet()) m_graph.exportToVTK(d_outVTKFilename.getFullPath()); - + if (d_outReportFilename.isSet()) m_graph.exportReportCSV(d_outReportFilename.getFullPath()); + + d_componentState.setValue(ComponentState::Valid); } template void SkeletonReader::draw(const sofa::core::visual::VisualParams* vparams) { + if (d_componentState.getValue() != ComponentState::Valid) + return; + using Color = sofa::type::RGBAColor; std::vector< type::Vec3 > dvec; - for (const SkeletonNode& node : m_graph.nodes()) { for (int childId : node.childrenIds()) @@ -104,7 +108,6 @@ void SkeletonReader::draw(const sofa::core::visual::VisualParams* vpa dvec.emplace_back(Coord(p0[0], p0[1], p0[2])); dvec.emplace_back(Coord(p1[0], p1[1], p1[2])); - vparams->drawTool()->drawLines(dvec, 2, Color::blue()); dvec.clear(); } diff --git a/src/MeshSkeletonizationPlugin/SkeletonResection/SkeletonResectionSimulator.cpp b/src/MeshSkeletonizationPlugin/SkeletonResection/SkeletonResectionSimulator.cpp new file mode 100644 index 0000000..88e59be --- /dev/null +++ b/src/MeshSkeletonizationPlugin/SkeletonResection/SkeletonResectionSimulator.cpp @@ -0,0 +1,22 @@ +#define SKELETONRESECTIONSIMULATOR_CPP +#include + +#include +#include + +namespace meshskeletonizationplugin +{ +using namespace sofa::defaulttype; + +void registerSkeletonResectionSimulator(sofa::core::ObjectFactory* factory) +{ + factory->registerObjects(sofa::core::ObjectRegistrationData( + "Simulates a vessel resection at one or more skeleton nodes and reports which " + "nodes and liver segments would lose blood supply, accounting for collateral " + "(loop/anastomosis) paths that bypass the cut") + .add< SkeletonResectionSimulator >()); +} + +template class SOFA_MESHSKELETONIZATIONPLUGIN_API SkeletonResectionSimulator; + +} // namespace meshskeletonizationplugin diff --git a/src/MeshSkeletonizationPlugin/SkeletonResection/SkeletonResectionSimulator.h b/src/MeshSkeletonizationPlugin/SkeletonResection/SkeletonResectionSimulator.h new file mode 100644 index 0000000..484a811 --- /dev/null +++ b/src/MeshSkeletonizationPlugin/SkeletonResection/SkeletonResectionSimulator.h @@ -0,0 +1,95 @@ +#pragma once + +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +using namespace sofa; +using namespace sofa::defaulttype; + +namespace meshskeletonizationplugin +{ + +template +class SkeletonResectionSimulator : public sofa::core::DataEngine +{ +public: + SOFA_CLASS(SOFA_TEMPLATE(SkeletonResectionSimulator, DataTypes), sofa::core::DataEngine); + + typedef typename DataTypes::Coord Coord; + typedef typename DataTypes::VecCoord VecCoord; + typedef typename Coord::value_type Real; + typedef type::Vec<3, Real> Vec3; + + /// The segmented skeleton to run resection scenarios against. If left + /// unset, the first SkeletonSegmentMapper found in the context is used. + sofa::core::objectmodel::SingleLink< + SkeletonResectionSimulator, + SkeletonSegmentMapper, + sofa::core::objectmodel::BaseLink::FLAG_STOREPATH | sofa::core::objectmodel::BaseLink::FLAG_STRONGLINK> + l_segmentMapper; + + /// Optional: the actual visual object (e.g. OglModel) for each segment, + /// in the same order as the linked mapper's segmentNames. If set, their + /// "color" Data is pushed directly every draw() call - more reliable + /// than a scene-level Data link, since most visual models don't re-read + /// a linked color every frame on their own. + sofa::core::objectmodel::MultiLink< + SkeletonResectionSimulator, + sofa::core::objectmodel::BaseObject, + sofa::core::objectmodel::BaseLink::FLAG_STOREPATH> + l_segmentVisualModels; + + // Inputs + sofa::core::objectmodel::Data> d_inCutNodeIds; ///< Node id(s) where the vessel is (candidate) severed + sofa::core::objectmodel::Data d_activationDelay; ///< Seconds after simulation start before affected nodes/segments actually switch color (default 5.0); everything shows as perfused (green) until then + sofa::core::objectmodel::DataFileName d_outReportFilename; ///< Optional CSV report: id, x, y, z, segment, for affected nodes only + + // Outputs + sofa::core::objectmodel::Data> d_outAffectedNodeIds; ///< Node ids that would lose blood supply (includes the cut nodes themselves) + sofa::core::objectmodel::Data d_outAffectedNodeCount; ///< Convenience: size of d_outAffectedNodeIds + sofa::core::objectmodel::Data> d_outAffectedSegmentIds; ///< Distinct segment indices touched by the affected nodes + sofa::core::objectmodel::Data> d_outAffectedSegmentNames; ///< Same, resolved to names via the linked SkeletonSegmentMapper's segmentNames + sofa::core::objectmodel::Data> d_outSegmentColors; ///< One color per segment, same order as the linked mapper's segmentNames - green until activationDelay elapses, then red for affected segments. Link an OglModel's "color" to e.g. "@resectionSim.segmentColors[0]" per segment (index syntax support depends on SOFA version). + + void init() override; + void doUpdate() override; + void draw(const sofa::core::visual::VisualParams* vparams) override; + +private: + SkeletonResectionSimulator(); + virtual ~SkeletonResectionSimulator() = default; + + /// Recomputes m_colorsActive from the current simulation time and + /// refreshes d_outSegmentColors accordingly. Called every draw() so the + /// activationDelay threshold takes effect without needing an event + /// listener (and the Sofa.Simulation.Core dependency that would add). + void updateSegmentColors(); + + /// Local working copy of the linked mapper's (already segmented) graph. + SkeletonGraph m_graph; + + /// Whether affected nodes/segments should currently render as + /// affected (red/black) rather than perfused (green) - false until + /// d_activationDelay seconds have elapsed since simulation start. + bool m_colorsActive{ false }; + + /// Guards so diagnostic messages below are logged once, not every frame. + bool m_loggedVisualModelLinkStatus{ false }; + std::vector m_loggedMissingColorData; +}; + +#if !defined(SKELETONRESECTIONSIMULATOR_CPP) +extern template class SOFA_MESHSKELETONIZATIONPLUGIN_API SkeletonResectionSimulator; +#endif + +} // namespace meshskeletonizationplugin \ No newline at end of file diff --git a/src/MeshSkeletonizationPlugin/SkeletonResection/SkeletonResectionSimulator.inl b/src/MeshSkeletonizationPlugin/SkeletonResection/SkeletonResectionSimulator.inl new file mode 100644 index 0000000..f7ac96b --- /dev/null +++ b/src/MeshSkeletonizationPlugin/SkeletonResection/SkeletonResectionSimulator.inl @@ -0,0 +1,269 @@ +#pragma once +#include + +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +namespace meshskeletonizationplugin +{ + +template +SkeletonResectionSimulator::SkeletonResectionSimulator() + : l_segmentMapper(initLink("segmentMapper", + "Segmented skeleton to run resection scenarios against; if empty, the first " + "SkeletonSegmentMapper found in the scene is used")) + , l_segmentVisualModels(initLink("segmentVisualModels", + "Optional: the visual object (e.g. OglModel) for each segment, same order as " + "segmentNames. Its 'color' Data is pushed directly every frame - more reliable " + "than linking a color attribute to segmentColors[i] in the scene.")) + , d_inCutNodeIds(initData(&d_inCutNodeIds, "cutNodeIds", + "Node id(s) where the vessel is (candidate) severed")) + , d_activationDelay(initData(&d_activationDelay, 5.0, "activationDelay", + "Seconds after simulation start before affected nodes/segments actually switch " + "color to red/black; everything renders as perfused (green) until then")) + , d_outReportFilename(initData(&d_outReportFilename, "outReportFilename", + "Optional path to export a CSV report (id, x, y, z, segment) of the affected nodes")) + , d_outAffectedNodeIds(initData(&d_outAffectedNodeIds, "affectedNodeIds", + "Node ids that would lose blood supply (includes the cut nodes themselves)")) + , d_outAffectedNodeCount(initData(&d_outAffectedNodeCount, 0, "affectedNodeCount", + "Convenience: size of affectedNodeIds")) + , d_outAffectedSegmentIds(initData(&d_outAffectedSegmentIds, "affectedSegmentIds", + "Distinct segment indices touched by the affected nodes")) + , d_outAffectedSegmentNames(initData(&d_outAffectedSegmentNames, "affectedSegmentNames", + "Same as affectedSegmentIds, resolved to names")) + , d_outSegmentColors(initData(&d_outSegmentColors, "segmentColors", + "One color per segment (same order as the linked mapper's segmentNames): green " + "until activationDelay elapses, then red for affected segments. Link an OglModel's " + "'color' Data to e.g. '@resectionSim.segmentColors[0]' for the first segment.")) +{ + addInput(&d_inCutNodeIds); + addOutput(&d_outAffectedNodeIds); + addOutput(&d_outAffectedNodeCount); + addOutput(&d_outAffectedSegmentIds); + addOutput(&d_outAffectedSegmentNames); + addOutput(&d_outSegmentColors); +} + +template +void SkeletonResectionSimulator::init() +{ + if (!l_segmentMapper) + { + l_segmentMapper.set(this->getContext()->template get>()); + } + if (!l_segmentMapper) + { + msg_error() << "No SkeletonSegmentMapper found. Set 'segmentMapper' to a " + "valid component path, or add one earlier in the scene."; + } + + setDirtyValue(); + update(); // computes each time without waiting for something to read an output Data +} + +template +void SkeletonResectionSimulator::doUpdate() +{ + if (!l_segmentMapper) + return; + // Working on copy of the already segmented graph so the candidate is re-run per candidate without touching the mapper + m_graph = l_segmentMapper->graph(); + + const auto& cutIds = d_inCutNodeIds.getValue(); + std::vector affected = m_graph.simulateResection(cutIds); + std::vector segments = m_graph.affectedSegments(affected); + + const auto& segmentNames = l_segmentMapper->d_inSegmentNames.getValue(); + std::vector segmentNameStrs; + segmentNameStrs.reserve(segments.size()); + for (int seg : segments) + segmentNameStrs.push_back(seg >= 0 && seg < static_cast(segmentNames.size()) + ? segmentNames[seg] : "unknown"); + + d_outAffectedNodeIds.setValue(affected); + d_outAffectedNodeCount.setValue(static_cast(affected.size())); + d_outAffectedSegmentIds.setValue(segments); + d_outAffectedSegmentNames.setValue(segmentNameStrs); + + updateSegmentColors(); // initial colors; draw() keeps these fresh every frame after this + + if (d_outReportFilename.isSet() && !d_outReportFilename.getFullPath().empty()) + { + std::ofstream out(d_outReportFilename.getFullPath(), std::ofstream::out | std::ofstream::trunc); + out << "# cutNodeIds=" << cutIds.size() + << ",affectedNodeCount=" << affected.size() + << ",affectedSegmentCount=" << segments.size() + << ",affectedSegments="; + for (std::size_t i = 0; i < segmentNameStrs.size(); ++i) + out << (i ? ";" : "") << segmentNameStrs[i]; + out << "\n"; + out << "id,x,y,z,segment\n"; + out << std::fixed << std::setprecision(6); + for (int id : affected) + { + const SkeletonNode* n = m_graph.node(id); + if (!n) + continue; + const auto& p = n->position(); + int seg = m_graph.segmentOf(id); + std::string segName = (seg >= 0 && seg < static_cast(segmentNames.size())) ? segmentNames[seg] : "unknown"; + out << id << "," << p[0] << "," << p[1] << "," << p[2] << "," << segName << "\n"; + } + } +} + +template +void SkeletonResectionSimulator::updateSegmentColors() +{ + // colors of the segments before and after (green --> red) + static const sofa::type::RGBAColor safeColor(0.20f, 0.80f, 0.20f, 1.0f); + static const sofa::type::RGBAColor affectedColor(0.90f, 0.10f, 0.10f, 1.0f); + + m_colorsActive = (this->getContext()->getTime() >= d_activationDelay.getValue()); + + if (!l_segmentMapper) + return; + + // --- Diagnostics: run once, so we can see exactly what's wrong instead of silently doing nothing. + if (!m_loggedVisualModelLinkStatus) + { + m_loggedVisualModelLinkStatus = true; + if (l_segmentVisualModels.empty()) + { + msg_warning() << "segmentVisualModels is empty: no visual objects linked, so no " + "mesh color will ever change. Check the 'segmentVisualModels' " + "attribute/paths in the scene."; + } + else + { + m_loggedMissingColorData.assign(l_segmentVisualModels.size(), false); + for (std::size_t i = 0; i < l_segmentVisualModels.size(); ++i) + { + sofa::core::objectmodel::BaseObject* obj = l_segmentVisualModels.get(i); + if (!obj) + { + msg_warning() << "segmentVisualModels[" << i << "] did not resolve to any object " + "(bad path?)."; + continue; + } + msg_info() << "segmentVisualModels[" << i << "] resolved to '" << obj->getName() + << "' (" << obj->getClassName() << ")."; + } + } + } + + const auto& segmentNames = l_segmentMapper->d_inSegmentNames.getValue(); + const auto& segments = d_outAffectedSegmentIds.getValue(); + std::set affectedSegSet(segments.begin(), segments.end()); + + sofa::helper::WriteAccessor>> colors(d_outSegmentColors); + colors.resize(segmentNames.size()); + for (std::size_t i = 0; i < segmentNames.size(); ++i) + colors[i] = (m_colorsActive && affectedSegSet.count(static_cast(i))) ? affectedColor : safeColor; + + // Push directly into each linked visual model's own "color" Data, if + // provided - this is what actually makes the mesh repaint, since a + // passive scene-level Data link to segmentColors[i] isn't guaranteed to + // be re-read every frame by the visual model itself. + for (std::size_t i = 0; i < l_segmentVisualModels.size() && i < colors.size(); ++i) + { + sofa::core::objectmodel::BaseObject* obj = l_segmentVisualModels.get(i); + if (!obj) + continue; + + // OglModel doesn't expose a plain "color" Data that is why this method was used + sofa::core::objectmodel::BaseData* materialData = obj->findData("material"); + if (!materialData) + { + if (i < m_loggedMissingColorData.size() && !m_loggedMissingColorData[i]) + { + m_loggedMissingColorData[i] = true; + std::ostringstream fields; + for (sofa::core::objectmodel::BaseData* d : obj->getDataFields()) + fields << d->getName() << " "; + msg_warning() << "'" << obj->getName() << "' (" << obj->getClassName() + << ") has no Data named 'material' either. Its actual Data fields are: " + << fields.str(); + } + continue; + } + + std::istringstream iss(materialData->getValueString()); + std::vector tokens{ std::istream_iterator(iss), std::istream_iterator() }; + + auto diffuseIt = std::find(tokens.begin(), tokens.end(), "Diffuse"); + if (diffuseIt == tokens.end() || std::distance(diffuseIt, tokens.end()) < 6) + { + if (i < m_loggedMissingColorData.size() && !m_loggedMissingColorData[i]) + { + m_loggedMissingColorData[i] = true; + msg_warning() << "Could not find a 'Diffuse r g b a' pattern in '" + << obj->getName() << "'s material string: " << materialData->getValueString(); + } + continue; + } + + const auto& c = colors[i]; + *(diffuseIt + 1) = "1"; // force useDiffuse on, so our color is actually applied + std::ostringstream rs, gs, bs, as; + rs << c[0]; gs << c[1]; bs << c[2]; as << c[3]; + *(diffuseIt + 2) = rs.str(); + *(diffuseIt + 3) = gs.str(); + *(diffuseIt + 4) = bs.str(); + *(diffuseIt + 5) = as.str(); + + std::ostringstream newStr; + for (std::size_t k = 0; k < tokens.size(); ++k) + newStr << (k ? " " : "") << tokens[k]; + materialData->read(newStr.str()); + } +} + +template +void SkeletonResectionSimulator::draw(const sofa::core::visual::VisualParams* vparams) +{ + updateSegmentColors(); // keeps segmentColors (and m_colorsActive) current every frame + + if (!vparams->displayFlags().getShowBehaviorModels()) + return; + + static const sofa::type::RGBAColor perfusedColor(0.20f, 0.80f, 0.20f, 1.0f); + static const sofa::type::RGBAColor affectedColor(0.0f, 0.0f, 0.0f, 1.0f); + static const sofa::type::RGBAColor cutColor(1.0f, 1.0f, 1.0f, 1.0f); + + const auto& affected = d_outAffectedNodeIds.getValue(); + std::set affectedSet(affected.begin(), affected.end()); + std::set cutSet(d_inCutNodeIds.getValue().begin(), d_inCutNodeIds.getValue().end()); + + std::vector points; + std::vector colors; + points.reserve(m_graph.nodes().size()); + colors.reserve(m_graph.nodes().size()); + + for (const SkeletonNode& n : m_graph.nodes()) + { + const auto& p = n.position(); + points.emplace_back(p[0], p[1], p[2]); + + if (!m_colorsActive) + colors.push_back(perfusedColor); // before activationDelay: everything shows as perfused + else if (cutSet.count(n.id())) + colors.push_back(cutColor); + else if (affectedSet.count(n.id())) + colors.push_back(affectedColor); + else + colors.push_back(perfusedColor); + } + + vparams->drawTool()->drawPoints(points, 7.0f, colors); +} + +} // namespace meshskeletonizationplugin diff --git a/src/MeshSkeletonizationPlugin/init.cpp b/src/MeshSkeletonizationPlugin/init.cpp index 6527428..4c994ef 100644 --- a/src/MeshSkeletonizationPlugin/init.cpp +++ b/src/MeshSkeletonizationPlugin/init.cpp @@ -30,6 +30,8 @@ namespace meshskeletonizationplugin extern void registerMeshSkeletonization(sofa::core::ObjectFactory* factory); //extern void registerSkeletonizationLoader(sofa::core::ObjectFactory* factory); extern void registerSkeletonReader(sofa::core::ObjectFactory* factory); +extern void registerSkeletonSegmentMapper(sofa::core::ObjectFactory* factory); +extern void registerSkeletonResectionSimulator(sofa::core::ObjectFactory* factory); //Here are just several convenient functions to help users know what the plugin contains extern "C" { @@ -82,7 +84,9 @@ void registerObjects(sofa::core::ObjectFactory* factory) { registerMeshSkeletonization(factory); // registerSkeletonizationLoader(factory); - registerSkeletonReader(factory); + registerSkeletonReader(factory); + registerSkeletonSegmentMapper(factory); + registerSkeletonResectionSimulator(factory); } }