-
Notifications
You must be signed in to change notification settings - Fork 0
/
ColiImage.m
336 lines (262 loc) · 12.6 KB
/
ColiImage.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
%COLIIMAGE
%
% This class handles a pair of images, one containing CheZ information
% and one containing CheY information. It contains a function to locate
% cells in the image, and a function to locate clusters in each cell,
% along with functions to plot all of these data.
%
% Class must inherit handle so it can be self modifying
classdef ColiImage < handle
properties
cellImage % location of the cell image
clusterImage % location of the cluster image
imageSize % size of the image
cells % cell array of each cell
pixelSize % pixel size of image, rather than camera
% i.e., 16µm camera w/ 100X objective =
% 160nm pixel size
end
methods
% Constructor function. Takes cell and cluster image locations
% and the pizel size
function self = ColiImage( cellImageLocation, clusterImageLocation,...
pixelSize )
% Calculate and store the important values
self.cellImage = cellImageLocation;
self.clusterImage = clusterImageLocation;
self.imageSize = size( imread( self.cellImage ) );
self.pixelSize = pixelSize;
self.getCells( 1000, 5000000, 15000000, 1.3, 0.85, [], 500 );
self.getClusters( 1.5, 600, 250, 750 );
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% getCells - function to find and categorise every cell in the %
% image. %
% %
% Inputs: %
% border - exclusion zone around edge of image to %
% prevent picking up half cells %
% minSize - minimum size of cell to be considered, in %
% square microns %
% maxSize - as above %
% aspectRatio - minimum ratio of major axis length to %
% minor axis length of cell %
% solidity - minimum ratio of cell area to convex area %
% of cell; removes bent/joined cells %
% maxBright - maximum mean brightness of cell %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function getCells( self, border, minSize, maxSize, aspectRatio,...
solidity, maxBright, expansion )
% load the error definition file
load( 'errors.mat' );
% read in the cell image
image = imread( self.cellImage );
% convert constraints from microns to pixels
minSize = minSize / ( self.pixelSize ^ 2 );
maxSize = maxSize / ( self.pixelSize ^ 2 );
border = ceil( border / self.pixelSize );
expansion = ceil( expansion / self.pixelSize );
% get the black level
blackLevel = self.getBlackLevel( 'cell', 16 );
% locate the cell edges using the canny filter. This can only
% tolerate a moderate level of noise, so run a quick gaussian
% over it first
cellEdges = edge(...
imfilter( image, fspecial( 'gaussian', 10, 15 ) ),...
'canny' );
% create a border mask to eliminate any cells touching the
% edge
borderMask = zeros( self.imageSize );
borderMask( border + 1:self.imageSize(1) - border - 1,...
border + 1:self.imageSize(2) - border - 1 ) = 1;
% bwmorph - close any small gaps in the edges of the cells
% imfill - fill in the cell outlines
% imopen - remove any small blobs that are clearly not cells
cellMask = imopen(...
imfill(...
bwmorph( cellEdges & borderMask, 'close' ),...
'holes' ),...
strel( 'disk', 4 ) );
% set up which checks to make on the cells
checks = { 'Image', 'BoundingBox', 'MeanIntensity', 'Area',...
'Solidity', 'MajorAxisLength', 'MinorAxisLength' };
% calculate the metrics of each cell
rp = regionprops( cellMask, image, checks );
% calculate mean intensity of CheZ
rpz = regionprops( cellMask, imread( self.clusterImage ),...
'MeanIntensity' );
% prepare the cellarray for filing
self.cells = cell( 1, numel( rp ) );
% loop through each cell
for i = 1:numel( rp )
% monitor the status of the cell
valid = 0;
% check the cell is not too small, if requested
if check_arg( minSize )
if rp(i).Area < minSize
valid = TOO_SMALL;
end
end
% check the cell is not too big, if requested
if check_arg( maxSize )
if rp(i).Area > maxSize
valid = TOO_BIG;
end
end
% check the cell is not too round, if requested
if check_arg( aspectRatio )
if rp(i).MajorAxisLength / rp(i).MinorAxisLength <...
aspectRatio
valid = TOO_ROUND;
end
end
% check the cell is solid enough, if requested
if check_arg( solidity )
if rp(i).Solidity < solidity
valid = TOO_BENT;
end
end
% check the cell is not too bright, if requested
if check_arg( maxBright )
if rp(i).MeanIntensity > maxBright
valid = TOO_BRIGHT;
end
end
% create a new cell object with all of the above
self.cells{i} = ColiCell( rp(i).Image, rp(i).BoundingBox(1),...
rp(i).BoundingBox(2),...
rp(i).MeanIntensity - blackLevel, valid,...
self.imageSize, rpz(i).MeanIntensity, expansion );
end
end
function validCells = validCells( self, varargin )
validCells = [];
for i = 1:numel( self.cells )
if self.cells{i}.validCell == 0
validCells = [ validCells, self.cells{i} ];
end
end
if nargin == 2
validCells = validCells(varargin{1});
end
end
function getClusters( self, threshold, searchRadius, minSize,...
maxSize )
ballRadius = floor( 1000 / self.pixelSize );
image = double( imread( self.clusterImage ) );
blackLevel = getBlackLevel( self, 'cluster', 8 );
cellBasalLevel = imopen( image, strel( 'ball', ballRadius,...
ballRadius, 0 ) );
statSet = statset( 'MaxIter', 25, 'TolFun', 1e-3,...
'TolX', 1e-3 );
minSize = minSize / self.pixelSize;
maxSize = maxSize / self.pixelSize;
searchRadius = floor( searchRadius / self.pixelSize );
pixels = meshgrid( -searchRadius:searchRadius );
xPix = reshape( pixels, 1, numel( pixels ) );
yPix = reshape( pixels', 1, numel( pixels ) );
validCells = self.validCells();
for i = 1:numel( validCells )
validCells(i).getClusters( image, cellBasalLevel,...
blackLevel, statSet, threshold, searchRadius,...
minSize, maxSize, xPix, yPix )
end
end
function blackLevel = getBlackLevel( self, im, blackLevelCutoff )
if strcmpi( im, 'cell' )
image = imread( self.cellImage );
elseif strcmpi( im, 'cluster' )
image = imread( self.clusterImage );
end
image = image(:);
blackLevel = mean( image(...
image < ( 2^16 / blackLevelCutoff ) ) );
end
function displayOverlay( self )
z = zeros( self.imageSize );
figure(1);
ceI = imread( self.cellImage );
clI = imread( self.clusterImage );
max(ceI(:))
max(clI(:))
im = cat( 3, ceI, clI, z ) ;
max(im(:))
size( im )
imshow( cat( 3, ceI, clI, z ) );
end
function displayCells( self, im )
load( 'errors.mat' );
if strcmpi( im, 'cell' )
image = imread( self.cellImage );
elseif strcmpi( im, 'cluster' )
image = imread( self.clusterImage );
end
cell_mask = zeros( self.imageSize );
size_mask = zeros( self.imageSize );
round_mask = zeros( self.imageSize );
bright_mask = zeros( self.imageSize );
bent_mask = zeros( self.imageSize );
for i = 1:numel( self.cells )
mask = self.cells{i}.expandedFullMask();
if self.cells{i}.validCell == 0
cell_mask = cell_mask | mask;
elseif self.cells{i}.validCell == TOO_BIG || self.cells{i}.validCell == TOO_SMALL
size_mask = size_mask | mask;
elseif self.cells{i}.validCell == TOO_ROUND
round_mask = round_mask | mask;
elseif self.cells{i}.validCell == TOO_BRIGHT
bright_mask = bright_mask | mask;
elseif self.cells{i}.validCell == TOO_BENT
bent_mask = bent_mask | mask;
end
end
o = ones( self.imageSize );
z = zeros( self.imageSize );
red = cat( 3, o, z, z );
green = cat( 3, z, o, z );
blue = cat( 3, z, z, o );
yellow = cat( 3, o, o, z );
orange = cat( 3, o, 0.5 * o, z );
figure(1);
imshow(image, [ min( image(:) ), max( image(:) ) ] );
hold on
valid = imshow( green );
set( valid, 'AlphaData', cell_mask * 192 );
round = imshow( red );
set( round, 'AlphaData', round_mask * 192 );
bent = imshow( blue );
set( bent, 'AlphaData', bent_mask * 192 );
bright = imshow( yellow );
set( bright, 'AlphaData', bright_mask * 128 );
sizes = imshow( orange );
set( sizes, 'AlphaData', size_mask * 128 );
hold off;
end
function cc = clusterCount( self )
validCells = self.validCells();
cc = zeros( 1, numel( validCells ) );
for i = 1:numel( validCells )
cc( i ) = validCells(i).clusterCount();
end
end
function varargout = clusterHistogram( self )
h = histc( self.clusterCount(), 0:4 );
varargout{1} = h;
if nargout == 0
bar( 0:4, h, 'histc' )
end
end
function cf = clusterFraction( self )
validCells = self.validCells();
cf = [];
for i = 1:numel( validCells )
cf = [ cf, validCells(i).clusterFraction ];
end
end
end
end
function ok = check_arg( arg )
ok = isa( arg, 'double' ) && numel( arg ) == 1 && arg > 0;
end