Description of makeTowerCons

0001 function  imData = makeTowerCons( data, CONST, xdim, ...
0002     disp_flag, skip, mag, fnum )
0003 % makeTowerCons : Makes normalized time and shape tower for a single cell.
0004 % It normalizes for time and cell shape.
0005 %
0006 % INPUT:
0007 %        data : (struct) Cell data to construct the Time-Lapse Tower
0008 %       CONST :  segmentation constants
0009 %        xdim : dimensions of image along x axis
0010 %   disp_flag : 0 or 1 value to show imagesa and progess
0011 %        skip : number of frames to be skipped when generating the tower
0012 %         mag : rescale image by mag
0013 %         fnum : flurescence channel to be used
0014 %
0015 % OUPUT:
0016 %       imData
0017 %           .tower : Grayscale masked consensus image for a cell
0018 %           .towerNorm : Normalized grayscale masked consensus image
0019 %           .towerC : Colorized masked cons image of a single cell
0020 %           .towerNormC : Normalized colorized masked cons image of a cell
0021 %           .towerRaw : Raw tower for a single cell, not scaled, or cropper
0022 %           .towerNormRaw : Normalized, not cropped (raw) tower for a cell
0023 %           .towerMask : binary mask of the cell for the tower
0024 %           .fmax : min and max intensity
0025 %           .imCell : cell array of the raw consenus images
0026 %           .imCellNorm : cell array of the raw normalized consenus images
0027 %           .maskCell  : cell array of binary cell masks
0028 %           .intWeight :
0029 %
0030 % Copyright (C) 2016 Wiggins Lab
0031 % Written by Paul Wiggins & Stella Stylianidou.
0032 % University of Washington, 2016
0033 % This file is part of SuperSegger.
0034 %
0035 % SuperSegger is free software: you can redistribute it and/or modify
0036 % it under the terms of the GNU General Public License as published by
0037 % the Free Software Foundation, either version 3 of the License, or
0038 % (at your option) any later version.
0039 %
0040 % SuperSegger is distributed in the hope that it will be useful,
0041 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0042 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0043 % GNU General Public License for more details.
0044 %
0045 % You should have received a copy of the GNU General Public License
0046 % along with SuperSegger.  If not, see <http://www.gnu.org/licenses/>.
0047 
0048 
0049 % set the color map here to generate the color images
0050 persistent colormap_;
0051 if isempty( colormap_ )
0052     colormap_ = colormap( 'jet' );
0053     %colormap_ = colormap( 'hsv' );
0054 end
0055 
0056 
0057 % fill unset input arguments
0058 if ~exist('skip','var') || isempty( skip )
0059     skip = 1;
0060 end
0061 
0062 if ~exist('disp_flag', 'var' ) || isempty( disp_flag )
0063     disp_flag = true;
0064 end
0065 
0066 if ~exist( 'xdim', 'var' ) || isempty( xdim )
0067     xdim = 1;
0068 end
0069 
0070 if isfield( CONST, 'view' ) && isfield( CONST.view, 'orientFlag' )
0071     orientFlag = CONST.view.orientFlag;
0072 else
0073     orientFlag = true;
0074 end
0075 
0076 
0077 if  ~isfield(data.CellA{1},'fluor1')
0078     disp ('no fluorescence field found');
0079     return;
0080 end
0081 
0082 % get the TimeStep variable for plotting
0083 TimeStep     = CONST.getLocusTracks.TimeStep;
0084 
0085 % get the number of frames
0086 numframe = numel( data.CellA );
0087 t0 = numframe;
0088 
0089 
0090 % set the consensus lengths
0091 L0 = 26*mag; % length of the cell in the first frame
0092 W0 =  9*mag; % width of the cells
0093 
0094 % T0 is the number of frames in the consensus
0095 if isfield( CONST.view, 'numCons' ) && ~isempty( CONST.view.numCons )
0096     T0 =  CONST.view.numCons;
0097 else
0098     T0 =  8;
0099 end
0100 
0101 % init variables
0102 imCell = cell( 1, T0 );
0103 ssCell = imCell;
0104 nii   = zeros(1, T0 );
0105 maskCell = imCell;
0106 imCellNorm = cell( 1, T0 );
0107 intWeight = zeros(1, T0 );
0108 
0109 
0110 f1mm = []; % init the max and min fluor values to empty.
0111 mm_name = ['fluor', num2str(fnum),'mm'];
0112 f_name  = ['fluor', num2str(fnum)];
0113 
0114 
0115 
0116 % figure out the mapping between the different time coords
0117 TT = (0:(T0-1))/(T0-1);
0118 tt = (0:(t0-1))/(t0-1);
0119 Td = abs( (0*tt'+1)*TT-tt'*(0*TT+1));
0120 [~,ord_tt] = min( Td' );
0121 [~,ord_TT] = min( Td );
0122 
0123 
0124 % straighten and align the cells, frame by frame
0125 for Tii = 1:T0
0126     
0127     ind = find( ord_tt == Tii );
0128     if isempty( ind )
0129         ind = ord_TT(Tii);
0130     end
0131     nii(Tii) = numel( ind );
0132     
0133     for ii = row(ind)
0134         % intDoCellOri: (i) rotate the image, (ii) adjust cell width to W0
0135         [imCell__, maskCell__, alpha__, ssCell__, xxCell__, yyCell__] = ...
0136             intDoCellOri( data.CellA{ii}, W0, L0, T0, Tii, mag, fnum );
0137         
0138         
0139         if isempty( imCell{Tii} )
0140             imCell{Tii}   = imCell__/nii(Tii);
0141             maskCell{Tii} = maskCell__;
0142             ssCell{Tii}   = ssCell__;
0143         else
0144             imCell{Tii}   = imCell{Tii} + imCell__/nii(Tii);
0145         end
0146         
0147         
0148         % set the max and min fluor values
0149         if isempty(f1mm)
0150             if isfield( data.CellA{ii}, mm_name )
0151                 f1mm = data.CellA{ii}.(mm_name) ;
0152             else
0153                 tmp_fluor = data.CellA{ii}.(f_name);
0154                 f1mm = [ min(tmp_fluor(:)), max(tmp_fluor(:))];
0155             end
0156         else
0157             if isfield( data.CellA{ii}, mm_name )
0158                 tmp_fluor = data.CellA{ii}.(mm_name) ;
0159                 f1mm = [ min([tmp_fluor,f1mm(1)]), max([tmp_fluor(2),f1mm(2)]) ];
0160             else
0161                 tmp_fluor = data.CellA{ii}.(f_name);
0162                 f1mm = [ min([f1mm(1),min( tmp_fluor(:))]), max([f1mm(2),max(tmp_fluor(:))])];
0163             end
0164         end
0165     end
0166     
0167     
0168     im    = imCell{Tii};
0169     mask  = logical(maskCell{Tii});
0170     
0171     im = im - min( im( mask ) );
0172     
0173     intWeight( Tii ) = sum( double(im(mask)))./sum( double(mask(:)));
0174     imCellNorm{Tii} = imCell{Tii}/intWeight( Tii ) ;
0175     
0176 end
0177 
0178 
0179 
0180 % used to calculate scaled images so that they all have dimensions L0, W0, and T0
0181 % [imCell, maskCell, ssCell, xxCell, yyCell, imCellS, maskCellS, ...
0182 %     intWeight, imCellNorm, imCellNormS, intWeightS] = ...
0183 %     intDoRescale( imCell, maskCell, ssCell, xxCell, yyCell, L0, W0, T0 );
0184 
0185 % Merge the images from the arrays into a single image
0186 [ imColorCons, imBWCons, towerIm, maskCons, nx, ny, max_x, max_y ] = ...
0187     towerMergeImages( imCell, maskCell, ssCell, xdim, skip, mag, CONST );
0188 
0189 % Merge the normalized images from the arrays into a single image
0190 [ imColorConsNorm, imBWConsNorm, towerImNorm, maskCons, nx, ny, max_x, max_y ] = ...
0191     towerMergeImages( imCellNorm, maskCell, ssCell, xdim, skip, mag, CONST );
0192 
0193 % if the disp flag is true, show the image.
0194 if disp_flag
0195     intDoDraw( imColorCons, T0, nx, ny, max_x, max_y, TimeStep, CONST );
0196     pause;
0197 end
0198 
0199 imData               = [];
0200 imData.tower         = imBWCons;
0201 imData.towerNorm     = imBWConsNorm;
0202 imData.towerC        = imColorCons;
0203 imData.towerNormC    = imColorConsNorm;
0204 imData.towerRaw      = towerIm;
0205 imData.towerNormRaw  = towerImNorm;
0206 imData.towerMask     = maskCons;
0207 imData.fmax          = f1mm;
0208 imData.imCell        = imCell;
0209 imData.imCellNorm    = imCellNorm;
0210 imData.maskCell      = maskCell;
0211 imData.intWeight     = intWeight;
0212 
0213 % .imCellScale   : Grayscale non-crop cons.
0214 % .imCellNormScale : Grayscale normalized non-crop cons.
0215 % .maskCellScale : double binary mask of the cell.
0216 %imData.imCellScale   = imCellScale;
0217 %imData.imCellNormScale  = imCellNormS;
0218 %imData.maskCellScale = maskCellScale;
0219 %imData.intWeightS    = intWeightS;
0220 
0221 
0222 end
0223 
0224 
0225 
0226 function intDoDraw( im, T0, nx, ny, max_x, max_y, TimeStep, CONST )
0227 % intDoDraw : draws the image.
0228 
0229 imshow( im );
0230 
0231 if CONST.view.falseColorFlag
0232     cc = 'w';
0233 else
0234     cc = 'b';
0235 end
0236 
0237 hold on;
0238 
0239 for ii = 1:T0
0240     yy = floor((ii-1)/nx);
0241     xx = ii-yy*nx-1;
0242     y = 1+yy*max_y;
0243     x = 1+xx*max_x;
0244     text( x+2, y+2, num2str((ii-1)*TimeStep),'Color',cc,'FontSize',12,'VerticalAlignment','Top');
0245 end
0246 
0247 dd = [1,ny*max_y+1];
0248 for xx = 1:(nx-1)
0249     plot( 0*dd + 1+xx*max_x, dd,[':',cc]);
0250 end
0251 
0252 dd = [1,nx*max_x+1];
0253 for yy = 1:(ny-1)
0254     plot( dd, 0*dd + 1+yy*max_y, [':',cc]);
0255 end
0256 end
0257 
0258 function [fluor1, mask_rot, alpha, ss, xx, yy] = intDoCellOri( ...
0259     celld, W0, L0, T0, Tii, mag, fnum )
0260 % intDoCellOri : orients the cell.
0261 
0262 f_name = ['fluor',num2str(fnum)];
0263 fl_name = ['fl',num2str(fnum)];
0264 debug_flag = false;
0265 Lii = L0*(1+(Tii-1)/(T0-1));
0266 
0267 persistent strel1;
0268 if isempty( strel1 )
0269     strel1 = strel('square',3);
0270 end
0271 
0272 if isfield( celld, 'pole' ) && ~isnan( celld.pole.op_ori ) && (celld.pole.op_ori ~= 0)
0273     ssign = sign(celld.pole.op_ori);
0274 else
0275     ssign = 1;
0276 end
0277 
0278 e1 = celld.coord.e1;
0279 alpha = 90-180/pi*atan2(e1(1),e1(2)) + 180*double(ssign==1);
0280 mask = celld.mask;
0281 
0282 % get the background fluorescence
0283 if isfield( celld, fl_name) && ...
0284         isfield( getfield( celld, fl_name), 'bg' ) && ...
0285         ~isnan( getfield( getfield(celld, fl_name), 'bg'))
0286     bg_fluor =  double( getfield( getfield( celld, fl_name), 'bg' ));
0287 else
0288     bg_fluor = 0;
0289 end
0290 
0291 
0292 % check to see if the cell runs up to the end
0293 xx_ = sum( mask, 1);
0294 xx_ = [find(logical(xx_), 1, 'first' ), find(logical(xx_), 1, 'last' )];
0295 
0296 yy_ = sum( mask, 2);
0297 yy_ = [find(logical(yy_), 1, 'first' ), find(logical(yy_), 1, 'last' )];
0298 
0299 ss_old = size( mask );
0300 fluor1  = double(getfield(celld,f_name))-bg_fluor;
0301 
0302 if xx_(1) < 3
0303     mask   = [ false( [ss_old(1),6] ), mask  ];
0304     fluor1 = [ zeros( [ss_old(1),6] ), fluor1];
0305 end
0306 
0307 ss_old = size( mask );
0308 
0309 if xx_(2) > ss_old(2)-2
0310     mask   = [ mask,   false( [ss_old(1),6] ) ];
0311     fluor1 = [ fluor1, zeros( [ss_old(1),6] ) ];
0312 end
0313 
0314 ss_old = size( mask );
0315 
0316 if yy_(1) < 3
0317     mask   = [ false( [6,ss_old(2)] ); mask  ];
0318     fluor1 = [ zeros( [6,ss_old(2)] ); fluor1];
0319 end
0320 
0321 ss_old = size( mask );
0322 
0323 if yy_(2) > ss_old(1)-2
0324     mask   = [ mask;   false( [6,ss_old(2)] ) ];
0325     fluor1 = [ fluor1; zeros( [6,ss_old(2)] ) ];
0326 end
0327 
0328 mask     = imrotate(       (mask),         alpha, 'bilinear' );
0329 mask = logical(imdilate(mask,strel1));
0330 fluor1   = imrotate( fluor1, alpha, 'bilinear');
0331 
0332 % scale the cell to the target size
0333 
0334 % compute the local width
0335 w = sum( mask,1 );
0336 
0337 w0 = median( w(logical(w)) );
0338 l0 = find( logical(w), 1, 'last' )-find( logical(w), 1, 'first' )+1;
0339 
0340 % magnify length and width independently to bring cell to approximately the
0341 % right size.
0342 mag_x = Lii/l0;
0343 mag_y = W0/w0;
0344 mag_v = [mag_y,mag_x];
0345 ss_old = size( mask );
0346 
0347 
0348 mask_rot = imresize( mask,   ss_old.*mag_v );
0349 fluor1   = imresize( fluor1, ss_old.*mag_v );
0350 
0351 % redefine width and length
0352 w = sum( mask_rot,1 );
0353 w0 = median( w(logical(w)) );
0354 x0_srt = find( logical(w), 1, 'first' ); % starting x coord
0355 x0_end = find( logical(w), 1, 'last' ); % ending x coord
0356 l0 = x0_end-x0_srt+1; % cell length
0357 
0358 ss = size( mask_rot );
0359 x = 1:ss(2);
0360 y = 1:ss(1);
0361 
0362 x_mid = 0.5*(x(1)+x(end));
0363 y_mid = 0.5*(y(1)+y(end));
0364 
0365 [X,Y] = meshgrid( x, y );
0366 
0367 % compute the center of intensity
0368 y_cen0 = sum(mask_rot.*Y)./w;
0369 
0370 x0_srt_ = x0_srt+floor(mag/2);
0371 x0_end_ = x0_end-floor(mag/2);
0372 dx0 = x0_end_-x0_srt_;
0373 
0374 ranger = floor( dx0*(0:.2:1)+x0_srt_);
0375 
0376 y_cen = interp1( ranger, y_cen0(ranger), x, 'linear','extrap' );
0377 
0378 % compute the model (theory) width
0379 RADIUS = W0/2;
0380 wt = intCellFit( x, x0_srt-1, x0_end+1, RADIUS );
0381 
0382 mask_rot_  = mask_rot;
0383 mask_rot__ = mask_rot;
0384 fluor1_    = fluor1;
0385 
0386 
0387 if debug_flag
0388     figure(1);
0389     clf;
0390     imshow( mask_rot_, []);
0391     hold on;
0392     plot( y_cen, 'r.-' );
0393 end
0394 
0395 
0396 for ii = x
0397     dy1 = (y - y_cen(ii));
0398     dy2 = (y - y_mid);
0399     
0400     mask_rot_(:,ii) = 0;
0401     if wt(ii)~=0
0402         mask_rot_(abs(dy2)<=wt(ii)/2,ii) = 1;
0403     end
0404     
0405     mask_rot__(:,ii)   = interp1( dy1, double(mask_rot__(:,ii)),   dy2,'nearest','extrap' );
0406     fluor1_(:,ii)   = interp1( dy1, fluor1(:,ii),   dy2,'nearest','extrap' );
0407 end
0408 
0409 if debug_flag
0410     figure(2);
0411     clf;
0412     imshow( fluor1_,[] );
0413     colormap jet
0414     
0415     figure(3);
0416     clf;
0417     imshow( mask_rot_,[] );
0418     colormap jet
0419     
0420     figure(4);
0421     clf;
0422     imshow( mask_rot__,[] );
0423     colormap jet
0424 end
0425 
0426 w  = sum(mask_rot_);
0427 xmin_ = max([1,find(w>0,1,'first')-1]);
0428 xmax_ = min([ss(2),find(w>0,1, 'last')+1]);
0429 
0430 ysum  = sum(mask_rot_');
0431 ymin_ = max([1,find(ysum>0,1,'first')-1]);
0432 ymax_ = min([ss(1),find(ysum>0,1, 'last')+1]);
0433 
0434 
0435 yy = ymin_-1+(1:(W0+2));
0436 xx = xmin_-1+(1:(Lii+2));
0437 
0438 
0439 mask_rot = mask_rot_(yy, xx);
0440 fluor1   = fluor1_(yy, xx);
0441 
0442 ss = size( mask_rot );
0443 
0444 end
0445 
0446 function [X] = intDoFit( x, ysum )
0447 % intDoFit : fits to the theoretical shape of the cell
0448 
0449 X(1) = find( ysum, 1, 'first');
0450 X(2) = find( ysum, 1, 'last');
0451 X(3) = max( ysum );
0452 
0453 X = fminsearch( @intDoFitInt, X );
0454 
0455     function err = intDoFitInt( X )
0456         y   = intCellFit( x, X(1), X(2), X(3) );
0457         
0458         err = sum( (ysum-y).^2 );
0459     end
0460 
0461 
0462 if 0
0463     clf;
0464     plot( x, ysum, 'y.-');
0465     hold on;
0466     plot( x, intCellFit( x, X(1), X(2), X(3) ), 'r.-');
0467     pause;
0468 end
0469 end
0470 
0471 function [imCell__, maskCell__, ssCell__, xxCell__, yyCell__, imCellScale, ...
0472     maskCellScale, intWeight, imCellNorm__, imCellNormS, intWeightS ] = ...
0473     intDoRescale( imCell, maskCell, ssCell, xxCell, yyCell, L0, W0, T0 )
0474 % intDoRescale : rescales the cell in both space and time and fixes
0475 % the normalization over time of the tower.
0476 
0477 
0478 nt = numel( imCell );
0479 
0480 % Renormalize the length of the cell to make it L = L0*(1+t/T)
0481 for ii = 1:nt
0482     tii = (ii-1)/(nt-1);
0483     ss = size( imCell{ii} );
0484     
0485     x = 1:ss(2);
0486     y = 1:ss(1);
0487     
0488     [ X, Y ] = meshgrid( x, y );
0489     
0490     xsum  = sum(maskCell{ii});
0491     
0492     xcom  = sum(maskCell{ii}(:).*X(:))/sum(maskCell{ii}(:));
0493     ycom  = sum(maskCell{ii}(:).*Y(:))/sum(maskCell{ii}(:));
0494     
0495     xi = (1:(2*L0+4))-(2*L0+4-1)/2;
0496     yi = (1:(W0+4))-(W0+4-1)/2;
0497     
0498     [Xi,Yi] = meshgrid( xi, yi );
0499     
0500     xmin_ = max([1,find(xsum>1,1,'first')-1]);
0501     xmax_ = min([ssCell{ii}(2),find(xsum>1,1, 'last')+1]);
0502     dx = xmax_ - xmin_ + 1;
0503     
0504     imCell_{ii}   = interp2( (X-xcom)/dx*2*L0, Y-ycom,   imCell{ii}, Xi-1, Yi-1);
0505     maskCell_{ii} = interp2( (X-xcom)/dx*2*L0, Y-ycom, double(maskCell{ii}), Xi-1, Yi-1);
0506     
0507 end
0508 
0509 %  Renormalize the length of the cell cycle to be length T0.
0510 for ii = 1:T0
0511     
0512     jj = (ii-1)/(T0-1)*(nt-1)+1;
0513     
0514     jjm = floor(jj);
0515     jjp = jjm + 1;
0516     djj = jj-jjm;
0517     
0518     if ii == 1
0519         imCell__{ii} = imCell_{ii};
0520         maskCell__{ii} = maskCell_{ii};
0521     elseif ii == T0
0522         imCell__{ii} = imCell_{end};
0523         maskCell__{ii} = maskCell_{end};
0524     else
0525         imCell__{ii} = (1-djj)*imCell_{jjm} + djj*imCell_{jjp};
0526         maskCell__{ii} = (1-djj)*maskCell_{jjm} + djj*maskCell_{jjp};
0527     end
0528     
0529 end
0530 
0531 % save the scaled versions
0532 imCellScale   = imCell__;
0533 imCellNormS   = imCell__;
0534 
0535 maskCellScale = maskCell__;
0536 intWeight = zeros(1,T0);
0537 intWeightS = zeros(1,T0);
0538 
0539 imCellNorm__ =  imCell__;
0540 
0541 % rescale the long axis dimension to scale uniformly
0542 for ii = 1:T0
0543     tii = (ii-1)/(T0-1);
0544     Lii = L0*(1+tii);
0545     
0546     xi = (1:(2*L0+4))-(2*L0+4-1)/2;
0547     yi = (1:(W0+4))-(W0+4-1)/2;
0548     
0549     [Xi,Yi] = meshgrid( xi, yi );
0550     
0551     imCell__{ii}   = interp2(  (Xi-1), Yi-1,   imCell__{ii}, (Xi-1)*(2*L0)/Lii, Yi-1);
0552     maskCell__{ii} = interp2(  (Xi-1), Yi-1, maskCell__{ii}, (Xi-1)*(2*L0)/Lii, Yi-1);
0553     
0554     imCell__{ii}(isnan(imCell__{ii}))     = 0;
0555     maskCell__{ii}(isnan(maskCell__{ii})) = 0;
0556     
0557     % normalizes fluor so that the mean fluor remains constant over time.
0558     % the abs value does not allow the sum to diverge. The sum can become 0
0559     % due to illumination non-unifority and real life, causing problems later.
0560     sum_ = sum(abs(imCell__{ii}(:)).*maskCell__{ii}(:))/sum(maskCell__{ii}(:));
0561     imCellNorm__{ii} = imCell__{ii}/sum_;
0562     intWeight(ii) = sum_;
0563     
0564     % normalize fluor so that the mean fluor remains constant over time.
0565     imCellScale{ii}(isnan(imCellScale{ii}))      = 0;
0566     maskCellScale{ii}(isnan(maskCellScale{ii}))  = 0;
0567     
0568     sumS_ = sum(imCellScale{ii}(:).*maskCellScale{ii}(:))/sum(maskCellScale{ii}(:));
0569     imCellNormS{ii} = imCellScale{ii}/sumS_;
0570     intWeightS(ii) = sumS_;
0571     ssCell__{ii} = size( imCell__{ii} );
0572     xxCell__{ii} = 1:ssCell__{ii}(2);
0573     yyCell__{ii} = 1:ssCell__{ii}(1);
0574     
0575 end
0576 end
0577 
0578 
0579
makeTowerCons

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
