[SOLVED] Modelling tissue deformation in needle insertion case

[fedev]

Dear all,

For an university project, I would like to simulate the case of needle insertion inside the body (a phantom model). For simplicity, my scene consists only on the outer most layer of the phantom together with a needle modeled as a syringe at the moment. The deformable tissue stay still on a workstation meanwhile a rigid needle moves towards that and should penetrate after an initial collision and deformation of the tissue.
My python scene is the following :

import Sofa
import os
from stlib.physics.rigid import Floor

class MainScene(Sofa.PythonScriptController):

    def createGraph(self, node):

        dt = 0.02 # In second
        displayFlags =  ['showForceFields', 'showCollisionModels', 'showBehavior']
        
        self.rootNode = node.getRoot()
        self.rootNode.dt = dt
        self.rootNode.gravity = [0, 0, 0]

        self.rootNode.createObject('RequiredPlugin', name='SofaMiscCollision')
        self.rootNode.createObject('RequiredPlugin', name='SofaPython')
        self.rootNode.createObject('RequiredPlugin', name='CImgPlugin')
        self.rootNode.createObject('RequiredPlugin', name='SofaOpenglVisual')

        self.rootNode.createObject('VisualStyle', displayFlags=displayFlags)
        self.rootNode.createObject('FreeMotionAnimationLoop')
        
        self.rootNode.createObject('GenericConstraintSolver', maxIterations = 1000, tolerance = 1e-07)
        self.rootNode.createObject('DefaultPipeline', verbose='0', depth="6", draw='1')
        self.rootNode.createObject('BruteForceDetection', name = 'N2')
        #self.rootNode.createObject('DiscreteIntersection', name='Intersection')
        self.rootNode.createObject("LocalMinDistance", name="Intersection", alarmDistance="0.3", contactDistance="0.1", useLMDFilters="0")
        self.rootNode.createObject('DefaultContactManager', name="Response", response="FrictionContact", responseParams='mu=0.8')

        #==========================
        # FAT NODE
        #==========================
                        
        meshFile = "data/phantom/Segmentation_Mar09_Fat.stl"
                    
        fatNode = self.rootNode.createChild("Fat")

        fatNode.createObject('EulerImplicitSolver', rayleighStiffness = 0.1, rayleighMass = 0.1)
        fatNode.createObject('CGLinearSolver', iterations=25, tolerance=1e-9, threshold=1e-9, printLog=True)
        fatNode.createObject('MechanicalObject', name = 'mechObject', template="Vec3d", dx="0", dy="0", dz="0", rx="0", ry="0", rz="0", scale="1.0")
        fatNode.createObject('UniformMass', template="Vec3d,double", name='mass', totalMass="1")
                    
        fatNode.createObject('RegularGridTopology', nx = 10, ny = 10, nz = 10, xmin=-155, xmax=125, ymin=-170, ymax=45, zmin=-120, zmax=10)
        #fatNode.createObject('RegularGridSpringForceField', name="Springs", stiffness="50", damping="1")
        #fatNode.createObject('HexahedronFEMForceField', template="Vec3d", name="FEM", poissonRatio="0.45", youngModulus="2000")
        fatNode.createObject('TetrahedronFEMForceField', template='Vec3d', name='FEM_tissue',
        method='polar', poissonRatio='0.495', youngModulus='200', computeVonMisesStress='1', showVonMisesStressPerNode='false', listening='1', updateStiffness='true', showStressColorMap='blue 1 0 0 1  1 0.5 0.5 1', isToPrint='1')
        fatNode.createObject('UncoupledConstraintCorrection')
        
        # Visual node
        fatVisNode = fatNode.createChild('Visual')
        fatVisNode.createObject('MeshSTLLoader', name='meshLoader', filename = meshFile)
        fatVisNode.createObject('OglModel', src='@meshLoader', color="#ecc854", name="Fat" )
        fatVisNode.createObject('BarycentricMapping', name="Visual Mapping", output="@Fat")

        # Collision node
        fatColNode = fatNode.createChild('Collision')
        fatColNode.createObject('MeshSTLLoader', name="meshLoader", filename = meshFile)
        #fatColNode.createObject('OglModel',name='Visual', src='@meshLoader', color="#ecc854")
        fatColNode.createObject('MeshTopology', src="@meshLoader")
        fatColNode.createObject('MechanicalObject', src="@meshLoader", name="CollisionObject", template="Vec3d", scale="1.0")
        fatColNode.createObject('TriangleCollisionModel', template="Vec3d")
        fatColNode.createObject('LineCollisionModel')
        fatColNode.createObject('PointCollisionModel')
        fatColNode.createObject('BarycentricMapping', name="Mechanical Mapping")

        # ==========================
        # NEEDLE NODE
        # ==========================

        meshFile = "data/syrette2.obj"
        scale=10
        
        needleNode = self.rootNode.createChild("Needle")

        needleNode.createObject('EulerImplicitSolver', name="ODE solver", rayleighStiffness="0.01", rayleighMass="1.0")
        needleNode.createObject('CGLinearSolver', name="linear solver", iterations="25", tolerance="1e-10", threshold="10e-10")
        needleNode.createObject('MechanicalObject', name="mechObject", template="Rigid3d",  dx = 0, dy = 150, dz = -50, rx = 0, ry = 0, rz = -90.0,  scale3d=[scale, scale, scale])
        needleNode.createObject('UniformMass', name="mass", totalMass="5")
        needleNode.createObject('UncoupledConstraintCorrection')

        # Visual node
        needleVisNode = needleNode.createChild("VisualModel")
        needleVisNode.createObject('MeshObjLoader', name='instrumentMeshLoader', filename=meshFile)
        needleVisNode.createObject('OglModel', name="InstrumentVisualModel", src='@instrumentMeshLoader', dy = -2*scale,  scale3d=[scale, scale, scale])
 
        needleVisNode.createObject('RigidMapping', name="MM-VM mapping", input="@../mechObject", output="@InstrumentVisualModel")

        # Collision node
        needleColNode = needleNode.createChild("CollisionModel")
        needleColNode.createObject('MeshObjLoader', filename=meshFile, name="loader")
        needleColNode.createObject('MeshTopology', src="@loader", name="InstrumentCollisionModel")
        needleColNode.createObject('MechanicalObject', src="@InstrumentCollisionModel", name="instrumentCollisionState", dy = -2*scale, scale3d=[scale, scale, scale])
        needleColNode.createObject('TriangleCollisionModel', name="instrumentTrinagle", contactStiffness="500", contactFriction="0.01")
        needleColNode.createObject('LineCollisionModel', name="instrumentLine", contactStiffness="500", contactFriction="0.01")
        needleColNode.createObject('PointCollisionModel', name="instrumentPoint",
        contactStiffness="500", contactFriction="0.01")
        needleColNode.createObject('RigidMapping', name="MM-CM mapping", input="@../mechObject", output="@instrumentCollisionState")
        
        needleNode.createObject('LinearMovementConstraint', template="Rigid3d",
                             indices = 0,
                             keyTimes=[0, 0.8, 1.6, 1.7],
                             movements= [[0, 0, 0, 0, 0, 0],
                                         [0, -80, 0, 0, 0, 0],
                                         [0, 0, 0, 0, 0, 0],
                                         [0, 0, 0, 0, 0, 0]])
                                         
        # Instantiate floor
        floor = Floor(self.rootNode, name="Floor", translation=[0.0, -180.0, -50.0], uniformScale=6, isAStaticObject=True, color=[0.3, 0.3, 0.4])
        
        return 0

def createScene(rootNode):
    obj = MainScene(rootNode)
    obj.createGraph(rootNode)

As said before, what I’m trying to achieve is to have an initial deformation of the tissue followed by a consequential rupture when the needle force is greater than an elastic response force of the tissue.
Sadly, I only achieved to move the needle inside the tissue without any deformation at all, as far.

Any help would be much much appreciated, thanks!

Federico

P.S. Here the objs files if needed.

[Cai]

Hi, I am new to sofa and Softrobot plugin. I have installed both successfully. They are located in /home/cmx/sofa/src/applications/plugins. I added the following command in the /home/cmx/sofa/src/applications/plugins/CMakeLists.txt to complie the plugins and succeeded.

set(SOFA_BUILD_METIS ON CACHE BOOL "Compile SOFA with Metis support by default as per dependency for the SoftRobots plugin." FORCE)
set(PLUGIN_SOFAPYTHON ON CACHE BOOL "Compile the SofaPython plugin by default as per dependency for the SoftRobots plugin." FORCE)
sofa_add_plugin(STLIB STLIB ON)
sofa_add_plugin(SoftRobots SoftRobots ON)

and I want to run :

mx@cmx-Lenovo-IdeaPad-Y510P:~/sofa/build/bin$ ./runSofa ../../src/applications/plugins/SoftRobots/docs/tutorials/CableGripper/details$ python grippercontroller.py
Traceback (most recent call last):
  File "grippercontroller.py", line 2, in <module>
    import Sofa
ImportError: No module named Sofa

I have some doubts:
1. how to run .py document with Sofa modules, it always can not find Sofa’s path?
2. How to import the correct path according to my file directory
Thanks

[fedev]

Dear Cai,

You should open a new post in the right forum section, “Getting started”, to get help on your issue.

[fedev]

Dear all,

I solved my issue mainly by just setting proper parameters to the contact manager, namely alarmDistance=”3″, contactDistance=”2.3″. Therefore I’m marking this thread as solved and leave the code for the community :).

Federico

[fabianv]

hi many thanks for your input, i’m new on sofa, how can run you scene?, i must be create another file?

thanks

[Author]

Posts