From 96941f59624f41064d87d031291f83b94b44d0fa Mon Sep 17 00:00:00 2001 From: Florian Vollrath Date: Fri, 29 Jul 2016 15:25:26 +0200 Subject: [PATCH] added a test folder for the medAlign function --- .../batchMedAlignWithTemporalBinning.m | 98 +++++++++++ medAlign_new/binning.m | 149 ++++++++++++++++ medAlign_new/medAlign.m | 163 ++++++++++++++++++ medAlign_new/testscript.m | 21 +++ 4 files changed, 431 insertions(+) create mode 100644 medAlign_new/batchMedAlignWithTemporalBinning.m create mode 100644 medAlign_new/binning.m create mode 100644 medAlign_new/medAlign.m create mode 100644 medAlign_new/testscript.m diff --git a/medAlign_new/batchMedAlignWithTemporalBinning.m b/medAlign_new/batchMedAlignWithTemporalBinning.m new file mode 100644 index 0000000..cb66f69 --- /dev/null +++ b/medAlign_new/batchMedAlignWithTemporalBinning.m @@ -0,0 +1,98 @@ + +nSubIm = 10; + +% set a query directory +dir_query = '/storage/letz/Data/Data/Leona_Enke/20160726_L59/imaging'; +dir_output = '/storage/letz/Data/Data/Leona_Enke/20160726_L59/imaging'; + +% level one +[dir_list, dir_count] = getDirsFromDir(dir_query); + +d = 1; + +for d = 1 : dir_count + + % sub dir name + sub_dir = [dir_query, filesep, dir_list{d}]; + subdir_list = getDirsFromDir(sub_dir); + disp(d); + + % remove finished + idx_rmv = strncmp('aligned',subdir_list,7); + subdir_list(idx_rmv) = []; + subdir_count = size(subdir_list,1); + + % loop through sub dirs + for k = 10 : subdir_count + + + + % current id + folder_id = sscanf(subdir_list{k},'Untitled%03d'); + + + + + % Compose image stack from filenames + + filePrefixGreen = 'ChanA_0001_0001_0001_'; %Gruen Kanal + filePrefixRed = 'ChanB_0001_0001_0001_'; %Rote Kanal + + m = dir([sub_dir,filesep,subdir_list{k},filesep,'ChanB_0001_0001_0001_*.tif']); + nImages = size(m,1); + + file_image = [sub_dir,filesep,subdir_list{k},filesep,'ChanB_0001_0001_0001_0001.tif']; + im = imread(file_image); + + + %Assembling Red Image stack + allIm = zeros(size(im,1),size(im,2),nImages); + for i = 1:nImages + filename= sprintf([sub_dir,filesep,subdir_list{k},filesep,filePrefixRed,'%04u.tif'],i); + allIm(:,:,i) = imread(filename); + end + + %temporal binning on red image stack + class_red = class(allIm); + allIm = cast(binning(allIm, [1 1 2], 'mean'), class_red); + + %Realignment + [rs, cs] = medAlign(allIm, nSubIm, 50, -inf); + + %Rotation and translation on red channel + filename= sprintf('%s%s%s%sred%03d.tif',dir_output,filesep,dir_list{d},filesep,folder_id); + for i = 1:nImages/2 + registered(:,:,i) = circshift(allIm(:,:,i),[rs(i),cs(i)]); + imwrite(uint16(registered(:,:,i)),filename, 'tif','WriteMode','append'); + end + + + %Assembling Green Image stack + allIm = zeros(size(im,1),size(im,2),nImages); + for i = 1:nImages + filename= sprintf([sub_dir,filesep,subdir_list{k},filesep,filePrefixGreen,'%04u.tif'],i); + allIm(:,:,i) = imread(filename); + end + + %temporal binning on green image stack + class_red = class(allIm); + allIm = cast(binning(allIm, [1 1 2], 'mean'), class_red); + + %Rotation and translation on green channel + filename= sprintf('%s%s%s%sgreen%03d.tif',dir_output,filesep,dir_list{d},filesep,folder_id); + %filename= sprintf([dir_output,filesep,filePrefixGreen,'_cc.tif'],i); + for i = 1:nImages/2 + registered(:,:,i) = circshift(allIm(:,:,i),[rs(i),cs(i)]); + imwrite(uint16(registered(:,:,i)),filename, 'tif','WriteMode','append'); + end + + + + end + + + + + +end + diff --git a/medAlign_new/binning.m b/medAlign_new/binning.m new file mode 100644 index 0000000..adcb241 --- /dev/null +++ b/medAlign_new/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 + + + + \ No newline at end of file diff --git a/medAlign_new/medAlign.m b/medAlign_new/medAlign.m new file mode 100644 index 0000000..0d7442b --- /dev/null +++ b/medAlign_new/medAlign.m @@ -0,0 +1,163 @@ +%MEDALIGN +% Function to register/realign calcium imaging data based +% on control point determination by cross correlation. +% +% Based on Daniel Castaño-Díez medianRegister.m script +% +% Copyright (C) 2016 Friedrich Kretschmer +% Max Planck Institute for Brain Research +% friedrich.kretschmer@brain.mpg.de +% +% [RS, CS] = MEDALIGN(ALLIM, NSUBIM, NPEAKS, LOCALTHRESHOLD); +% +% ALLIM, the imagestack to be alligned +% NSUBIM, the number of images to be computed in one iteration +% higher is better but needs to be adjusted to available memory +% (RAM or GPU) +% NPEAKS, the number of peaks in the cross correlation +% to use for allignement (Daniel used 50) +% THRESHOLD, above which CC maxima need to lie. Use -inf if you want to +% include all values +% +% GPU computation can be activated by setting GPU=TRUE in the +% first line of the function +% +function [rs, cs] = medAlign(allIm, nSubIm, nPeaks, localThreshold) + %run this function on gpu + gpu = false; + + nImages = size(allIm,3); + nRows = size(allIm,1); + nCols = size(allIm,2); + + rs = nan(nImages,1); + cs = nan(nImages,1); + + %Calculate template + if gpu + tmp = gpuArray(allIm); + template = gpuArray(median(tmp,3)); + clear tmp + else + template = median(allIm,3); + end + template = template-mean(template(:)); + + %Calculate crosscorrelation portion for template + v2_square = repmat((template).^2,1,1,nSubIm); + F = repmat(fft2(template),1,1,nSubIm); + clear template templateMat + + for imIdx = 1:nSubIm:nImages + %Get until end if subset is > nImages when grabbing last chunk + %Instead of resizing v2 and F one could change indexing below + %but in this way the code below stays more clear + if imIdx+nSubIm-1 > nImages + nSubIm = nSubIm+(nImages-(imIdx+nSubIm-1)); + v2_square = v2_square(:,:,1:nSubIm); + F = F(:,:,1:nSubIm); + end + + if gpu + imMat = gpuArray(allIm(:,:,imIdx:imIdx+nSubIm-1)); + else + imMat = double(allIm(:,:,imIdx:imIdx+nSubIm-1)); + end + + %Substract mean intensity from each image + meanMat = mean(mean(imMat)); + meansMat = repmat(meanMat,size(imMat,1),size(imMat,2)); + clear meanMat + imMat = imMat-meansMat; + clear meansMat + + %Calculate cross correlation + M = fft2(imMat); + + v1_square = (imMat).^2; + clear imMat + tmp = sqrt(sum(sum(v1_square,1),2) .* sum(sum(v2_square,1),2)); + clear v1_square + denominator = repmat(tmp,nRows,nCols); + clear tmp + numerator = real(ifft2((M).*conj(F))); + clear M + cc = numerator./denominator; + clear numerator denominator + % Unfortunately imregionalmax is only 2D when ran on GPU + if gpu + BW = gpuArray(false(nRows,nCols,nSubIm)); + end + for i = 1:nSubIm + BW(:,:,i) = imregionalmax(cc(:,:,i)); + end + + % Get intensities at maxima + intensityMap = cc .* BW; + clear cc BW + + % Discard maxima whose value is too low + intensityMap(intensityMapnRows/2) = (rMat(rMat>nRows/2) - nRows); + rsMat(rMat<-nRows/2) = (rMat(rMat<-nRows/2) + nRows); + clear rMat + csMat(cMat>nCols/2) = (cMat(cMat>nCols/2) - nCols); + csMat(cMat<-nCols/2) = (cMat(cMat<-nCols/2) + nCols); + clear cMat + + rsMat = rsMat-1; + csMat = csMat-1; + + %Shift in y direction + rsMat = round(-rsMat); + + %Shift in x direction + csMat = round(-csMat); + + %Calculate mean shift and gather from gpu if necessary + if gpu + %rs(imIdx:imIdx+nSubIm-1) = gather(mean(rsMat)); + %cs(imIdx:imIdx+nSubIm-1) = gather(mean(csMat)); + rs(imIdx:imIdx+nSubIm-1) = gather(rsMat(1,:)); + cs(imIdx:imIdx+nSubIm-1) = gather(csMat(1,:)); + else + %rs(imIdx:imIdx+nSubIm-1) = mean(rsMat); + %cs(imIdx:imIdx+nSubIm-1) = mean(csMat); + rs(imIdx:imIdx+nSubIm-1) = rsMat(1,:); + cs(imIdx:imIdx+nSubIm-1) = csMat(1,:); + end + clear rsMat csMat + end + clear v2_square F +end \ No newline at end of file diff --git a/medAlign_new/testscript.m b/medAlign_new/testscript.m new file mode 100644 index 0000000..ed5abca --- /dev/null +++ b/medAlign_new/testscript.m @@ -0,0 +1,21 @@ +%% load and prepare all single tiffs +%importing +[data_g,data_r] = CellSortLoadSingleFiles('ChanA*.tif', 'ChanB*.tif'); + +%Realignment +nSubIm = 10; +[rs, cs] = medAlign(data_r, nSubIm, 50, -inf); + +%Rotation and translation on red channel +for i = 1:size(data_r,3) + data_r_corr(:,:,i) = circshift(data_r(:,:,i),[rs(i),cs(i)]); +end + +%Rotation and translation on red channel +for i = 1:size(data_g,3) + data_g_corr(:,:,i) = circshift(data_g(:,:,i),[rs(i),cs(i)]); +end + +binningFactor = [1 1 2]; %x, y, t +[data_g_corr, data_r_corr] = CellSortTemporalBinning(data_g_corr, data_r_corr, binningFactor); +