-
Notifications
You must be signed in to change notification settings - Fork 6
/
readGEFlow.py
169 lines (115 loc) · 7.66 KB
/
readGEFlow.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
import pydicom,os, glob, scipy.io, numpy, vtk, sys, saveVTK, math, eddyNoise
from clint.textui import colored
from vtk.util import numpy_support
def readGEFlow(inputFlags, PatientDataStruc):
print( colored.green("\nLooking for flow data... \n"))
print(colored.cyan("Mag folder is: " + PatientDataStruc.MagPath))
print(colored.cyan("R/L folder is: " + PatientDataStruc.FlowPathRL))
print(colored.cyan("A/P folder is: " + PatientDataStruc.FlowPathAP))
print(colored.cyan("S/I folder is: " + PatientDataStruc.FlowPathSI))
for folderNumber in range(0,4):
sliceLocation = []
# flow files list
lstFilesDCM = []
triggerTime = []
if folderNumber == 0:
folderPath = PatientDataStruc.MagPath
print(colored.cyan("Reading the Magnitude files."))
if folderNumber == 1:
if inputFlags.segmentation:
break
folderPath = PatientDataStruc.FlowPathRL
print(colored.cyan("Reading the flow files (R/L)."))
if folderNumber == 2:
if inputFlags.segmentation:
break
folderPath = PatientDataStruc.FlowPathAP
print(colored.cyan("Reading the flow files (A/P)."))
if folderNumber == 3:
if inputFlags.segmentation:
break
folderPath = PatientDataStruc.FlowPathSI
print(colored.cyan("Reading the flow files (S/I)."))
################## Reading time of flight files
# listing magnitude files
for dirName, subdirList, fileList in os.walk( folderPath + "/"):
for filename in fileList:
# if ".dcm" in filename.lower(): # check whether the file's DICOM
lstFilesDCM.append(os.path.join(dirName,filename))
ds = pydicom.read_file(lstFilesDCM[-1])
sliceLocation.append(ds.SliceLocation)
triggerTime.append(ds.TriggerTime)
# Get ref file
#RefDs = dicom.read_file(lstFilesDCM[0])
triggerTimeTemp = sorted(set(triggerTime), key=float)
#print(triggerTimeTemp)
#sliceLocationTemp = set(sliceLocation)
sliceLocationTemp = sorted(set(sliceLocation), key=float)
#print(sliceLocationTemp)
if folderNumber == 0:
#ConstPixelDims = (int(RefDs.Rows), int(RefDs.Columns),66, int(ds.CardiacNumberOfImages))
#print(ConstPixelDims)
ReadData = numpy.zeros(PatientDataStruc.MagVecSize, dtype=numpy.double)
for iFile in lstFilesDCM:
dsTemp = pydicom.read_file (iFile)
ReadData[:,:,sliceLocationTemp.index(dsTemp.SliceLocation),triggerTimeTemp.index(dsTemp.TriggerTime)]= dsTemp.pixel_array.astype('float')
magDataTemp = ReadData.mean(3)
if inputFlags.mat:
scipy.io.savemat(inputFlags.output + "/mag.mat", mdict={'magDataTemp': magDataTemp})
# numpy.save(args.output +"/mag", magDataTemp)
else:
if folderNumber == 1:
# ConstFlowPixelDims = (int(RefDs.Rows), int(RefDs.Columns), 66, 3, int(ds.CardiacNumberOfImages))
flowData = numpy.zeros(PatientDataStruc.FlowVecSize, dtype=numpy.double)
# print(PatientDataStruc.FlowVecSize)
for iFile in lstFilesDCM:
dsTemp = pydicom.read_file (iFile)
flowData[:,:,sliceLocationTemp.index(dsTemp.SliceLocation), folderNumber-1,triggerTimeTemp.index(dsTemp.TriggerTime)]= dsTemp.pixel_array.astype('float')
if inputFlags.mat:
scipy.io.savemat(inputFlags.output + "/vel.mat", mdict={'flowData': flowData})
#print(flowData.shape)
if inputFlags.segmentation is False:
### The combination of -x +y and -z and permuted x and y is working. Ali Aug24 2017
UOrg = inputFlags.velocitysign[0] * (flowData[:, :, :, inputFlags.velocityorder[0]].squeeze())
VOrg = inputFlags.velocitysign[1] * (flowData[:, :, :, inputFlags.velocityorder[1]].squeeze())
WOrg = inputFlags.velocitysign[2] * (flowData[:, :, :, inputFlags.velocityorder[2]].squeeze())
flowCorrected = numpy.zeros([flowData.shape[0], flowData.shape[1], flowData.shape[2],3,flowData.shape[4]])
if (inputFlags.eddycurrent is not None and inputFlags.randomnoise is not None):
# flowCorrectedtemp2 = eddyNoise.randNoiseV2(magDataTemp, UOrg, VOrg, WOrg, int(inputFlags.randomnoise), 0)
flowCorrectedtemp = eddyNoise.randNoise(UOrg, VOrg, WOrg, inputFlags.randomnoise)
flowCorrected = eddyNoise.eddyCurrentCorrection(inputFlags, flowCorrectedtemp[:,:,:,0], flowCorrectedtemp[:,:,:,1], flowCorrectedtemp[:,:,:,2], magDataTemp, inputFlags.eddythreshold, inputFlags.eddyplane, STDPower=2, plotEddyPlane=0, plotPlain=20)
if (inputFlags.eddythreshold is not None and inputFlags.randomnoise is None):
flowCorrected = eddyNoise.eddyCurrentCorrection(inputFlags, UOrg, VOrg, WOrg, magDataTemp, inputFlags.eddythreshold, inputFlags.eddyplane, STDPower=2, plotEddyPlane=0, plotPlain=20)
if (inputFlags.randomnoise is not None and inputFlags.eddythreshold is None):
flowCorrected = eddyNoise.randNoise(UOrg, VOrg, WOrg, inputFlags.randomnoise)
if (inputFlags.eddythreshold is None and inputFlags.randomnoise is None):
flowCorrected[:, :, :,0, :] = UOrg
flowCorrected[:, :, :,1, :] = VOrg
flowCorrected[:, :, :,2, :] = WOrg
vTemp = numpy.amax(flowCorrected, axis=4)
#flowCorrected.mean(4)
vTemp2 = numpy.sqrt(vTemp[:,:,:,0]**2 + vTemp[:,:,:,1]**2+ vTemp[:,:,:,2]**2)
vCMRA = numpy.multiply(vTemp2, magDataTemp)
print(colored.green("\nGetting ready to write files... This takes a little bit of time"))
magSize = magDataTemp.shape
totalNodes = magSize[0] * magSize[1] * magSize[2]
if (inputFlags.vtk == False and inputFlags.mat == False):
print(colored.yellow("We will ONLY save in npy format, since you didnt select your preference! (VTK or MAT)"))
numpy.save(inputFlags.output +"/FlowData", flowCorrected)
# if not inputFlags.segmentation:
# numpy.save(inputFlags.output +"/FlowData", flowCorrected)
if inputFlags.vtk:
if inputFlags.segmentation:
saveVTK.saveVTKSeg(magDataTemp,False,False, PatientDataStruc.PixelSize, totalNodes, inputFlags.output)
else:
saveVTK.saveVTK(vCMRA, flowCorrected, PatientDataStruc.PixelSize, totalNodes, inputFlags.output)
saveVTK.saveVTKSeg(vCMRA,True,False, PatientDataStruc.PixelSize, totalNodes, inputFlags.output)
saveVTK.saveVTKSeg(magDataTemp,False,False, PatientDataStruc.PixelSize, totalNodes, inputFlags.output)
if inputFlags.mat:
if inputFlags.segmentation:
with open(inputFlags.output + "/FlowData.mat", 'wb') as matlabFile:
scipy.io.savemat(matlabFile, mdict={'magnitude': magDataTemp})
else:
with open(inputFlags.output + "/FlowData.mat", 'wb') as matlabFile:
scipy.io.savemat(matlabFile, mdict={'velocity': flowData})
scipy.io.savemat(matlabFile, mdict={'magnitude': magDataTemp})