-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathCoronaryPerfusionTerritories.py
175 lines (144 loc) · 6.2 KB
/
CoronaryPerfusionTerritories.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
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
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 CoronaryPerfusionTerritories():
def __init__(self,Args):
#Input argumenets
self.Args=Args
#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 and CenterlineFileNames[i].find("left")>=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):
#Get the areas for all of the surface caps
Centroids=self.Get_Centroids(self.CenterlineFileNamesPaths,self.CenterlineFileNames)
#Read the LV mesh file
print ("--- Reading Left Ventricle Myocardial Blood Flow Data: %s"%self.VentricleFileName)
VentricleMesh=self.Read_Vtu(self.VentricleFileName)
#Get the Nodes and Scalars of the LV Mesh
print ("--- Getting Coordinates from the Left Ventricle Mesh")
VentricleCoordinates,ImageScalars=self.Get_Point_Data(VentricleMesh)
#Read the Centerline Coordinates
print ("--- Reading Centerline Files")
CenterlineCoordinates=self.Read_CL_Coords(self.CenterlineFileNamesPaths)
#For each LV node, find the shortest distance to the CL
print ("--- Getting the shortest distance for each Ventricle Node to closest centerline")
ClosestCenterLines=self.Get_Shortest_Distance(CenterlineCoordinates,VentricleCoordinates)
#Now write the data with a new territory division
print ("--- Writing the territory maps")
self.Write_Territories(VentricleMesh,ClosestCenterLines)
def Write_Territories(self,LV_mesh,Closest_CL):
#Create New Array in LV mesh with CL names
territories=vtk.vtkIntArray()
territories.SetNumberOfComponents(1)
territories.SetName("TerritoryMaps")
N=LV_mesh.GetNumberOfPoints()
for i in range(N):
index=self.filenames.index(Closest_CL[i])
territories.InsertNextValue(index)
LV_mesh.GetPointData().AddArray(territories)
writer=vtk.vtkXMLUnstructuredGridWriter()
writer.SetInputData(LV_mesh)
writer.SetFileName(self.Args.OutputFileName)
writer.Update()
def Get_Shortest_Distance(self,CL_coords,LV_coords):
#For each Node, Loop over all of the centerlines
N=len(LV_coords)
Closest_CL=[]
count_p=0.05
progress_old=-1
for i,LV_coord in enumerate(LV_coords):
progress_=self.PRINT_PROGRESS(i,N,progress_old)
progress_old=progress_
distance={}
#Loop over all of the centerlines
for Key in CL_coords.keys():
distance[Key]=np.min(DISTANCE.cdist(CL_coords[Key],[LV_coords[i]]))
#Store the Centerline that is closest to the point
Closest_CL.append(min(distance,key=distance.get))
return Closest_CL
def Read_CL_Coords(self,filenames):
CL_coords={}
for filename in filenames:
print ("-----%s"%filename)
filename_short=filename.split("/")[-1]
reader=vtk.vtkXMLPolyDataReader()
reader.SetFileName(filename)
reader.Update()
CL_data=reader.GetOutput()
N=CL_data.GetNumberOfPoints()
points_=np.zeros(shape=(N,3))
for i in range(N):
points_[i]=CL_data.GetPoint(i)
CL_coords[filename_short]=points_
return CL_coords
def Get_Point_Data(self,LV_mesh):
print ("------ Getting LV Coordinates and MBF")
N=LV_mesh.GetNumberOfPoints()
Coords=np.zeros(shape=(N,3))
imageScalars=np.zeros(N)
for i in range(N):
coord_=LV_mesh.GetPoint(i)
Coords[i,0]=coord_[0]
Coords[i,1]=coord_[1]
Coords[i,2]=coord_[2]
imageScalars[i]=LV_mesh.GetPointData().GetArray(self.Args.ArrayName).GetValue(i)
return Coords,imageScalars
def Read_Vtu(self,filename):
reader=vtk.vtkXMLUnstructuredGridReader()
reader.SetFileName(filename)
reader.Update()
return reader.GetOutput()
def Get_Centroids(self,surface_files,filenames):
print ("-"*50)
print ("------ Computing Centroids of Surface Outlets")
Centroids={}
count=0
for surface_file in surface_files:
print ("----Processing %s"%filenames[count])
#Read the Surface File
reader=vtk.vtkXMLPolyDataReader()
reader.SetFileName(surface_file)
reader.Update()
#Get the Centroid
centerOfMassFilter=vtk.vtkCenterOfMass()
centerOfMassFilter.SetInputData(reader.GetOutput())
centerOfMassFilter.SetUseScalarsAsWeights(False)
centerOfMassFilter.Update()
center_=centerOfMassFilter.GetCenter()
Centroids[filenames[count]]=[center_[0],center_[1],center_[2]]
del reader,centerOfMassFilter
count+=1
return Centroids
def PRINT_PROGRESS(self,i,N,progress_old):
progress_=(int((float(i)/N*100+0.5)))
if progress_%10==0 and progress_%10!=progress_old: print (" Progress: %d%%"%progress_)
return progress_%10
if __name__=="__main__":
#Description
parser = argparse.ArgumentParser(description="This script will computed terriorites of the myocardium using coronary vessels.")
#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")
#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")
#Output argumenets
parser.add_argument('-ofile', '--OutputFileName', type=str, required=True, dest="OutputFileName",help="The output filename of the projected surface")
args=parser.parse_args()
CoronaryPerfusionTerritories(args).main()