[Skeletonization+ SOFA] Liver resection simulation: segment mapping + candidate-cut analysis

.gitignore

@@ -5,6 +5,5 @@
 *.vtk
 *.obj
 *.sxn
-*.scn
 *.csv
 

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)

TestScene.scn

@@ -0,0 +1,80 @@
+<?xml version="1.0"?>
+<Node name="root" dt="0.01" gravity="0 0 0">
+    <RequiredPlugin pluginName="MeshSkeletonizationPlugin"/>
+    <RequiredPlugin pluginName="Sofa.GL.Component.Rendering3D"/>
+    <RequiredPlugin pluginName="Sofa.Component.SceneUtility"/>
+    <RequiredPlugin pluginName="Sofa.Component.IO.Mesh"/>
+    <RequiredPlugin pluginName="Sofa.Component.Visual"/>
+
+    <!--
+      ONE VisualStyle for the whole scene. showBehaviorModels is required for
+      SkeletonReader/SkeletonSegmentMapper/SkeletonResectionSimulator's draw()
+      to render anything - without it their colored skeleton points are silently
+      skipped even though the mesh itself (showVisualModels) still shows up.
+    -->
+    <VisualStyle displayFlags="showVisualModels showBehaviorModels showWireframe"/>
+
+    <AddDataRepository path="../data/"/>
+    <AddDataRepository path="./data/"/>
+
+    <DefaultAnimationLoop/>
+    <DefaultVisualManagerLoop/>
+
+    <Node name="Mesh">
+        <MeshOBJLoader name="loader" filename="patient_08\visual/Patient_08_venoussystem_visual.obj" />
+        <MeshSkeletonization template="Vec3d" name="skel"
+            inputVertices="@loader.position"
+            inputTriangles="@loader.triangles"
+            outputSkeleton="\patient_08\Patient_08_venoussystem_skeleton.txt" />
+
+        <Node name="visual">
+            <OglModel src="@../loader"/>
+        </Node>
+
+        <SkeletonReader name="reader"
+            filename="./patient_08\Patient_08_venoussystem_skeleton.txt"
+            inputVertices="@loader.position"
+            outputVTK="./patient_08\Patient_08_venoussystem_skeleton.vtk"
+            outputReport="Results\Patient_08_venoussystem_report_skeleton_mapper.csv" />
+    </Node>
+
+    <Node name="LiverSegments">
+        <!--
+          segmentNames must stay in the SAME order as these loaders, since
+          that order is what maps index 0/1/2... to "II"/"III"/"IVa"... below.
+        -->
+        <MeshOBJLoader name="segII"   filename="patient_08\collision/couinaud_Segment_II____Lat__Sup_collision.obj"/>
+        <MeshOBJLoader name="segIII"  filename="patient_08\collision/couinaud_Segment_III___Lat__Inf_collision.obj"/>
+        <MeshOBJLoader name="segIVa"  filename="patient_08\collision/couinaud_Segment_IVa___Med__Sup_collision.obj"/>
+        <MeshOBJLoader name="segIVb"  filename="patient_08\collision/couinaud_Segment_IVb___Med__Inf_collision.obj"/>
+        <MeshOBJLoader name="segV"    filename="patient_08\collision/couinaud_Segment_V_____Ant__Inf__R_collision.obj"/>
+        <MeshOBJLoader name="segVI"   filename="patient_08\collision/couinaud_Segment_VI____Post__Inf__R_collision.obj"/>
+        <MeshOBJLoader name="segVII"  filename="patient_08\collision/couinaud_Segment_VII___Post__Sup__R_collision.obj"/>
+        <MeshOBJLoader name="segVIII" filename="patient_08\collision/couinaud_Segment_VIII___Ant__Sup__R_collision.obj"/>
+
+
+        <Node name="visII"><OglModel name="oglII" src="@../segII"/></Node>
+        <Node name="visIII"><OglModel name="oglIII" src="@../segIII"/></Node>
+        <Node name="visIVa"><OglModel name="oglIVa" src="@../segIVa"/></Node>
+        <Node name="visIVb"><OglModel name="oglIVb" src="@../segIVb"/></Node>
+        <Node name="visV"><OglModel name="oglV" src="@../segV"/></Node>
+        <Node name="visVI"><OglModel name="oglVI" src="@../segVI"/></Node>
+        <Node name="visVII"><OglModel name="oglVII" src="@../segVII"/></Node>
+        <Node name="visVIII"><OglModel name="oglVIII" src="@../segVIII"/></Node>
+
+        <SkeletonSegmentMapper name="segmentMapper" template="Vec3d"
+            skeletonReader="@../Mesh/reader"
+            segmentMeshes="@segII @segIII @segIVa @segIVb @segV @segVI @segVII @segVIII"
+            segmentNames="II III IVa IVb V VI VII VIII"
+            outSegmentReportFilename="\Results\Patient_08_liver_segment_report_before.csv"/>
+
+
+        <SkeletonResectionSimulator name="resectionSim" template="Vec3d"
+            segmentMapper="@segmentMapper"
+            cutNodeIds="135 100 150 1"
+            activationDelay="5.0"
+            segmentVisualModels="@visII/oglII @visIII/oglIII @visIVa/oglIVa @visIVb/oglIVb @visV/oglV @visVI/oglVI @visVII/oglVII @visVIII/oglVIII"
+            outReportFilename="\Results\Patient_08_resection_report_after.csv"/>
+    </Node>
+
+</Node>

scenes/BeamFEMSkeleton.scn

@@ -12,12 +12,7 @@
 		<MeshSkeletonization template="Vec3d" name="skel" inputVertices="@myLoader.position" inputTriangles="@myLoader.triangles" />
 	</Node>	
 
-	<!-- NOTE: 
-	The position and lines should be passed directly from MeshSkeletonization to BeamFEMForceField.	
-	I will implement a quick Engine to do that.
-	In the mean time use generated files and copy/paste them.
-	You may have to order the points and edges before feeding them to BeamFEMForceField.
-	 -->
+
 	<Node name="Beams" activated="1">    		
 	        <MechanicalObject template="Rigid3d" name="BeamDOFs" position="0.00671732 0.00198818 0.200274 0 0 0 1 0.00525121 0.00274144 0.200697 0 0 0 1 0.0144205 0.000200659 0.196766 0 0 0 1 0.00381526 0.00325765 0.201011 0 0 0 1 -0.000413926 0.00384828 0.201621 0 0 0 1 -0.00269336 0.00367566 0.201811 0 0 0 1 -0.0162832 0.0218957 0.194973 0 0 0 1 -0.0157561 0.0209032 0.195621 0 0 0 1 -0.0148986 0.0193463 0.196598 0 0 0 1 -0.013763 0.01737 0.197801 0 0 0 1 -0.0118151 0.0142045 0.199543 0 0 0 1 -0.0127952 -0.00277541 0.202707 0 0 0 1 -0.0179481 -0.00390045 0.201449 0 0 0 1 -0.0100881 -0.000453375 0.20268 0 0 0 1 -0.0191585 -0.00378408 0.201018 0 0 0 1 -0.0198925 -0.00365268 0.200709 0 0 0 1 -0.00592703 0.00263035 0.20215 0 0 0 1 -0.00484675 0.00326465 0.201986 0 0 0 1 -0.0110013 0.0129638 0.200084 0 0 0 1 -0.00679096 0.00660121 0.201484 0 0 0 1 -0.00584271 0.00502494 0.201721 0 0 0 1 0.00970073 0.000678009 0.199255 0 0 0 1 0.0084043 0.00112596 0.199779 0 0 0 1 0.011668 0.000427765 0.198275 0 0 0 1 0.00241641 0.00359541 0.201263 0 0 0 1 -0.00201755 0.00384477 0.201749 0 0 0 1 -0.0143672 -0.00409901 0.202465 0 0 0 1 -0.0114475 -0.00163015 0.202736 0 0 0 1 -0.018478 -0.00390527 0.201304 0 0 0 1 -0.00707019 0.00210355 0.202297 0 0 0 1 -0.00991798 0.011348 0.200619 0 0 0 1 0.0129828 0.000301091 0.197578 0 0 0 1 0.0138834 0.000232424 0.197076 0 0 0 1 0.00247815 0.00358254 0.201252 0 0 0 1 0.00104453 0.00377688 0.20146 0 0 0 1 -0.0148784 0.0193105 0.19662 0 0 0 1 -0.012762 0.0157058 0.198763 0 0 0 1 -0.0114838 -0.00166135 0.202737 0 0 0 1 -0.00850305 0.000904422 0.202519 0 0 0 1 -0.020055 -0.00362431 0.200637 0 0 0 1 -0.00790542 0.00833984 0.20123 0 0 0 1 0.0139482 0.000228183 0.197039 0 0 0 1 0.00394671 0.00321437 0.200984 0 0 0 1 -0.00846776 0.000934929 0.202514 0 0 0 1 -0.00784623 0.00824833 0.201244 0 0 0 1 -0.00194938 0.00384693 0.201744 0 0 0 1 -0.0117845 0.0141575 0.199564 0 0 0 1 -0.0159144 -0.00385712 0.202004 0 0 0 1 -0.014472 -0.00557562 0.202643 0 0 0 1 -0.0176483 -0.00389706 0.201531 0 0 0 1 -0.00895798 0.00991337 0.200912 0 0 0 1 0.00925139 0.000830593 0.199438 0 0 0 1 -0.0173011 -0.00389325 0.201627 0 0 0 1 " />
 	        <MeshTopology name="lines" points="@BeamDOFs.position" lines="51 22 22 0 0 1 51 21 42 3 1 42 21 23 3 33 23 31 41 32 41 2 33 24 24 34 34 4 4 45 45 25 25 5 5 17 6 7 7 8 8 35 35 9 9 36 36 10 46 10 26 47 26 48 26 11 52 49 47 52 49 12 37 27 11 37 28 12 28 14 15 14 43 38 13 38 43 29 16 17 18 30 30 50 44 19 44 40 19 20 20 17 31 32 18 46 13 27 29 16 15 39 50 40 "/>

src/MeshSkeletonizationPlugin/SegmentMapping/SkeletonSegmentMapper.cpp

@@ -0,0 +1,22 @@
+#define SKELETONSEGMENTMAPPER_CPP
+#include <MeshSkeletonizationPlugin/SegmentMapping/SkeletonSegmentMapper.inl>
+
+#include <sofa/core/ObjectFactory.h>
+#include <sofa/defaulttype/VecTypes.h>
+
+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<Vec3Types> >());
+}
+
+template class SOFA_MESHSKELETONIZATIONPLUGIN_API SkeletonSegmentMapper<Vec3Types>;
+
+} // namespace meshskeletonizationplugin

src/MeshSkeletonizationPlugin/SegmentMapping/SkeletonSegmentMapper.h

@@ -0,0 +1,79 @@
+#pragma once
+
+#include <MeshSkeletonizationPlugin/config.h>
+
+#include <sofa/type/Vec.h>
+#include <sofa/defaulttype/VecTypes.h>
+#include <sofa/core/DataEngine.h>
+#include <sofa/core/objectmodel/DataFileName.h>
+#include <sofa/core/objectmodel/Link.h>
+#include <sofa/core/visual/VisualParams.h>
+#include <sofa/core/loader/MeshLoader.h>
+
+#include <MeshSkeletonizationPlugin/SkeletonGraph/SkeletonGraph.h>
+#include <MeshSkeletonizationPlugin/SkeletonGraph/SkeletonReader.h>
+
+using namespace sofa;
+using namespace sofa::defaulttype;
+
+namespace meshskeletonizationplugin
+{
+
+template <class DataTypes>
+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<DataTypes> found in the context is used.
+    sofa::core::objectmodel::SingleLink<
+        SkeletonSegmentMapper<DataTypes>,
+        SkeletonReader<DataTypes>,
+        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<DataTypes>,
+        sofa::core::loader::MeshLoader,
+        sofa::core::objectmodel::BaseLink::FLAG_STOREPATH>
+        l_segmentMeshes;
+
+    //Component parameters
+    // Inputs
+    sofa::core::objectmodel::Data<sofa::type::vector<std::string>> 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<sofa::type::vector<int>> d_outNodeSegments;         ///< Segment index per skeleton node (into l_segmentMeshes/d_inSegmentNames, -1 = unknown), indexed like graph().nodes()
+    sofa::core::objectmodel::Data<sofa::type::vector<std::string>> d_outNodeSegmentNames; ///< Same as d_outNodeSegments, resolved to names ("unknown" if -1)
+    sofa::core::objectmodel::Data<int> 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<SkeletonSegmentMapper<DataTypes>>()->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<defaulttype::Vec3Types>;
+#endif
+
+} // namespace meshskeletonizationplugin