diff --git a/functions/trialAlignmentScript_Belen.m b/functions/trialAlignmentScript_Belen.m index 1bf96fe..86e9443 100644 --- a/functions/trialAlignmentScript_Belen.m +++ b/functions/trialAlignmentScript_Belen.m @@ -88,59 +88,59 @@ %%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'); + 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 -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)]); + %%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 - 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'); + %%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 -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'); + %%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 -end %%visualize data numFrames = 0; numFrameList = zeros(1,numberGreen); diff --git a/functions/trialAlignmentScript_newMicroscope.m b/functions/trialAlignmentScript_newMicroscope.m new file mode 100644 index 0000000..62a7566 --- /dev/null +++ b/functions/trialAlignmentScript_newMicroscope.m @@ -0,0 +1,190 @@ +% 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!!!'); + + diff --git a/functions/trialAlignmentScript_test.m b/functions/trialAlignmentScript_test.m new file mode 100644 index 0000000..5930bfb --- /dev/null +++ b/functions/trialAlignmentScript_test.m @@ -0,0 +1,106 @@ +% trialAlignmentScript_test +clc; +close all; +% +clear variables; + +%% user input +%pathInput = '/Users/tushevg/Desktop/new_2P_test/TSeries-06072017-1138-048'; +%pathInput = '/Users/tushevg/Desktop/002'; +pathInput = '/Users/tushevg/Desktop/new_2P_test'; +pathOutput = '/Users/tushevg/Desktop/testout'; + +numberOfShiftCorrections = 4; +binningFactor = [1, 1, 2]; +nSubIm = 2; + +labelTrial = 'Cycle'; +labelChannel = 'Ch'; + +indexRed = 1; +indexGreen = 2; + + +%% explore parent directory +dirInfo = dir(pathInput); +subDirs = {dirInfo.name}'; +subDirs(strncmp('.',subDirs,1) | ~[dirInfo.isdir]') = []; +subDirsCount = length(subDirs); + +%% explore child direcotries +d = 1; +%for d = 1 : subDirsCount + + %%% extract files from sub directory + subDirInfo = dir([pathInput, filesep, subDirs{d}, filesep, '*.tif']); + fileNames = {subDirInfo.name}'; + + %%% extract lists of trial/channel names + trialNames = regexp(fileNames, [labelTrial,'[0-9]+'], 'once', 'match'); + [trialUnqNames,~,idxunq_trials] = unique(trialNames); + trialCount = max(idxunq_trials); + + channelNames = regexp(fileNames, [labelChannel,'([12]|[AB])'], 'once', 'match'); + [~, ~, idxunq_channels] = unique(channelNames); + + %%% allocate images + imgInfo = imfinfo([pathInput, filesep, subDirs{d}, filesep, fileNames{1}]); + imgStack = sum((idxunq_trials == 1) & (idxunq_channels == indexRed)); + imgRed = zeros(imgInfo.Height,... + imgInfo.Width,... + imgStack,... + trialCount,... + sprintf('uint%d',imgInfo.BitDepth)); + imgGreen = zeros(imgInfo.Height,... + imgInfo.Width,... + imgStack,... + trialCount,... + sprintf('uint%d',imgInfo.BitDepth)); + + + %% process trials + t = 1; + %for t = 1 : trialCount + + %%% processing trial + fprintf('### processing trial %s\n', trialUnqNames{t}); + + %%% current trial and channels per trial + idx_trial = (idxunq_trials == t); + idx_red = idx_trial & (idxunq_channels == indexRed); + idx_green = idx_trial & (idxunq_channels == indexGreen); + + %%% full file lists + list_red = strcat(pathInput, filesep, subDirs{d}, filesep, fileNames(idx_red)); + list_green = strcat(pathInput, filesep, subDirs{d}, filesep, fileNames(idx_green)); + + %%% read images + fprintf('# read image stacks %d ... ', imgStack); + tic + for s = 1 : imgStack + imgRed(:, :, s, t) = imread(list_red{s}); + imgGreen(:, :, s, t) = imread(list_green{s}); + end + fprintf('done. %.4f sec\n', toc); + + %%% do shift correction on Red channel + fprintf('# shift correction on Red channel ... '); + tic + [rs, cs] = medAlign(imgRed(:,:,:,t), nSubIm, 50, -inf); + fprintf('done. %.4f sec\n', toc); + + %%% re-align channels to shift + fprintf('# re-align channels ... '); + tic + for s = 1 : imgStack + imgRed(:, :, s, t) = circshift(imgRed(:, :, s, t), [rs(s), cs(s)]); + imgGreen(:, :, s, t) = circshift(imgGreen(:, :, s, t), [rs(s), cs(s)]); + end + fprintf('done. %.4f sec\n', toc); + + %end + + %% translate to old version + + + \ No newline at end of file