import os
import sys
sys.path.insert(0, os.path.abspath('../source/py/myAceplaque'))
sys.path.insert(0, os.path.abspath('../source/py/myModulesV2'))
sys.path.insert(0, os.path.abspath('../source/py/mesh-importer_v2'))
sys.path.insert(0, os.path.abspath('../source/py/myPatientsModulesV2'))
#import myModulesV2.importMeshFromTXT as imesh
#import myModulesV2.postResults as pres
import pandas as pd
#import myPatientsModulesV2.nrrdAux as nr
import collections
from tabulate import tabulate
import pyvista as pv
import glob
import os
import json
import numpy as np
import matplotlib.pyplot as plt
[docs]
class myPatient:
def __init__(self, name, customPath=None):
"""
Represents a patient loading both FEM and clinical data.
Automations once the name is defined:
- loading the mesh
- compute the average radii (internel, external, wall thickness)
- load the materials
- load the JSON with the FEM information and update it
- load `matsSegmented`, a segmeneted version for the material: 0 for `CALC`, 1 for `MATX` and 2 for `LRNC-IPH`
Attributes:
`name` (str): The name of the patient.
`mesh`: The 3D mesh representing the patient's anatomy.
`radii`: The computed radii based on the patient's anatomy.
`materials`: The materials associated with the patient.
`fields` (pd.DataFrame): The fields associated with the patient.
`nrrdHeader` (collections.OrderedDict): The header information for the patient's NRRD (Nearly Raw Raster Data) file.
`info` (collections.OrderedDict): Additional information about the patient.
`FEMinfo`: Information related to the Finite Element Method (FEM) analysis for the patient.
`outputPath` (str): The path where the patient's results are stored.
Methods:
-----
`showMesh()`: Displays the patient's mesh.
`loadClinicalData()`: Loads the clinical data for the patient from NRRD files.
`showClinicalInfo()`: Displays the clinical information associated with the patient NRRD files.
`showInfo()`: Displays the additional information associated with the patient.
`showMaterials()`: Displays the patient's materials calling `getMaterials()`.
`showField()`: Displays the field associated with the patient calling `getField()`.
`showClinicalData()`: Displays the clinical data associated with the patient calling `getClinicalData()`.
`showFEMJSON()`: Displays the FEM information associated with the patient.
`getMaterials()`: Retrieves the materials associated with the patient.
`getField()`: Retrieves the field associated with the patient.
`getFieldValues()`: Retrieves the numerical values associated with the `fieldName` field.
`getClinicalData()`: Retrieves the clinical volumetric data from Pyvista volumes.
`showHist()`: Displays the histogram of the `fieldName` results for the `step` step.
`getHist()`: Retrieves the histogram of the `fieldName` results for the `step` step.
`getCompareHist()`: Retrieves the histogram of the `fieldName` results between the classic and growth results respectively for the `classicStep` and `growthStep` steps.
`compareFields()`: Compares the `fieldName` results between the classic and growth results respectively for the `classicStep` and `growthStep` steps.
`compareHist()`: Compares the histogram of the `fieldName` results between the classic and growth results respectively for the `classicStep` and `growthStep` steps.
"""
pv.global_theme.allow_empty_mesh = True
self.name = name
self.customPath = customPath
self.loadMesh()
self.radii=self.computeRadii()
self.loadMaterials()# it load materials based on c2 coefficient
self.fields = pd.DataFrame()
self.nrrdHeader = collections.OrderedDict()
self.info = collections.OrderedDict()
self.FEMinfo = self.loadFEMJSON()
self.updateFEMJSON()
self.outputPath = pres.pathForResults(self)+self.name
self.updateMaterialsFromNodesInfo()# it load field 'matsSegmented' based on numerical id
[docs]
def loadMesh(self):
"""
`myPatient.loadMesh`
Loads the mesh for the patient based on the name provided.
"""
try:
# Call the method
#self.mesh = imesh.loadMeshFromTXT(self.name, self.customPath)
self.mesh = imesh.loadMeshFromHDF5(self.name, self.customPath)
except ValueError as e:
# If a ValueError is raised, this block will be executed
print(f"An error occurred loading the mesh: {e}")
return
[docs]
def loadMaterials(self):
"""
`myPatient.loadMaterials` use the c2 coefficient from the FEM results to load the materials.
It works using the `loadMaterialFromTXT` function from `importMeshFromTXT.py`.
"""
#self.mesh=imesh.loadMaterialFromTXT(self.mesh,self.name, self.customPath)
self.mesh = imesh.loadMaterialFromHDF5(self.mesh,self.name, self.customPath)
return
[docs]
def showMesh(self, plotter=None):
"""
`myPatient.showMesh` shows the 3D mesh for the patient.
"""
p = pres.get3Dmesh(self.mesh,plotter=plotter)
p.show()
[docs]
def showMaterials(self,slices=False, plotter=None):
"""
Visualizes the materials in the mesh calling the `getMaterials` function.
## Parameters:
`slices` (bool): If True, the materials will be displayed including principal slices. Default is False.
"""
p=self.getMaterials(slices=slices,plotter=plotter)
p.show()
[docs]
def getMaterials(self,slices=False,plotter=None, normal=None, center = None, addBar=False, addMesh= False, addLegend = True):
"""
`myPatient.getMaterials` retrieves the materials associated with the mesh field.
Args:
`slices` (bool): If True, the materials will be displayed in slices. Default is False. If only `slices=True` is specified, the principal slices will be displayed.
`plotter` (pyvista.plotting.plotting.Plotter): The plotter object. Default is None.
`normal` (list nx, nx, nz): Normal vector for slicing. Default is None.
`center` (list cx, cy, cz): Center point for slicing. Default is None which leads to myPatient.PlaqueSliceCenter.
"""
if slices and normal is None and center is None:
center = self.plaqueSliceCenter
p = pres.get3DslicesWithMaterials(self.mesh, plotter=plotter, addBar = True, center = center)
elif slices and normal or center is not None:
if center is None:
center = self.plaqueSliceCenter
if not addMesh:
p = pres.getSlicemeshWithMaterials(self.mesh, normal = normal, plotter=plotter, center=center, addBar=addBar)
else:
p = pres.get3DslicesWithMaterials(self.mesh, plotter=plotter, normal=normal, center=center, addBar=addBar)
else:
p = pres.get3DmeshWithMaterials(self.mesh, plotter=plotter, addLegend=addLegend)
return p
[docs]
def showField(self,fieldName='u',step=1,GRLab='Classic',clim=[-10, 80], scaling=None, plotter=None, slices=False, normal=None, center = None):
"""
`myPatient.showField` shows the field `fieldName` for the patient by visualizing the Pyvista results from the `getField` function.
## Parameters:
`fieldName` (str): Name of the field to show. Default is 'u'.
`step` (int): Load step value for the field. Default is 1.
`GRLab` (str): The label for the field. Default is 'Classic' otherwise you can specify 'G&R'.
`clim` (list): List of two values representing the color limit. Default is [-10, 80].
`scaling` (string): Scaling factor. Default is None.
`plotter` (pyvista.plotting.plotting.Plotter): The plotter object. Default is None.
`slices` (bool): If True, the field will be displayed in slices. Default is False.
`normal` (list nx, nx, nz): Normal vector for slicing. Default is None.
`center` (list cx, cy, cz): Center point for slicing. Default is None which leads to myPatient.PlaqueSliceCenter.
"""
p = self.getField(fieldName=fieldName,step=step,GRLab=GRLab,clim=clim, scaling=scaling, plotter=plotter,slices=slices, normal=normal, center = center)
p.show()
[docs]
def compareFields(self,fieldName='u',classicStep=1,growthStep=1, clim=[-10, 80], scaling=None, normal=None, center = None):
"""
`myPatient.compareFields`
Compares the `fieldName` results between the classic and growth results respectively for the `classicStep` and `growthStep` steps.
You can also specify the color limit `clim`, the scaling factor `scaling`, the normal vector for slicing `normal` and the center point for slicing 'center' for the slice.
If you specify only the normal the my_data.PlaqueSliceCenter will be used as center point.
If you specify only the center point the normal will be set to [0,0,1].
## Parameters:
`fieldName` (str): Name of the field to compare. Default is 'u'.
`classicStep` (int): Step value for classic field. Default is 1.
`growthStep` (int): Step value for growth field. Default is 1.
`clim` (list): List of two values representing the color limit. Default is [-10, 80].
`scaling` (string): Scaling factor. Default is None.
`normal` (list nx, nx, nz): Normal vector for slicing. Default is None.
`center` (list cx, cy, cz): Center point for slicing. Default is None.
Returns:
None
"""
if center is None:
center = self.plaqueSliceCenter
if isinstance(classicStep, (int, float)):
classicStepValue=classicStep
else:
if classicStep == 'MaxLoad':
classicStepValue=self.classicStep['MaxLoad']
elif classicStep == 'MAP':
classicStepValue=self.classicStep['MAP']
if isinstance(growthStep, (int, float)):
growthStepValue=growthStep
else:
if classicStep == 'MaxLoad':
growthStepValue=self.growthStep['MaxLoad']
elif classicStep == 'MAP':
growthStepValue=self.growthStep['MAP']
p= pv.Plotter(shape=(1, 3))
p.subplot(0, 0)
self.getField(fieldName=fieldName,slices=True,step=classicStepValue, GRLab='Classic', clim=clim, scaling=scaling, plotter=p, normal=normal, center = center)
p.subplot(0, 1)
self.getField(fieldName=fieldName,slices=True,step=growthStepValue, GRLab='G&R', clim=clim, scaling=scaling, plotter=p, normal=normal, center = center)
p.subplot(0, 2)
pres.getSlicemeshWithMaterials(self.mesh, normal = normal, plotter=p, center=center, addBar=False)
p.show()
[docs]
def getField(self,fieldName='u', step=1, GRLab='Classic', clim=[-10, 80], scaling=None, plotter=None,slices=False, normal=None, center = None, threeshold=None, textLabel = True, scalarBar = True, cmap = None):
"""
`myPatient.getField` retrieves the field associated with the patient using the Pyvista mesh data.
## Parameters:
`fieldName` (str): Name of the field to retrieve. Default is 'u'.
`step` (int): Load step value for the field. Default is 1.
`GRLab` (str): The label for the field. Default is 'Classic' otherwise you can specify 'G&R'.
`clim` (list): List of two values representing the color limit. Default is [-10, 80].
`scaling` (string): Scaling factor. Default is None.
`plotter` (pyvista.plotting.plotting.Plotter): The plotter object. Default is None.
`slices` (bool): If True, the field will be displayed in slices. Default is False. If you specifiy only `slices=True` the principal slices will be displayed.
`normal` (list nx, nx, nz): Normal vector for slicing. Default is None.
`center` (list cx, cy, cz): Center point for slicing. Default is None which leads to myPatient.PlaqueSliceCenter.
"""
if not (step,fieldName,GRLab) in self.fields.columns:
self.loadFields(fieldsName=[fieldName],step=step,GR=GRLab=='G&R')
mesh=self.mesh
field=str(fieldName)+'_step_'+str(step)+'_'+GRLab
#print(field)
scalingFac = computeScaleFactor(scaling)
mesh[field]=scalingFac*self.fields[(step,fieldName,GRLab)][0:len(mesh.points)]
if threeshold is not None:
if len(threeshold)==1:
mesh=mesh.threshold(value=threeshold, scalars=field)
if len(threeshold)==2:
if threeshold[0]=='>':
mesh=mesh.threshold(value=threeshold[1], scalars=field)
if threeshold[0]=='<':
mesh=mesh.threshold(value=threeshold[1], scalars=field, invert=True)
if not slices:
if hasattr(self,'plaqueSliceCenter') and center is None:
center = self.plaqueSliceCenter
p = pres.get3DmeshWithFields(mesh,field=field,clim=clim,scaling=scaling, plotter=plotter, center = center, cmap=cmap)
else:
if center is None:
center = self.plaqueSliceCenter
p = pres.get3DSlicesWithFields(mesh,field=field,clim=clim,scaling=scaling, plotter=plotter, normal=normal, center=center, textLabel = textLabel, scalarBar = scalarBar, cmap=cmap)
return p
[docs]
def getFieldValues(self,fieldName='u',step=1,GRLab='Classic', scaling = None, limitRegions=None):
"""
`myPatient.getFieldValues` retrieves the array of values associated with the specified `fieldName`.
## Parameters:
`fieldName` (str): Name of the field to retrieve. Default is 'u'.
`step` (int): Load step value for the field. Default is 1.
`GRLab` (str): The label for the field. Default is 'Classic' otherwise you can specify 'G&R'.
`scaling` (string): Scaling factor. Default is None.
`limitRegions` (str): The region to limit the field values. Default is None but you can specify:
- 'VascularTissue' to filter data for all the vascular wall
- 'PlaqueRegion' to filter data for al tissues in the plaque region
- 'PlaqueRegionVascularTissue' to filter data for the vascular wall in plaque region
"""
scalingFac = computeScaleFactor(scaling)
if not (step,fieldName,GRLab) in self.fields.columns:
self.loadFields(fieldsName=[fieldName],step=step,GR=GRLab=='G&R')
#σ1
out = scalingFac*self.fields[(step,fieldName,GRLab)].values[0:len(self.mesh.points)]
if limitRegions is not None:
if limitRegions == 'VascularTissue':
out = out[self.nodesIdx.wall]
elif limitRegions == 'PlaqueRegion':
out = out[self.nodesPlaqueRegionIdx.all]
elif limitRegions == 'PlaqueRegionVascularTissue':
out = out[self.nodesPlaqueRegionIdx.wall]
else:
print("Region not recognized")
return out
[docs]
def loadFields(self,fieldsName=['u','v'],step=1,GR=False):
"""
`myPatient.loadFields` loads the fields associated with the patient's data.
"""
if GR:
data=pres.importGrowthResults(self,self.name,step,customPath=self.customPath)
GRLab='G&R'
else:
data=pres.importClassicResults(self,self.name,step,customPath=self.customPath)
GRLab='Classic'
for field in fieldsName:
self.fields[(step,field,GRLab)]=data[field]
return
[docs]
def loadClinicalData(self,allData=False):
"""
`myPatient.loadClinicalData` loads the clinical data for the patient.
If `allData=True`, it will load all the available fields. Otherwise it only loads the 'lumenSegmentation' and 'wallSegmentation' fields.
It also computes the geometrical features like volume ratios and bounding boxes.
"""
if self.name in ['A_643', 'B_710', 'C_753', 'D_779', 'E_735', 'F_740', 'G_757', 'H_759', 'I_771', 'J_CP20-009', 'K_860', 'L_1081', 'M_1055', 'N_785', 'O_809', 'P_817', 'Q_874', 'R_780']:
if not hasattr(self, 'clinicalData'):
self.clinicalData = nr.setUpNRRDField(self)
if allData:
fieldsToLoad=nr.findAllPtFields(self.name)
else:
fieldsToLoad=['lumenSegmentation','wallSegmentation']
for field in fieldsToLoad:
self.clinicalData = nr.loadNRRD(self,mat=field)
# self.clinicalData = nr.loadNRRD(self,mat='lumenSegmentation')
# self.clinicalData = nr.loadNRRD(self,mat='wallSegmentation')
self.clinicalData = nr.postProcessNRRD(self)#compute the wall as difference between totalWall - lumen
nr.loadNRRDcomposition(self)
# compute the geometrical features
self.computeGeometricalFeatures()
return
[docs]
def getClinicalData(self,mat='general',materials=False,addLume=False, plotter=None):
"""
`myPatient.getClinicalData` retrieves the clinical data associated with the patient from NRRD volumetric data.
## Parameters
`mat` (str): The material to display. Default is 'general'. You can specify:
- 'lumen' to display the lumen.
- 'wall' to display the wall.
- 'LRNC' or 'IPH' to display the Lipid Rich Necrotic Core or Intraplaque Hemorrhage.
- 'CALC' to display the calcification.
`materials` (bool): If True, all materials will be displayed. Default is False.
"""
#if mat=='wall':
if plotter is None:
p = pv.Plotter(window_size=[600, 800])
else:
p = plotter
color='red'
if mat=='lumen':
mats=['lumen']
colors=['red']
opac=[0.1]
elif mat=='wall':
mats=['wall']
colors=['red']
opac=[0.1]
elif (mat=='LRNC'or mat=='IPH'):
mats=[mat,'wall']
colors=['orange','gray']
opac=[0.4,0.1]
elif mat=='CALC':
mats=['CALC','wall']
colors=['navy','gray']
opac=[0.4,0.1]
if materials:
mats=['wall','CALC','LRNC','IPH','MATX']
colors=['gray','navy','orange','orange','red']
opac=[0.1,0.4,0.4,0.4,0.1]
if addLume:
mats.append('lumen')
colors.append('cyan')
opac.append(0.05)
for i in range(len(mats)):
thresholded=nr.computeVolumetricData(self.clinicalData,mat=mats[i])
p.add_mesh(thresholded, color=colors[i], show_scalar_bar=True, opacity=opac[i],label=mats[i])
p.title = 'Volumetric data:'
p.add_legend(loc='lower right',bcolor='white', face='None')
# Display the plot
return p
[docs]
def showClinicalData(self,mat='general',materials=False,addLume=False, plotter=None):
"""
`myPatient.showClinicalData` shows the clinical data for the patient. It calls the `getClinicalData` function.
## Parameters
`mat` (str): The material to display. Default is 'general'. You can specify:
- 'lumen' to display the lumen.
- 'wall' to display the wall.
- 'LRNC' or 'IPH' to display the Lipid Rich Necrotic Core or Intraplaque Hemorrhage.
- 'CALC' to display the calcification.
`materials` (bool): If True, all materials will be displayed. Default is False.
"""
p = self.getClinicalData(self,mat=mat,materials=materials,addLume=addLume, plotter=plotter)
p.show()
[docs]
def showClinicalInfo(self):
"""
`myPatient.showClinicalInfo` shows the clinical information for the patient in a tabular form.
"""
ids=list(self.nrrdHeader.keys())
keys=list(self.nrrdHeader[ids[0]].keys())
keys.sort(key=str.lower, reverse=True)
od=[]
for i in range(len(ids)):
od.append(self.nrrdHeader[ids[i]])
rows=[]
for i in range(len(keys)):
rows.append([od[j].get(keys[i]) for j in range(len(od))])
rows[i].insert(0,keys[i])
# Print the rows in a tabular form
print(tabulate(rows, headers=[od[j].get('region') for j in range(len(od))]))
[docs]
def computeGeometricalFeatures(self):
"""
`myPatient.computeGeometricalFeatures` computes the geometrical features for the patient:
- Volume ratios
- Bounding boxes
- Sizes
"""
#compute the geometrical features
volumeList=[]
boundsList=[]
sizesList=[]
# here we can append wall and lumen ????
regions = self.info['AvailableRegions']
regions.append('wall')
regions.append('lumen')
for reg in regions:
volume, bounds,sizes = nr.computeVolumetricData(self.clinicalData,mat=reg,returnBounds=True)
volumeList.append((reg, volume))
boundsList.append((reg, bounds))
sizesList.append((reg, sizes))
volumeDict=collections.OrderedDict(volumeList)
volumeDict['LRNC+IPH']=volumeDict['LRNC']+volumeDict['IPH'] #sum the two volumes
self.info['Volume']=volumeDict
boundDict=collections.OrderedDict(boundsList)
self.info['Bounds']=boundDict
sizeDict=collections.OrderedDict(sizesList)
self.info['Sizes']=sizeDict
sizeDict['LRNC+IPH']=[max(sizeDict['LRNC'][i], sizeDict['IPH'][i]) for i in range(len(sizeDict['LRNC']))] #sum the two volumes
return
[docs]
def showInfo(self):
"""
`myPatient.showInfo` shows the additional information for the patient in a tabular form.
"""
for mainkey in self.info.keys():
print("\n")
od = self.info[mainkey]
if isinstance(od, collections.OrderedDict):
rows = [[key, od[key]] for key in od.keys()]
print(tabulate(rows, headers=['Key', mainkey]))
else:
print(mainkey,od)
print("\n")
[docs]
def loadFEMJSON(self):
"""
`myPatient.loadFEMJSON` loads the FEM JSON file for the patient.
In particular it loads:
- PlaqueSliceCenter
- NodesGlobal
- NodesPlaqueRegion
- Classic FEA steps
- Growth FEA steps
"""
# load the json file
json_files = glob.glob(os.path.join(pres.pathForResults(self)+self.name, '*.json'))
#print(json_files)
with open(json_files[0], 'r') as f:
data = json.load(f)
if 'PlaqueSliceCenter' in data.keys():
self.plaqueSliceCenter=data['PlaqueSliceCenter']
if data['PlaqueSliceCenter'] == [0,0,0]:
print("PlaqueSliceCenter may not be defined: "+ str(data['PlaqueSliceCenter']) )
if 'NodesGlobal' in data.keys():
# shift by 1 the nodes index
self.nodesIdx=collections.OrderedDict()
self.nodesIdx.calc=np.array(data['NodesGlobal']['CALC'])-1
self.nodesIdx.wall=np.array(data['NodesGlobal']['Wall'])-1
if 'FIBR' in data['NodesGlobal']:
self.nodesIdx.fibr=np.array(data['NodesGlobal']['FIBR'])-1
else :
self.nodesIdx.fibr=[]
self.nodesIdx.lrnciph=np.array(data['NodesGlobal']['LRNC+IPH'])-1
if 'NodesPlaqueRegion' in data.keys():
# shift by 1 the nodes index
self.nodesPlaqueRegionIdx=collections.OrderedDict()
calcNodes=np.array(data['NodesPlaqueRegion']['CALC'])-1
wallNodes=np.array(data['NodesPlaqueRegion']['Wall'])-1
if 'FIBR' in data['NodesGlobal']:
fibrNodes=np.array(data['NodesPlaqueRegion']['FIBR'])-1
else :
fibrNodes=[]
lrncNodes=np.array(data['NodesPlaqueRegion']['LRNC+IPH'])-1
self.nodesPlaqueRegionIdx.calc=calcNodes
self.nodesPlaqueRegionIdx.wall=wallNodes
self.nodesPlaqueRegionIdx.lrnciph=lrncNodes
self.nodesPlaqueRegionIdx.fibr=fibrNodes
self.nodesPlaqueRegionIdx.all=np.concatenate((calcNodes,wallNodes,lrncNodes,fibrNodes))
if 'Classic' in data.keys():
classicSteps={'MaxLoad':data['Classic']['MaxLoad'], 'MAP':data['Classic']['MAP']+1}#this +1 is a long story, have a look to the simualtion step
self.classicStep = collections.OrderedDict(classicSteps)
if 'Growth' in data.keys():
#growthSteps={'MaxLoad':data['Growth']['MaxLoad']+1, 'MAP':data['Growth']['MAP']+1, 'GrowthFirstPartialStep':data['Growth']['GrowthFirstPartialStep'], 'GrowthFirstCompleteStep':data['Growth']['GrowthFirstCompleteStep'], 'GrowthLastCompleteStep':data['Growth']['GrowthLastCompleteStep']}
growthSteps = {
'MaxLoad': safe_increment(data['Growth']['MaxLoad']),
'MAP': safe_increment(data['Growth']['MAP']),
'GrowthFirstPartialStep': safe_increment(data['Growth']['GrowthFirstPartialStep']),
'GrowthFirstCompleteStep': safe_increment(data['Growth']['GrowthFirstCompleteStep']),
'GrowthLastCompleteStep': safe_increment(data['Growth']['GrowthLastCompleteStep'])
}
self.growthStep = collections.OrderedDict(growthSteps)
return data
[docs]
def showFEMJSON(self):
"""
`myPatient.showFEMJSON` shows the FEM JSON file for the patient in a tabular form."""
# Check if data is a dictionary
data=self.FEMinfo
if isinstance(data, dict):
# Create a list of rows where each row is a list of a key and its corresponding value in data
rows = [[key, data[key]] for key in data.keys()]
noDictRow=[]
for row in rows:
rowDict=row[-1]
if isinstance(rowDict, dict):
singleRow = [[key, rowDict[key]] for key in rowDict.keys()]
for i in range(len(singleRow)):
line = singleRow[i]
if isinstance(line[-1], list) and len(line[-1])>5:
singleRow[i][-1]=line[-1][0:5] + ['...'] + [str(len(line[-1]))+ ' elements']
print(tabulate(singleRow, headers=['\033[1m' +row[0]+'\033[0m', '']))
print('\n')
# Print the rows in a tabular form
else:
noDictRow.append(row)
print(tabulate(noDictRow, headers=['Key', 'Value']))
else:
print("Data is not a dictionary")
[docs]
def updateFEMJSON(self):
"""
`myPatient.updateFEMJSON` updates the FEM JSON file for the patient adding:
- ClinicalVolumes
- ClinicalSize
- ClinicalBoundingBoxVolume
- ClinicalDispersionRatio
- Radii
- PlaqueSliceCenter
- MinLumenRadii
"""
if self.name in ['A_643', 'B_710', 'C_753', 'D_779', 'E_735', 'F_740', 'G_757', 'H_759', 'I_771', 'J_CP20-009', 'K_860', 'L_1081', 'M_1055', 'N_785', 'O_809', 'P_817', 'Q_874', 'R_780']:
if not self.info:
self.loadClinicalData(allData=False) #load lumen, wall and compisition
# probably htey need to be rescaled to mm
## make function to extract volume from wall, CALC, LRNC+IPH
# for the wall use the main volumetric mesh. For the others use the composition
self.FEMinfo['ClinicalVolumes']=self.info['Volume']
## make a function to extract max and min bounds from wall, CALC, LRNC+IPH
self.FEMinfo['ClinicalSize']=self.info['Sizes']
vols = [[key, self.info['Sizes'][key]] for key in self.info['Sizes'].keys()]
for k in range(len(vols)):
vols[k]=[vols[k][0], vols[k][1][0]*vols[k][1][1]*vols[k][1][2]]
self.FEMinfo['ClinicalBoundingBoxVolume']=collections.OrderedDict(vols)
for k in range(len(vols)):
actualVol = self.FEMinfo['ClinicalVolumes'][vols[k][0]]
boundVol = vols[k][-1] #!!! this need the previous computation also
vols[k]=[vols[k][0], actualVol/boundVol]
self.FEMinfo['ClinicalDispersionRatio']=collections.OrderedDict(vols)
self.FEMinfo['Radii']=self.radii
if not hasattr(self,'plaqueSliceCenter'):
self.plaqueSliceCenter=[0,0,0]
self.FEMinfo['PlaqueSliceCenter']=self.plaqueSliceCenter
self.FEMinfo['MinLumenRadii']=self.computeRadii(maxObstruction=True)
json_files = glob.glob(os.path.join(pres.pathForResults(self)+self.name, '*.json'))
with open(json_files[0], 'w') as f:
json.dump(self.FEMinfo,f,indent=True)
[docs]
def computeRadii(self, maxObstruction = False):
"""
`myPatient.computeRadii`
Computes the average radii (internel, external, wall thickness) based on the mesh area.
`maxObstruction=True` is used to compute the lumen in the plaque region (minimal region) instead of the average.
"""
extRads=[]
intRads=[]
wallThicks=[]
x1, x2, y1, y2, z1, z2 = self.mesh.bounds
dataset={}
steps=15
center = self.mesh.center
if maxObstruction:
center[2] = self.plaqueSliceCenter[2]
slice = self.mesh.slice(normal=[0,0,1], origin=center)
filled=slice.fill_holes(hole_size=1000)
annulus=slice.area
externalArea=filled.area
lumenArea=externalArea-annulus
extRad=np.sqrt(externalArea/np.pi)
lumRad=np.sqrt(lumenArea/np.pi)
wallThickness=extRad-lumRad
#print("External radius: ", extRad)
dataset['InternalRadius']=lumRad
dataset['ExternalRadius']=extRad
dataset['WallThickness']=wallThickness
else:
for i in range(steps):
center[2] = z1 + (z2-z1)*(0.05+i/steps)
slice = self.mesh.slice(normal=[0,0,1], origin=center)
filled=slice.fill_holes(hole_size=1000)
annulus=slice.area
externalArea=filled.area
lumenArea=externalArea-annulus
extRad=np.sqrt(externalArea/np.pi)
extRads.append(extRad)
lumRad=np.sqrt(lumenArea/np.pi)
intRads.append(lumRad)
wallThickness=extRad-lumRad
wallThicks.append(wallThickness)
#print("External radius: ", extRad)
dataset['InternalRadius']=np.mean(intRads)
dataset['ExternalRadius']=np.mean(extRads)
dataset['WallThickness']=np.mean(wallThicks)
return collections.OrderedDict(dataset)
[docs]
def updateMaterialsFromNodesInfo(self):
"""
It load materials as a segmented field based on the numerical id of the nodes."""
tmp = np.ones(len(self.mesh.points))
calcIdx=self.nodesIdx.calc
lrncIdx=self.nodesIdx.lrnciph
fibrIdx=self.nodesIdx.fibr
for i in calcIdx:
tmp[i] = 0
for i in lrncIdx:
tmp[i] = 2
for i in fibrIdx:
tmp[i] = 3
self.mesh.point_data['matsSegmented'] = tmp
[docs]
def showHist(self,fieldName='u',step=1, GRLab='Classic', scaling=None,limitRegions=None,bins=150):
"""
`myPatient.showHist` shows the histogram of the `fieldName` results for the patient. It works by calling the `getHist` function.
## Parameters:
`fieldName` (str): Name of the field to show. Default is 'u'.
`step` (int): Load step value for the field. Default is 1.
`GRLab` (str): The label for the field. Default is 'Classic' otherwise you can specify 'G&R'.
`scaling` (string): Scaling factor. Default is None.
`limitRegions` (str): The region to limit the field values. Default is None but you can specify:
- 'VascularTissue' to filter data for all the vascular wall
- 'PlaqueRegion' to filter data for al tissues in the plaque region
- 'PlaqueRegionVascularTissue' to filter data for the vascular wall in plaque region
`bins` (int): Number of bins for the histogram. Default is 150. """
ax = plt.subplot()
ax = self.getHist(fieldName=fieldName,step=step, GRLab=GRLab, scaling=scaling,limitRegions=limitRegions,ax=ax,bins=bins)
plt.show()
[docs]
def compareHist(self,fieldName='u',classicStep=1,growthStep=1, scaling=None,limitRegions=None,bins=200, clim=None, density=True):
"""
`myPatient.compareHist` compares the histogram of the `fieldName` results between the classic and growth results respectively for the `classicStep` and `growthStep` steps.
It works by calling the `getCompareHist` function.
## Parameters:
`fieldName` (str): Name of the field to compare. Default is 'u'.
`classicStep` (int): Step value for classic field. Default is 1.
`growthStep` (int): Step value for growth field. Default is 1.
`scaling` (string): Scaling factor. Default is None.
`limitRegions` (str): The region to limit the field values. Default is None but you can specify:
- 'VascularTissue' to filter data for all the vascular wall
- 'PlaqueRegion' to filter data for al tissues in the plaque region
- 'PlaqueRegionVascularTissue' to filter data for the vascular wall in plaque region
`bins` (int): Number of bins for the histogram. Default is 200.
`clim` (list): List of two values representing the color limit. Default is None.
`density` (bool): If True, the histogram will be displayed as a density plot. Default is True."""
ax = plt.subplot()
self.getCompareHist(fieldName=fieldName,classicStep=classicStep,growthStep=growthStep, scaling=scaling,limitRegions=limitRegions,bins=bins, clim=clim, density=density, ax = ax)
plt.show()
[docs]
def getCompareHist(self,fieldName='u',classicStep=1,growthStep=1, scaling=None,limitRegions=None,bins=200, clim=None, density=True, ax = None, colors=['blue','red'], opacity=0.5):
"""
`myPatient.getCompareHist` retrieves the histogram of the `fieldName` results between the classic and growth results respectively for the `classicStep` and `growthStep` steps.
## Parameters:
`fieldName` (str): Name of the field to compare. Default is 'u'.
`classicStep` (int): Step value for classic field. Default is 1.
`growthStep` (int): Step value for growth field. Default is 1.
`scaling` (string): Scaling factor. Default is None.
`limitRegions` (str): The region to limit the field values. Default is None but you can specify:
- 'VascularTissue' to filter data for all the vascular wall
- 'PlaqueRegion' to filter data for al tissues in the plaque region
- 'PlaqueRegionVascularTissue' to filter data for the vascular wall in plaque region
`bins` (int): Number of bins for the histogram. Default is 200.
`clim` (list): List of two values representing the color limit. Default is None.
`density` (bool): If True, the histogram will be displayed as a density plot. Default is True."""
if isinstance(classicStep, (int, float)):
classicStepValue=classicStep
else:
if classicStep == 'MaxLoad':
classicStepValue=self.classicStep['MaxLoad']
elif classicStep == 'MAP':
classicStepValue=self.classicStep['MAP']
if isinstance(growthStep, (int, float)):
growthStepValue=growthStep
else:
if growthStep == 'MaxLoad':
growthStepValue=self.growthStep['MaxLoad']
elif growthStep == 'MAP':
growthStepValue=self.growthStep['MAP']
self.getHist(fieldName=fieldName,step=growthStepValue, GRLab='G&R', scaling=scaling,limitRegions=limitRegions,ax=ax, bins=bins, label='Growth',alpha=1, density=density, color=colors[1])
self.getHist(fieldName=fieldName,step=classicStepValue, GRLab='Classic', scaling=scaling,limitRegions=limitRegions,ax=ax, bins=bins, label='Classic', alpha=1, density=density, color=colors[0])
self.getHist(fieldName=fieldName,step=growthStepValue, GRLab='G&R', scaling=scaling,limitRegions=limitRegions,ax=ax, bins=bins, alpha=opacity, density=density, color=colors[1])
ax.legend(loc='upper right')
if clim is not None:
ax.set_xlim(clim)
return
[docs]
def getHist(self,fieldName='u',step=1, GRLab='Classic', scaling=None,limitRegions=None,ax=None,bins=30, label='', alpha=1, density = False, color='blue'):
"""
`myPatient.showHist` shows the histogram of the `fieldName` results for the patient. It works by calling the `getHist` function.
## Parameters:
`fieldName` (str): Name of the field to show. Default is 'u'.
`step` (int): Load step value for the field. Default is 1.
`GRLab` (str): The label for the field. Default is 'Classic' otherwise you can specify 'G&R'.
`scaling` (string): Scaling factor. Default is None.
`limitRegions` (str): The region to limit the field values. Default is None but you can specify:
- 'VascularTissue' to filter data for all the vascular wall
- 'PlaqueRegion' to filter data for al tissues in the plaque region
- 'PlaqueRegionVascularTissue' to filter data for the vascular wall in plaque region
`bins` (int): Number of bins for the histogram. Default is 150. """
if ax is None:
fig, ax = plt.subplots()
data=self.getFieldValues(fieldName=fieldName,step=step, GRLab=GRLab, scaling=scaling,limitRegions=limitRegions)
ax.hist(data, bins=bins, label=label, alpha=alpha, density=density, color=color)
return
[docs]
def getIndexes(self,isMap=True, printBool=False, limitRegions=None):
if limitRegions is None:
limitRegions='PlaqueRegionVascularTissue'
else:
limitRegions=limitRegions
if isMap:
classicStep='MAP'
growthStep='MAP'
mapString='SYS'
clim=[-10,60]
else:
classicStep='MaxLoad'
growthStep='MaxLoad'
mapString='MaxLoad'
clim=[-10,120]
if isinstance(classicStep, (int, float)):
classicStepValue=classicStep
else:
if classicStep == 'MaxLoad':
classicStepValue=self.classicStep['MaxLoad']
elif classicStep == 'MAP':
classicStepValue=self.classicStep['MAP']
if isinstance(growthStep, (int, float)):
growthStepValue=growthStep
else:
if classicStep == 'MaxLoad':
growthStepValue=self.growthStep['MaxLoad']
elif classicStep == 'MAP':
growthStepValue=self.growthStep['MAP']
stressClassic=self.getFieldValues(fieldName='σ1',step=classicStepValue, GRLab='Classic', scaling='kPa',limitRegions=limitRegions)
stressGR=self.getFieldValues(fieldName='σ1',step=growthStepValue, GRLab='G&R', scaling='kPa',limitRegions=limitRegions)
max90Pctile=np.percentile(stressClassic, 95)
ClassicPercentileVolume=len([x for x in stressClassic if x>max90Pctile])/len(stressClassic)
GRPercentileVolume=len([x for x in stressGR if x>max90Pctile])/len(stressGR)
indexes={}
indexes['MaxStressVolumeRedution']=100*(ClassicPercentileVolume-GRPercentileVolume)/ClassicPercentileVolume
indexes['ClassicVolumeHigher90Percentile']=100*ClassicPercentileVolume
indexes['GRVolumeHigher90Percentile']=100*GRPercentileVolume
indexes['PercentileValues']=max90Pctile
indexes['MaxStressReduction']=100*(max(stressClassic)-max(stressGR))/max(stressClassic)
indexes['99PcMaxStressReduction']=100*(np.percentile(stressClassic, 99)-np.percentile(stressGR, 99))/np.percentile(stressClassic, 99)
indexes['MeanStressReduction']=100*(np.mean(stressClassic)-np.mean(stressGR))/np.mean(stressClassic)
indexes['MaxStressGrowth']=max(stressGR)
indexes['MaxStressClassic']=max(stressClassic)
indexes['MaxStress99PCGrowth']=np.percentile(stressGR, 99)
indexes['MaxStress99PCClassic']=np.percentile(stressClassic, 99)
indexes['MeanStressGrowth']=np.mean(stressGR)
indexes['MeanStressClassic']=np.mean(stressClassic)
indexes['STDStressGrowth']=np.std(stressGR)
indexes['STDStressClassic']=np.std(stressClassic)
if printBool:
print(str(self.name).replace("_"," "),end='& ')
print("%.2f" % indexes['MaxStressClassic'],end='& ')
print("%.2f" %indexes['MaxStressGrowth'],end='& ')
print("%.2f" %indexes['MaxStressReduction'],end='& ')
print("%.2f" % indexes['MeanStressClassic'],end='& ')
print("%.2f" %indexes['MeanStressGrowth'],end='& ')
print("%.2f" %indexes['MeanStressReduction'],end='& ')
print("%.2f" % indexes['ClassicVolumeHigher90Percentile'],end='& ')
print("%.2f" %indexes['GRVolumeHigher90Percentile'],end='& ')
print("%.2f" %indexes['MaxStressVolumeRedution'],end='\\\ \hline ')
print("")
return indexes
[docs]
def computeScaleFactor(scalingString=None):
if scalingString is None:
return 1
elif scalingString=='kPa':
return 1e-3
else:
return 0
# Function to safely increment fields if they are not empty
[docs]
def safe_increment(value):
if isinstance(value, list) and not value: # Check if it's an empty list
return None # or return value if you want to keep the empty list
elif isinstance(value, (int, float)):
return value + 1 # Increment numeric values
else:
return value