-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathImageAnalysisMyocardiumPerfusionTerritoriesPlot.py
132 lines (106 loc) · 5.23 KB
/
ImageAnalysisMyocardiumPerfusionTerritoriesPlot.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
#This script will read the territories file
#and write a tecplot file with histogram plots
import vtk
import numpy as np
import scipy as sp
import os
from glob import glob
from scipy.spatial import distance as DISTANCE
import argparse
class ImageAnalysisMyocardiumPerfusionTerritoriesPlot():
def __init__(self,Args):
#Input argumenets
self.Args=Args
#If output filename is not defined, create your own
if self.Args.OutputFileName is None:
InputFileNameStripped=self.Args.InputFileName.split("/")[-1]
OutputFolder =self.Args.InputFileName.replace(InputFileNameStripped,"")
OutputFileName=OutputFolder+"MyocardiumPerfusionTerritoriesPlot.dat"
self.Args.OutputFileName=OutputFileName
#Name of the Coronary Artery Tags
self.MBF_data={"RCA":[], "LAD":[], "Diag":[], "Septal":[], "LCx":[], "LAD_Diag_Sep":[]}
#The LV filename that contains MBF
self.VentricleFileName=self.Args.InputFileName
#Coronary Centerline Files
CenterlineFileNamesPaths=sorted(glob("%s/wall_*cl.vtp"%self.Args.CenterlinesFolder))
CenterlineFileNames=[filename.split("/")[-1] for filename in CenterlineFileNamesPaths]
#Separate the Coronary Territories into Left and Right side
self.CenterlineFileNamesPaths=[]
self.CenterlineFileNames=[]
for i in range(len(CenterlineFileNames)):
if CenterlineFileNames[i].find("RCA")>=0:
self.CenterlineFileNames.append(CenterlineFileNames[i])
self.CenterlineFileNamesPaths.append(CenterlineFileNamesPaths[i])
if CenterlineFileNames[i].find("LCA")>=0:
self.CenterlineFileNames.append(CenterlineFileNames[i])
self.CenterlineFileNamesPaths.append(CenterlineFileNamesPaths[i])
def main(self):
#Read the LV mesh file that contains the territory mappings and MBF
print ("--- Reading Left Ventricle Myocardial Blood Flow Data and Territory Maps: %s"%self.VentricleFileName)
VentricleMesh=self.Read_Vtu(self.VentricleFileName)
#Compute the vessel-specific MBF
MBF_vessel=self.compute_mbf(VentricleMesh)
#Write the tecplot file with the Vessel S
self.write_MBF_data(MBF_vessel)
def compute_mbf(self,VentricleMesh):
#Separate the vessels into various territories
MBF_data=self.MBF_data
#Get the number of points in the array
N=VentricleMesh.GetNumberOfPoints()
#Loop over all points and append to the correct array
for i in range(N):
value_=VentricleMesh.GetPointData().GetArray(self.Args.ArrayName).GetValue(i)
if value_==0: continue
territory_idx_=VentricleMesh.GetPointData().GetArray("TerritoryMaps").GetValue(i)
filename_=self.CenterlineFileNames[territory_idx_]
if filename_.find("RCA")>=0: MBF_data["RCA"].append(value_)
elif filename_.find("lad")>=0: MBF_data["LAD"].append(value_)
elif filename_.find("diag")>=0:MBF_data["Diag"].append(value_)
elif filename_.find("septal")>=0: MBF_data["Septal"].append(value_)
elif filename_.find("lcx")>=0: MBF_data["LCx"].append(value_)
elif filename_.find("LCA")>=0: MBF_data["LAD"].append(value_)
else:
print ("A value was not found for any territory")
exit(1)
if filename_.find("LCA")>=0 and filename_.find("lcx")<0:
MBF_data["LAD_Diag_Sep"].append(value_)
return MBF_data
def write_MBF_data(self,MBF_data):
outfile=open(self.Args.OutputFileName,'w')
outfile.write('TITLE="MBF(mL/min/100g)"\n')
outfile.write('VARIABLES = "Region","MBF"\n')
if self.Args.PostCABG==0:
count=1
tag="pre"
else:
count=1.25
tag="post"
for key,elem in MBF_data.items():
outfile.write('Zone T= "%s_%s", I=3, F=POINT\n'%(key,tag))
mean_=np.mean(elem)
stdev_=np.std(elem)
outfile.write("%.02f %.05f\n"%(count,mean_+stdev_))
outfile.write("%.02f %.05f\n"%(count,mean_))
outfile.write("%.02f %.05f\n"%(count,mean_-stdev_))
count+=1
outfile.close()
def Read_Vtu(self,filename):
reader=vtk.vtkXMLUnstructuredGridReader()
reader.SetFileName(filename)
reader.Update()
return reader.GetOutput()
if __name__=="__main__":
#Description
parser = argparse.ArgumentParser(description="This script will generate tecplot files from perfusion territory vtu files.")
#Input filename of the perfusion map
parser.add_argument('-ifile', '--InputFileName', type=str, required=True, dest="InputFileName",help="Volumetric Mesh that contains the Myocardial Blood Flow Data and Territory Maps")
#Input folder that contains the coronary centerlines
parser.add_argument('-centerlinesfolder', '--CenterlinesFolder', type=str, required=True, dest="CenterlinesFolder",help="Folder that contains the centerline files")
#Array Name of the Data
parser.add_argument('-arrayname', '--ArrayName', type=str, required=False,default="ImageScalars", dest="ArrayName",help="The name of the array containing the MBF values")
#Pre- or Post-CABG
parser.add_argument('-postCABG', '-PostCABG', type=int, required=True, default=0, dest="PostCABG", help="0 for pre-cabg (default) and 1 for post-cabg")
#Output argumenets
parser.add_argument('-ofile', '--OutputFileName', type=str, required=False, dest="OutputFileName",help="The output filename to store the perfusion plots and histograms")
args=parser.parse_args()
ImageAnalysisMyocardiumPerfusionTerritoriesPlot(args).main()