trackOptiMakeCell : generates the CellA field indexed by the region number
which contains information about each cell in each region in each frame.
It goes through the dirname/*err.mat files and computes the characteristics of the
cells in each frame and puts them in the CellA structure that is indexed
by the region number, not the cell ID. This code figures out the pole
age and also computes statistics on the fluorescence channels as
well as fitting the locus positions. These last two features are
controlled by parameters set in the loadConstants.m or
loadConstantsMine.m files.
data.CellA{1}.
mask : logical (1 and 0) cell mask, unoriented
r_offset : location of cell in full image from top left corner
BB : coordinates of cell bounding box, containing the
cell and the pad
edgeFlag : if the cell is on the edge of the image (bad)
phase : the cropped phase image of the cell
coord : coordinates, area, and orientation, see below
length : (1) length and (2) width of the cell
pole : orientation of the cell pole, see below
fluor1 : the 1st fluor channel image
fuor1mm : min and max of fluor1
fl1 : statistics of fluor1, e.g. background level
fluor2 : the 2nd fluor channel image
fluor2mm : min and max of fluor2
fl2 : statistics of fluor2
cell_dist : distance to the edge of the colony
gray : average phase gray value in cell region
locus1 : If focus fitting was run, data on the fit
(locations, score..), see below
locus2 : Same as above for channel 2
r : global coordinates of cell centroid (mid-point of cell)
error : segmentation error list
ehist : ehist is the sum of all errors in the region?s history
contactHist : ?
stat0 : stat0 flag that is true if the cell is born without error
The coord field contains a lot of cell specific info:
data.CellA{1}.coord =
A: Area of cell mask
r_center: geometrical center of the cell
box: coordinates of box surrounding cell
xaxis: coordinates of major axis
yaxis: coordinates of minor axis
I: Moment of inertia of cell mask
e1: priniple axis (major) unit vector
e2: priniple axis (minor) unit vector
rcm: center of mass position of mask
The pole field contains info pertaining to the cell pole and pole ages:
data.CellA{1}.pole =
e1: major axis direction.
op_ori: 1 if old pole is in direction of e1
-1 if old pole is in opposite e1
op_age: age of old pole in cell cycles
NaN if no birth is observed
np_age: age of new pole in cell cycles
The locus field contains info from locus fitting - if it ran.
data.CellA{1}.locus1(1) =
r: Spot position in global coords
score: score from spot fit.
intensity: raw intensity
b: spot width from fit
shortaxis: Locus position in local coords (short axis)
longaxis: Spot position in local coords (long axis)
INPUT :
dirname: seg folder eg. maindirectory/xy1/seg
CONST: are the segmentation constants.
header : string displayed with information
Copyright (C) 2016 Wiggins Lab
Written by Stella Stylianidou & Paul Wiggins.
University of Washington, 2016
This file is part of SuperSegger.
SuperSegger is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
SuperSegger is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with SuperSegger. If not, see <http://www.gnu.org/licenses/>.0001 function trackOptiMakeCell(dirname,CONST,header) 0002 % trackOptiMakeCell : generates the CellA field indexed by the region number 0003 % which contains information about each cell in each region in each frame. 0004 % 0005 % It goes through the dirname/*err.mat files and computes the characteristics of the 0006 % cells in each frame and puts them in the CellA structure that is indexed 0007 % by the region number, not the cell ID. This code figures out the pole 0008 % age and also computes statistics on the fluorescence channels as 0009 % well as fitting the locus positions. These last two features are 0010 % controlled by parameters set in the loadConstants.m or 0011 % loadConstantsMine.m files. 0012 % 0013 % data.CellA{1}. 0014 % mask : logical (1 and 0) cell mask, unoriented 0015 % r_offset : location of cell in full image from top left corner 0016 % BB : coordinates of cell bounding box, containing the 0017 % cell and the pad 0018 % edgeFlag : if the cell is on the edge of the image (bad) 0019 % phase : the cropped phase image of the cell 0020 % coord : coordinates, area, and orientation, see below 0021 % length : (1) length and (2) width of the cell 0022 % pole : orientation of the cell pole, see below 0023 % fluor1 : the 1st fluor channel image 0024 % fuor1mm : min and max of fluor1 0025 % fl1 : statistics of fluor1, e.g. background level 0026 % fluor2 : the 2nd fluor channel image 0027 % fluor2mm : min and max of fluor2 0028 % fl2 : statistics of fluor2 0029 % cell_dist : distance to the edge of the colony 0030 % gray : average phase gray value in cell region 0031 % locus1 : If focus fitting was run, data on the fit 0032 % (locations, score..), see below 0033 % locus2 : Same as above for channel 2 0034 % r : global coordinates of cell centroid (mid-point of cell) 0035 % error : segmentation error list 0036 % ehist : ehist is the sum of all errors in the region?s history 0037 % contactHist : ? 0038 % stat0 : stat0 flag that is true if the cell is born without error 0039 % 0040 % The coord field contains a lot of cell specific info: 0041 % data.CellA{1}.coord = 0042 % A: Area of cell mask 0043 % r_center: geometrical center of the cell 0044 % box: coordinates of box surrounding cell 0045 % xaxis: coordinates of major axis 0046 % yaxis: coordinates of minor axis 0047 % I: Moment of inertia of cell mask 0048 % e1: priniple axis (major) unit vector 0049 % e2: priniple axis (minor) unit vector 0050 % rcm: center of mass position of mask 0051 % 0052 % The pole field contains info pertaining to the cell pole and pole ages: 0053 % data.CellA{1}.pole = 0054 % e1: major axis direction. 0055 % op_ori: 1 if old pole is in direction of e1 0056 % -1 if old pole is in opposite e1 0057 % op_age: age of old pole in cell cycles 0058 % NaN if no birth is observed 0059 % np_age: age of new pole in cell cycles 0060 % 0061 % The locus field contains info from locus fitting - if it ran. 0062 % data.CellA{1}.locus1(1) = 0063 % r: Spot position in global coords 0064 % score: score from spot fit. 0065 % intensity: raw intensity 0066 % b: spot width from fit 0067 % shortaxis: Locus position in local coords (short axis) 0068 % longaxis: Spot position in local coords (long axis) 0069 % 0070 % INPUT : 0071 % dirname: seg folder eg. maindirectory/xy1/seg 0072 % CONST: are the segmentation constants. 0073 % header : string displayed with information 0074 % 0075 % Copyright (C) 2016 Wiggins Lab 0076 % Written by Stella Stylianidou & Paul Wiggins. 0077 % University of Washington, 2016 0078 % This file is part of SuperSegger. 0079 % 0080 % SuperSegger is free software: you can redistribute it and/or modify 0081 % it under the terms of the GNU General Public License as published by 0082 % the Free Software Foundation, either version 3 of the License, or 0083 % (at your option) any later version. 0084 % 0085 % SuperSegger is distributed in the hope that it will be useful, 0086 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0087 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0088 % GNU General Public License for more details. 0089 % 0090 % You should have received a copy of the GNU General Public License 0091 % along with SuperSegger. If not, see <http://www.gnu.org/licenses/>. 0092 0093 if ~exist('header','var') 0094 header = []; 0095 end 0096 0097 if(nargin<1 || isempty(dirname)) 0098 dirname = '.'; 0099 end 0100 dirname = fixDir(dirname); 0101 0102 verbose = CONST.parallel.verbose; 0103 0104 % Get the track/error file names 0105 contents=dir([dirname '*_err.mat']); 0106 num_im = numel(contents); 0107 0108 if CONST.parallel.show_status 0109 h = waitbar( 0, 'Making Cells.'); 0110 cleanup = onCleanup( @()( delete( h ) ) ); 0111 else 0112 h = []; 0113 end 0114 0115 % loop through all the cells. 0116 for i = 1:num_im; 0117 0118 if (i ==1) && (1 == num_im) % snapshots 0119 data_r = []; 0120 data_c = loaderInternal([dirname,contents(i).name]); 0121 data_f = []; 0122 elseif i == 1; % first frame 0123 data_r = []; 0124 data_c = loaderInternal([dirname,contents(i).name]); 0125 data_f = loaderInternal([dirname,contents(i+1).name]); 0126 elseif i==num_im; % last frame 0127 data_r = loaderInternal([dirname,contents(i-1).name]); 0128 data_c = loaderInternal([dirname,contents(i).name]); 0129 data_f = []; 0130 else 0131 data_r = loaderInternal([dirname,contents(i-1).name]); 0132 data_c = loaderInternal([dirname,contents(i).name]); 0133 data_f = loaderInternal([dirname,contents(i+1).name]); 0134 end 0135 0136 % first frame, find the total number of channels 0137 if i==1 0138 nc = 0; 0139 tmp_fn = fieldnames( data_c ); 0140 nf = numel( tmp_fn ); 0141 for j = 1:nf; 0142 if numel(strfind(tmp_fn{j},'fluor')==1)&& ~numel((strfind(tmp_fn{j},'fluor0'))) 0143 nc = nc+1; 0144 end 0145 end 0146 end 0147 0148 % set min max of fluorescence 0149 for j = 1:nc 0150 tmp = getfield( data_c,['fluor',num2str(j)]); 0151 ff(j,:) = [min(tmp(:)),max(tmp(:))]; 0152 end 0153 0154 % make the cell array holding the cells that has one element for each 0155 % region. 0156 data_c.CellA = cell(1,data_c.regs.num_regs); 0157 dist_mask = makeColonyDist( data_c.mask_bg ); 0158 0159 for ii = 1:data_c.regs.num_regs 0160 % Cut out a region with pixel pad size PAD_SIZE. 0161 celld = struct(); 0162 PAD_SIZE = 5; 0163 ss = size(data_c.phase); 0164 0165 celld.cellLength = [data_c.regs.L1(ii),data_c.regs.L2(ii)]; 0166 [xx,yy] = getBBpad( data_c.regs.props(ii).BoundingBox, ss, PAD_SIZE); 0167 celld.xx = xx; 0168 celld.yy = yy; 0169 celld.mask = logical(data_c.regs.regs_label(yy,xx)==ii); 0170 celld.r_offset = [xx(1),yy(1)]; 0171 celld.BB = [xx(1),yy(1),xx(end)-xx(1),yy(end)-yy(1)]; 0172 tmpEdge = [ ceil(data_c.regs.props(ii).BoundingBox(1)),... 0173 ceil(data_c.regs.props(ii).BoundingBox(2)),... 0174 floor(sum(data_c.regs.props(ii).BoundingBox([1,3]))),... 0175 floor(sum(data_c.regs.props(ii).BoundingBox([2,4])))]; 0176 celld.edgeFlag = any( tmpEdge == [1,1,ss(2),ss(1)]); 0177 0178 % copy information from regions 0179 if CONST.trackOpti.NEIGHBOR_FLAG 0180 celld.contactHist = data_c.regs.contactHist(ii); 0181 end 0182 0183 0184 % record the number of cell neighbors 0185 if CONST.trackOpti.NEIGHBOR_FLAG 0186 nei_ = numel(trackOptiNeighbors(data_c,ii)); 0187 data_c.regs.numNeighbors{ii} = nei_ ; 0188 celld.numNeighbors = nei_ ; 0189 end 0190 0191 celld.phase = data_c.phase(yy,xx); 0192 0193 % Keep track of pole age and what direction the old pole is in. 0194 if isempty(data_r) || isempty(data_c.regs.map.r{ii}) 0195 celld = toMakeCell(celld,[],data_c.regs.props(ii)); 0196 celld.pole.e1 = celld.coord.e1; 0197 celld.pole.op_ori = 0; 0198 celld.pole.op_age = NaN; 0199 celld.pole.np_age = NaN; 0200 else 0201 if data_c.regs.error.r(ii) 0202 try 0203 celld = toMakeCell(celld, data_r.CellA{data_c.regs.map.r{ii}(1)}.pole.e1,data_c.regs.props(ii)); 0204 catch ME 0205 printError(ME); 0206 end 0207 celld.pole.e1 = celld.coord.e1; 0208 celld.pole.op_ori = 0; 0209 celld.pole.op_age = NaN; 0210 celld.pole.np_age = NaN; 0211 0212 elseif data_c.regs.birthF(ii) && ( data_c.regs.sisterID(ii) ) && ~isempty(find( data_c.regs.sisterID(ii) == data_c.regs.ID )) 0213 0214 cell_old = data_r.CellA{data_c.regs.map.r{ii}(1)}; 0215 celld = toMakeCell(celld, cell_old.pole.e1,data_c.regs.props(ii)); 0216 celld.pole.e1 = celld.coord.e1; 0217 e1 = celld.pole.e1; 0218 op_ori = cell_old.pole.op_ori; 0219 0220 jj = find( data_c.regs.sisterID(ii) == data_c.regs.ID ); 0221 jj = jj(1); 0222 rs = data_c.regs.props(jj).Centroid; 0223 r0 = data_c.regs.props(ii).Centroid; 0224 0225 dr = (r0-rs)*e1; 0226 celld.pole.op_ori = sign(dr); 0227 0228 if celld.pole.op_ori == cell_old.pole.op_ori 0229 celld.pole.op_age = cell_old.pole.op_age+1; 0230 else 0231 celld.pole.op_age = cell_old.pole.np_age+1; 0232 end 0233 0234 celld.pole.np_age = 1; 0235 0236 debug_flag = 0; 0237 if debug_flag 0238 clf; 0239 kk = find( data_c.regs.motherID(ii) == data_r.ID ); 0240 imshow( 0.5*cat(3, ag(data_c.regs.regs_label==jj), ag(data_c.regs.regs_label==ii),ag(data_r.regs_label==kk))); 0241 hold on; 0242 0243 tmp = celld; 0244 r = tmp.coord.r_center; 0245 xaxisx = r(1) + [0,tmp.length(1)*tmp.coord.e1(1)]/2; 0246 xaxisy = r(2) + [0,tmp.length(1)*tmp.coord.e1(2)]/2; 0247 yaxisx = r(1) + [0,tmp.length(2)*tmp.coord.e2(1)]/2; 0248 yaxisy = r(2) + [0,tmp.length(2)*tmp.coord.e2(2)]/2; 0249 old_pole = r + tmp.length(1)*tmp.coord.e1'*tmp.pole.op_ori/2; 0250 new_pole = r - tmp.length(1)*tmp.coord.e1'*tmp.pole.op_ori/2; 0251 0252 plot( r(1), r(2), 'g.'); 0253 plot( xaxisx, xaxisy, 'g-'); 0254 plot( yaxisx, yaxisy, 'g-'); 0255 plot( old_pole(1), old_pole(2), 'r.'); 0256 plot( new_pole(1), new_pole(2), 'w.'); 0257 0258 tmp = data_r.CellA{kk}; 0259 r = tmp.coord.r_center; 0260 r = tmp.coord.r_center; 0261 xaxisx = r(1) + [0,tmp.length(1)*tmp.coord.e1(1)]/2; 0262 xaxisy = r(2) + [0,tmp.length(1)*tmp.coord.e1(2)]/2; 0263 yaxisx = r(1) + [0,tmp.length(2)*tmp.coord.e2(1)]/2; 0264 yaxisy = r(2) + [0,tmp.length(2)*tmp.coord.e2(2)]/2; 0265 old_pole = r + tmp.length(1)*tmp.coord.e1'*tmp.pole.op_ori/2; 0266 new_pole = r - tmp.length(1)*tmp.coord.e1'*tmp.pole.op_ori/2; 0267 0268 plot( r(1), r(2), 'g.'); 0269 plot( xaxisx, xaxisy, 'g-'); 0270 plot( yaxisx, yaxisy, 'g-'); 0271 if ~tmp.pole.op_ori 0272 plot( old_pole(1), old_pole(2), 'r.'); 0273 plot( new_pole(1), new_pole(2), 'w.'); 0274 end 0275 end 0276 else 0277 if (data_c.regs.map.r{ii}(1) == 0) 0278 data_c.regs.map.r{ii}(1) = []; 0279 end 0280 celld = toMakeCell(celld, data_r.CellA{data_c.regs.map.r{ii}(1)}.pole.e1,data_c.regs.props(ii)); 0281 celld.pole = data_r.CellA{data_c.regs.map.r{ii}(1)}.pole; 0282 celld.pole.e1 = celld.coord.e1; 0283 end 0284 end 0285 0286 0287 % copy fluorescence fields, compute fluorescence statistics and find loci 0288 for j = 1:nc 0289 tmp = data_c.(['fluor',num2str(j)]); 0290 celld.(['fluor',num2str(j)]) = tmp(yy,xx); 0291 celld.(['fluor',num2str(j),'mm']) = ff(j,:) ; 0292 0293 if isfield( CONST.trackLoci, 'fluorFlag' ) && CONST.trackLoci.fluorFlag 0294 tmp = trackOptiCellFluor( tmp(yy,xx), celld.mask, celld.r_offset); 0295 else 0296 tmp = []; 0297 end 0298 tmp.bg = data_c.(['fl',num2str(j),'bg']); 0299 celld.(['fl',num2str(j)]) = tmp; 0300 0301 end 0302 0303 % compute the distance to the edge of the colony 0304 if data_c.regs.ID(ii) 0305 0306 [xx,yy] = getBBpad( data_c.regs.props(ii).BoundingBox, ss, 0); 0307 mask_cell = logical(data_c.regs.regs_label(yy,xx)==ii); 0308 dist_mask_crop = dist_mask(yy,xx); 0309 cell_dist = min(dist_mask_crop(mask_cell)); 0310 celld.cell_dist = cell_dist; 0311 0312 % calculate average phase gray value in cell region 0313 celld.gray = mean(double(celld.phase(celld.mask))); 0314 0315 end 0316 0317 data_c.CellA{ii} = celld; 0318 0319 end 0320 0321 % save the updated err files. 0322 dataname = [dirname,contents(i).name]; 0323 save(dataname,'-STRUCT','data_c'); 0324 0325 if CONST.parallel.show_status 0326 waitbar(i/num_im,h,['Making Cells--Frame: ',num2str(i),'/',num2str(num_im)]); 0327 elseif verbose 0328 disp([header, 'MakeCell frame: ',num2str(i),' of ',num2str(num_im)]); 0329 end 0330 end 0331 0332 if CONST.parallel.show_status 0333 close(h); 0334 end 0335 0336 end 0337 0338 function data = loaderInternal(filename) 0339 data = load( filename ); 0340 end