Modélisation d'un fémur pour déterminer la probabilité de fracture

Par Jeanne GUILBAUD et Julie BOUSSARD, le 8 décembre 2022.

Introduction

Les fractures ostéoporotiques de la hanche provoquent une perte de la qualité de vie et peuvent même être mortelles chez les personnes âgées. Il est donc nécessaire de pouvoir analyser le fémur afin d'en déterminer le risque de fracture.

Pour atteindre cet objectif nous avons dans un premier temps segmenter un fémur, avant de le mailler avec un maillage d'abord surfacique puis volumique. En superposant le fémur segmenté avec des images issues de CT-scan nous avons pu déterminer le module d'Young de l'os dans les différents points du maillage. Par méthode des éléments finis nous avons pu déterminer la déformation élastique au sein du fémur. A partir de la valeur de cette déformation élastique, nous pouvons déterminer la probabilité de fracture du fémur étudié.

Import des bibliothèques nécessaires au projet

import pyvista as pv import math import sys import gmsh import pydicom import numpy as np import matplotlib.pyplot as plt import glob from fenics import * from ufl import nabla_div

Séparation surfacique du fémur en trois parties: head, body et base

Dans cette partie, nous souhaitons séparer la surface du fémur segmenté en 3 parties limites afin d'appliquer sur chacune d'entre elles les meilleures conditions limites possibles. On considèrera que la tête subira l'entièreté du poids de l'individu tandis que la base sera fixe.

filename

On lit le fichier .stl et on affiche les coordonnées limites pour pouvoir ensuite séparer la surface en trois.

femur = pv.read(filename) pv.start_xvfb() bounds = femur.bounds print(bounds)

On sépare ensuite dans un premier temps la tête du reste. Puis, on coupe le reste en deux : le corps et la base.

Il y a deux manière de faire le premier découpage en fonction de latéralité du fémur. En effet si il s'agit d'un fémur droit ou d'un fémur gauche, le plan de séparation ne sera pas le même (il fera soit un angle de 20° soit de -20° avec la verticale). Pour pouvoir choisir quelle méthode utiliser, on ajoute un argument d'entrée: latéralité.

lateralite

Calcul des différentes surfaces

if lateralite == 'gauche': head,rest = femur.clip(normal =(-math.sin(math.pi/9),0.0,-math.cos(math.pi/9)),origin=(bounds[1]-0.15*(bounds[1]-bounds[0]),0.5*(bounds[3]-bounds[2]),bounds[5]-0.1*(bounds[5]-bounds[4])), return_clipped=True) else: head,rest = femur.clip(normal =(math.sin(math.pi/9),0.0,-math.cos(math.pi/9)),origin=(bounds[1]-0.6*(bounds[1]-bounds[0]),0.2*(bounds[3]-bounds[2]),bounds[5]-0.1*(bounds[5]-bounds[4])), return_clipped=True) body,base = rest.clip(origin = (0.0,0.0,(bounds[5])-0.90*(bounds[5]-bounds[4])),normal ="-z",return_clipped=True)

Affichage des trois parties

p = pv.Plotter() p.add_mesh(head, color='blue', label='Femur Head') p.add_mesh(base, color='red', label='Femur Base') p.add_mesh(body, label='Body') p.add_legend() p.add_title('accidentally french femur') p.show()

Enregistrement des trois surfaces dans trois fichiers stl différents

head.save('./Output/head.stl') body.save('./Output/body.stl') base.save('./Output/base.stl')

Génération d'un maillage 3D

geofile

gmsh.initialize(sys.argv) #change the verbosity level to display more or less output : # (0: silent except for fatal errors, 1: +errors, 2: +warnings, 3: +direct, 4: +information, 5: +status, 99: +debug) gmsh.option.setNumber("General.Verbosity", 5) gmsh.option.setNumber("Mesh.LabelType",2) gmsh.option.setNumber("Mesh.SaveAll",1) gmsh.option.setNumber("Mesh.SaveElementTagType",2) gmsh.open(geofile) gmsh.option.setNumber("Mesh.Algorithm", 6) gmsh.model.occ.synchronize() gmsh.model.mesh.generate(3) gmsh.write("femur.mesh") # MEDIT gmsh.write("femur.inp") # ABAQUS gmsh.write("femur.vtk") # VTK gmsh.finalize()

# read the vtu mesh with pyvista pv.start_xvfb() mesh = pv.read("/work/Data/femur.stl") cpos = mesh.plot()

mesh2 = pv.read("femur.vtk") cpos = mesh2.plot()

Cartographie des propriétés matériaux

Nous avons construit un model surfacique 3D du fémur, il nous faut désormais modéliser les propriétés matériaux de l'os, notamment son module d'Young aux différents points du maillage. Pour déterminer le module d'Young, nous nous appuyons sur des images médicales réalisées par CT-scan que nous superposons au fémur. Nous relevons les valeurs de l'Hounsfield Unit pour chaque pixel de l'image. Grâce à une relation déjà existant, nous sommes alors capables de relier l'Hounsfield Unit avec la densité minérale osseuse. En utilisant une relation élasticité/densité nous pouvons dès lors remonter à la valeur du module d'Young.

Affichage de l'image médical et du fémur segmenté

reader = pv.DICOMReader("/work/CT_dataset") mesh = pv.read("/work/Data/femur.stl") dataset = reader.read() pv.start_xvfb() plotter = pv.Plotter(off_screen=True) plotter.add_mesh(mesh) plotter.add_mesh(dataset) plotter.show(screenshot="myscreenshot.png") print(dataset)

Nous remarquons que le fémur et les images ne sont pas superposés. Ceci s'explique par le fait que les deux objets ne sont pas représentés dans le même système de coordonnées. Les images médicales sont dans un système de coordonnées (i,j,k) où i et j vont de 1 au nombre de pixel de l'image et k va de 1 au nombre d'image, tandis que le fémur segmenté utilise les coordonnées physiques. Pour permettre la superposition et des deux objets et donc le calcul du module d'Young, nous allons utiliser le format DICOM.

Le Format DICOM

Le format DICOM (Digital Imaging and Communications in Medicine) est le format principalement utilisé pour traiter des images médicales et leurs données. Pour superposer les images médicales avec notre fémur, nous allons calculer les coordonnées de chaque voxel des fichiers dicom grâce à la formule ci-dessous.

Où:

Pxyz sont les coordonnées en mm du voxel (i,j) dans le plan image de la trame, ce sont nos inconnues.

Sxyz position en mm à partir de l'origine du RCS, les valeurs sont données par le vecteur ImagePositionPatient.

Xxyz: données par les trois premières valeurs ImageOrientationPatient.

Yxyz: données par les trois premières valeurs ImageOrientationPatient.

i est l'indice de colonne de chaque image, il commence à 0.

Δi résolution en mm des pixels de la colonne de Pixel Spacing.

j est l'indice de ligne de chaque image, il commence à 0

Δj résolution en mm des pixels de la ligne de Pixel Spacing.

Toutes ces variables sont contenues dans les informations de l'image, il suffit d'utiliser le bon argument pour obtenir leur valeur.

Calcul des coordonnées des voxels de chaque image

Nous allons utiliser pydicom pour manipuler la série d'images enregistrées au format dicom.

Nous allons désormais appliquer la formule décrite ci-dessus afin de calculer les coordonnées des voxels.

def readDicomWithCoords(pathname): #THIS FUNCTION RETURN THE HOUNSFIELD UNIT IN THE 3D IMAGE files = [] if pathname[-1] == "/": newpath = pathname + "*.dcm" else: newpath = pathname + "/*.dcm" for fname in glob.glob(newpath, recursive=False): print("loading: {}".format(fname)) files.append(pydicom.dcmread(fname)) print("file count: {}".format(len(files))) # skip files with no SliceLocation (eg scout views) slices = [] skipcount = 0 for f in files: if hasattr(f, 'SliceLocation'): slices.append(f) else: skipcount = skipcount + 1 print("skipped, no SliceLocation: {}".format(skipcount)) # ensure they are in the correct order slices = sorted(slices, key=lambda s: s.SliceLocation) # pixel aspects, assuming all slices are the same ps = slices[0].PixelSpacing ss = slices[0].SliceThickness print(ps) print(ss) # create 3D array img_shape = list(slices[0].pixel_array.shape) img_shape.append(len(slices)) img3d = np.zeros(img_shape) NumberOfVoxels = img3d.size # fill 3D array with the images from the files for i, s in enumerate(slices): img2d = s.pixel_array img3d[:, :, i] = img2d.T x = np.zeros(img_shape) y = np.zeros(img_shape) z = np.zeros(img_shape) for k, s in enumerate(slices): # remplissage de la matrice M décrite ci-dessus delta_i = s.PixelSpacing[0] delta_j = s.PixelSpacing[1] M14 = s.ImagePositionPatient[0] M24 = s.ImagePositionPatient[1] M34 = s.ImagePositionPatient[2]-s.SliceThickness/2 M11 = s.ImageOrientationPatient[0]*delta_i M21 = s.ImageOrientationPatient[1]*delta_i M31 = s.ImageOrientationPatient[2]*delta_i M12 = s.ImageOrientationPatient[3]*delta_j M22 = s.ImageOrientationPatient[4]*delta_j M32 = s.ImageOrientationPatient[5]*delta_j for i in range (img_shape[0]): for j in range (img_shape[1]): # calcul des cordonnées x[i,j,k] = M11*i+M12*j+M14 y[i,j,k] = M21*i+M22*j+M24 z[i,j,k] = M31*i+M32*j+M34 return img3d,x,y,z

Maintenant que nous disposons des coordonnées des voxels nous pouvons interpoler l'intensité des images médicales directement sur notre maillage en utilisant des fonctions d'interpolation de Scipy.

img, x,y,z = readDicomWithCoords("/work/CT_dataset") #create a new structured grid with the coordinates of the image data_to_probe = pv.StructuredGrid(x, y, z) #add the image intensity as a point data field named HU data_to_probe.point_data['HU'] = img.ravel(order='F') # Interpolate the image data onto the stl mesh interpHU = mesh.sample(data_to_probe) #visualize the scene pv.start_xvfb() plotter = pv.Plotter(off_screen=True) plotter.add_mesh(data_to_probe,opacity=0.5,color=True,scalars='HU') plotter.add_mesh(interpHU,color=True,scalars='HU') plotter.show(screenshot="myscreenshot.png")

On interpole ensuite l'Hounsfield Unit sur le fémur uniquement pour ne garder que la partie qui nous intéresse.

interpHU.plot(scalars="HU")

Nous devons désormais échantillonner l'Hounsfield Unit sur un nuage de points.

# extract the points from mesh pts = mesh.points print(pts.shape) # this gets a list of points (x,y,z) as a pyvista ndarray which is equivalent to numpy ndarray point_cloud = pv.PolyData(pts) # Interpolate the image data onto the point cloud results = point_cloud.sample(data_to_probe) results.plot(scalars="HU") extract_val=results.point_data['HU'] print(extract_val)

Conversion de l'Hounsfield Unit en module d'Young

Grâce aux deux articles suivants, nous avons pu identifier une relation entre HU et E. Il s'agit d'une relation linéaire élasticité/densité que l'on corrige à cause des incertitudes de mesure.

#Le module d'Young calculé ici est un tableau de toutes les valeurs de E en chaque voxel. rhoQCT=-0.0039357+0.000791701*extract_val rhoAsh = 0.079+ 0.877*rhoQCT E=(14664*rhoAsh**1.49)*1000000 #On met le module d'Young en Pascal print('E=',E.max(),'Pa')

Pour vérifier la cohérence de nos résultats on affiche le maximum des valeurs de E on obtient un maximum de 18408 MPa, ce qui correspond bien à l'ordre de grandeur du module d'Young de l'os.

Nous avons donc jusqu'ici réussi à segmenter, délimiter les surfaces puis mailler le fémur que nous avons superposer avec les images médicales. En interpolant la valeur de HU sur le fémur segmenté nous avons alors pu en déterminer le module d'Young. Toutes ces étapes sont pour le moment uniquement appliquées à la surface du fémur. Dans la suite du projet, nous allons considérer le fémur dans son volume entier affin de pouvoir en déterminer la probabilité de fracture.

Obtention de la probabilité de fracture

geofile

Dans cette dernière partie, nous cherchons à obtenir la probabilité de fracture. Pour cela, nous allons générer un nouveau maillage en volume sur lequel nous appliquerons les différentes contraintes auxquelles l'os est soumis. Une fois ces contraintes définies, nous utiliserons la contrainte maximale pour en déduire la probabilité de fracture maximale.

Initialisation et mise en place des éléments nécessaires à la création du maillage .msh

gmsh.initialize(sys.argv) #change the verbosity level to display more or less output : # (0: silent except for fatal errors, 1: +errors, 2: +warnings, 3: +direct, 4: +information, 5: +status, 99: +debug) gmsh.option.setNumber("General.Verbosity", 5) gmsh.option.setNumber("Mesh.LabelType",2) #gmsh.option.setNumber("Mesh.SaveAll",1) gmsh.option.setNumber("Mesh.SaveElementTagType",2) gmsh.open(geofile) gmsh.option.setNumber("Mesh.Algorithm", 6) gmsh.model.occ.synchronize() gmsh.model.mesh.generate(3) gmsh.write("Femur.msh") # GMSH MESH FILE FORMAT gmsh.finalize()

Création du maillage

import meshio msh = meshio.read("Femur.msh") print(msh) def create_mesh(mesh, cell_type, prune_z=False): cells = mesh.get_cells_type(cell_type) cell_data = mesh.get_cell_data("gmsh:physical", cell_type) field_data = mesh.field_data out_mesh = meshio.Mesh(points=mesh.points, cells={ cell_type: cells}, cell_data={"physical_tags": [cell_data]}, field_data=field_data) if prune_z: out_mesh.prune_z_0() return out_mesh tetra_mesh = create_mesh(msh, "tetra") triangle_mesh = create_mesh(msh, "triangle") meshio.write("femur.xdmf", tetra_mesh) meshio.write("femur_bc.xdmf", triangle_mesh)

Calcul du module d'Young E et des différentes variables

Dans cette partie, il a été vérifié que tous les fichiers .stl et toutes les distances considérées sont en mm.

La force qui s'exerce sur le fémur est le poids du corps, on suppose qu'elle s'applique uniquement sur la tête du fémur, sur sa surface. On considère une personne appliquant 55kg sur son fémur.

Poids = 9.81*55/1000 ## poids en kg.mm.s^-2

On redéfinit maintenant différentes variables variables présentes dans notre maillage comme étant les parties ou caractéristique du fémur.

mesh = Mesh() with XDMFFile("femur.xdmf") as infile: infile.read(mesh) mvc = MeshValueCollection("size_t", mesh, 2) with XDMFFile("femur_bc.xdmf") as infile: infile.read(mvc, "physical_tags") mf = cpp.mesh.MeshFunctionSizet(mesh, mvc) physical_names = {"base": 9 , "body": 8 , "head" : 7}

Pour calculer la force surfacique, il nous faut calculer la surface : celle de la tête du fémur. L'aire est en mm*mm. Enfin, on peut calculer la force surfacique appliquée.

ds_custom = Measure("ds", domain=mesh, subdomain_data=mf, subdomain_id=physical_names["head"]) area = assemble(1*ds_custom) ##calcul de l'aire de head area = area ## aire en mm^2 print("L'aire est de ", area,"mm^2") Load = Poids/area ## calcul de Fs V = VectorFunctionSpace(mesh, 'P',1)

Ainsi, la tête de notre fémur est soumise au poids du corps.

Reste maintenant à exploiter notre deuxième condition aux limites, la base fixée :

bc = DirichletBC(V, Constant((0, 0, 0)),mf,physical_names["base"])

Nous pouvons maintenant définir les contraintes et tensions qui s'appliquent sur notre fémur. Puis calculer le module d'Young, comme nous l'avons fait précédemment.

U = FunctionSpace(mesh,'P',1) n = U.dim() di = mesh.geometry().dim() dof_coordinates = U.tabulate_dof_coordinates() z = dof_coordinates[:, 2] nu = 0.3 def calcul_E(HU): rhoQCT=-0.0039357+0.000791701*HU rhoAsh = 0.079+ 0.877*rhoQCT E=14664*rhoAsh**1.49 *1000000 ## calcul du module d'Young en Pascal return E

point_cloud = pv.PolyData(dof_coordinates) # Interpolate the image data onto the point cloud results = point_cloud.sample(data_to_probe) results.plot(scalars="HU") extract_val=results.point_data['HU']

Ci dessus, nous pouvons visualiser HU réparti sur le maillage de notre fémur. Les valeurs de HU nous servent maintenant à en déduire E grâce à la fonction calcul_E . Les autres variables et constantes sont calculées telles que :

E = Function(U,name="Young Modulus") E.vector().set_local(calcul_E(extract_val)) mu = E/(2*(1+nu)) lambda_ = E*nu/((1+nu)*(1-2*nu)) def epsilon(u): return 0.5*(nabla_grad(u) + nabla_grad(u).T) #return sym(nabla_grad(u)) def sigma(u): return lambda_*nabla_div(u)*Identity(d) + 2*mu*epsilon(u) u = TrialFunction(V) d = u.geometric_dimension() # space dimension v = TestFunction(V) f = Constant((0, 0, 0)) T = Constant((-Load*np.sin(np.pi/9), 0, -Load*np.cos(np.pi/9))) a = inner(sigma(u), epsilon(v))*dx L = dot(f, v)*dx + dot(T, v)*ds_custom

On peut maintenant calculer la solution.

# Compute solution u = Function(V,name="Displacement") solve(a == L, u, bc)

Calcul des contraintes et de la déformation

Par soucis de simplification, on considère une contrainte de Von Mises s'appliquant sur notre fémur.

s = sigma(u) - (1./3)*tr(sigma(u))*Identity(d) # deviatoric stress von_Mises = sqrt(3./2*inner(s, s)) TAU = FunctionSpace(mesh, 'P', 1) # same as V von_Mises = project(von_Mises, TAU) eps= sqrt(3./2*inner(s, s)) eps = project(eps, TAU) von_Mises.rename('Mises','Von Mises Stress') with XDMFFile('FemurSolution.xdmf') as file: file.parameters.update( { "functions_share_mesh": True, "rewrite_function_mesh": False }) file.write(u, 0.) file.write(von_Mises, 0.) file.write(E,0.) file.write(eps,0.)

Calcul de la probabilité maximale à partir de la déformation

La probabilité de fracture étant proportionnelle à la déformation subie par l'os, nous n'utilisons que la valeur max de epsilon pour calculer une probabilité maximale de fracture. Nous avons :

def Proba(eps): epsilon0=0.0015 D_cri = 0.38 beta=5 epsilonc=0.0174 gamma=np.log(1/(1-D_cri))/(np.log(epsilon0/epsilonc)) kappa=1/(epsilon0**gamma) if eps<=epsilon0: D=0 elif eps >=epsilonc: D=D_cri else : D=1-kappa*(eps**gamma) P=(D/D_cri)**(beta/2-1) return P print('la probabilité de facture est ',Proba(eps.vector().max()))

Finalement, nous obtenons une probabilité de fracture de 0,484.

Conclusion

Pour conclure ce projet, les différentes étapes nous ont permis d'évaluer les contraintes s'appliquant sur le fémur à partir de sa segmentation. Il s'agit ici d'un exemple d'application de modélisation au service de la santé, la donnée en entrée étant un fémur et un jeu de données d'imagerie médicale en sortie une probabilité de risque de fracture.

Toutefois, la probabilité maximale est d'autant plus utile qu'elle est accompagnée de sa localisation et il serait même encore mieux d'établir une carte des probabilité. Cela permettrait de pouvoir détailler les zones à risque et ainsi d'aider au mieux le patient à prévenir une potentielle fracture.

Critiques

Plus précisément, la probabilité finale trouvée n'est pas en accord avec la probabilité attendue, même si elle reste physiquement possible. En effet, la probabilité de fracture devrait être proche de 0 pour un fémur de cette taille subissant un poids de 55kg. Il doit donc il y avoir une erreur d'ordre de grandeur sur l'une de nos variables : nous n'avons pas réussi à l'identifier.

Bibliographie

Helgason, Benedikt, et al. "Mathematical relationships between bone density and mechanical properties: a literature review." Clinical biomechanics 23.2 (2008): 135-146.

Schileo, Enrico, et al. "Subject-specific finite element models can accurately predict strain levels in long bones." Journal of biomechanics 40.13 (2007): 2982-2989.

Enrique López, "A mechanical model for predicting the probability of osteoporotic hip fractures based in DXA measurements and finite element simulation." BioMedCentral (2012)