trialAlignmentScript_newMicroscope.m

% trialAlignmentScript_newMicroscope
clc;
close all;
clear variables;
warning('off','all');

%% user input
pathInput = '/Users/tushevg/Desktop/new_2P_test/TSeries-06072017-1138-048';

numberOfShiftCorrections = 4;
binningFactor = [1, 1, 2];
nSubIm = 2;

labelTrial = 'Cycle';
labelChannel = 'Ch';

indexRed = 1;
indexGreen = 2;

%% explore parent directory
dirInfo = dir([pathInput, filesep, '*.tif']);
fileNames = {dirInfo.name}';

%% extract lists of trial/channel names
trialNames = regexp(fileNames, [labelTrial,'[0-9]+'], 'once', 'match');
[trialUnqNames,~,indexTrial] = unique(trialNames);
trialCount = max(indexTrial);

channelNames = regexp(fileNames, [labelChannel,'([12]|[AB])'], 'once', 'match');
[~, ~, indexChannel] = unique(channelNames);

%% allocate containers
data_r = repmat({}, trialCount, 1);
data_g = repmat({}, trialCount, 1);

%% read all frames and trials
for t = 1 : trialCount
    
    %%% processing trial
    fprintf('### processing trial %s\n', trialUnqNames{t});

    %%% current trial and channels per trial
    idx_trial = (indexTrial == t);
    idx_red = idx_trial & (indexChannel == indexRed);
    idx_green = idx_trial & (indexChannel == indexGreen);
        
    %%% full file lists
    list_red = strcat(pathInput, filesep, fileNames(idx_red));
    list_green = strcat(pathInput, filesep, fileNames(idx_green));

    %%% allocate image stack
    imgInfo = imfinfo(list_red{1});
    imgStack = min([length(list_red),length(list_green)]);
    imgRed = zeros(imgInfo.Height,...
                   imgInfo.Width,...
                   imgStack,...
                   sprintf('uint%d',imgInfo.BitDepth));
    imgGreen = zeros(imgInfo.Height,...
                   imgInfo.Width,...
                   imgStack,...
                   sprintf('uint%d',imgInfo.BitDepth));
               
    %%% read images
    fprintf('# read image stacks %d ... ', imgStack);
    tic
    for s = 1 : imgStack
        imgRed(:, :, s) = readTiffNoStack(list_red{s});
        imgGreen(:, :, s) = readTiffNoStack(list_green{s});
    end
    fprintf('done. %.4f sec\n', toc);
    
    %%% fill containers
    data_r(t) = {imgRed};
    data_g(t) = {imgGreen};
    
end

%% ------------------------- %%
%% --- SHIFT CORRECTIONS --- %%
%% ------------------------- %%

%% shift correct trials to each other
%do shift correction
disp('Performing Shift Correction On All Trials..');
reference = median(data_g{1},3);
numberGreen = size(data_g,2);
for i = 1:numberGreen
    disp(['Working on Trial ' num2str(i)]);
    medianProj = median(data_g{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_g{j}, nSubIm, 50, -inf);
    temp_r = data_r{j};
    temp_g = data_g{j};
    %shift red and green channel
    parfor f = 1:size(data_g{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_g{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_g{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_g{1},3);
numberGreen = size(data_g,2);
for i = 1:numberGreen
    disp(['Working on Trial ' num2str(i)]);
    medianProj = median(data_g{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
%% save tiff files
for 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));
for 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 folders(i).name  '_complete_r.tif'], options);
saveastiff(data_g_array, [outFolder filesep folders(i).name  '_complete_g.tif'], options);

%% Finished
clear variables
disp('Finished Complete Script!!!');