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 
0038 %
0039 % June-August 2016 (Cellsort 1.41): Fixed an inconsistency between create_tcov and
0040 % creat_xcov. In the current version, both methods agree with the previous
0041 % version of create_tcov.
0042 %
0043 
0044 tic
0045 fprintf('-------------- CellsortPCA %s: %s -------------- \n', date, fn)
0046 
0047 %-----------------------
0048 % Check inputs
0049 if isempty(dir(fn))
0050     error('Invalid input file name.')
0051 end
0052 if (nargin<2)||(isempty(flims))
0053     nt_full = tiff_frames(fn);
0054     flims = [1,nt_full];
0055 end
0056 
0057 useframes = setdiff((flims(1):flims(2)), badframes);
0058 nt = length(useframes);
0059 
0060 if nargin<3 || isempty(nPCs)
0061     nPCs = min(150, nt);
0062 end
0063 if nargin<4 || isempty(dsamp)
0064     dsamp = [1,1];
0065 end
0066 if nargin<5 || isempty(outputdir)
0067     outputdir = [pwd,'/cellsort_preprocessed_data/'];
0068 end
0069 if nargin<6
0070     badframes = [];
0071 end
0072 if isempty(dir(outputdir))
0073     mkdir(pwd, '/cellsort_preprocessed_data/')
0074 end
0075 if outputdir(end)~='/';
0076     outputdir = [outputdir, '/'];
0077 end
0078 
0079 [~, fname] = fileparts(fn);
0080 if isempty(badframes)
0081     fnmat = [outputdir, fname, '_',num2str(flims(1)),',',num2str(flims(2)), '_', date,'.mat'];
0082 else
0083     fnmat = [outputdir, fname, '_',num2str(flims(1)),',',num2str(flims(2)),'_selframes_', date,'.mat'];
0084 end
0085 if ~isempty(dir(fnmat))
0086     fprintf('CELLSORT: Movie %s already processed;', ...
0087         fn)
0088     forceload = input(' Re-load data? [0-no/1-yes] ');
0089     if isempty(forceload) || forceload==0
0090         load(fnmat)
0091         return
0092     end
0093 end
0094 
0095 [pixw,pixh] = size(imread(fn,1));
0096 npix = pixw*pixh;
0097 
0098 fprintf('   %d pixels x %d time frames;', npix, nt)
0099 
0100 % Create covariance matrix
0101 if nt < npix 
0102     fprintf(' using temporal covariance matrix.\n')
0103     [covmat, mov, movm, movtm] = create_tcov(fn, pixw, pixh, useframes, nt, dsamp);
0104 else
0105     fprintf(' using spatial covariance matrix.\n')
0106     [covmat, mov, movm, movtm] = create_xcov(fn, pixw, pixh, useframes, nt, dsamp);
0107 end
0108 
0109 covtrace = trace(covmat) / npix;
0110 movm = reshape(movm, pixw, pixh);
0111 
0112 if nt < npix 
0113     % Perform SVD on temporal covariance
0114     [mixedsig, CovEvals] = cellsort_svd(covmat, nPCs, nt, npix);
0115 
0116     % Load the other set of principal components
0117     [mixedfilters] = reload_moviedata(pixw*pixh, mov, mixedsig, CovEvals);
0118 else
0119     % Perform SVD on spatial components
0120     [mixedfilters, CovEvals] = cellsort_svd(covmat, nPCs, nt, npix);
0121     mixedfilters = mixedfilters' * npix;
0122     
0123     % Load the other set of principal components
0124     [mixedsig] = reload_moviedata(nt, mov', mixedfilters', CovEvals);
0125     mixedsig = mixedsig' / npix^2;
0126 end
0127 mixedfilters = reshape(mixedfilters, pixw,pixh,nPCs);
0128 
0129 %------------
0130 % Save the output data
0131 save(fnmat,'mixedfilters','CovEvals','mixedsig', ...
0132     'movm','movtm','covtrace')
0133 fprintf(' CellsortPCA: saving data and exiting; ')
0134 toc
0135 
0136     function [covmat, mov, movm, movtm] = create_xcov(fn, pixw, pixh, useframes, nt, dsamp)
0137         %-----------------------
0138         % Load movie data to compute the spatial covariance matrix
0139 
0140         npix1 = pixw*pixh;
0141 
0142         % Downsampling
0143         if length(dsamp)==1
0144             dsamp_time = dsamp(1);
0145             dsamp_space = 1;
0146         else
0147             dsamp_time = dsamp(1);
0148             dsamp_space = dsamp(2); % Spatial downsample
0149         end
0150 
0151         if (dsamp_space==1)
0152             mov = zeros(pixw, pixh, nt);
0153             for jjind=1:length(useframes)
0154                 jj = useframes(jjind);
0155                 mov(:,:,jjind) = imread(fn,jj);
0156                 if mod(jjind,500)==1
0157                     fprintf(' Read frame %4.0f out of %4.0f; ', jjind, nt)
0158                     toc
0159                 end
0160             end
0161         else
0162             [pixw_dsamp,pixh_dsamp] = size(imresize( imread(fn,1), 1/dsamp_space, 'bilinear' ));
0163             mov = zeros(pixw_dsamp, pixh_dsamp, nt);
0164             for jjind=1:length(useframes)
0165                 jj = useframes(jjind);
0166                 mov(:,:,jjind) = imresize( imread(fn,jj), 1/dsamp_space, 'bilinear' );
0167                 if mod(jjind,500)==1
0168                     fprintf(' Read frame %4.0f out of %4.0f; ', jjind, nt)
0169                     toc
0170                 end
0171             end
0172         end
0173 
0174         fprintf(' Read frame %4.0f out of %4.0f; ', jjind, nt)
0175         toc
0176         mov = reshape(mov, npix1, nt);
0177 
0178         % DFoF normalization of each pixel
0179         movm = mean(mov,2); % Average over time
0180         movmzero = (movm==0);
0181         movm(movmzero) = 1;
0182         mov = mov ./ (movm * ones(1,nt)) - 1; % Compute Delta F/F
0183         mov(movmzero, :) = 0;
0184 
0185         if dsamp_time>1
0186             mov = filter(ones(dsamp,1)/dsamp, 1, mov, [], 2);
0187             mov = downsample(mov', dsamp)';
0188         end
0189 
0190         movtm = mean(mov,1); % Average over space
0191         mov = mov - ones(size(mov,1),1)*movtm; 
0192 
0193         covmat = (mov*mov')/size(mov,2);
0194         covmat = covmat * size(mov,2)/size(mov,1); % Rescale to gree with create_tcov
0195         toc
0196     end
0197 
0198     function [covmat, mov, movm, movtm] = create_tcov(fn, pixw, pixh, useframes, nt, dsamp)
0199         %-----------------------
0200         % Load movie data to compute the temporal covariance matrix
0201         npix1 = pixw*pixh;
0202 
0203         % Downsampling
0204         if length(dsamp)==1
0205             dsamp_time = dsamp(1);
0206             dsamp_space = 1;
0207         else
0208             dsamp_time = dsamp(1);
0209             dsamp_space = dsamp(2); % Spatial downsample
0210         end
0211 
0212         if (dsamp_space==1)
0213             mov = zeros(pixw, pixh, nt);
0214             for jjind=1:length(useframes)
0215                 jj = useframes(jjind);
0216                 mov(:,:,jjind) = imread(fn,jj);
0217                 if mod(jjind,500)==1
0218                     fprintf(' Read frame %4.0f out of %4.0f; ', jjind, nt)
0219                     toc
0220                 end
0221             end
0222         else
0223             [pixw_dsamp,pixh_dsamp] = size(imresize( imread(fn,1), 1/dsamp_space, 'bilinear' ));
0224             mov = zeros(pixw_dsamp, pixh_dsamp, nt);
0225             for jjind=1:length(useframes)
0226                 jj = useframes(jjind);
0227                 mov(:,:,jjind) = imresize( imread(fn,jj), 1/dsamp_space, 'bilinear' );
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         end
0234 
0235         fprintf(' Read frame %4.0f out of %4.0f; ', jjind, nt)
0236         toc
0237         mov = reshape(mov, npix1, nt);
0238 
0239         % DFoF normalization of each pixel
0240         movm = mean(mov,2); % Average over time
0241         movmzero = (movm==0); % Avoid dividing by zero
0242         movm(movmzero) = 1;
0243         mov = mov ./ (movm * ones(1,nt)) - 1;
0244         mov(movmzero, :) = 0;
0245 
0246         if dsamp_time>1
0247             mov = filter(ones(dsamp,1)/dsamp, 1, mov, [], 2);
0248             mov = downsample(mov', dsamp)';
0249         end
0250 
0251         c1 = (mov'*mov)/npix1;
0252         movtm = mean(mov,1); % Average over space
0253         covmat = c1 - movtm'*movtm;
0254         clear c1
0255     end
0256 
0257     function [mixedsig, CovEvals, percentvar] = cellsort_svd(covmat, nPCs, nt, npix1)
0258         %-----------------------
0259         % Perform SVD
0260 
0261         covtrace1 = trace(covmat) / npix1;
0262 
0263         opts.disp = 0;
0264         opts.issym = 'true';
0265         if nPCs<size(covmat,1)
0266             [mixedsig, CovEvals] = eigs(covmat, nPCs, 'LM', opts);  % pca_mixedsig are the temporal signals, mixedsig
0267         else
0268             [mixedsig, CovEvals] = eig(covmat);
0269             CovEvals = diag( sort(diag(CovEvals), 'descend'));
0270             nPCs = size(CovEvals,1);
0271         end
0272         CovEvals = diag(CovEvals);
0273         if nnz(CovEvals<=0)
0274             nPCs = nPCs - nnz(CovEvals<=0);
0275             fprintf(['Throwing out ',num2str(nnz(CovEvals<0)),' negative eigenvalues; new # of PCs = ',num2str(nPCs),'. \n']);
0276             mixedsig = mixedsig(:,CovEvals>0);
0277             CovEvals = CovEvals(CovEvals>0);
0278         end
0279 
0280         mixedsig = mixedsig' * nt;
0281         CovEvals = CovEvals / npix1;
0282 
0283         percentvar = 100*sum(CovEvals)/covtrace1;
0284         fprintf([' First ',num2str(nPCs),' PCs contain ',num2str(percentvar,3),'%% of the variance.\n'])
0285     end
0286 
0287     function [mixedfilters] = reload_moviedata(npix1, mov, mixedsig, CovEvals)
0288         %-----------------------
0289         % Re-load movie data
0290         nPCs1 = size(mixedsig,1);
0291 
0292         Sinv = inv(diag(CovEvals.^(1/2)));
0293 
0294         movtm1 = mean(mov,1); % Average over space
0295         movuse = mov - ones(npix1,1) * movtm1;
0296         mixedfilters = reshape(movuse * mixedsig' * Sinv, npix1, nPCs1);
0297     end
0298 
0299     function j = tiff_frames(fn)
0300         %
0301         % n = tiff_frames(filename)
0302         %
0303         % Returns the number of slices in a TIFF stack.
0304         %
0305         % Modified April 9, 2013 for compatibility with MATLAB 2012b
0306 
0307         j = length(imfinfo(fn));
0308     end
0309 end