%% 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

