-
Notifications
You must be signed in to change notification settings - Fork 0
/
videoFrame.m
546 lines (470 loc) · 25.4 KB
/
videoFrame.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
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
classdef videoFrame
% videoFrame V1.61
% removed unnecessary code from detectBrightfieldDropletsFULL and
% adjusted radius correction
% Comes with videoMaker to analyze droplets moving in image sequence
properties
% might not be worth it... stats % regionprops of the image after image processing
analysis % results based on shape request
time % currentTime in video of frame
end
properties (Dependent)
Positions % Locations of objects in image
Radii % radius dimension of the objects (in px)
end
properties (SetAccess = protected)
type % Image type (fluorescence or brightfield)
rect % Image rectangle
Image % original image given
rescaled % true if image was rescaled for analysis (by 2x)
end
methods
function obj = videoFrame(Im, rect, time, varargin) % finds approximate region properties in image for further analysis
%---------------------------------------------------------------
% Step 1 Initialize immutable properties
%---------------------------------------------------------------
if nargin == 0
return
end
if nargin > 1 && ~isempty(rect)
obj.rect = rect;
if any(cellfun(@(str)strcmp(str,'rescale'),varargin))
obj.rescaled = true;
Im = imcrop(imresize(Im,2),[2*rect(1) 2*rect(2)+rect(4)/2 rect(3)*2 rect(4)]);
else
Im = imcrop(Im,rect);
end
end
if nargin > 2
obj.time = time;
end
if size(Im,3) == 3
obj.Image = rgb2gray(Im);
else
obj.Image = Im;
end
if any(cellfun(@(str)strcmp(str,'brightfield'),varargin))
obj.type = 'brightfield';
else
obj.type = 'fluorescence';
end
%---------------------------------------------------------------
% Step 2 Aproximate properties of image
% (skipping this for now since it is slow)
%---------------------------------------------------------------
% obj.stats = approximateFrameProperties(obj,'Area','WeightedCentroid','Eccentricity','MajorAxisLength','MinorAxisLength');
%---------------------------------------------------------------
% Step 3 Run analysis
%---------------------------------------------------------------
implementedAnalysis = {'circle','full','combined'}; % potentially ellipse, droplet, rectangle
performAnalysisFor = find(ismember(implementedAnalysis,lower(varargin)));
for n = performAnalysisFor
forType = implementedAnalysis{n};
obj.analysis.(forType) = runAnalysis(obj,forType);
end
end
function obj = clean(obj)
obj.Image = [];
end
function h = plot(obj,ax,varargin)
% if isempty(obj.analysis)
% h = plot(obj.stats.WeightedCentroid(1), obj.stats.WeightedCentroid(2),'ro');
% return
% end
if nargin < 2 || isempty(ax)
ax = axes(figure);
imshow(obj.Image)
end
if isempty(obj.analysis)
return
elseif any(cellfun(@(field)isa(field,'numeric'),varargin))
analysisToPlot = varargin{find(cellfun(@(field)isa(field,'numeric'),varargin),1)};
else
analysisToPlot = length(fieldnames(obj.analysis));
end
% CHANGE if varargin specifies anything, else plot first item
hold(ax,'on')
analysisFor = fieldnames(obj.analysis);
switch analysisFor{analysisToPlot}
case 'circle'
if any(cellfun(@(str)strcmp(str,'noOutline'),varargin))
h = plot(obj.analysis.circle.centers(:,1),obj.analysis.circle.centers(:,2),'r+');
else
h = viscircles(ax,obj.analysis.circle.centers,obj.analysis.circle.radii);
end
case 'full'
if any(cellfun(@(str)strcmp(str,'noOutline'),varargin))
h = plot(obj.analysis.full.centers(:,1),obj.analysis.full.centers(:,2),'r+');
else
h = viscircles(ax,obj.analysis.full.centers,obj.analysis.full.radii);
end
case 'combined'
otherwise
return
end
hold(ax,'off')
end
function analysis = runAnalysis(obj,type)
switch type
case 'circle'
if strcmp(obj.type,'brightfield')
analysis = detectBrightfieldDroplets(obj);
else
analysis = detectDroplets(obj); % previously detectDroplets(obj)
end
case 'full'
if ~strcmp(obj.type,'brightfield')
analysis = detectDropletsFULL(obj);
else
analysis = detectBrightfieldDropletsFULL(obj);
end
case 'combined'
analysis = combineAnalysis(obj);
otherwise
analysis = [];
return
end
end
function report = reportData(obj,type)
analysisFor = fieldnames(obj.analysis);
if nargin < 2
type = analysisFor{end};
elseif ~ismember(analysisFor, type)
report = []; return
end
if obj.rescaled
report = [(obj.analysis.(type).centers(:,1) + obj.rect(1)*0)/2, (obj.analysis.(type).centers(:,2) + obj.rect(4)/2)/2, obj.analysis.(type).radii/2];
else
report = [obj.analysis.(type).centers, obj.analysis.(type).radii];
end
report(:,4) = obj.time;
report = double(report);
end
function positions = positions(obj,type)
analysisFor = fieldnames(obj.analysis);
if nargin < 2
type = analysisFor{1};
elseif ~ismember(analysisFor, type)
positions = []; return
end
if obj.rescaled
positions = [(obj.analysis.(type).centers(:,1) + obj.rect(1)*0)/2, (obj.analysis.(type).centers(:,2) + obj.rect(4)/2)/2];
else
positions = obj.analysis.(type).centers;
end
positions(:,3) = obj.time;
positions = double(positions);
end
end
methods (Hidden = true)
function analysis = detectDroplets(obj)
% setJump = 3;
warning('off')
% TRY AGAIN ISOLATING DROPLETS AND FITTING ONLY RADII THAT FIT
% INTO SIZE OF REGION
%---------------------------------------------------------------
% Step 1 Isolate droplet regions
%---------------------------------------------------------------
% A = obj.Image > 30; % thresholding original image
% B = bwareaopen(A,15); % excluding small regions
% C = regionprops(A,'BoundingBox'); % isolating droplet regions
%
% image = cell(uint8(zeros(size(obj.Image))));
% for n = 1:numel(C)
% image{n} = imcrop(obj.Image,C(n).BoundingBox);
% end
%images{1} = obj.Image; images{2} = imresize(obj.Image,2);
rescale = table([0.93; 0.93; 0.93; 0.93],[30 60;10 30;4 10;1 4],'VariableNames',{'Scale','inRange'});
if obj.rescaled
rescale = table(rescale{:,1},rescale{:,2}*2,'VariableNames',{'Scale','inRange'});
end
%---------------------------------------------------------------
% Step 1 Set radius relations (Deprecated: Approximate radii from stats)
%---------------------------------------------------------------
% circularObjects = [obj.stats.Eccentricity] <= 0.75;
% radiusRange = [min([1 floor([obj.stats(circularObjects).MajorAxisLength]/4 + [obj.stats(circularObjects).MinorAxisLength]/4)]) ...
% max(ceil([obj.stats(circularObjects).MajorAxisLength]/4 + [obj.stats(circularObjects).MinorAxisLength]/4))];
% range = radiusRange(1)*[1 setJump];
%---------------------------------------------------------------
% Step 2 Find circles in ranges
%---------------------------------------------------------------
centers = nan(100,2); radii = nan(100,1); metric = nan(100,1);
% while max(range) < radiusRange(2)
for n = 1:height(rescale)
sensitivity = rescale{n,1}; inRange = rescale{n,2};
[C, R, M] = imfindcircles(obj.Image,inRange,'Sensitivity',sensitivity);%0.92);
% rescale center and radius
%C = C./scale; R = R./scale;
currentIndex = max([find(radii(~isnan(radii)),1,'last'),0])+1;
for m = 1:numel(R)
overlap = (C(m,1) - centers(:,1)).^2 + (C(m,2) - centers(:,2)).^2 < radii.^2;
if any(overlap)
% replace better fits
centers(overlap & metric < M(m),:) = repmat(C(m,:),sum(overlap & metric < M(m)),1);
radii(overlap & metric < M(m),:) = R(m);
metric(overlap & metric < M(m),:) = M(m);
else
% add new circle
centers(currentIndex,:) = C(m,:);
radii(currentIndex) = R(m);
metric(currentIndex) = M(m);
currentIndex = currentIndex + 1;
end
end
% range = range(2)*[1 setJump];
end
warning('on')
% Clean up centers and radii (remove identical doubles)
results = unique([centers(metric > 0.1,:), radii(metric > 0.1)],'rows');
analysis = struct('for','circles',...
'centers',results(~isnan(results(:,3)),1:2),...
'radii',results(~isnan(results(:,3)),3),...
'quality',struct('absolute',sum(metric(~isnan(metric))),'relative',median(metric(~isnan(metric)))));
end
function analysis = detectFluorescentDroplets(obj)
% threshold image
A = obj.Image;
B = adaptthresh(A,0.3); B(B>0.75) = 0.75; % B = A; B(A<25) = 0; B(A>150) = 255; % potentially do without
C = imbinarize(A,B); % imextendedmax(B,5);
% find coarse estimate of droplets and filter bad shapes
% (conjoined droplets)
imageStats = regionprops(C,A,'BoundingBox','Extent','Image');
badEstimates = [imageStats.Extent] < 0.75; blank = zeros(size(A)); % start blank guess for droplet locations
% Note: perfect circles should have extend around 0.78, larger
% values usually occur for very small droplets
% filter out conjoined droplets and find a better estimate
% based on distance maxima
localRadii = bwdist(~C);
localMaxima = logical(imregionalmax(localRadii).*bwmorph(imregionalmax(medfilt2(localRadii)),'thicken',2)); % finds maxima and roughly filters local maxima between two droplets
for i = find(~badEstimates) % remove regions of good estimates from image
n = imageStats(i);
Ys = ceil(n.BoundingBox(1)):floor(n.BoundingBox(1)+n.BoundingBox(3));
Xs = ceil(n.BoundingBox(2)):floor(n.BoundingBox(2)+n.BoundingBox(4));
localMaxima(Xs,Ys) = zeros(size(n.Image));
blank(Xs,Ys) = n.Image; % fill up blank guess with good estimates for droplets
end
% find located maxima and draw into blank
localStats = regionprops(localMaxima,'Centroid','PixelIdxList');
centers = vertcat(localStats.Centroid);
radii = arrayfun(@(x)mean(localRadii(x.PixelIdxList))-3,localStats);
[colums, rows] = meshgrid(1:size(A,2),1:size(A,1));
inCirc = @(c,r)(colums-c(1)).^2 + (rows-c(2)).^2 < r.^2;
for n = 1:length(radii)
blank(inCirc(centers(n,:),radii(n))) = 1;
end
% extend the size of droplets to overlap with original image
finalGuess = bwmorph(blank,'thicken',25);
F = double(A).*finalGuess; F = F > 255/exp(1); % overlap and threshold guess with original image
finalStats = regionprops(F, A,'WeightedCentroid','Eccentricity','EquivDiameter');
finalStats([finalStats.Eccentricity] >= 0.75) = [];
analysis = struct('for','circles',...
'centers',vertcat(finalStats.WeightedCentroid),...
'radii',vertcat(finalStats.EquivDiameter)/2);
end
function analysis = detectBrightfieldDroplets(obj)
% thresh = double(max(obj.Image(:)) - min(obj.Image(:)))/255 +0.1;
A = imbinarize(imcomplement(obj.Image)); %,thresh); %0.7); % binarize the inverse brightfield image
B = imfill(A,'holes'); % fill any fully enclosed white regions
C = B; C(A) = 0; % remove outlines of said regions
stats = regionprops(C, obj.Image, 'Area','WeightedCentroid','Eccentricity','MajorAxisLength','MinorAxisLength');
stats([stats.Eccentricity] >= 0.75) = [];
D = xor(bwareaopen(A,1), bwareaopen(A,min([10*max([stats.Area]),300000]))); % isolate dark droplet rim in original image
d = bwdist(imcomplement(D)); d = median(d(d>0)); % d is the half-width of the dark rim of the droplet
analysis = struct('for','circles',...
'centers',vertcat(stats.WeightedCentroid),...
'radii',mean([vertcat(stats.MajorAxisLength),vertcat(stats.MinorAxisLength)]/2,2) +1.75*d);
end
function analysis = detectDropletsFULL(obj)
warning('off')
% imfindcircle parameters
inRange = [20 50]; resizeFactor = [1 2 5 10];
%inRange = [10 25]; resizeFactor = [1 2 5 0.5];
% noarmalize and rescale images
objImage = uint8(double(obj.Image)*255/max(double(obj.Image(:))));
images = arrayfun(@(n)imresize(objImage,n),resizeFactor,'Uni',0);
%---------------------------------------------------------------
% Step 1 Find circles using imfindcircles
%---------------------------------------------------------------
C(1:numel(resizeFactor)+1,1) = {zeros(0,2)}; R = cell(numel(resizeFactor),1); M = R;
n = 1;
for m = resizeFactor
[C{n}, R{n}, M{n}] = imfindcircles(images{n},inRange,'Sensitivity',0.85);
% rescale center and radius
C{n} = C{n}/m; R{n} = R{n}/m;
n = n+1;
end
% clear circles with low metric
centers = cell2mat(C);
radii = cell2mat(R);
metric = cell2mat(M);
centers(metric < 0.35,:) = [];
radii(metric < 0.35) = [];
metric(metric < 0.35) = [];
% clear overlapping circles
[m, n] = meshgrid(1:numel(radii), 1:numel(radii));
X = centers(:,1); Y = centers(:,2);
overlap = any( ((X(m) - X(n)).^2 + (Y(m) - Y(n)).^2 <= radii(m).^2) & metric(m) > metric(n) ,2);
centers(overlap,:) = [];
radii(overlap) = [];
% metric(overlap) = [];
%---------------------------------------------------------------
% Step 2 Remove found circles and try again
%---------------------------------------------------------------
% remove cirlces already found
Im = images{1};
[colums, rows] = meshgrid(1:size(Im,2),1:size(Im,1));
inCirc = @(c,r)(colums-c(1)).^2 + (rows-c(2)).^2 <= (r+2.5).^2;
for i = 1:length(radii)
Im(inCirc(centers(i,:),radii(i))) = 0;
end
A = double(Im); A = A/median(A(A>0));
B = A; B(Im<25) = 0; % B(I>200) = 255; % potentially do without
D = imopen(logical(B),strel('disk',5)); % imextendedmax(B,5);
% find coarse estimate of droplets and filter bad shapes
% (conjoined droplets)
imageStats = regionprops(D,A,'BoundingBox','Extent','Image');
badEstimates = [imageStats.Extent] < 0.75; blank = zeros(size(A)); % start blank guess for droplet locations
% Note: perfect circles should have extend around 0.78, larger
% values usually occur for very small droplets
% filter out conjoined droplets and find a better estimate
% based on distance maxima
localRadii = bwdist(~D);
localMaxima = logical(imregionalmax(localRadii).*bwmorph(imregionalmax(medfilt2(localRadii)),'thicken',2)); % finds maxima and roughly filters local maxima between two droplets
for i = find(~badEstimates) % remove regions of good estimates from image
n = imageStats(i);
Ys = ceil(n.BoundingBox(1)):floor(n.BoundingBox(1)+n.BoundingBox(3));
Xs = ceil(n.BoundingBox(2)):floor(n.BoundingBox(2)+n.BoundingBox(4));
localMaxima(Xs,Ys) = zeros(size(n.Image));
blank(Xs,Ys) = n.Image; % fill up blank guess with good estimates for droplets
end
% find located maxima and draw into blank
localStats = regionprops(localMaxima,'Centroid','PixelIdxList');
centroid = vertcat(localStats.Centroid);
radius = arrayfun(@(x)mean(localRadii(x.PixelIdxList))-3,localStats);
[colums, rows] = meshgrid(1:size(A,2),1:size(A,1));
inCirc = @(c,r)(colums-c(1)).^2 + (rows-c(2)).^2 < r.^2;
for n = 1:length(radius)
blank(inCirc(centroid(n,:),radius(n))) = 1;
end
% get sizes of shapes not connected to inlet/outlet
finalGuess = bwmorph(blank,'thicken',2);
F = double(A).*finalGuess; F = F/max(F(:)) > exp(-1.5);
finalStats = regionprops(logical(imclearborder(F)), obj.Image,'WeightedCentroid','Eccentricity','EquivDiameter','EulerNumber','Image');
finalStats([finalStats.Eccentricity] >= 0.75) = [];
finalStats([finalStats.EulerNumber] < 1) = [];
finalStats([finalStats.EquivDiameter]/2 < 2) = [];
%---------------------------------------------------------------
% Step 3 Consolidate and report results
%---------------------------------------------------------------
warning('on')
centers = [centers; vertcat(finalStats.WeightedCentroid)];
radii = [radii; vertcat(finalStats.EquivDiameter)/2];
analysis = struct('for','circles',...
'centers',centers,...
'radii',radii,...
'quality',struct('absolute',sum(metric(~isnan(metric))),'relative',median(metric(~isnan(metric)))));
end
function analysis = detectBrightfieldDropletsFULL(obj)
warning('off')
% imfindcircle parameters
inRange = [20 50]; resizeFactor = [1 2 5 10];
% noarmalize and rescale images
objImage = uint8(double(obj.Image)*255/max(double(obj.Image(:))));
images = arrayfun(@(n)imresize(objImage,n),resizeFactor,'Uni',0);
%---------------------------------------------------------------
% Step 1 Find circles using imfindcircles
%---------------------------------------------------------------
C(1:numel(resizeFactor)+1,1) = {zeros(0,2)}; R = cell(numel(resizeFactor),1); M = R;
n = 1;
for m = resizeFactor
[C{n}, R{n}, M{n}] = imfindcircles(images{n},inRange,'Sensitivity',0.95,'ObjectPolarity','dark');
% rescale center and radius
C{n} = C{n}/m; R{n} = R{n}/m;
n = n+1;
end
% clear circles with low metric
centers = cell2mat(C);
radii = cell2mat(R);
metric = cell2mat(M);
centers(metric < 0.35,:) = [];
radii(metric < 0.35) = [];
metric(metric < 0.35) = [];
% clear overlapping circles
[m, n] = meshgrid(1:numel(radii), 1:numel(radii));
X = centers(:,1); Y = centers(:,2);
overlap = any( ((X(m) - X(n)).^2 + (Y(m) - Y(n)).^2 <= radii(m).^2) & metric(m) > metric(n) ,2);
centers(overlap,:) = [];
radii(overlap) = [];
metric(overlap) = [];
radii = 0.9*radii; % adjust radii
analysis = struct('for','circles',...
'centers',centers,...
'radii',radii,...
'quality',struct('absolute',sum(metric(~isnan(metric))),'relative',median(metric(~isnan(metric)))));
end
function combined = combineAnalysis(obj)
% currently combineAnalysis only combines full and circle
% analysis to avoid evaluating circle quality
if isempty(obj.analysis) || ~all(ismember({'full','circle'},fieldnames(obj.analysis)))
return
end
if strcmp(obj.type,'fluorescence')
full = obj.analysis.full;
circle = obj.analysis.circle;
circle.radii(circle.radii > 2) = [];
circle.centers(circle.radii > 2,:) = [];
if isempty(circle.radii)
combined = full;
elseif isempty(full.radii)
combined = circle;
else
[X,Y] = meshgrid(1:numel(circle.radii),1:numel(full.radii));
overlap = any(reshape(sqrt((circle.centers(X,1) - full.centers(Y,1)).^2 ...
+ (circle.centers(X,2) - full.centers(Y,2)).^2),size(X)) < circle.radii(X));
combined.centers = vertcat(full.centers, circle.centers(~overlap,:));
combined.radii = vertcat(full.radii, circle.radii(~overlap));
combined.for = full.for;
end
else
return
end
end
function stats = approximateFrameProperties(obj,varargin) % WORK IN PROGRESS (TOO SLOW)
A = obj.Image;
A(A < 30) = 0; % thresholds low intensity parts of image
% B = imbinarize(A,0.6);
% C = imcomplement(B);
% tic; D = C; for n = 1:9; D = imclose(D,strel('line',7,20*n)); end; toc
% F = bwmorph(B,'hbreak',inf);
B = imtophat(imgaussfilt(A),strel('disk',15)); % smooths image and removes non-uniform brightness
C = imopen(localcontrast(B,1,1),strel('disk',1)); % increases contrast locally and smooth single pixel objects
D = imbinarize(C,0.6); % creates black-and-white image
E = bwmorph(D,'hbreak',inf); % removes small connections between regions
F = imfill(E,'holes'); % fills in holes in regions (from local contrast algorithm)
stats = regionprops(F,A,varargin{:});
end
end
methods (Static, Hidden = true)
function [centers, radii, metric] = detectCircles(Im,radius,varargin)
% finds circles in Im
% validate input
validateattributes(Im, {'numeric'},{'3d'})
validateattributes(radius, {'numeric'},{'nondecreasing','2d'})
if length(radius) == 1
radius = [radius*0.9 radius*1.1];
else
validateattributes(radius, {'numeric'},{'size',[1,2]})
end
% check if Image is RGB and needs to be converted still
if size(Im,3) == 3
Im = rgb2gray(Im);
end
% check input parameters
if strcmp('dark',varargin); polarity = 'dark'; else; polarity = 'bright'; end
if strcmp('TwoStage',varargin); method = 'TwoStage'; else; method = 'PhaseCode'; end
sensitivity = [varargin{cellfun(@(y)isa(y,'double'),varargin)}, 0.85];
% find circles
[centers, radii, metric] = imfindcircles(Im,radius,'ObjectPolarity',polarity,'Sensitivity',sensitivity(1), 'Method', method);
end
end
end