-
Notifications
You must be signed in to change notification settings - Fork 7
/
volReader.py
executable file
·302 lines (255 loc) · 12.2 KB
/
volReader.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
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
#!/usr/bin/env python
#
# Aaron Y. Lee MD MSCI (University of Washington) Copyright 2019
#
# Code ported from Markus Mayer's excellent work (https://www5.cs.fau.de/research/software/octseg/)
#
# Also thanks to who contributed to the original openVol.m in Markus's project
# Radim Kolar, Brno University, Czech Republic
# Kris Sheets, Retinal Cell Biology Lab, Neuroscience Center of Excellence, LSU Health Sciences Center, New Orleans
import struct, array, datetime, codecs
import numpy as np
from collections import OrderedDict
class volFile():
def __init__(self, filename):
"""
Parses Heyex Spectralis *.vol files.
Args:
filename (str): Path to vol file
Returns:
volFile class
"""
self.__parseVolFile(filename)
@property
def oct(self):
"""
Retrieve OCT volume as a 3D numpy array.
Returns:
3D numpy array with OCT intensities as 'uint8' array
"""
return self.wholefile["cScan"]
@property
def irslo(self):
"""
Retrieve IR SLO image as 2D numpy array
Returns:
2D numpy array with IR reflectance SLO image as 'uint8' array.
"""
return self.wholefile["sloImage"]
@property
def grid(self):
"""
Retrieve the IR SLO pixel coordinates for the B scan OCT slices
Returns:
2D numpy array with the number of b scan images in the first dimension
and x_0, y_0, x_1, y_1 defining the line of the B scan on the pixel
coordinates of the IR SLO image.
"""
wf = self.wholefile
grid = []
for bi in range(len(wf["slice-headers"])):
bscanHead = wf["slice-headers"][bi]
x_0 = int(bscanHead["startX"] / wf["header"]["scaleXSlo"])
x_1 = int(bscanHead["endX"] / wf["header"]["scaleXSlo"])
y_0 = int(bscanHead["startY"] / wf["header"]["scaleYSlo"])
y_1 = int(bscanHead["endY"] / wf["header"]["scaleYSlo"])
grid.append([x_0, y_0, x_1, y_1])
return grid
def renderIRslo(self, filename, renderGrid=False):
"""
Renders IR SLO image as a PNG file and optionally overlays grid of B scans
Args:
filename (str): filename to save IR SLO image
renderGrid (bool): True will render red lines for the location of the B scans.
Returns:
None
"""
from PIL import Image, ImageDraw
wf = self.wholefile
a = np.copy(wf["sloImage"])
if renderGrid:
a = np.stack((a,)*3, axis=-1)
a = Image.fromarray(a)
draw = ImageDraw.Draw(a)
grid = self.grid
for (x_0, y_0, x_1, y_1) in grid:
draw.line((x_0,y_0, x_1, y_1), fill=(255,0,0), width=3)
a.save(filename)
else:
Image.fromarray(a).save(filepre)
def renderOCTscans(self, filepre = "oct", renderSeg=False):
"""
Renders OCT images a PNG file and optionally overlays segmentation lines
Args:
filepre (str): filename prefix. OCT Images will be named as "<prefix>-001.png"
renderSeg (bool): True will render colored lines for the segmentation of the RPE, ILM, and NFL on the B scans.
Returns:
None
"""
from PIL import Image
wf = self.wholefile
for i in range(wf["cScan"].shape[0]):
a = np.copy(wf["cScan"][i])
if renderSeg:
a = np.stack((a,)*3, axis=-1)
for li in range(wf["segmentations"].shape[0]):
for x in range(wf["segmentations"].shape[2]):
a[int(wf["segmentations"][li,i,x]),x, li] = 255
Image.fromarray(a).save("%s-%03d.png" % (filepre, i))
def __parseVolFile(self, fn):
wholefile = OrderedDict()
decode_hex = codecs.getdecoder("hex_codec")
with open(fn, "rb") as fin:
header = OrderedDict()
header["version"] = fin.read(12)
header["octSizeX"] = struct.unpack("I", fin.read(4))[0] # lateral resolution
header["numBscan"] = struct.unpack("I", fin.read(4))[0]
header["octSizeZ"] = struct.unpack("I", fin.read(4))[0] # OCT depth
header["scaleX"] = struct.unpack("d", fin.read(8))[0]
header["distance"] = struct.unpack("d", fin.read(8))[0]
header["scaleZ"] = struct.unpack("d", fin.read(8))[0]
header["sizeXSlo"] = struct.unpack("I", fin.read(4))[0]
header["sizeYSlo"] = struct.unpack("I", fin.read(4))[0]
header["scaleXSlo"] = struct.unpack("d", fin.read(8))[0]
header["scaleYSlo"] = struct.unpack("d", fin.read(8))[0]
header["fieldSizeSlo"] = struct.unpack("I", fin.read(4))[0] # FOV in degrees
header["scanFocus"] = struct.unpack("d", fin.read(8))[0]
header["scanPos"] = fin.read(4)
header["examTime"] = struct.unpack("l", fin.read(8))[0] / 1e7 # (second) from 01.01.1601 => divide by 10^7 becasue the value is originally in the unit of 100nsec.
header["examTime"] = datetime.datetime.utcfromtimestamp(header["examTime"])
header["examTime"] = header["examTime"] - (datetime.date(1970,1,1) - datetime.date(1601,1,1))
header["scanPattern"] = struct.unpack("I", fin.read(4))[0]
header["BscanHdrSize"] = struct.unpack("I", fin.read(4))[0]
header["ID"] = fin.read(16)
header["ReferenceID"] = fin.read(16)
header["PID"] = struct.unpack("I", fin.read(4))[0]
header["PatientID"] = fin.read(21)
header["unknown2"] = fin.read(3)
header["DOB"] = struct.unpack("d", fin.read(8))[0] - 25569
header["DOB"] = datetime.datetime.utcfromtimestamp(header["DOB"] * 24 * 60 * 60 ) # needs to be checked
header["VID"] = struct.unpack("I", fin.read(4))[0]
header["VisitID"] = fin.read(24)
header["VisitDate"] = struct.unpack("d", fin.read(8))[0] - 25569
header["VisitDate"] = datetime.datetime.utcfromtimestamp(header["VisitDate"] * 24 * 60 * 60 ) # needs to be checked
header["GridType"] = struct.unpack("I", fin.read(4))[0]
header["GridOffset"] = struct.unpack("I", fin.read(4))[0]
wholefile["header"] = header
fin.seek(2048)
U = array.array("B")
U.fromstring(fin.read(header["sizeXSlo"] * header["sizeYSlo"]))
U = np.array(U).astype("uint8").reshape((header["sizeXSlo"],header["sizeYSlo"]))
wholefile["sloImage"] = U
sloOffset = 2048 + header["sizeXSlo"] * header["sizeYSlo"]
octOffset = header["BscanHdrSize"] + header["octSizeX"] * header["octSizeZ"] * 4
bscans = []
bscanheaders = []
segmentations = None
for i in range(header["numBscan"]):
fin.seek(16 + sloOffset + i * octOffset)
bscanHead = OrderedDict()
bscanHead["startX"] = struct.unpack("d", fin.read(8))[0]
bscanHead["startY"] = struct.unpack("d", fin.read(8))[0]
bscanHead["endX"] = struct.unpack("d", fin.read(8))[0]
bscanHead["endY"] = struct.unpack("d", fin.read(8))[0]
bscanHead["numSeg"] = struct.unpack("I", fin.read(4))[0]
bscanHead["offSeg"] = struct.unpack("I", fin.read(4))[0]
bscanHead["quality"] = struct.unpack("f", fin.read(4))[0]
bscanHead["shift"] = struct.unpack("I", fin.read(4))[0]
bscanheaders.append(bscanHead)
# extract OCT B scan data
fin.seek(header["BscanHdrSize"] + sloOffset + i * octOffset)
U = array.array("f")
U.fromstring(fin.read(4 * header["octSizeX"] * header["octSizeZ"]))
U = np.array(U).reshape((header["octSizeZ"],header["octSizeX"]))
# remove out of boundary
v = struct.unpack("f", decode_hex('FFFF7F7F')[0])
U[U == v] = 0
# log normalize
U = np.log(10000 * U + 1)
U = (255. * (np.clip(U, 0, np.max(U)) / np.max(U))).astype("uint8")
bscans.append(U)
# extract OCT segmentations data
fin.seek(256 + sloOffset + i * octOffset)
U = array.array("f")
U.fromstring(fin.read(4 * header["octSizeX"] * bscanHead["numSeg"]))
U = np.array(U)
U[U == v] = 0.
if segmentations == None:
segmentations = []
for j in range(bscanHead["numSeg"]):
segmentations.append([])
for j in range(bscanHead["numSeg"]):
segmentations[j].append(U[j*header["octSizeX"]:(j+1) * header["octSizeX"]].tolist())
wholefile["cScan"] = np.array(bscans)
wholefile["segmentations"] = np.array(segmentations)
wholefile["slice-headers"] = bscanheaders
self.wholefile = wholefile
@property
def fileHeader(self):
"""
Retrieve vol header fields
Returns:
Dictionary with the following keys
- version: version number of vol file definition
- numBscan: number of B scan images in the volume
- octSizeX: number of pixels in the width of the OCT B scan
- octSizeZ: number of pixels in the height of the OCT B scan
- distance: unknown
- scaleX: resolution scaling factor of the width of the OCT B scan
- scaleZ: resolution scaling factor of the height of the OCT B scan
- sizeXSlo: number of pixels in the width of the IR SLO image
- sizeYSlo: number of pixels in the height of the IR SLO image
- scaleXSlo: resolution scaling factor of the width of the IR SLO image
- scaleYSlo: resolution scaling factor of the height of the IR SLO image
- fieldSizeSlo: field of view (FOV) of the retina in degrees
- scanFocus: unknown
- scanPos: Left or Right eye scanned
- examTime: Datetime of the scan (needs to be checked)
- scanPattern: unknown
- BscanHdrSize: size of B scan header in bytes
- ID: unknown
- ReferenceID
- PID: unknown
- PatientID: Patient ID string
- DOB: Date of birth
- VID: unknown
- VisitID: Visit ID string
- VisitDate: Datetime of visit (needs to be checked)
- GridType: unknown
- GridOffset: unknown
"""
return self.wholefile["header"]
def bScanHeader(self, slicei):
"""
Retrieve the B Scan header information per slice.
Args:
slicei (int): index of B scan
Returns:
Dictionary with the following keys
- startX: x-coordinate for B scan on IR. (see getGrid)
- startY: y-coordinate for B scan on IR. (see getGrid)
- endX: x-coordinate for B scan on IR. (see getGrid)
- endY: y-coordinate for B scan on IR. (see getGrid)
- numSeg: 2 or 3 segmentation lines for the B scan
- quality: OCT signal quality
- shift: unknown
"""
return self.wholefile["slice-headers"][slicei]
def saveGrid(self, outfn):
"""
Saves the grid coordinates mapping OCT Bscans to the IR SLO image to a text file. The text file
will be a tab-delimited file with 5 columns: The bscan number, x_0, y_0, x_1, y_1 in pixel space
scaled to the resolution of the IR SLO image.
Args:
outfn (str): location of where to output the file
Returns:
None
"""
grid = self.grid
with open(outfn, "w") as fout:
fout.write("bscan\tx_0\ty_0\tx_1\ty_1\n")
ri = 0
for r in grid:
r = [ri] + r
fout.write("%s\n" % "\t".join(map(str, r)))
ri += 1