CellsortPCA

 [mixedsig, mixedfilters, CovEvals, covtrace, movm, movtm] = CellsortPCA(fn, flims, nPCs, dsamp, outputdir, badframes)

 CELLSORT
 Read TIFF movie data and perform singular-value decomposition (SVD)
 dimensional reduction.

 Inputs:
   fn - movie file name. Must be in TIFF format.
   flims - 2-element vector specifying the endpoints of the range of
   frames to be analyzed. If empty, default is to analyze all movie
   frames.
   nPCs - number of principal components to be returned
   dsamp - optional downsampling factor. If scalar, specifies temporal
   downsampling factor. If two-element vector, entries specify temporal
   and spatial downsampling, respectively.
   outputdir - directory in which to store output .mat files
   badframes - optional list of indices of movie frames to be excluded
   from analysis

 Outputs:
   mixedsig - N x T matrix of N temporal signal mixtures sampled at T
   points.
   mixedfilters - N x X x Y array of N spatial signal mixtures sampled at
   X x Y spatial points.
   CovEvals - largest eigenvalues of the covariance matrix
   covtrace - trace of covariance matrix, corresponding to the sum of all
   eigenvalues (not just the largest few)
   movm - average of all movie time frames at each pixel
   movtm - average of all movie pixels at each time frame, after
   normalizing each pixel deltaF/F

 Eran Mukamel, Axel Nimmerjahn and Mark Schnitzer, 2009
 Email: eran@post.harvard.edu, mschnitz@stanford.edu

0001 function [mixedsig, mixedfilters, CovEvals, covtrace, movm, ...
0002     movtm] = CellsortPCA(fn, flims, nPCs, dsamp, outputdir, badframes)
0003 % [mixedsig, mixedfilters, CovEvals, covtrace, movm, movtm] = CellsortPCA(fn, flims, nPCs, dsamp, outputdir, badframes)
0004 %
0005 % CELLSORT
0006 % Read TIFF movie data and perform singular-value decomposition (SVD)
0007 % dimensional reduction.
0008 %
0009 % Inputs:
0010 %   fn - movie file name. Must be in TIFF format.
0011 %   flims - 2-element vector specifying the endpoints of the range of
0012 %   frames to be analyzed. If empty, default is to analyze all movie
0013 %   frames.
0014 %   nPCs - number of principal components to be returned
0015 %   dsamp - optional downsampling factor. If scalar, specifies temporal
0016 %   downsampling factor. If two-element vector, entries specify temporal
0017 %   and spatial downsampling, respectively.
0018 %   outputdir - directory in which to store output .mat files
0019 %   badframes - optional list of indices of movie frames to be excluded
0020 %   from analysis
0021 %
0022 % Outputs:
0023 %   mixedsig - N x T matrix of N temporal signal mixtures sampled at T
0024 %   points.
0025 %   mixedfilters - N x X x Y array of N spatial signal mixtures sampled at
0026 %   X x Y spatial points.
0027 %   CovEvals - largest eigenvalues of the covariance matrix
0028 %   covtrace - trace of covariance matrix, corresponding to the sum of all
0029 %   eigenvalues (not just the largest few)
0030 %   movm - average of all movie time frames at each pixel
0031 %   movtm - average of all movie pixels at each time frame, after
0032 %   normalizing each pixel deltaF/F
0033 %
0034 % Eran Mukamel, Axel Nimmerjahn and Mark Schnitzer, 2009
0035 % Email: eran@post.harvard.edu, mschnitz@stanford.edu
0036 
0037 tic
0038 fprintf('-------------- CellsortPCA %s: %s -------------- \n', date, fn)
0039 
0040 %-----------------------
0041 % Check inputs
0042 if isempty(dir(fn))
0043     error('Invalid input file name.')
0044 end
0045 if (nargin<2)||(isempty(flims))
0046     nt_full = tiff_frames(fn);
0047     flims = [1,nt_full];
0048 end
0049 
0050 useframes = setdiff((flims(1):flims(2)), badframes);
0051 nt = length(useframes);
0052 
0053 if nargin<3 || isempty(nPCs)
0054     nPCs = min(150, nt);
0055 end
0056 if nargin<4 || isempty(dsamp)
0057     dsamp = [1,1];
0058 end
0059 if nargin<5 || isempty(outputdir)
0060     outputdir = [pwd,'/cellsort_preprocessed_data/'];
0061 end
0062 if nargin<6
0063     badframes = [];
0064 end
0065 if isempty(dir(outputdir))
0066     mkdir(pwd, '/cellsort_preprocessed_data/')
0067 end
0068 if outputdir(end)~='/';
0069     outputdir = [outputdir, '/'];
0070 end
0071 
0072 % Downsampling
0073 if length(dsamp)==1
0074     dsamp_time = dsamp(1);
0075     dsamp_space = 1;
0076 else
0077     dsamp_time = dsamp(1);
0078     dsamp_space = dsamp(2); % Spatial downsample
0079 end
0080 
0081 [fpath, fname] = fileparts(fn);
0082 if isempty(badframes)
0083     fnmat = [outputdir, fname, '_',num2str(flims(1)),',',num2str(flims(2)), '_', date,'.mat'];
0084 else
0085     fnmat = [outputdir, fname, '_',num2str(flims(1)),',',num2str(flims(2)),'_selframes_', date,'.mat'];
0086 end
0087 if ~isempty(dir(fnmat))
0088     fprintf('CELLSORT: Movie %s already processed;', ...
0089         fn)
0090     forceload = input(' Re-load data? [0-no/1-yes] ');
0091     if isempty(forceload) || forceload==0
0092         load(fnmat)
0093         return
0094     end
0095 end
0096 
0097 fncovmat = [outputdir, fname, '_cov_', num2str(flims(1)), ',', num2str(flims(2)), '_', date,'.mat'];
0098 
0099 [pixw,pixh] = size(imread(fn,1));
0100 npix = pixw*pixh;
0101 
0102 fprintf('   %d pixels x %d time frames;', npix, nt)
0103 if nt<npix
0104     fprintf(' using temporal covariance matrix.\n')
0105 else
0106     fprintf(' using spatial covariance matrix.\n')
0107 end
0108 
0109 % Create covariance matrix
0110 if nt < npix
0111     [covmat, mov, movm, movtm] = create_tcov(fn, pixw, pixh, useframes, nt, dsamp);
0112 else
0113     [covmat, mov, movm, movtm] = create_xcov(fn, pixw, pixh, useframes, nt, dsamp);
0114 end
0115 covtrace = trace(covmat) / npix;
0116 movm = reshape(movm, pixw, pixh);
0117 
0118 if nt < npix
0119     % Perform SVD on temporal covariance
0120     [mixedsig, CovEvals, percentvar] = cellsort_svd(covmat, nPCs, nt, npix);
0121 
0122     % Load the other set of principal components
0123     [mixedfilters] = reload_moviedata(pixw*pixh, mov, mixedsig, CovEvals);
0124 else
0125     % Perform SVD on spatial components
0126     [mixedfilters, CovEvals, percentvar] = cellsort_svd(covmat, nPCs, nt, npix);
0127 
0128     % Load the other set of principal components
0129     [mixedsig] = reload_moviedata(nt, mov', mixedfilters, CovEvals);
0130 end
0131 mixedfilters = reshape(mixedfilters, pixw,pixh,nPCs);
0132 
0133 firstframe_full = imread(fn,1);
0134 firstframe = firstframe_full;
0135 if dsamp_space>1
0136     firstframe = imresize(firstframe, size(mov(:,:,1)),'bilinear');
0137 end
0138 
0139 %------------
0140 % Save the output data
0141 save(fnmat,'mixedfilters','CovEvals','mixedsig', ...
0142     'movm','movtm','covtrace')
0143 fprintf(' CellsortPCA: saving data and exiting; ')
0144 toc
0145 
0146     function [covmat, mov, movm, movtm] = create_xcov(fn, pixw, pixh, useframes, nt, dsamp)
0147         %-----------------------
0148         % Load movie data to compute the spatial covariance matrix
0149 
0150         npix = pixw*pixh;
0151 
0152         % Downsampling
0153         if length(dsamp)==1
0154             dsamp_time = dsamp(1);
0155             dsamp_space = 1;
0156         else
0157             dsamp_time = dsamp(1);
0158             dsamp_space = dsamp(2); % Spatial downsample
0159         end
0160 
0161         if (dsamp_space==1)
0162             mov = zeros(pixw, pixh, nt);
0163             for jjind=1:length(useframes)
0164                 jj = useframes(jjind);
0165                 mov(:,:,jjind) = imread(fn,jj);
0166                 if mod(jjind,500)==1
0167                     fprintf(' Read frame %4.0f out of %4.0f; ', jjind, nt)
0168                     toc
0169                 end
0170             end
0171         else
0172             [pixw_dsamp,pixh_dsamp] = size(imresize( imread(fn,1), 1/dsamp_space, 'bilinear' ));
0173             mov = zeros(pixw_dsamp, pixh_dsamp, nt);
0174             for jjind=1:length(useframes)
0175                 jj = useframes(jjind);
0176                 mov(:,:,jjind) = imresize( imread(fn,jj), 1/dsamp_space, 'bilinear' );
0177                 if mod(jjind,500)==1
0178                     fprintf(' Read frame %4.0f out of %4.0f; ', jjind, nt)
0179                     toc
0180                 end
0181             end
0182         end
0183 
0184         fprintf(' Read frame %4.0f out of %4.0f; ', jjind, nt)
0185         toc
0186         mov = reshape(mov, npix, nt);
0187 
0188         % DFoF normalization of each pixel
0189         movm = mean(mov,2); % Average over time
0190         movmzero = (movm==0);
0191         movm(movmzero) = 1;
0192         mov = mov ./ (movm * ones(1,nt)) - 1; % Compute Delta F/F
0193         mov(movmzero, :) = 0;
0194 
0195         if dsamp_time>1
0196             mov = filter(ones(dsamp,1)/dsamp, 1, mov, [], 2);
0197             mov = downsample(mov', dsamp)';
0198         end
0199 
0200         movtm = mean(mov,2); % Average over space
0201         clear movmzeros
0202 
0203         c1 = (mov*mov')/size(mov,2);
0204         toc
0205         covmat = c1 - movtm*movtm';
0206         clear c1
0207     end
0208 
0209     function [covmat, mov, movm, movtm] = create_tcov(fn, pixw, pixh, useframes, nt, dsamp)
0210         %-----------------------
0211         % Load movie data to compute the temporal covariance matrix
0212         npix = pixw*pixh;
0213 
0214         % Downsampling
0215         if length(dsamp)==1
0216             dsamp_time = dsamp(1);
0217             dsamp_space = 1;
0218         else
0219             dsamp_time = dsamp(1);
0220             dsamp_space = dsamp(2); % Spatial downsample
0221         end
0222 
0223         if (dsamp_space==1)
0224             mov = zeros(pixw, pixh, nt);
0225             for jjind=1:length(useframes)
0226                 jj = useframes(jjind);
0227                 mov(:,:,jjind) = imread(fn,jj);
0228                 if mod(jjind,500)==1
0229                     fprintf(' Read frame %4.0f out of %4.0f; ', jjind, nt)
0230                     toc
0231                 end
0232             end
0233         else
0234             [pixw_dsamp,pixh_dsamp] = size(imresize( imread(fn,1), 1/dsamp_space, 'bilinear' ));
0235             mov = zeros(pixw_dsamp, pixh_dsamp, nt);
0236             for jjind=1:length(useframes)
0237                 jj = useframes(jjind);
0238                 mov(:,:,jjind) = imresize( imread(fn,jj), 1/dsamp_space, 'bilinear' );
0239                 if mod(jjind,500)==1
0240                     fprintf(' Read frame %4.0f out of %4.0f; ', jjind, nt)
0241                     toc
0242                 end
0243             end
0244         end
0245 
0246         fprintf(' Read frame %4.0f out of %4.0f; ', jjind, nt)
0247         toc
0248         mov = reshape(mov, npix, nt);
0249 
0250         % DFoF normalization of each pixel
0251         movm = mean(mov,2); % Average over time
0252         movmzero = (movm==0); % Avoid dividing by zero
0253         movm(movmzero) = 1;
0254         mov = mov ./ (movm * ones(1,nt)) - 1;
0255         mov(movmzero, :) = 0;
0256 
0257         if dsamp_time>1
0258             mov = filter(ones(dsamp,1)/dsamp, 1, mov, [], 2);
0259             mov = downsample(mov', dsamp)';
0260         end
0261 
0262         c1 = (mov'*mov)/npix;
0263         movtm = mean(mov,1); % Average over space
0264         covmat = c1 - movtm'*movtm;
0265         clear c1
0266     end
0267 
0268     function [mixedsig, CovEvals, percentvar] = cellsort_svd(covmat, nPCs, nt, npix)
0269         %-----------------------
0270         % Perform SVD
0271 
0272         covtrace = trace(covmat) / npix;
0273 
0274         opts.disp = 0;
0275         opts.issym = 'true';
0276         if nPCs<size(covmat,1)
0277             [mixedsig, CovEvals] = eigs(covmat, nPCs, 'LM', opts);  % pca_mixedsig are the temporal signals, mixedsig
0278         else
0279             [mixedsig, CovEvals] = eig(covmat);
0280             CovEvals = diag( sort(diag(CovEvals), 'descend'));
0281             nPCs = size(CovEvals,1);
0282         end
0283         CovEvals = diag(CovEvals);
0284         if nnz(CovEvals<=0)
0285             nPCs = nPCs - nnz(CovEvals<=0);
0286             fprintf(['Throwing out ',num2str(nnz(CovEvals<0)),' negative eigenvalues; new # of PCs = ',num2str(nPCs),'. \n']);
0287             mixedsig = mixedsig(:,CovEvals>0);
0288             CovEvals = CovEvals(CovEvals>0);
0289         end
0290 
0291         mixedsig = mixedsig' * nt;
0292         CovEvals = CovEvals / npix;
0293 
0294         percentvar = 100*sum(CovEvals)/covtrace;
0295         fprintf([' First ',num2str(nPCs),' PCs contain ',num2str(percentvar,3),'%% of the variance.\n'])
0296     end
0297 
0298     function [mixedfilters] = reload_moviedata(npix, mov, mixedsig, CovEvals)
0299         %-----------------------
0300         % Re-load movie data
0301         nPCs = size(mixedsig,1);
0302 
0303         Sinv = inv(diag(CovEvals.^(1/2)));
0304 
0305         movtm = mean(mov,1); % Average over space
0306         movuse = mov - ones(npix,1) * movtm;
0307         mixedfilters = reshape(movuse * mixedsig' * Sinv, npix, nPCs);
0308     end
0309
0310 
0311     function j = tiff_frames(fn)
0312         %
0313         % n = tiff_frames(filename)
0314         %
0315         % Returns the number of slices in a TIFF stack.
0316         %
0317         % Modified April 9, 2013 for compatibility with MATLAB 2012b
0318 
0319         j = length(imfinfo(fn));
0320     end
0321 end