Method to generate rotation of ultrasound probe for various collision points

[HaydenGowing]

Hello,

My project consists of deforming a breast containing a tumor inside of it, with an ultrasound probe as seen in the image below.

The issue we have been running into is that we need to deform this breast from a couple thousand different points on the breast mesh. I have all of the points (x,y,z) stored in a file which I know I can cycle through using the LinearMovementConstraint class. The problem is how to rotate/transform the probe to be perpendicular to the desired point on the breast mesh. Is there a class or method that I can use in SOFA in order to generate the desired rotation of my probe to align with its target position for a large amount of positions?

Thanks,
Hayden

import Sofa
import numpy as np
import os,time 
import pywavefront



class controlFrame(Sofa.Core.Controller):

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

        root.addObject('AnimationLoopParallelScheduler', name='mainLoop', threadNumber='0')
        root.addObject('DefaultVisualManagerLoop')
        root.addObject('DefaultAnimationLoop')
        root.addObject('VisualStyle', displayFlags=" showVisualModels showBehaviorModels")
        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=100, 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/BreastTum.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.3615 -1.1784 -0.9338', radii='0.75', drawSphere='1', drawTetrahedra='1')
        breast.addObject('SphereROI', name='breastmodel', centers='2.8 -1.5 -2.3', radii='5', drawSphere='1', drawTetrahedra='1')

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

        tumor_FEM_indices = np.load('tumorFEMindices.npy')
        heterogeneous_FEM = create_list(73258, tumor_FEM_indices)
        breast.addObject('TetrahedronFEMForceField', name="FEM", youngModulus=heterogeneous_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/BreastTum.obj')
        breast_collision.addObject("MeshTopology", src="@../loader")
        breast_collision.addObject('MechanicalObject', name='Collision', scale3d=[0.1, 0.1, 0.1])
        breast_collision.addObject("TriangleCollisionModel",simulated='1', moving='1')
        breast_collision.addObject("LineCollisionModel",simulated='1', moving='1')
        breast_collision.addObject("PointCollisionModel",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="@../stl_load", scale3d=[0.1, 0.1, 0.1])
        # breast_visu.addObject("MeshExporter", filename='DeformedBreast', position='@VisualModel.position', edges='@VisualModel.edges', triangles='@VisualModel.triangles', exportEveryNumberOfSteps='1', format='obj' )
        breast_visu.addObject("OBJExporter", filename='DeformedBreast', exportEveryNumberOfSteps='1')#, format='obj')
        breast_visu.addObject("BarycentricMapping", name="VisualMapping", input="@../Collision", output="@VisualModel")

        pos = probePosition()
        probe = root.addChild('Probe')
        probe.addObject('TransformEngine', name='T', template='Rigid3d', rotation=pos)
        probe.addObject("MechanicalObject", name='USprobe', template="Rigid3d", scale3d=[0.1, 0.1, 0.1])
        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)


        target_position = [2, -0.4, -.75]
        probe.addObject('LinearMovementConstraint', keyTimes='0 7 8 12', movements=[[0, 0, 0], [0, 0, 0], 
        target_position, [0, 0, 0], 
        [-2, 0.4, 0.75],  [0, 0, 0], 
        [0, 0, 0], [0, 0, 0]])

        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", template="Vec3", name='Probemodel', scale3d=[0.1, 0.1, 0.1])
        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')
        
    
    def onAnimateBeginEvent(self, __):

        pass

    def onAnimateEndEvent(self, __):
        pass

[bakpaul]

Hello Hayden.

I could think of several tricks that would do what you are looking for using weak constraints and the solvers to automatically converge to the right orientation, but those would be cumbersome. To answer you question precisely, no component doing exactly that currently exists in the api in my knowledge. But you can compute the normals of the surface mesh using the component NormalsFromPoints, then, by finding the id of closest triangle to the contact point, you can define the orientation of the probe to be along its normal in the opposite direction. This is easily doable in python. The only remaining unknown is the rotation of the probe along this direction, which I suppose is a parameter of your method.

[hugtalbot]

Hi @HaydenGowing

I guess you are aiming at training a NN from your simulation
Why not applying various solicitations to the breast tissue by using a "force" approach using the SurfacePressureForceField ?
You could define iteratively various regions over which to apply the force on the surface and deform the breast.
Would this suit your need?

Cheers

[zhanglepy]

Hello,
I am very interested in the research you have done, and I would like to learn more about the details of your research so that I can better cite your work in the future. So, should we establish a connection? My email is [email protected] .

[hugtalbot]

Hi @zhanglepy
Thanks a lot for your interest.
How to cite us is described here: https://www.sofa-framework.org/applications/publications/#CiteSOFA
How to get involved is described here: https://www.sofa-framework.org/community/get-involved/

For more, I would be glad to establish contact per email as you suggested
Best wishes,

[github-actions[bot]]

⚠️ ⚠️ ⚠️
@HaydenGowing
Feedback has been given to you by the project reviewers, however we have not received a response from you. Without further news in the coming weeks, this discussion will be automatically closed in order to keep this forum clean and fresh 🌱 Thank you for your understanding

[zhanglepy]

@HaydenGowing Hello, have you solved the rotation/transformation problem of the probe? Or real-time control, looking forward to your reply! If you have any new research findings, could you share them so that I can easily cite your latest paper.