diff --git a/inPhase_Segmentation_v11.m b/inPhase_Segmentation_v11.m
new file mode 100644
index 0000000000000000000000000000000000000000..c21701bb869f7223f0efeea8c98df7edbfc6eadc
--- /dev/null
+++ b/inPhase_Segmentation_v11.m
@@ -0,0 +1,374 @@
+%% inPhase segemntation routine
+% It assumes pairs of BF images and fluorescent maximum projections in .tif format
+% These pairs should have the same name with "BF" and "Max" as labels.
+% these pairs are supposed to have the used temperature in the name of the file (_xxC_)
+% Example:
+% BF_proteinName_20C
+% Max_proteinName_20C
+% This routine runs iteratively through all subfolders in the selected folder.
+
+% User input for the folder to start evaluation
+pathi = uigetdir();
+filesTiff = dir([pathi '/**/*.tif']);
+folderExp = unique({filesTiff.folder});
+
+% Save images of segmentation for all individual emulsion droplets
+saveImage = true;
+
+% Experimental setup related constants
+pixelSize = 0.1625;
+darkPixel = 100;
+
+% Iterates trough all the unique folders and searches for tiff pairs for
+% analysis
+
+for k = 1 : length(folderExp)
+
+ currFolderExp = folderExp(k);
+ currFilesTiff = dir([currFolderExp{1} '/**/*.tif']);
+ NameExp = {currFilesTiff.name};
+ IndxBF = contains(NameExp, 'BF');
+ IndxFluo = contains(NameExp, 'Max');
+ NameBF = natsort(NameExp(IndxBF));
+ NameFluo = natsort(NameExp(IndxFluo));
+ NExp = numel(NameFluo);
+
+ % Read conctrations of old evaluations and ask for user input for new ones
+ filesResults = dir([currFolderExp{1} '/**/Results*.mat']);
+
+ if isempty(filesResults)
+ % assumes no previous evaluation
+ newEval = true;
+ prompt = {'Enter Protein [c]:',...
+ 'Enter protein name (e.g. PGL3',...
+ 'Enter labeled fraction:',...
+ 'Enter pH:',...
+ 'Enter Salt:',...
+ 'Enter Drug [c]:',...
+ 'Enter additional Info:'};
+ dlgTitle = 'User input PD analysis';
+ dims = [1 35];
+ defInput = {'1 2 3 4 5 6 7 8', 'PGL3', '0.05', '7.4', '150', '0', ''};
+ answer = inputdlg(prompt, dlgTitle, dims, defInput);
+ proteinConcentrations = str2num(answer{1})'; %#ok<*ST2NM>
+ proteinConcentrations = repmat(proteinConcentrations, NExp/length(proteinConcentrations), 1);
+ proteinName = answer{2};
+ labelingFraction = str2num(answer{3})';
+ labelingFraction = repmat(labelingFraction, NExp/length(labelingFraction), 1);
+ pH = str2num(answer{4})';
+ pH = repmat(pH, NExp/length(pH), 1);
+ saltConcentration = str2num(answer{5})';
+ saltConcentration = repmat(saltConcentration, NExp/length(saltConcentration), 1);
+ drugConcentration = str2num(answer{6})';
+ drugConcentration = repmat(drugConcentration, NExp/length(drugConcentration), 1);
+ additionalInfo = answer{7};
+ else contains([filesResults.name], 'ResultsNew.mat')
+ % assumes to have the new version of ResultsNew.mat
+ newEval = false;
+ indx = contains({filesResults.name}, 'ResultsNew.mat');
+ load(fullfile(filesResults(indx).folder, filesResults(indx).name))
+ proteinConcentrations = cell2mat(getDataAllFields(Results, 'proteinConcentration'))';
+ proteinConcentrations = repmat(proteinConcentrations, NExp/length(proteinConcentrations), 1);
+ proteinName = getDataAllFields(Results, 'proteinName');
+ proteinName = cell2mat(unique(proteinName));
+ labelingFraction = cell2mat(getDataAllFields(Results, 'labelingFraction'))';
+ labelingFraction = repmat(labelingFraction, NExp/length(labelingFraction), 1);
+ try
+ pH = cell2mat(getDataAllFields(Results, 'pH'))';
+ pH = repmat(pH, NExp/length(pH), 1);
+ catch
+ prompt = {'Enter pH:'};
+ dlgTitle = 'User input PD analysis';
+ defInput = {'7.4'};
+ dims = [1 35];
+ answer = inputdlg(prompt, dlgTitle, dims, defInput);
+ pH = str2num(answer{1})';
+ pH = repmat(pH, NExp/length(pH), 1);
+ end
+ try
+ saltConcentration = cell2mat(getDataAllFields(Results, 'saltConcentration'))';
+ saltConcentration = repmat(saltConcentration, NExp/length(saltConcentration), 1);
+ catch
+ prompt = {'Enter Salt:'};
+ dlgTitle = 'User input PD analysis';
+ defInput = {'150'};
+ dims = [1 35];
+ answer = inputdlg(prompt, dlgTitle, dims, defInput);
+ saltConcentration = str2num(answer{1})';
+ saltConcentration = repmat(saltConcentration, NExp/length(saltConcentration), 1);
+ end
+ try
+ additionalInfo = getDataAllFields(Results, 'additionalInfo');
+ additionalInfo = cell2mat(unique(additionalInfo));
+ catch
+ additionalInfo = [];
+ end
+ try
+ drugConcentration = cell2mat(getDataAllFields(Results, 'saltConcentration'))';
+ drugConcentration = repmat(drugConcentration, NExp/length(drugConcentration), 1);
+ catch
+ drugConcentration = [];
+ end
+ end
+
+ % Segmentation of emulsion droplets in the specified folder
+
+ % Loads a pair of BF and Fluo data and segements condensates in
+ % individual emulsion droplets
+ Results = struct;
+ for kk = 1 : NExp
+
+ tic
+ % Display number of experiments in command window
+ disp(['Folder: ', num2str(k) ,' / ', num2str(length(folderExp)), ' | Experiment: ', num2str(kk) ,' / ', num2str(NExp)])
+
+ % Reset data container
+ DataIndividual = struct;
+ ExperimentInfo = struct;
+
+ % get temperature and date from experiment name
+ try
+ Indx = regexp(NameFluo{kk}, '_\d\dC_');
+ temperature = str2double(NameFluo{kk}(Indx+1 : Indx+2));
+ catch
+ temperature = NaN;
+ end
+ try
+ experimentDate = fullfile(currFolderExp{1},NameFluo{kk});
+ indx = regexp(experimentDate, '\d\d\d\d\d\d\_');
+ experimentDate = datetime(str2double(experimentDate(indx : indx+5)), 'ConvertFrom', 'YYYYMMDD', 'Format','yy-MM-dd');
+ catch
+ experimentDate = NaN;
+ end
+
+ ExperimentInfo.saltConcentration = saltConcentration(kk);
+ ExperimentInfo.temperature = temperature;
+ ExperimentInfo.proteinConcentration = proteinConcentrations(kk);
+ ExperimentInfo.drugConcentration = drugConcentration(kk);
+ ExperimentInfo.labelingFraction = labelingFraction(kk);
+ ExperimentInfo.pH = pH(kk);
+ ExperimentInfo.proteinName = proteinName;
+ ExperimentInfo.experimentDate = experimentDate;
+ ExperimentInfo.fileNameFluo = NameFluo{kk};
+ ExperimentInfo.fileNameBF = NameBF{kk};
+ ExperimentInfo.additionalInfo = additionalInfo;
+ ExperimentInfo.folderExperiment = currFolderExp{1};
+ ExperimentInfo.pixelSize = pixelSize;
+ ExperimentInfo.darkPixel = darkPixel;
+
+ % load images
+ currImage = im2double(imread(fullfile(currFolderExp{1}, NameFluo{kk})) - darkPixel); % changed 21.06.16 (- darkPixel)
+ currBF = im2double(imread(fullfile(currFolderExp{1}, NameBF{kk})));
+
+ % find water-in -oil emulsions in reduced size image
+ % first image filtering
+ scaleFac = 0.10;
+ appPixel = pixelSize / scaleFac;
+ sens_Circle = 0.85; % tricky parameter
+ currImg_adj = imadjust(imresize(currImage, scaleFac));
+ currBF_adj2 = imadjust(imresize(currBF, scaleFac));
+ currBF_adj = imadjust(imresize(currBF, scaleFac));
+ currBF_adj = imgaussfilt(currBF_adj, 3);
+ currBF_adj = imbinarize(currBF_adj);
+ currmaskDist = - bwdist(~currBF_adj);
+ maskDist2 = imimposemin(currmaskDist, imextendedmin(currmaskDist, 2));
+ maskWatershed = watershed(maskDist2);
+ currBF_adj(maskWatershed == 0) = 0;
+ currBF_adj = logical(currBF_adj);
+ currBF_adj = imfill(currBF_adj, 'holes');
+ currBF_adj = bwareaopen(currBF_adj, round(4*pi^-1*(2*30/appPixel)^2));
+ % second detection of circular structures of the right size
+ [centers1, radii1, metric1] = imfindcircles(currBF_adj, [round(50/appPixel) round(150/appPixel)], 'Sensitivity', sens_Circle, 'Method', 'TwoStage', 'ObjectPolarity', 'bright');
+ [centers2, radii2, metric2] = imfindcircles(currBF_adj, [round(151/appPixel) round(450/appPixel)], 'Sensitivity', sens_Circle, 'Method', 'TwoStage', 'ObjectPolarity', 'bright');
+ emulsionCenters = [centers1; centers2] ./ scaleFac;
+ emulsionRadii = [radii1; radii2] ./ scaleFac;
+
+ % first quality control
+ metric = [metric1; metric2];
+ Indx = metric > 0.15;
+ emulsionCenters = emulsionCenters(Indx,:);
+ emulsionRadii = emulsionRadii(Indx);
+
+ % quality control of circles (find intersecting/inside circles and remove the ones that intersect 2 other)
+ if ~isempty(emulsionRadii)
+ circleIntersect = zeros(numel(emulsionRadii));
+ circleInside = zeros(numel(emulsionRadii));
+ for kkk = 1 : numel(emulsionRadii)
+ circleIntersect(:,kkk) = sqrt((emulsionCenters(:, 1) - emulsionCenters(kkk, 1)).^2 +...
+ (emulsionCenters(:, 2) - emulsionCenters(kkk, 2)).^2) <...
+ emulsionRadii + emulsionRadii(kkk);
+ circleInside(:,kkk) = sqrt((emulsionCenters(:, 1) - emulsionCenters(kkk, 1)).^2 +...
+ (emulsionCenters(:, 2) - emulsionCenters(kkk, 2)).^2) <...
+ min(emulsionRadii, emulsionRadii(kkk));
+ end
+ circleIntersect = abs(sum(tril(circleIntersect),2) - sum(tril(circleIntersect),1)');
+ circleInside = sum(tril(circleInside),2);
+ Indx = or(circleIntersect >= 2, circleInside >= 2);
+ emulsionCenters = emulsionCenters(~Indx, :);
+ emulsionRadii = emulsionRadii(~Indx);
+ % and circles that cut the margins of the image
+ Indx = or(or((emulsionCenters(:,1) - emulsionRadii) < 0, (emulsionCenters(:,2) - emulsionRadii) < 0),...
+ or((emulsionCenters(:,1) + emulsionRadii) > size(currBF,2), (emulsionCenters(:,2) + emulsionRadii) > size(currBF,1)));
+ emulsionCenters = emulsionCenters(~Indx, :);
+ emulsionRadii = emulsionRadii(~Indx);
+ end
+
+ % save the BF and Fluo image with the circles overlayed
+ if saveImage
+ hf1 = figure('visible','off');
+ imshow(currBF_adj2)
+ viscircles(gca(hf1), emulsionCenters .* scaleFac, emulsionRadii .* scaleFac,'Color','r');
+ hf2 = figure('visible','off');
+ imshow(currImg_adj)
+ viscircles(gca(hf2), emulsionCenters .* scaleFac, emulsionRadii .* scaleFac,'Color','r');
+ end
+
+ % generate mask from the circles found
+ bboxes = [(emulsionCenters - emulsionRadii + 1), (2*emulsionRadii -1), (2*emulsionRadii -1)];
+ emulsionRadius = emulsionRadii .* pixelSize;
+ Indx = emulsionRadius <= 50;
+ emulsionVolume = zeros(size(emulsionRadius));
+ emulsionVolume(Indx) = 4/3 .* pi .* emulsionRadius(Indx) .^3;
+ emulsionVolume(~Indx) = 4/3 .* pi .* emulsionRadius(~Indx) .^2 .* 50; % parafilm is ~100um in height!
+ DataTable = table(emulsionRadius, emulsionVolume);
+
+ % Export figures of the current parafilm strips
+ if saveImage
+ if ~exist([currFolderExp{1}, filesep, 'EvalPics'], 'dir')
+ mkdir([currFolderExp{1}, filesep, 'EvalPics'])
+ end
+ try
+ matlab.graphics.internal.export.exportTo(get(hf1,'CurrentAxes'), [currFolderExp{1}, filesep, 'EvalPics', filesep, NameFluo{kk}(1:end-4), '_OverviewBF', '.pdf'])
+ matlab.graphics.internal.export.exportTo(get(hf2,'CurrentAxes'), [currFolderExp{1}, filesep, 'EvalPics', filesep, NameFluo{kk}(1:end-4), '_OverviewFluo', '.pdf'])
+ close(hf1)
+ close(hf2)
+ catch
+ saveas(get(hf1,'CurrentAxes'), [currFolderExp{1}, filesep, 'EvalPics', filesep, NameFluo{kk}(1:end-4), '_OverviewBF', '.pdf'])
+ saveas(get(hf2,'CurrentAxes'), [currFolderExp{1}, filesep, 'EvalPics', filesep, NameFluo{kk}(1:end-4), '_OverviewFluo', '.pdf'])
+ close(hf1)
+ close(hf2)
+ end
+ end
+
+ % Display number of emulsion droplets in command window
+ disp(['Number of emulsion droplets: ', num2str(length(emulsionRadii))])
+
+ % Analyse the individual emulsion droplets
+ for kkk = 1 : length(emulsionRadii)
+
+ % generate a mask for the current emulsion droplet
+ currRadius = emulsionRadii(kkk);
+ [currX, currY] = meshgrid(1 : 2*currRadius, 1 : 2*currRadius);
+ emulsionMask = (currX-currRadius).^2 + (currY-currRadius).^2 < emulsionRadii(kkk)^2;
+
+ % cut out fluo image with boundingbox (using actual size of currImage)
+ % and catch problems with crop at the edges of the image
+ currBBox = round(bboxes(kkk,:));
+ if currBBox(1) == 0
+ currBBox(3:4) = [size(emulsionMask,1), size(emulsionMask,2)-1];
+ elseif currBBox(2) == 0
+ currBBox(3:4) = [size(emulsionMask,1)-1, size(emulsionMask,2)];
+ elseif currBBox(2) + currBBox(3) >= size(currImage, 1)
+ currBBox(3:4) = [size(emulsionMask,1)-1, size(emulsionMask,2)];
+ offSet = currBBox(2) + currBBox(3) - size(currImage, 1);
+ if offSet < 0
+ offSet = 0;
+ end
+ emulsionMask = emulsionMask(1 : end - offSet, :);
+ elseif currBBox(1) + currBBox(4) >= size(currImage, 2)
+ currBBox(3:4) = [size(emulsionMask,1), size(emulsionMask,2)-1];
+ offSet = currBBox(1) + currBBox(4) - size(currImage, 2);
+ emulsionMask = emulsionMask(:, 1 : end - offSet);
+ else
+ currBBox(3:4) = size(emulsionMask) - 1;
+ end
+ currImageCrop = imcrop(currImage, currBBox);
+ currImageCrop(~emulsionMask) = 0;
+
+ % last resort if all the fun before does not work...
+ if any(size(currImageCrop) ~= size(emulsionMask))
+ try
+ currImageCrop = currImageCrop(1:size(emulsionMask,1), 1:size(emulsionMask,2));
+ catch
+ currImageCrop = [];
+ emulsionMask = [];
+ end
+ end
+
+ % simple BF analysis
+ currImageCropBF = imadjust(imcrop(currBF, currBBox));
+ currImageCropBF = imgaussfilt(currImageCropBF, 3);
+ currImageCropBF2 = imgradient(currImageCropBF);
+ currImageCropBF2(~emulsionMask) = mean(currImageCropBF2(emulsionMask));
+ [centerBF, radiusBF, metricBF] = imfindcircles(imadjust(currImageCropBF2), [10, round(size(currImageCropBF2, 1)/5)], 'Sensitivity', 0.6, 'ObjectPolarity', 'dark');
+ if saveImage
+ hf1 = figure('visible', 'off');
+ imshow(currImageCropBF)
+ viscircles(gca(hf1), centerBF, radiusBF, 'Color', 'r');
+ % save
+ matlab.graphics.internal.export.exportTo(get(hf1,'CurrentAxes'), [currFolderExp{1}, filesep, 'EvalPics', filesep, NameFluo{kk}(1:end-4), 'SegBF_', num2str(kkk), '.png'])
+ close(hf1)
+ end
+
+ % image segementation and data handling
+ [DataOut, maskDroplets] = condensate_Segmentation_V4(currImageCrop, emulsionMask, pixelSize);
+
+ DataTable.dropletTotalVolumeBF(kkk) = sum(4/3*pi .* (radiusBF.*pixelSize) .^3);
+ DataTable.emulsionVolume(kkk) = emulsionVolume(kkk);
+ DataTable.emulsionCenters(kkk,:) = emulsionCenters(kkk,:);
+ DataTable.soluteIntensity_mean(kkk) = DataOut.soluteIntensity_mean;
+ DataTable.soluteIntensity_std(kkk) = DataOut.soluteIntensity_std;
+ DataTable.dropletTotalVolume(kkk) = DataOut.dropletTotalVolume;
+ DataTable.volumeFraction(kkk) = DataOut.dropletTotalVolume / emulsionVolume(kkk);
+ DataTable.dropletIntensities_mean(kkk) = DataOut.dropletIntensities_mean;
+ DataTable.dropletIntensities_std(kkk) = DataOut.dropletIntensities_std;
+ DataTable.nDroplets(kkk) = length(DataOut.dropletIntensities);
+ DataTable.saltConcentration(kkk) = saltConcentration(kk);
+ DataTable.temperautre(kkk) = temperature;
+ DataTable.proteinConcentration(kkk) = proteinConcentrations(kk);
+ DataTable.pH(kkk) = pH(kk);
+ DataTable.labelingFraction(kkk) = labelingFraction(kk);
+ DataTable.proteinName(kkk,:) = proteinName;
+ DataTable.experimentDate(kkk) = experimentDate;
+
+ DataIndividual.dropletRadii{kkk} = DataOut.dropletRadii;
+ DataIndividual.dropletIntensities{kkk} = DataOut.dropletIntensities;
+ DataIndividual.dropletIntensitiesMax{kkk} = DataOut.dropletIntensitiesMax;
+ DataIndividual.dropletVolumes{kkk} = DataOut.dropletVolumes;
+ DataIndividual.dropletEccentricity{kkk} = DataOut.dropletEccentricity;
+ DataIndividual.dropletCircularity{kkk} = DataOut.dropletCircularity;
+ DataIndividual.maskDroplets{kkk} = maskDroplets;
+
+ % Overlay segmentation on croped image and save
+ if saveImage
+ try
+ hf1 = figure('visible', 'off');
+ imshow(currImageCrop, 'DisplayRange', [min(min(currImageCrop)) max(max(currImageCrop))])
+ saveas(hf1,[currFolderExp{1}, filesep, 'EvalPics', filesep, NameFluo{kk}(1:end-4),'FluoCrop_', num2str(kkk), '.png']);
+ close(hf1)
+ hf2 = figure('visible', 'off');
+ maskDroplets(maskDroplets>1) = 1;
+ Lrgb = label2rgb(maskDroplets, 'jet', 'k', 'shuffle');
+ himage = imshow(Lrgb);
+ saveas(hf2,[currFolderExp{1}, filesep, 'EvalPics', filesep, NameFluo{kk}(1:end-4),'SegFluo_', num2str(kkk), '.png']);
+ close(hf2)
+ catch
+ end
+ end
+
+ end
+
+ % Output data
+ currFileName = ['Experiment_', num2str(temperature, '%02d'), 'C_', num2str(kk, '%03d')];
+ Results.(currFileName).DataTable = DataTable;
+ Results.(currFileName).DataIndividual = DataIndividual;
+ Results.(currFileName).ExperimentInfo = ExperimentInfo;
+ toc
+
+ end
+
+ % save Output data as .mat
+ save([currFolderExp{1}, filesep, 'ResultsNew.mat'], 'Results')
+
+end
+