changed structure for script for rawData processing, added loading of…

… single files, binning, shift correction

CellSortExportTifStack.m

@@ -0,0 +1,5 @@
+function CellSortExportTifStack(data, fn)
+for i=1:size(data,3)
+    imwrite(data(:,:,i),fn, 'WriteMode', 'append');
+end
+end

CellSortLoadSingleFiles.m

@@ -0,0 +1,24 @@
+function [data_g, data_r] = CellSortLoadSingleFiles(green, red)
+fileListGreen = dir(green);
+fileListRed = dir(red);
+
+numFilesGreen = size(fileListGreen,1);
+numFilesRed = size(fileListRed,1);
+if numFilesGreen ~= numFilesRed
+disp('Green and Red Channel have different number of Images, take the smaller value..');
+end
+numFiles = min([numFilesRed, numFilesGreen]);
+
+%read first frame to get resolution fily type etc..
+temp = readTiff(fileListGreen(1).name);
+data_g = zeros(size(temp,1), size(temp,2), numFiles, class(temp));
+temp = readTiff(fileListRed(1).name);
+data_r = zeros(size(temp,1), size(temp,2), numFiles, class(temp));
+
+for i=1:numFiles
+    temp = readTiff(fileListGreen(i).name);
+    data_g(:,:,i) = temp;
+    temp = readTiff(fileListRed(i).name);
+    data_r(:,:,i) = temp;
+end
+end

CellSortShiftCorrection.m

@@ -0,0 +1,13 @@
+function [data_g, data_r] = CellSortShiftCorrection(data_g, data_r)
+%% shiftcorrection
+%import green channel
+disp('Performing Shift Correction..');
+%shift correction
+firstRedFrame = data_r(:,:,1);
+parfor r = 2:size(data_g, 3)
+    shiftVector = findshift(firstRedFrame, data_r(:,:,r), 'iter');
+    data_r(:,:,r)  = shift(data_r(:,:,r),shiftVector);
+    data_g(:,:,r)  = shift(data_g(:,:,r),shiftVector);
+end
+disp('finished');
+end

CellSortTemporalBinning.m

@@ -0,0 +1,4 @@
+function [data_g, data_r] = CellSortTemporalBinning(data_g, data_r, binningFactor)
+data_g = binning(data_g, bininngFactor, 'mean');
+data_r = binning(data_r, bininngFactor, 'mean');
+end

binning.m

@@ -0,0 +1,149 @@
+function A = binning(A,bin,fun_str)
+%
+% A=binning(A,bin,fun)
+%
+% This function rearranges the array A according to the specifications
+% of the vector bin. In the first dimension of the generated temporary 
+% array the bin elements from all dimension are arranged. By applying 
+% functions like 'sum', 'mean' or 'std' to the first dimension, a real
+% binning ('sum'), a averaging ('mean') and error calculations can be
+% done. 
+%
+% A=binning(A,bin) The input array 'A' is binned by the vector 'bin'.
+% If 'bin' is a scalar, 'bin' elements in each dimenion of 'A' are
+% summarized.
+%
+% Note: The input array 'A' is cut to the maxial integer factor of
+% 'bin'. 
+%
+% A=binning(A,bin,fun) The function 'fun' is applied to the first
+% dimension of the rearranged array. Beside the real binning with the
+% function 'sum' a number of additional functions can be applied. Up
+% to now the following functions are implemented:
+% 'sum', 'mean', 'prod'
+% 
+% Special comment:
+% The functionality can be easily extended by implementing additional
+% functions, or due to the time consuming permutations by the
+% combination of appropriate functions like 'mean' and 'std'.
+%
+% Author: A. Zeug (last modified 19.07.2007)
+% Version: 2.0.0
+
+% comments to: azeug@gwdg.de
+
+tic; t1=toc;
+%% Input validation
+if nargin < 3
+  fun_str='sum';
+  fun_type = 1;
+elseif sum(strcmp(fun_str,{'sum';'mean';'prod'}))>0
+  fun_type = 1;
+elseif sum(strcmp(fun_str,{'std';'var'}))>0
+  fun_type = 2;
+else
+  error('Function not implemented yet')
+end
+fun=str2func(fun_str); % defines function type
+
+if nargin < 1
+  help binning
+  return
+elseif nargin < 2
+  error('Please specify binning vector.')
+end
+
+%% bin vector dimension
+trans = false;
+dims  = ndims(A);
+if isvector(A)
+  dims = 2;
+  if size(A,2)==1; trans = true; end
+end
+if length(bin)==1; 
+  bin = bin * ones(1, dims); 
+elseif length(bin)~=dims;
+  error('Wrong length of binning vector.')
+end
+
+%% resizing of A to fit bin
+sz=size(A);
+if dims==1; sz=length(A); end
+bin=min([sz;bin]);
+for i=1:dims
+  index{i}=1:(bin(i)*floor(sz(i)/bin(i)));
+end
+A = A(index{:});
+
+%% reshaping A
+sz = size(A); if dims==1; sz = length(A); end
+v = []; sz_new = sz./bin;
+for i=1:dims
+  v=[v, bin(i), sz_new(i)];
+end
+A=reshape(A,v);
+
+%% Calculate binning
+if fun_type==1
+  % standard type 'sum', 'mean','prod' 
+  for i=1:dims
+      if isequal(fun, @sum)
+          A=fun(A,(2*i-1),'native');
+      else
+          A=fun(A,(2*i-1));
+      end
+  end
+else
+  % type 'std', 'var'
+  % need to reshape again to calc function at ones
+  v = [];
+  for i=1:dims
+    v=[v, i, dims+i];
+  end
+  A = permute(A, v);
+  A = reshape(A, [prod(bin) sz_new]);
+  try
+    A = fun( double(A) );
+  catch
+    t=toc; fprintf('(%07.3fs) binning: direct calculation failed, start subsequential method (%0.5fs)\n',t,t-t1); t1=t;
+    try
+      A1 = zeros(sz_new,'double');
+      fun1 = @double;
+    catch
+      t=toc; fprintf('(%07.3fs) binning: calculating %s with double precession failed, try single (%0.5fs)\n',t,func2str(fun),t-t1); t1=t;
+      A1 = zeros(sz_new, 'single');
+      fun1 = @single;
+    end
+    h = waitbar(0,'Please wait...','Name','binning: Plese Wait!');
+    % subsequential calculation
+    sub_dim = find(sz_new == max(sz_new));
+    for i = 1:dims; cc{i}=':'; end
+    step = floor(1E6 / (prod(sz_new(sz_new < max(sz_new)))*prod(bin)));
+    for ii = 0:step:sz_new(sub_dim) % subsequentially...
+      waitbar(ii/sz_new(sub_dim),h,['Calculated binning: ' num2str(100*ii/sz_new(sub_dim),'%0.2f') '%']);
+      ii_end = min(sz_new(sub_dim), ii+step) - ii;
+      cc{sub_dim}=(1:ii_end)+ii;
+      A1(cc{:}) = fun(fun1(A(:,cc{:})));
+    end
+    A = A1;
+    delete(h); drawnow;
+    t=toc; fprintf('(%07.3fs) binning: calculation with %s precession susseed (%0.5fs)\n',t,func2str(fun1),t-t1); t1=t;
+  end
+end
+
+%
+%A = reshape(A, sz_new);
+
+if dims==1 
+  A=squeeze(A);
+  if trans==true
+    A=A';
+  end
+else
+  A = reshape(A, sz_new);
+
+end
+
+
+
+

main.m

@@ -1,18 +1,25 @@
-%% filename
-fn_green = 'green_13.tiff';
-fn_red = 'red_13.tiff';
+%% load and prepare all single tiffs
+%importing
+[data_g,data_r] = CellSortLoadSingleFiles('ChanA*.tif', 'ChanB*.tif');
+%% shiftCorrection
+%generates shift corrected files for the green and red channels, needs
+%DIPlib (available for free)
+[data_g, data_r] = CellSortShiftCorrection(data_g,data_r);
+%temporal binning
+binningFactor = [1 1 2]; %x, y, t
+[data_g, data_r] = CellSortTemporalBinning(data_g, data_r, binningFactor);
+%% export
+fileName_g = 'green.tif';
+fileName_r = 'red.tif';
+CellSortExportTifStack(data_g, fileName_g);
+CellSortExportTifStack(data_r, fileName_r);
 %% merge files
 fn_green = 'merged_green.tiff';
 fn_red = 'merged_red.tiff';
 fileList = {'green_13.tiff', 'green_14.tiff','green_15.tiff','green_1.tiff','green_2.tiff','green_3.tiff'};
 CellSortMergeTiff(fileList, fn_green);
 fileList = {'red_13.tiff', 'red_14.tiff','red_15.tiff','red_1.tiff','red_2.tiff','red_3.tiff'};
 CellSortMergeTiff(fileList, fn_red);
-
-%% shiftCorrection
-%generates shift corrected files for the green and red channels, needs
-%DIPlib (available for free)
-[fn_green, fn_red, data_g, data_r] = shiftCorrection(fn_green,fn_red, true);
 %% load files
 %if you haven't used the shiftCorrection you have to load the files extra
 fn_green = 'green_13_corr.tif';

main_old.m

@@ -0,0 +1,134 @@
+%% filename
+fn_green = 'green_13.tiff';
+fn_red = 'red_13.tiff';
+%% merge files
+fn_green = 'merged_green.tiff';
+fn_red = 'merged_red.tiff';
+fileList = {'green_13.tiff', 'green_14.tiff','green_15.tiff','green_1.tiff','green_2.tiff','green_3.tiff'};
+CellSortMergeTiff(fileList, fn_green);
+fileList = {'red_13.tiff', 'red_14.tiff','red_15.tiff','red_1.tiff','red_2.tiff','red_3.tiff'};
+CellSortMergeTiff(fileList, fn_red);
+
+%% shiftCorrection
+%generates shift corrected files for the green and red channels, needs
+%DIPlib (available for free)
+[fn_green, fn_red, data_g, data_r] = shiftCorrection(fn_green,fn_red, true);
+%% load files
+%if you haven't used the shiftCorrection you have to load the files extra
+fn_green = 'green_13_corr.tif';
+fn_red = 'red_13_corr.tif';
+data_g = readTiff(fn_green);
+data_r = readTiff(fn_red);
+%% 1. PCA
+%calculate PCAs
+nPCs = 100;
+flims = [];
+dSamp = [1 1]; %temporal, spatial downsampling
+[mixedsig, mixedfilters, CovEvals, covtrace, movm, movtm] = ...
+    CellsortPCA(fn_green, flims, nPCs, dSamp, [], [], data_g);
+
+%% 2a. Choose PCs
+%choose PCAs, type f/b in command line to go back or forward. type two number
+%to select the range "1 ENTER 25 ENTER" would select the PCAs from 1 to 25
+[PCuse] = CellsortChoosePCs(fn_green, mixedfilters);
+
+%% 2b. Plot PC spectrum
+%optionally, you can display the PC spectrum..
+CellsortPlotPCspectrum(fn_green, CovEvals, PCuse)
+
+%% 3a. ICA
+%calculate ICAs, set mu to a value between 0 and 1, where 1 means pure
+%temporal identification and 0 pure spatial. A value of 0.1-0.2 seems to be
+%best.. make sure that the algorithm converges (see command line). If not
+%run it again
+nIC = length(PCuse);
+mu = 0.2;
+
+[ica_sig, ica_filters, ica_A, numiter] = CellsortICA(mixedsig, mixedfilters, CovEvals, PCuse, mu, nIC);
+
+%% 3b. Plot ICs
+%optionally, plot ICs (you have to use at least 20.. or change the number
+%on the buttom
+tlims = [];
+dt = 0.1;
+
+figure(2)
+CellsortICAplot('series', ica_filters, ica_sig, movm, tlims, dt, [], [], [1:20]);
+
+%% 4a. Segment contiguous regions within ICs
+%search for indepent segments, you have to adjust these values well:
+%smwidth: width of a smoothing filter (1 should be okay)
+%thresh: threshold for normalized standard deviation in ICs (verey
+%        sensitive, standard: 2)
+%arealims: limit for the detected areas in pixeks: [min max] or [min]
+%plotting: 1 = yes or 0 = no.. shall a plot be shown or not?
+
+smwidth = 1;
+thresh = 2;
+arealims = [250];
+plotting = 1;
+
+[ica_segments, segmentlabel, segcentroid] = CellsortSegmentation(ica_filters, smwidth, thresh, arealims, plotting);
+
+%% 4b. CellsortApplyFilter
+%filter signal and extract regions signals
+subtractmean = 1;
+
+cell_sig = CellsortApplyFilter(fn_green, ica_segments, flims, movm, subtractmean, data_g);
+
+%% 5. CellsortFindspikes
+%optionally, find spikes in signal
+deconvtau = 0;
+spike_thresh = 2;
+normalization = 1;
+
+[spmat, spt, spc] = CellsortFindspikes(ica_sig, spike_thresh, dt, deconvtau, normalization);
+
+%% Show results
+%optionally, plot signals
+figure(2)
+CellsortICAplot('series', ica_filters, ica_sig, movm, tlims, dt, 1, 2, [1:20], spt, spc);
+%% prepare georgis data for visualization
+%optionally, if you want to use the alternative evaluation georgi was using, import the
+%mask and run this section
+fn='green_13.tiff';
+info = imfinfo(fn);
+num_images = numel(info);
+data_g = readTiff(fn);
+%make mask logical
+mask1 = logical(mask.*(mask > 0));
+ica_segments = permute(mask, [3 1 2]);
+mask1 = repmat(mask, [1 1 1 size(data_g,3)]);
+%generate cell_sig for every feature
+cellsig = zeros(size(mask1,3), size(data_g,3));
+for i=1:size(mask1,3)
+    cell_sig(i,:) = squeeze(mean(mean(squeeze(mask1(:,:,i,:)).*data_g(:,:,:),1),2));
+end
+clear mask
+%% show signal of cells (cell_sig)
+%caluclate correlation matrix for cell_sig
+[ica_segments, cluster, subPlotCols, subPlotWidth, sortedIdx] = calcCorrEff(cell_sig, ica_segments);
+alternativePlot(ica_segments, cluster, cell_sig, subPlotCols, sortedIdx, data_r);
+
+%% let user deselect and merge rois
+%merge -> change cluster ->recereate plots
+merge = {[1,2,3], [5,6]};
+delete = [8,4,11];
+figure;
+cluster_new = cluster;
+%delete
+for i = 1:size(delete,2)
+    cluster_new{delete(i)} = [];
+end
+%merge
+for i = 1:size(merge, 2)
+    for j = 2:size(merge{i},2)
+        cluster_new{merge{i}(1)} = [cluster_new{merge{i}(1)}, cluster_new{merge{i}(j)}];
+        cluster_new{merge{i}(j)} = [];
+    end
+end
+cluster_new = cluster_new(~cellfun(@isempty, cluster_new));
+subPlotCols = size(cluster_new, 2);
+
+%redraw
+alternativePlot(ica_segments, cluster_new, cell_sig, subPlotCols, sortedIdx, data_r);

readTiff.m

@@ -1,7 +1,15 @@
 function data = readTiff(filename)
 info = imfinfo(filename);
+
+if info.BitDepth == 16 || info.BitDepth == 12
+    filetype = 'uint16';
+elseif info.BitDepth == 8
+    filetype = 'uint8';
+else
+    filetype = 'uint16';
+end
 num_images = numel(info);
-data = zeros(info(1).Width, info(1).Height, num_images);
+data = zeros(info(1).Width, info(1).Height, num_images, filetype);
 t_g = Tiff(filename, 'r');
 for k = 1:num_images
     t_g.setDirectory(k);