Unable to Visualize Heterogeneous material.

[HaydenGowing]

Hello all,

I am working on a project to simulate the behavior of a firm tumor within a breast model that will be deformed by an ultrasound probe hundreds of times. The position data of all the nodes of the breast/tumor will be recorded and used to train a machine learning model which will then be used to visualize real tumor movement on a HoloLens for surgeons.

I am using the sphereROI to define my size and shape of tumor, I collected the tetrahedron indices from this ROI in order to set the proper FEM of the breast model. I added drawHetergeneousTetra = 1, unfortunately there is no difference in the visualization. In the image attached you can see the blue sphere representing my tumour from the sphereROI but that obviously does not move/deform when use the ultrasound probe. Additionally I am aware that the tetrahedrons contained within my sphereROI do not encapsulate a spherical shaped object and just the arbitrary tetrahedrons inside of it, would it be better to use triangles or points within the sphereROI instead?

If anyone has any suggestions to help my visualize and/or any improvements to how I should go about this project in general I would greatly appreciate it!

(sphereROI visualized and tetrahedrons within sphereROI visualized)

class ControlFrame(Sofa.Core.Controller):

    def __init__(self, node):
        Sofa.Core.Controller.__init__(self)
        self.root = self.CreateGraph(node)

    def CreateGraph(self, root):
        root.addObject('AnimationLoopParallelScheduler', name='mainLoop', threadNumber='0')
        root.addObject('DefaultVisualManagerLoop')
        root.addObject('DefaultAnimationLoop')
        root.addObject('VisualStyle', displayFlags=" showVisualModels")
        root.addObject('RequiredPlugin',
                       pluginName="MultiThreading Sofa.Component.Collision.Detection.Algorithm CImgPlugin SofaMiscCollision SofaGeneralSimpleFem SofaOpenglVisual SofaBaseMechanics SofaBaseTopology SofaSparseSolver SofaImplicitOdeSolver SofaTopologyMapping SofaBoundaryCondition SofaEngine SofaGeneralDeformable SofaMiscFem SofaMeshCollision SofaGeneralLoader")

        root.gravity = [0, 0, 0]
        root.dt = 0.1
        root.addObject("EulerImplicitSolver", name="cg_odesolver", printLog=False, rayleighStiffness=0.1,
                       rayleighMass=1)
        root.addObject("CGLinearSolver", iterations=25, name="linear solver", tolerance=1.0e-9, threshold=1.0e-9)
        root.addObject("DefaultPipeline")
        root.addObject("BruteForceBroadPhase")
        root.addObject("BVHNarrowPhase")
        # root.addObject("GenericConstraintSolver")
        root.addObject("DefaultContactManager", name="Response", response="PenalityContactForceField") #FrictionContactConstraint used initially PenalityContactForceField
        root.addObject("MinProximityIntersection", name="Proximity", alarmDistance=0.1, contactDistance=0.05)

        breast = root.addChild('Breast')
        breast.addObject('VisualStyle', displayFlags=" showWireFrame")
        breast.addObject("EulerImplicitSolver", name="cg_odesolver", printLog=False, rayleighStiffness=0.1, rayleighMass=1)
        breast.addObject("CGLinearSolver", iterations=25, name="linear solver", tolerance=1.0e-9, threshold=1.0e-9)
        # breast.addObject("SparseLDLSolver", name="linear solver")
        breast.addObject("MeshGmshLoader", name="loader", filename="models/BreastModelFull.msh", createSubelements='true')
        breast.addObject('MechanicalObject', name="mo", src="@loader", scale3d=[0.1, 0.1, 0.1])
        breast.addObject('TetrahedronSetTopologyContainer', name='topology', src="@loader")
        breast.addObject('TetrahedronSetGeometryAlgorithms', name='TetraAlg', template='Vec3d')
        breast.addObject('DiagonalMass', totalMass="5")


        breast.addObject('SphereROI', name='tumor', centers='2 -0.5 -1', radii='0.6', drawSphere='1', drawTetrahedra='1')
        breast.addObject('SphereROI', name='breastmodel', centers='2 -1.5 -1', radii='5', drawSphere='0') #, drawTetrahedra='1')


        breast.addObject('IndexValueMapper', name='breastindex', indices='@breastmodel.tetrahedronIndices', value = '2000')#, input='@tumorindex.outputValues')
        breast.addObject('IndexValueMapper', name='tumorindex', indices='@tumor.tetrahedronIndices', value='10000', inputValues='@breastindex.outputValues')

        tumor_indices = [1669, 1723, 2128, 2693, 2958, 3971, 3979, 5851, 6093, 6266, 6357, 6867, 7331, 8541]
        breast_FEM = create_list(15921, tumor_indices)
        

        breast.addObject('TetrahedronFEMForceField', name="FEM", youngModulus=breast_FEM, poissonRatio="0.4", drawHeterogeneousTetra='1') #@tumorindex.outputValues

        breast_indices = list(range(1,647)) #All points along bottom of breast model
        breast.addObject('FixedConstraint', name="FixedConstraint", indices=breast_indices)

        self.breastPos = root.Breast.getObject('mo')
        self.breast = root.Breast

        breast_collision = breast.addChild('collision')
        breast_collision.addObject("MeshOBJLoader", name='stl_load', filename='models/BreastModelTest.obj')
        breast_collision.addObject("MeshTopology", src="@../loader")
        breast_collision.addObject('MechanicalObject', name='Collision', scale3d=[0.1, 0.1, 0.1])
        breast_collision.addObject("TriangleCollisionModel", selfCollision='0',simulated='1', moving='1')
        breast_collision.addObject("LineCollisionModel", selfCollision='0',simulated='1', moving='1')
        breast_collision.addObject("PointCollisionModel", selfCollision='0',simulated='1', moving='1')
        breast_collision.addObject('BarycentricMapping', name='CollisionMapping', input='@../mo', output='@Collision')
        
        breast_visu = breast_collision.addChild('visu')
        # breast_visu.addObject("MeshSTLLoader", name='visu_load', filename='models/BreastModelNew.stl')
        breast_visu.addObject("OglModel", name="VisualModel", src="@../../loader", scale3d=[0.1, 0.1, 0.1])
        breast_visu.addObject("BarycentricMapping", name="VisualMapping", input="@../Collision", output="@VisualModel")


        probe = root.addChild('Probe')
        probe.addObject("MechanicalObject", name='USprobe', scale3d=[0.1, 0.1, 0.1], template="Rigid3d")
        probe.addObject("UniformMass", totalMass=3)
        probe.addObject("EulerImplicitSolver", name="cg_odesolver", printLog=False, rayleighStiffness=0.1, rayleighMass=1)
        probe.addObject("CGLinearSolver", iterations=25, name="linear solver", tolerance=1.0e-9, threshold=1.0e-9)
        probe.addObject('LinearMovementConstraint', keyTimes='0 7 8 12', movements=[[0, 0, 0], [0, 0, 0],
        [2, -0.4, -0.75], [0, 0, 0], 
        [-2, 0.4, 0.75],  [0, 0, 0], 
        [0, 0, 0], [0, 0, 0]])
        # probe.addObject("SparseLDLSolver", name="linear solver")
        # probe.addObject('UncoupledConstraintCorrection')

        probe_collision = probe.addChild('ProbeCollision')
        probe_collision.addObject("MeshOBJLoader", name="loader", filename="models/ProbeCollide.obj")
        probe_collision.addObject("MeshTopology", src="@loader")
        probe_collision.addObject("MechanicalObject", scale3d=[0.1, 0.1, 0.1], template="Vec3", name='Probemodel')
        probe_collision.addObject("RigidMapping", name='collisionmap', input='@../USprobe', output='@Probemodel')
        probe_collision.addObject("TriangleCollisionModel", name="TriangleModel",simulated='1', moving='1')
        probe_collision.addObject("LineCollisionModel", name="LineModel",simulated='1', moving='1')
        probe_collision.addObject("PointCollisionModel", name="PointModel",simulated='1', moving='1')

        probe_visual = probe.addChild('ProbeVisual')
        probe_visual.addObject('MeshOBJLoader', name='visloader', filename='models/VisualProbe.obj')
        probe_visual.addObject('OglModel', name='probevis', src='@visloader', scale3d=[0.1, 0.1, 0.1])
        probe_visual.addObject('RigidMapping', name='ProbeMapping', input='@../USprobe', output='@probevis')

        self.UsProbe = root.Probe.getObject('USprobe')


        return root

[hugtalbot]

Hi @HaydenGowing

Interesting work. To define an heterogeneous elasticity of your breast, you should define a vector with a size equal to the number of tetrahedron in your whole mesh. This vector would be a vector corresponding the Young's modulus per tetrahedron. Then modify the vector entries corresponding to the tetrahedra within the sphere ROI by setting a different Young's modulus (I guess a higher stiffness).

Then, for the visualization, you should indeed to activate the booleans : drawHeterogeneousTetra, computeVonMisesStress and vonMisesPerElement

Hope this helps.
An example XML scene is given here : https://github.com/sofa-framework/sofa/blob/master/examples/Component/SolidMechanics/FEM/Heterogeneous-TetrahedronFEMForceField.scn