trackOptiAlignPad

 trackOptiAlignPad : aligns phase images to correct for microscope drift.
 To keep as much data as possible, instead of cropping the resulting
 images it builds a larger image that encompases all drift positions.
 It saves the alignment information in an array called crop_box and saves
 the aligned images in a directory called dirname/align.
 If the images are out of focused or have a high error they are skipped
 in alignment. Images can be aligned to any channel (fluorescence or phase)
 by setting CONST.imAlign.AlignChannel (default is 1).

 INPUT :
       dirname_ : folder were .tif image files in NIS format are contained
       workers : number of workers for parallel computation
       CONST : segmentation constants
 OUTPUT :
       crop_box : information about alignement


 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 [crop_box] = trackOptiAlignPad(dirname_, workers, CONST, targetd)
0002 % trackOptiAlignPad : aligns phase images to correct for microscope drift.
0003 % To keep as much data as possible, instead of cropping the resulting
0004 % images it builds a larger image that encompases all drift positions.
0005 % It saves the alignment information in an array called crop_box and saves
0006 % the aligned images in a directory called dirname/align.
0007 % If the images are out of focused or have a high error they are skipped
0008 % in alignment. Images can be aligned to any channel (fluorescence or phase)
0009 % by setting CONST.imAlign.AlignChannel (default is 1).
0010 %
0011 % INPUT :
0012 %       dirname_ : folder were .tif image files in NIS format are contained
0013 %       workers : number of workers for parallel computation
0014 %       CONST : segmentation constants
0015 % OUTPUT :
0016 %       crop_box : information about alignement
0017 %
0018 %
0019 % Copyright (C) 2016 Wiggins Lab
0020 % Written by Stella Stylianidou & Paul Wiggins.
0021 % University of Washington, 2016
0022 % This file is part of SuperSegger.
0023 %
0024 % SuperSegger is free software: you can redistribute it and/or modify
0025 % it under the terms of the GNU General Public License as published by
0026 % the Free Software Foundation, either version 3 of the License, or
0027 % (at your option) any later version.
0028 %
0029 % SuperSegger is distributed in the hope that it will be useful,
0030 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0031 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0032 % GNU General Public License for more details.
0033 %
0034 % You should have received a copy of the GNU General Public License
0035 % along with SuperSegger.  If not, see <http://www.gnu.org/licenses/>.
0036 
0037 
0038 if exist( dirname_, 'dir' )
0039 
0040 if ~exist('workers', 'var') || isempty( workers )
0041     workers = 0;
0042 end
0043 
0044 % Upsampling factor used for image alignment.
0045 % Images will be registered to within 1/precision of a pixel
0046 precision = 100;
0047 
0048 if dirname_ == '.'
0049     dirname_ = pwd;
0050 end
0051 dirname_ = fixDir(dirname_);
0052 
0053 verbose = CONST.parallel.verbose;
0054 file_filter = '*.tif';
0055 contents=dir([dirname_ file_filter]);
0056 num_im = numel( contents );
0057 
0058 nt  = [];
0059 nc  = [];
0060 nxy = [];
0061 nz  = [];
0062 
0063 % extract name information for each image
0064 for i = 1:num_im
0065     nameInfo = ReadFileName(contents(i).name);
0066     nt  = [nt, nameInfo.npos(1,1)];
0067     nc  = [nc, nameInfo.npos(2,1)];
0068     nxy = [nxy,nameInfo.npos(3,1)];
0069     nz  = [nz, nameInfo.npos(4,1)];
0070 end
0071 
0072 nt  = sort(unique(nt));
0073 nc  = sort(unique(nc));
0074 nxy = sort(unique(nxy));
0075 nz  = sort(unique(nz));
0076 
0077 
0078 
0079 % set the align channel. Default is 1 (bright field)
0080 if isfield(CONST, 'imAlign') && isfield(CONST.imAlign, 'AlignChannel') && ...
0081         any(CONST.imAlign.AlignChannel == nc)
0082     nc = [CONST.imAlign.AlignChannel, ...
0083         nc(nc~=CONST.imAlign.AlignChannel)];
0084 else
0085     CONST.imAlign.AlignChannel = 1;
0086     CONST.imAlign.medFilt      = false;
0087     CONST.imAlign.AlignToFirst = false;
0088 end
0089 
0090 
0091 if ~exist( 'targetd', 'var' ) || isempty( targetd )
0092     targetd = [dirname_,'align',filesep];
0093 else
0094      targetd = fixDir(targetd);
0095     
0096 end
0097 
0098 
0099 mkdir(targetd);
0100 num_xy = numel(nxy);
0101 
0102 if (workers>0) && (num_xy>1)
0103     SHOW_FLAG = false;
0104 else
0105     SHOW_FLAG = true;
0106     workers = 0;
0107 end
0108 
0109 
0110 crop_box = cell(1, num_xy);
0111 
0112 % parallelized alignment for each xy
0113 parfor(jj=1:num_xy, workers)
0114     crop_box{jj} = intFrameAlignXY( SHOW_FLAG, nt, nz, nc, nxy(jj), ...
0115         dirname_, targetd, nameInfo, precision, CONST);
0116 end
0117 
0118 end
0119 
0120 end
0121 
0122 
0123 
0124 function crop_box = intFrameAlignXY( SHOW_FLAG, nt, nz, nc, ...
0125     nnxy, dirname, targetd, nameInfo, precision, CONST)
0126 % intFrameAlignXY : Frame alignment for one xy directory
0127 %
0128 % INPUT :
0129 %       SHOW_FLAG : to display the waitbar
0130 %       nt : array with numbers of time frames
0131 %       nz : array with numbers of z frames
0132 %       nc : array with numbers of channels
0133 %       nnxy : xy directory number
0134 %       dirname : dirname where images are at
0135 %       targetd : align directory where aligned images are temporarily placed
0136 %       nameInfo : name information
0137 %       precision : images registered to within 1/precision of a pixel
0138 %       CONST : segmentation constants
0139 %
0140 % OUTPUT :
0141 %       crop_box : information about alignement
0142 %
0143 
0144 FOCUS_NUM_LIM = 0;
0145 ERR_LIM = 1000000;
0146 
0147 if SHOW_FLAG
0148     h = waitbar(0,'aligning frames: ');
0149     cleanup = onCleanup( @()( delete( h ) ) );
0150 end
0151 
0152 verbose = CONST.parallel.verbose;
0153 
0154 nnt = numel( nt );
0155 OutArray   = zeros( nnt, 2);
0156 FocusArray = zeros( nnt, 1);
0157 
0158 countt   = 0;
0159 nz = sort(nz);
0160 
0161 nnz = numel(nz);
0162 nnz2 = ceil(nnz/2); % half of z axis
0163 nz = [nz(nnz2), nz([1:(nnz2-1),(nnz2+1):nnz])];
0164 initFlag = false; % sets the second image the first time through
0165 % or if the focus or error are not within the limits
0166 
0167 % computing the alignment values
0168 for it = nt;
0169     if SHOW_FLAG
0170         waitbar(countt/nnt,h);
0171     end
0172     
0173     countt = countt+1;
0174     
0175     for iz = nz
0176         for ic = nc
0177             nameInfo.npos(:,1) = [it; ic; nnxy; iz];
0178             in_name =  [dirname, MakeFileName(nameInfo)];
0179             if verbose
0180                 disp(['trackOptiAlignPad : Image name: ',in_name]);
0181             end
0182             im = imread(in_name);
0183             
0184             if numel(size(im)) > 2
0185                 disp('Images are in color - attempting to convert to monochromatic.');
0186                 im = convertToMonochromatic(im);
0187             end
0188             
0189             if (ic == nc(1)) && (iz == nnz2 | iz == -1)
0190                 % align phase image at half the z axis
0191                 if ~initFlag % first time through, or high error
0192                     phaseBef = im;
0193                 end
0194                 
0195                 % applies a median filter to the phase image
0196                 if isfield(CONST.imAlign,'medFilt') && CONST.imAlign.medFilt
0197                     im_     = medfilt2( im, [3,3], 'symmetric' );
0198                     phaseBef_ = medfilt2( phaseBef, [3,3], 'symmetric' );
0199                 else
0200                     im_     = im;
0201                     phaseBef_ = phaseBef;
0202                 end
0203                 
0204                 [out,errNum,focusNum] = intAlignIm(im_, phaseBef_, precision );
0205                 if verbose
0206                     disp(['focusNum: ',num2str(focusNum),' errNum: ',num2str(errNum)]);
0207                 end
0208                 FOCUS_FLAG = (focusNum > FOCUS_NUM_LIM) & (errNum < ERR_LIM);
0209                 
0210                 if FOCUS_FLAG % focused, with low alignment error
0211                     initFlag = true;
0212                     OutArray(countt,:) = out(3:4);
0213                     FocusArray(countt) = true;
0214                 end
0215             end
0216             
0217             im = intShiftIm(im, out);
0218             
0219             if SHOW_FLAG
0220                 if ic == nc(1) % phase image
0221                     if ic>0
0222                         try
0223                         figure(ic)
0224                         clf;
0225                         imshow( cat( 3, ag(im), ...
0226                             double(FOCUS_FLAG)*(0.5*ag(im)+0.5*ag(phaseBef)), ...
0227                             double(FOCUS_FLAG)*ag(phaseBef) ));
0228                         end
0229                     end
0230                 else % fluorescence channel
0231                     backer = .8*ag(phaseBef);
0232                     backer0 = backer;
0233                     fluor  = ag(uint16(im-median(im(:))));
0234                     
0235                     figure(ic)
0236                     clf;
0237                     if mod(ic,2) == 0
0238                         imshow( cat( 3, backer, backer0+fluor, backer0));
0239                     else
0240                         imshow( cat( 3, backer+fluor, backer0, backer0));
0241                     end
0242                 end
0243                 drawnow;
0244             end
0245             
0246             if FOCUS_FLAG
0247                 if ic == nc(1) && ~CONST.imAlign.AlignToFirst
0248                     phaseBef = im; % set the previous image to current
0249                 end
0250                 out_name = [targetd, MakeFileName(nameInfo)];
0251             else % non focused or high alignment error
0252                 disp( ['Skipping frame: ', in_name] );
0253             end
0254         end
0255     end
0256 end
0257 
0258 % instead of cropping the image, we add a pad region outside the image.
0259 
0260 ss = size(phaseBef);
0261 
0262 maxy = ceil( max(-OutArray(:,1)));
0263 maxx = ceil( max(-OutArray(:,2)));
0264 miny = floor(min(-OutArray(:,1)));
0265 minx = floor(min(-OutArray(:,2)));
0266 
0267 sspad = [ss(1)+maxy-miny,ss(2)+maxx-minx];
0268 countt   = 0;
0269 initFlag = false;
0270 
0271 crop_box = [OutArray,OutArray];
0272 crop_box(:,1) = -crop_box(:,1) + 1-miny;
0273 crop_box(:,2) = -crop_box(:,2) + 1-minx;
0274 crop_box(:,3) = crop_box(:,1) + ss(1);
0275 crop_box(:,4) = crop_box(:,2) + ss(2);
0276 
0277 for it = nt;
0278     if SHOW_FLAG
0279         waitbar(countt/nnt,h);
0280     end
0281     
0282     countt = countt+1;
0283     
0284     for iz = nz % go through z axis
0285         for ic = nc % go through each channel
0286             nameInfo.npos(:,1) = [it; ic; nnxy; iz];
0287             in_name =  [dirname, MakeFileName(nameInfo)];
0288             if verbose
0289                 disp(['Image name: ',in_name]);
0290             end
0291             im_ = imread( in_name );
0292             if numel(size(im_)) > 2
0293                 disp('Images are in color - attempting to convert to monochromatic.');
0294                 im_ = convertToMonochromatic(im_);
0295             end
0296             im = zeros(sspad, class(im_) ) + mean( im_(:));
0297             im((1-miny):(ss(1)-miny),(1-minx):(ss(2)-minx)) = im_;
0298             out = [0,0,OutArray(countt,:)];
0299             
0300             im = intShiftIm(im, out + CONST.imAlign.out{ic} - ...
0301                 CONST.imAlign.out{1});
0302             
0303             if ~initFlag % first time through, set previous image
0304                 phaseBef = im;
0305             end
0306             
0307             if SHOW_FLAG
0308                 if ic == nc(1)
0309                     figure(ic)
0310                     clf;
0311                     imshow(cat( 3, ag(im), ...
0312                         double(FOCUS_FLAG)*(0.5*ag(im)+0.5*ag(phaseBef)), ...
0313                         double(FOCUS_FLAG)*ag(phaseBef) ));
0314                 else
0315                     backer = .8*ag(phaseBef);
0316                     backer0 = backer;
0317                     fluor  = ag(uint16(im-median(im(:))));
0318                     
0319                     figure(ic)
0320                     clf;
0321                     if mod(ic,2) == 0
0322                         imshow( cat( 3, backer, backer0+fluor, backer0));
0323                     else
0324                         imshow( cat( 3, backer+fluor, backer0, backer0));
0325                     end
0326                 end
0327                 drawnow;
0328                 hold on;
0329                 xxxx = [crop_box(countt,2),crop_box(countt,2),...
0330                     crop_box(countt,4),crop_box(countt,4),...
0331                     crop_box(countt,2)];
0332                 yyyy = [crop_box(countt,1),crop_box(countt,3),...
0333                     crop_box(countt,3), crop_box(countt,1),...
0334                     crop_box(countt,1)];
0335                 plot( xxxx,yyyy,'y');
0336             end
0337             
0338             if FocusArray(countt) % image focused, with low error
0339                 % save shifted image to target directory
0340                 out_name = [targetd, MakeFileName(nameInfo)];
0341                 imwrite(uint16(im), out_name ,'tif','Compression', 'none');
0342                 initFlag = true;
0343                 
0344                 if ic == nc(1)
0345                     phaseBef = im;
0346                 end
0347             else
0348                 disp(['Skipping frame: ', in_name]);
0349             end
0350         end
0351     end
0352     
0353 end
0354 
0355 % done shifting the frames.
0356 if SHOW_FLAG
0357     close(h);
0358 end
0359 
0360 end