Skip to content
Snippets Groups Projects
Commit ea855d27 authored by fritsch's avatar fritsch
Browse files

Upload New File

parent d6491eaf
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
inPhase_Segmentation_v11.m 0 → 100644
+ 374
− 0
View file @ ea855d27
%% 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
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment