condensate_Segmentation_V4.m

function [DataOut, maskDroplets] = condensate_Segmentation_V4(imageCrop, emulsionMask, pixelSize)
% in version 4 the watersheding was improved an switched on again

% calibrated for Olympus 40xUPLXAPO air objective  at Olympus IX83 microscope 
seSize = 5;
se = strel('disk', seSize);
se2 = strel('disk', 20);

% start image analysis
if std2(imageCrop) <= 5 / 2^16
    % no droplets in image
    maskDroplets = 0;
else
    % there are droplets start segmentation
    % gradient based particle detection
    imCropFilterd = roifilt2(imageCrop, emulsionMask, 'imgaussfilt', 2);
    imCropFilterd = roifilt2(imCropFilterd, emulsionMask, 'imgradient');
    imCropFilterd = imCropFilterd ./ (mean2(imCropFilterd) + 2 * std2(imCropFilterd));
    imCropFilterd(imCropFilterd <= 1) = 0;
    
    % bw treshholding step
    thresh = adaptthresh(imCropFilterd, 0.1) ;
    maskDroplets = imbinarize(imCropFilterd, thresh);
    
    % BW filtering
    maskDroplets = imfill(maskDroplets, 'holes');
    maskDroplets = imerode(maskDroplets, se); % removes small objects (careful with se)
    maskDroplets = bwareaopen(maskDroplets, seSize);
    
    % watershed to separate connected particles
    % inspired by https://blogs.mathworks.com/steve/2013/11/19/watershed-transform-question-from-tech-support/
    maskDist = - bwdist(~maskDroplets);
    maskDist2 = imimposemin(maskDist, imextendedmin(maskDist, 2));
    maskWatershed = watershed(maskDist2);
    maskDroplets(maskWatershed == 0) = 0;
    maskDroplets = logical(maskDroplets);
    
    % derive solute mask and mean intensity
    maskSolute = emulsionMask - imdilate(maskDroplets, se2);
    maskSolute(maskSolute < 0) = 0;
    maskSolute = logical(imerode(maskSolute, se2));
end

if ~sum(sum(maskDroplets)) == 0
    dropletDiameters = regionprops(maskDroplets, 'EquivDiameter');
    dropletRadii = cat(1, dropletDiameters.EquivDiameter) ./2 .* pixelSize;
    dropletVolumes = 4/3 .* pi .* dropletRadii.^3;
    dropletRegProps = regionprops(maskDroplets, imageCrop, 'MeanIntensity', 'MaxIntensity', 'Eccentricity', 'Circularity', 'Centroid');
    
    % Data output in case of condensates present
    DataOut.dropletRadii = dropletRadii;
    DataOut.dropletVolumes = dropletVolumes;
    DataOut.dropletTotalVolume = sum(dropletVolumes);
    DataOut.dropletIntensities = [dropletRegProps.MeanIntensity]';
    DataOut.dropletIntensities_mean = mean([dropletRegProps.MeanIntensity]);
    DataOut.dropletIntensities_std = std([dropletRegProps.MeanIntensity]);
    DataOut.dropletIntensitiesMax = [dropletRegProps.MaxIntensity]';
    DataOut.dropletEccentricity = [dropletRegProps.Eccentricity]';
    DataOut.dropletCircularity = [dropletRegProps.Circularity]';
    DataOut.dropletCentroid = [dropletRegProps.Circularity]';
    DataOut.soluteIntensity_mean = mean(imageCrop(maskSolute));
    DataOut.soluteIntensity_std = std(imageCrop(maskSolute));
else
    % Data output in case of no condensates present
    DataOut.dropletRadii = NaN;
    DataOut.dropletVolumes = NaN;
    DataOut.dropletTotalVolume =  NaN;
    DataOut.dropletIntensities =  NaN;
    DataOut.dropletIntensities_mean = NaN;
    DataOut.dropletIntensities_std = NaN;
    DataOut.dropletIntensitiesMax = NaN;
    DataOut.dropletEccentricity = NaN;
    DataOut.dropletCircularity = NaN;
    DataOut.dropletCentroid = NaN;
    DataOut.soluteIntensity_mean = mean(imageCrop(emulsionMask));
    DataOut.soluteIntensity_std = std(imageCrop(emulsionMask));
end

end