alignmentScript.m

%%
clear all
close all
clc
%% Set Paramters
inFolder = 'D:\data\Leona';
outFolder = 'D:\data\Out1234';
binningFactor = [1, 1, 2];
nSubIm = 200;%number of subImages for med align
numImages = 50; %number of images that will be meaned for global shift correction
referenceChannel = 1; %0: green, 1: red; reference channel for shift corrections
numShiftCorrections = 1; %number of shift corrections
bScope = 0; %set 1 if b scope, 0 if not

%% Load Files
tic;
currentFolder = pwd;
cd(inFolder);
folders = dir;
binningCounter = [1,1,1];
for k = length(folders):-1:1
    % remove non-folders
    if ~folders(k).isdir
        folders(k) = [ ];
        continue
    end
    % remove folders starting with .
    fname = folders(k).name;
    if fname(1) == '.'
        folders(k) = [ ];
    end
end

numFolders = length(folders);

%iterate through all folders
for f = 1:numFolders
    %go into folder
    cd(inFolder);
    cd(folders(f).name);
    if bscope
        greenFiles = dir('ChanA_*_*.tif');
        redFiles = dir('ChanB_*_*.tif');
    else
        greenFiles = dir('*Ch2_*.tif');
        redFiles = dir('*Ch1_*.tif');
    end
    %check if fileNumber is the same
    if size(greenFiles, 1) ~= size(redFiles, 1)
        error(['Unequal number of files in trial folder "', folders(f).name, '"!!']);
    end
    numFiles = size(greenFiles, 1);
    
    %load first red frame to get resolution and being able to preallocate
    temp = readTiffNoStack(redFiles(1).name);
    res = size(temp);
    
    %preallocate data_r
    data_r = zeros(res(1), res(2), numFiles, class(temp));
    data_g = zeros(res(1), res(2), numFiles, class(temp));
    %load data
    disp(['Loading Files In Folder ', folders(f).name]);
    parfor i=1:numFiles
        data_r(:,:,i) = readTiffNoStack(redFiles(i).name);
        data_g(:,:,i) = readTiffNoStack(greenFiles(i).name);
    end
    
    %save stacks in memory
    data_r_complete{f} = data_r;
    data_g_complete{f} = data_g;
end
clear data_r data_g

%go back to old current folder
cd(currentFolder);
data_r = data_r_complete;
data_g = data_g_complete;
clear data_r_complete data_g_complete
%create big array
numFrames = 0;
numFrameList = zeros(1,numFolders);
for i = 1:numFolders
    numFrames = numFrames + size(data_g{i},3);
    numFrameList(i+1) = numFrames;
end

%copy data in big matrix

%red channel
data_r_array = zeros(size(data_g{i},1),size(data_g{i},2),numFrames, class(data_r{1}));
for i = 1:numFolders
    data_r_array(:,:,(numFrameList(i)+1):numFrameList(i+1)) = data_r{i};
    data_r{i} = [];
end

%green channel
data_g_array = zeros(size(data_g{i},1),size(data_g{i},2),numFrames, class(data_g{1}));
for i = 1:numFolders
    data_g_array(:,:,(numFrameList(i)+1):numFrameList(i+1)) = data_g{i};
    data_g{i} = [];
end
clear data_r_complete data_g_complete data_r data_g
data_r = data_r_array;
data_g = data_g_array;
clear data_g_array data_r_array
%% Do Shift Correction
%wb = waitbar(0, 'Performing Shift Correction');
disp('Performing first shift correction');
for j = 1:numShiftCorrections
    if referenceChannel
        referenceChannelData = data_r;
    else
        referenceChannelData = data_g;
    end
    parfor f=1:size(referenceChannelData,3)
        referenceChannelData(:,:,f) = imgaussfilt(referenceChannelData(:,:,f),2);
    end
    
    %Med Align
    [rs, cs] = medAlign(referenceChannelData, nSubIm, 50, -inf);
    %shift red and green channel
    %waitbar(j/numShiftCorrections, wb, 'Performing Med Align..');
    parfor i = 1:size(data_g,3)
        data_r(:,:,i) = circshift(data_r(:,:,i),[rs(i),cs(i)]);
        data_g(:,:,i) = circshift(data_g(:,:,i),[rs(i),cs(i)]);
    end
end
%matVis(data_r,data_g);
%% Do binning
disp('Binning');
data_r = cast((binning(data_r, binningFactor, 'mean')), class(temp));
data_g = cast((binning(data_g, binningFactor, 'mean')), class(temp));
binningCounter = binningCounter .* binningFactor;
numFrameList = numFrameList/binningFactor(3);
%% Do Shift Correction 2
%wb = waitbar(0, 'Performing Shift Correction');
disp('Performing more shift corrections');
for j = 1:numShiftCorrections
    if referenceChannel
        referenceChannelData = data_r;
    else
        referenceChannelData = data_g;
    end
    %Med Align
    parfor f=1:size(referenceChannelData,3)
        referenceChannelData(:,:,f) = imgaussfilt(referenceChannelData(:,:,f),2);
    end
    [rs, cs] = medAlign(referenceChannelData, nSubIm, 50, -inf);
    %shift red and green channel
    %waitbar(j/numShiftCorrections, wb, 'Performing Med Align..');
    parfor i = 1:size(data_g,3)
        data_r(:,:,i) = circshift(data_r(:,:,i),[rs(i),cs(i)]);
        data_g(:,:,i) = circshift(data_g(:,:,i),[rs(i),cs(i)]);
    end
end
%matVis(data_r,data_g);
%% Do binning 2
disp('Second Binning');
data_r = cast(binning(data_r, binningFactor, 'mean'),class(temp));
data_g = cast(binning(data_g, binningFactor, 'mean'), class(temp));
binningCounter = binningCounter .* binningFactor;
numFrameList = numFrameList/binningFactor(3);

%% Correct for slow global drifts
disp('Correction for slow global drifts');
% waitbar(0, wb, 'Performing Global Shift Correction');
if referenceChannel
    referenceChannelData = data_r;
else
    referenceChannelData = data_g;
end


parfor f=1:size(referenceChannelData,3)
    referenceChannelData(:,:,f) = imgaussfilt(referenceChannelData(:,:,f),2);
end
%generate mean images dataset
imageCounter = 1;
i=1;
while imageCounter <= size(referenceChannelData,3)
    %     waitbar(imageCounter/size(referenceChannelData,3), wb,'Create Chunked Data Set');
    if imageCounter+numImages-1 > size(referenceChannelData,3)
        shiftData(:,:,i) = mean(referenceChannelData(:,:,imageCounter:end),3);
    else
        shiftData(:,:,i) = mean(referenceChannelData(:,:,imageCounter:imageCounter+numImages-1),3);
    end
    lastIntervalSize = size(referenceChannelData,3) - imageCounter+1;
    imageCounter = imageCounter + numImages;
    i=i+1;
end

%Med Align
[rs, cs] = medAlign(shiftData, size(shiftData,3), 50, -inf);

%shift red and green channel
for i = 1:size(shiftData,3)
    %     waitbar(i/size(shiftData,3), wb,'Perform Global Shift Correction');
    if i == size(shiftData,3)
        imagesToCorrect = lastIntervalSize;
    else
        imagesToCorrect = numImages;
    end
    for j=1:imagesToCorrect
        data_r(:,:,(i-1)*numImages+j) = circshift(data_r(:,:,(i-1)*numImages+j),[rs(i),cs(i)]);
        data_g(:,:,(i-1)*numImages+j) = circshift(data_g(:,:,(i-1)*numImages+j),[rs(i),cs(i)]);
    end
end
data_r(:,:,(i-1)*numImages+j:end) = circshift(data_r(:,:,(i-1)*numImages+j:end),[rs(i),cs(i)]);
data_g(:,:,(i-1)*numImages+j:end) = circshift(data_g(:,:,(i-1)*numImages+j:end),[rs(i),cs(i)]);
clear shiftData
%% Export Files
%save trial files
disp('Saving single trials in files');
for i=1:numFolders
    %     waitbar(i/numFolders, wb, 'Saving Trial Tiff Files..');
    data_r_save = (data_r(:,:,floor(numFrameList(i)+1):floor(numFrameList(i+1))));
    data_g_save = (data_g(:,:,floor(numFrameList(i)+1):floor(numFrameList(i+1))));
    
    saveastiff(data_r_save, [outFolder filesep folders(i).name  '_r.tif']);
    saveastiff(data_g_save, [outFolder filesep folders(i).name  '_g.tif']);
end
clear data_r_save data_g_save

%save big tiff file
disp('Saving Big Data Set');
%waitbar(0, wb, 'Saving Big Tiff Files..');
options.big = true;
saveastiff(data_r, [outFolder filesep 'complete_r.tif'], options);
%waitbar(0.5, wb, 'Saving Big Tiff Files..');
saveastiff(data_g, [outFolder filesep 'complete_g.tif'], options);
%%
% close(wb);
disp('Finished Script :)!!!');
toc