trialAlignmentScript_Belen.m

clear all
close all
clc
%% Load Trial Files
numberShiftCorrections = 4;
inFolder = 'D:\data\Leona';
outFolder = 'D:\data\out';
binningFactor = [1, 1, 2];
nSubIm = 10;



%% load single tif files of trials to ram
%operations
currentFolder = pwd;
cd(inFolder);
folders = dir;

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);
    greenFiles = dir('ChanA_*_*.tif');
    redFiles = dir('ChanB_*_*.tif');
    %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);
    data_g = zeros(res(1), res(2), numFiles);
    %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
    
    disp(['Shift Correction Of Files In Folder ', folders(f).name]);
    %do the shift correction
    [rs, cs] = medAlign(data_r, nSubIm, 50, -inf);
    
    %shift red and green channel
    disp(['Shift Data Arrays Of Files In Folder ', folders(f).name]);
    parfor i = 1:numFiles
        data_r(:,:,i) = circshift(data_r(:,:,i),[rs(i),cs(i)]);
        data_g(:,:,i) = circshift(data_g(:,:,i),[rs(i),cs(i)]);
    end
        
    %save stacks in memory
    cd(outFolder);
    disp(['Write Files Of Folder ', folders(f).name ' in ' outFolder ' to RAM']);
    data_r_complete{f} = data_r;
    data_g_complete{f} = data_g;
end

disp('Finished Shift Correction And Fusing Of Data!');
%go back to old current folder
cd(currentFolder);
clear data_r data_g
data_r = data_r_complete;
data_g = data_g_complete;
clear data_r_complete data_g_complete
numberGreen = size(data_g,2);
%% do shift correction
%%shift correct trials to each other
%do shift correction
for shiftCorrectionCounter = 1:numberShiftCorrections
    disp('Performing Shift Correction On All Trials..');
    reference = median(data_r{1},3);
    for i = 1:numberGreen
        disp(['Working on Trial ' num2str(i)]);
        medianProj = median(data_r{i},3);
        shiftVector = findshift(reference, medianProj, 'iter', 0, 0.25*size(reference,1));
        %shift frames
        for f = 1:size(data_r{i},3)
            data_r{i}(:,:,f)  = imtranslate(data_r{i}(:,:,f),shiftVector', 'linear');
            data_g{i}(:,:,f)  = imtranslate(data_g{i}(:,:,f),shiftVector', 'linear');
        end
    end
    %%do another med align
    for j =1:numberGreen
        disp(['Working on Trial ' num2str(j)]);
        [rs, cs] = medAlign(data_r{j}, nSubIm, 50, -inf);
        temp_r = data_r{j};
        temp_g = data_g{j};
        %shift red and green channel
        parfor f = 1:size(data_r{j},3)
            temp_r(:,:,f) = circshift(temp_r(:,:,f),[rs(f),cs(f)]);
            temp_g(:,:,f) = circshift(temp_g(:,:,f),[rs(f),cs(f)]);
        end
        data_r{j} = temp_r;
        data_g{j} = temp_g;
    end
    %%do another shift correction for every frame, takes some time
    reference = median(data_r{1},3);
    numberGreen = size(data_g,2);
    for i = 1:numberGreen
        disp(['Working on Trial ' num2str(i)]);
        for f = 1:size(data_r{i},3)
        shiftVector = findshift(reference, data_r{j}(:,:,f), 'iter',0, 0.25*size(reference,1));
        %shift frames

            data_r{i}(:,:,f)  = imtranslate(data_r{i}(:,:,f),shiftVector', 'linear');
            data_g{i}(:,:,f)  = imtranslate(data_g{i}(:,:,f),shiftVector', 'linear');
        end
    end
    %%do another shift correctionto align trials to each other
    reference = median(data_r{1},3);
    numberGreen = size(data_g,2);
    for i = 1:numberGreen
        disp(['Working on Trial ' num2str(i)]);
        medianProj = median(data_r{i},3);
        shiftVector = findshift(reference, medianProj, 'iter',0,  0.25*size(reference,1));
        %shift frames
        for f = 1:size(data_r{i},3)
            data_r{i}(:,:,f)  = imtranslate(data_r{i}(:,:,f),shiftVector', 'linear');
            data_g{i}(:,:,f)  = imtranslate(data_g{i}(:,:,f),shiftVector', 'linear');
        end
    end
end
%%visualize data
numFrames = 0;
numFrameList = zeros(1,numberGreen);
for i = 1:numberGreen
    numFrames = numFrames + size(data_g{i},3);
    numFrameList(i+1) = numFrames;
end

%preallocate data
data_r_array = zeros(size(data_g{i},1),size(data_g{i},2),numFrames, class(data_r{1}));
data_g_array = zeros(size(data_g{i},1),size(data_g{i},2),numFrames, class(data_g{1}));
disp('Creating Big Data Set..');
for i = 1:numberGreen
    data_r_array(:,:,(numFrameList(i)+1):numFrameList(i+1)) = data_r{i};
    data_g_array(:,:,(numFrameList(i)+1):numFrameList(i+1)) = data_g{i};
end
%% save tiff files
parfor i=1:numberGreen    
    %binning
    data_r{i} = uint16(binning(data_r{i}, binningFactor, 'mean'));
    data_g{i} = uint16(binning(data_g{i}, binningFactor, 'mean'));

    saveastiff(data_r{i}, [outFolder filesep folders(i).name  '_r.tif']);
    saveastiff(data_g{i}, [outFolder filesep folders(i).name  '_g.tif']);
end
%% generate quality plot
%put them in single arrays
numFrames = 0;
numFrameList = zeros(1,numberGreen);
for i = 1:numberGreen
    numFrames = numFrames + size(data_g{i},3);
    numFrameList(i+1) = numFrames;
end

%preallocate data
data_r_array = zeros(size(data_g{i},1),size(data_g{i},2),numFrames, class(data_r{1}));
data_g_array = zeros(size(data_g{i},1),size(data_g{i},2),numFrames, class(data_g{1}));
disp('Creating Big Data Set..');
for i = 1:numberGreen
    data_r_array(:,:,(numFrameList(i)+1):numFrameList(i+1)) = data_r{i};
    data_r{i} = [];
    data_g_array(:,:,(numFrameList(i)+1):numFrameList(i+1)) = data_g{i};
    data_g{i} = [];
end

%calculate quality plot
reference = mean(data_r_array, 3);
qualityPlot = zeros(1,size(data_r_array, 3));
parfor i = 1:size(data_r_array,3)
    qualityPlot(i) = sqrt(squeeze(mean(mean((cast(data_r_array(:,:,i),'like', reference) - reference).^2,1),2)));
end
figure('Name','Quality Plot');
plot(qualityPlot);
for i=1:numberGreen
    vline(numFrameList(i), 'black', greenFiles(i).name);
end
title('Quality Of Shift Correction');
xlabel('frame');
ylabel('std of frame');
%% export complete files
disp('Saving Big Data Set..');
options.big = true;
saveastiff(data_r_array, [outFolder filesep  'complete_r.tif'], options);
saveastiff(data_g_array, [outFolder filesep  'complete_g.tif'], options);
%% Finished
clear all
disp('Finished Complete Script!!!');