-
Notifications
You must be signed in to change notification settings - Fork 5
/
Copy pathOCT.m
505 lines (445 loc) · 16.9 KB
/
OCT.m
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
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
classdef OCT < handle
% OCT
%
% Description:
% Represents features and analyses of a single OCT b-scan
%
% Constructor:
% x = OCT(imageID, imagePath)
%
% Inputs:
% imageID Image number (also image title)
% imagePath Folder where image and data is stored, optional
% If not provided, opens up UI to choose a folder
%
% Methods:
% obj.update()
% obj.crop()
% obj.doAnalysis()
% Visualization methods:
% obj.show()
% obj.plotRatio()
% obj.plotSizes()
% obj.showSegmentation()
%
% History:
% 7Aug2018 - SSP - working version compiled from previous functions
% 11Aug2018 - SSP - added crop function
% 14Aug2018 - SSP - added reload option
% 2Jan2019 - SSP - changed edges to control points
% 3Jan2019 - SSP - added scale property
% 4Jan2019 - SSP - crop now happens before rotation and scaling
% ---------------------------------------------------------------------
% Identifiers
properties (SetAccess = private)
imagePath
imageID
refID
end
% Segmentation
properties (SetAccess = public)
RPE = [];
ILM = [];
Choroid = [];
ChoroidParams = [];
ControlPoints = [];
CropValues = [];
Shift = [];
Theta = [];
Scale = [];
end
% Analysis
properties (SetAccess = private)
choroidSize
retinaSize
choroidRatio
end
% Lazy loading of images and analysis
properties (Dependent = true, Hidden = true)
octImage
rawImage
imageName
end
% Transient copy of the full original image to support lazy loading
properties (Hidden = true, SetAccess = private)
originalImage
end
% Features which may have been extracted from the image
properties (Constant = true, Hidden = true)
FEATURES = {'rpe', 'ilm', 'choroid', 'parabola', 'controlpoints',...
'crop', 'shift', 'theta', 'scale'};
end
methods
function obj = OCT(imageID, imagePath)
% OCT Constructor
assert(isnumeric(imageID), 'Image ID must be numeric');
obj.imageID = imageID;
if nargin < 2 || ~isfolder(imagePath)
obj.imagePath = uigetdir();
else
obj.imagePath = imagePath;
end
fprintf('image path: %s\n', obj.imagePath);
obj.pull();
if ~isempty(obj.Choroid) && ~isempty(obj.RPE)
obj.doAnalysis();
end
end
function addRef(obj, refID)
obj.refID = refID;
end
function update(obj)
obj.pull();
end
function im = getSemiProcessedImage(obj, varargin)
ip = inputParser();
addParameter(ip, 'Crop', false, @islogical);
addParameter(ip, 'Scale', false, @islogical);
addParameter(ip, 'Rotate', false, @islogical);
parse(ip, varargin{:});
im = obj.fetchImage(true);
if ip.Results.Crop
im = obj.doCrop(im);
end
if ip.Results.Scale
im = obj.doScale(im);
end
if ip.Results.Rotate
im = obj.doRotate(im);
end
end
function cropValues = crop(obj, saveValues)
if nargin < 2
saveValues = false;
else
assert(islogical(saveValues), 'Save values should be t/f');
end
fh = figure();
[~, cropValues] = imcrop(obj.octImage);
% Return if no crop specified.
if isempty(cropValues)
obj.CropValues = [];
obj.saveJSON();
return;
end
cropValues = [ceil(cropValues(1:2)), floor(cropValues(3:4))];
if saveValues
obj.CropValues = cropValues;
obj.saveJSON();
end
delete(fh);
end
end
% Analysis methods
methods
function doAnalysis(obj)
% Determine the full x-axis
xpts = obj.getXPts();
% Evaluate choroid parabola at the x-axis points
CHOROID = parabola(xpts, obj.ChoroidParams);
% Interpolate the RPE and ILM boundaries to the new x-axis
iRPE = interp1(obj.RPE(:, 1), obj.RPE(:, 2), xpts);
iILM = interp1(obj.ILM(:, 1), obj.ILM(:, 2), xpts);
% Choroid to RPE-Choroid boundary
obj.choroidSize = abs(CHOROID - iRPE);
% RPE-Choroid boundary to ILM
obj.retinaSize = abs(iRPE - iILM);
% Ratio of choroid size to retina size
obj.choroidRatio = obj.choroidSize./obj.retinaSize;
% Shift, if necessary
% if ~isempty(obj.Shift)
% obj.choroidRatio = circshift(obj.choroidRatio, -obj.Shift);
% if obj.Shift > 0
% obj.choroidRatio(end-abs(obj.Shift):end) = NaN;
% elseif obj.Shift < 0
% obj.choroidRatio(1:obj.Shift) = NaN;
% end
% end
end
function xpts = getXPts(obj, doShift)
% GETXPTS Get image x-axis, shifted or unshifted
if nargin < 2
doShift = false;
end
if ~isempty(obj.RPE) && ~isempty(obj.ILM)
x_min = min([obj.RPE(:, 1); obj.ILM(:, 1)]);
x_max = max([obj.RPE(:, 1); obj.ILM(:, 1)]);
xpts = x_min:x_max;
else
xpts = [];
end
if doShift && ~isempty(obj.Shift)
xpts = xpts + obj.Shift;
fprintf('Applying shift of %.3g pixels\n', obj.Shift);
end
end
function saveJSON(obj)
% SAVEJSON Store features as a single .json file
warning('off', 'MATLAB:structOnObject');
S = struct(obj);
S = rmfield(S, {'choroidSize', 'retinaSize', 'choroidRatio'});
S = rmfield(S, {'imageName', 'FEATURES'});
S = rmfield(S, {'octImage', 'originalImage', 'rawImage'});
savejson('', S, [obj.imagePath, filesep, obj.imageName, '.json']);
warning('on', 'MATLAB:structOnObject');
fprintf('Saved %s.json\n', obj.imageName);
end
function igor = igor(obj, smoothFac)
if nargin < 2
smoothFac = 12;
end
if isempty(obj.choroidRatio)
obj.doAnalysis();
end
igor = [obj.getXPts(true)', smooth(obj.choroidRatio, smoothFac)];
openvar('igor');
end
end
% Visualization methods
methods
function show(obj, ax)
% SHOW Plot the OCT image (with rotation and crop)
if nargin < 2
ax = axes('Parent', figure());
end
imagesc(ax, obj.octImage); colormap(gray);
title(ax, sprintf('Image #%u', obj.imageID));
axis image off
tightfig(gcf);
end
function ax = plotRatio(obj, varargin)
% PLOTRATIO Plot the choroid:retina ratio
ip = inputParser();
ip.CaseSensitive = false;
ip.KeepUnmatched = true;
addParameter(ip, 'Relative', false, @islogical);
addParameter(ip, 'Smooth', 10, @isnumeric);
addParameter(ip, 'ShowData', false, @islogical);
addParameter(ip, 'Color', rgb('light blue'),...
@(x) isvector(x) || ischar(x));
addParameter(ip, 'LineWidth', 1, @isnumeric);
addParameter(ip, 'Axes', [], @ishandle);
parse(ip, varargin{:});
smoothFac = ip.Results.Smooth;
if isempty(obj.choroidRatio)
obj.doAnalysis();
end
xpts = obj.getXPts(true);
if isempty(ip.Results.Axes)
fh = figure('Name', [obj.imageName, ' Choroid'],...
'Renderer', 'painters');
ax = axes('Parent', fh);
h = plot(ax, [min(xpts), max(xpts)], [1, 1], '--',...
'Color', [0.3, 0.3, 0.3], 'LineWidth', 0.75,...
'Tag', 'NormLine');
set(get(get(h, 'Annotation'), 'LegendInformation'),...
'IconDisplayStyle', 'off');
set(ax, 'Box', 'off'); grid(ax, 'on');
ylabel(ax, 'choroid:retina ratio');
xlabel(ax, 'x-axis (pixels)');
figPos(fh, 0.7, 0.7);
else
ax = ip.Results.Axes;
end
hold(ax, 'on');
if ip.Results.Relative && ~isempty(obj.refID)
refOCT = OCT(obj.refID, obj.imagePath);
refRatio = refOCT.choroidRatio;
thisRatio = obj.choroidRatio;
if numel(refRatio) > numel(xpts)
refRatio(numel(xpts)+1:end) = [];
elseif numel(xpts) > numel(refRatio)
thisRatio(numel(refRatio)+1:end) = [];
% xpts = refOCT.getXPts();
xpts(numel(refRatio)+1:end) = [];
end
plot(ax, xpts, smooth(thisRatio - refRatio, smoothFac),...
'Color', ip.Results.Color, 'LineWidth', 1.5);
set(findobj(ax, 'Tag', 'NormLine'), 'YData', [0, 0]);
ylim(ax, [-0.5, 0.5]);
else
plot(ax, xpts, smooth(obj.choroidRatio, smoothFac),...
'Color', ip.Results.Color,...
'LineWidth', ip.Results.LineWidth,...
'Tag', obj.imageName, 'DisplayName', obj.imageName);
ylim(ax, [0, 2]);
end
% Raw data
if ip.Results.ShowData && ~ip.Results.Relative
h = plot(ax, xpts, obj.choroidRatio, '.',...
'Color', hex2rgb('334de6'), 'MarkerSize', 4,...
'Tag', 'Data');
set(get(get(h, 'Annotation'), 'LegendInformation'),...
'IconDisplayStyle', 'off');
end
xlim(ax, [0, max(xpts)]);
end
function plotSizes(obj, ax)
% PLOTRATIO Plot raw sizes of choroid and retina
if isempty(obj.choroidRatio)
obj.doAnalysis();
end
xpts = obj.getXPts(true);
if nargin < 2
ax = axes('Parent', figure('Renderer', 'painters'),...
'Box', 'off', 'XLim', [0, max(xpts)]);
hold(ax, 'on');
grid(ax, 'on');
xlabel(ax, 'x-axis (pixels)');
ylabel(ax, 'Width (pixels)');
title(ax, [obj.imageName ' - Choroid and Retina Sizes']);
end
p1 = plot(ax, xpts, smooth(obj.choroidSize, 10),...
'b', 'LineWidth', 1, 'Tag', [obj.imageName, '_Choroid']);
p2 = plot(ax, xpts, smooth(obj.retinaSize, 10),...
'r', 'LineWidth', 1, 'Tag', [obj.imageName, '_RPE']);
if nargin < 2
legend([p1, p2], {'Choroid', 'Retina'});
end
end
function showSegmentation(obj)
% SHOWSEGMENTATION Show image with overlaid segmentation
if isempty(obj.RPE)
return
end
obj.show(); hold on;
plot(obj.RPE(:, 1), obj.RPE(:, 2), 'b');
plot(obj.ILM(:, 1), obj.ILM(:, 2), 'b');
if ~isempty(obj.Choroid)
plot(obj.Choroid(:, 1), obj.Choroid(:, 2), 'r');
end
end
end
% Dependent set/get methods
methods
function imageName = get.imageName(obj)
imageName = ['im', num2str(obj.imageID)];
end
function octImage = get.octImage(obj)
if isempty(obj.originalImage)
octImage = obj.fetchImage();
else
octImage = obj.originalImage;
end
% Crop/rotate/transform, if necessary
octImage = obj.doCrop(octImage);
octImage = obj.doScale(octImage);
octImage = obj.doRotate(octImage);
end
function rawImage = get.rawImage(obj)
rawImage = obj.fetchImage(true);
end
end
% Private image processing methods
methods (Access = private)
function im = fetchImage(obj, getRaw)
% FETCHIMAGE Read image from file
if nargin < 2
getRaw = false;
end
if getRaw
str = [obj.imagePath, 'raw', filesep, num2str(obj.imageID), '.png'];
else
str = [obj.imagePath, num2str(obj.imageID), '.png'];
end
im = imread(str);
% Convert to grayscale if necessary
if numel(size(im)) == 3
im = rgb2gray(im);
end
end
function im = doCrop(obj, im)
% DOCROP Crop the image, if necessary
if ~isempty(obj.CropValues)
im = imcrop(im, obj.CropValues);
fprintf('Cropped by %u, %u\n',...
abs(round([size(im, 1) - obj.CropValues(3),...
size(im, 2) - obj.CropValues(4)])));
end
end
function im = doScale(obj, im)
% DOSCALE Scale the image, if necessary
if ~isempty(obj.Scale)
im = imresize(im, obj.Scale);
fprintf('Scaled by %.4g%%\n', 100 * obj.Scale);
end
end
function im = doRotate(obj, im)
% DOROTATE Rotate the image, if necessary
if ~isempty(obj.Theta)
im = imrotate(im, obj.Theta);
fprintf('Applied rotation of %.2f degrees\n', obj.Theta);
end
end
end
% Private IO methods
methods (Access = private)
function pull(obj)
% RELOAD Reload image and data from .json files
obj.originalImage = obj.fetchImage();
hasJSON = obj.loadJSON();
if ~hasJSON
obj.saveJSON();
% disp('Loading .txt files instead');
% obj.loadTXT();
end
if ~isempty(obj.RPE) && ~isempty(obj.Choroid)
obj.doAnalysis();
end
end
function tf = loadJSON(obj)
% LOADJSON Load a .json file
jsonPath = [obj.imagePath, filesep, obj.imageName, '.json'];
if ~exist(jsonPath, 'file')
str = strsplit(jsonPath, filesep);
fprintf('Did not find file %s\n', str{end});
tf = false;
else
S = loadjson(jsonPath);
tf = true;
if isfield(S, 'imageID') && ~isempty(S.imageID)
assert(obj.imageID == S.imageID,...
'Image ID numbers do not match!');
end
obj.refID = S.refID;
obj.RPE = S.RPE; obj.ILM = S.ILM;
obj.Shift = S.Shift; obj.Theta = S.Theta; obj.Scale = S.Scale;
obj.CropValues = S.CropValues;
obj.ControlPoints = S.ControlPoints;
obj.ChoroidParams = S.ChoroidParams;
obj.Choroid = S.Choroid;
end
end
function loadTXT(obj)
% LOAD Load parameters from file, if they exist
obj.Choroid = obj.fetch(obj.getPath('choroid'));
obj.ChoroidParams = obj.fetch(obj.getPath('parabola'));
obj.RPE = obj.fetch(obj.getPath('rpe'));
obj.ILM = obj.fetch(obj.getPath('ilm'));
obj.ControlPoints = obj.fetch(obj.getPath('controlpoints'));
obj.CropValues = obj.fetch(obj.getPath('crop'));
obj.Shift = obj.fetch(obj.getPath('shift'));
obj.Theta = obj.fetch(obj.getPath('theta'));
obj.Scale = obj.fetch(obj.getPath('scale'));
end
function str = getPath(obj, x)
% GETPATH Returns path of saved feature .txt file
x = lower(x);
assert(ismember(x, obj.FEATURES),...
'Path value not in features list!');
str = [obj.imagePath, filesep, obj.imageName, '_', x, '.txt'];
end
end
methods (Static)
function x = fetch(filePath)
% FETCH Reads a .txt file if existing, returns message if not
if ~exist(filePath, 'file')
str = strsplit(filePath, filesep);
fprintf('Did not find file %s\n', str{end});
x = [];
else
x = dlmread(filePath);
end
end
end
end