%% Init
clc,clear all, close all
disp('### Semi-automatic image analysis')
% preselectVector = [1 3];
preselectVector = [0 0];
primaryLimitStruct.convexity = 0.90; % 0.9 before
primaryLimitStruct.diameter = 100; % 50 before
primaryLimitStruct.sphericity = 1; % off
primaryLimitStruct.sphericity2 = 1; % off
primaryLimitStruct.sphericity3 = 1; % off
primaryLimitStruct.angleGap = 36; % 36 before

% secondaryLimitStruct.convexity = 0.95;
% secondaryLimitStruct.sphericity = 0.1;
secondaryLimitStruct.sphericity2 = 0.96;
secondaryLimitStruct.sphericity3 = 0.995;
secondaryLimitStruct.minSize = 0;
% secondaryLimitStruct.maxSize = 2000;
secondaryLimitStruct.relRoughness = 1;
secondaryLimitStruct.roughness = 50;
% secondaryLimitStruct.maxRoundness = 0.95;
secondaryLimitStruct.minBrightFraction = 0.015;

flagViewPrimaryIgnored =  1;
flagSavePrimaryIgnored = 1;
flagViewSecondaryIgnored = 1;
flagSaveSecondaryIgnored = 1;

%% Select Experiment
disp('# Select experiment folder:')
% generate folder list from dir-command
fileList = dir;
filterVector = (1:length(fileList)) == 0;
for iFile = 1:length(fileList)
    if fileList(iFile).isdir && ~(fileList(iFile).name(1) == '.')
        filterVector(iFile) = 1;
    end
end
fileList = {fileList(filterVector).name};
% ask user or use preselection
if preselectVector(1) == 0
    selection = selectFromList(fileList, 'Exit');
    if selection == 0
        return
    end
else
    selection = preselectVector(1);
end
experimentFolder = fileList{selection};

%% Select subfolder
disp('# Select subfolder:')
fileList = dir(experimentFolder);
filterVector = (1:length(fileList)) == 0;
for iFile = 1:length(fileList)
    if fileList(iFile).isdir && ~(fileList(iFile).name(1) == '.')
        filterVector(iFile) = 1;
    end
end
fileList = {fileList(filterVector).name};
% ask user or use preselection
if preselectVector(1) == 0
    selection = selectFromList(fileList, 'Exit');
    if selection == 0
        return
    end
else
    selection = preselectVector(2);
end
experimentFolder = [experimentFolder filesep fileList{selection}];

%% create output folder
[a,b,c] = mkdir(experimentFolder, 'analyse');
analyseFolder = [experimentFolder filesep 'analyse'];

%% Find files to use
disp(['# Finding files in folder: ' experimentFolder])
fileList = dir(experimentFolder);
filterVector = (1:length(fileList)) == 0;
for iFile = 1:length(fileList)
    if ~fileList(iFile).isdir && strcmp(fileList(iFile).name((end-3):end), '.csv')
        filterVector(iFile) = 1;
    end
end
fileList = {fileList(filterVector).name};
disp(['### ' num2str(length(fileList)) ' files found, now analyzing']);

%% Figure
hFigure = figure(42);
set(42, 'InvertHardCopy', 'off');
hAxisImage = subplot(1,2,1);%axis off
hAxisFit = subplot(1,2,2);
set(hFigure, 'Position', [0 800 1280 600])

%% main loop
for iFile = 1:length(fileList)
    %% read data
    fileID = fopen([experimentFolder filesep fileList{iFile}]);
    scanData = textscan(fileID, '%u', 'delimiter', ',');
    scanData = scanData{1}';
    fclose(fileID);
    %% transform to seperate particles
    zeroIndexVector = [0 find(scanData==0)];
    nZero = length(zeroIndexVector);
    clear boundaryData
    if nZero > 2
        boundaryData(nZero-2).cartX = [];
        boundaryData(nZero-2).cartY = [];
        for iZero = 1:(nZero-2)
            boundaryData(iZero).cartX = double(scanData(...
                (zeroIndexVector(iZero)+1+2):2:(zeroIndexVector(iZero+1)-1)));
            boundaryData(iZero).cartY = double(scanData(...
                (zeroIndexVector(iZero)+1+3):2:(zeroIndexVector(iZero+1)-1)));
        end
        disp(['# ' num2str(length(boundaryData)) ' curves found in ' fileList{iFile}])
    else
        disp('# No boundary curve found')
        continue
    end
    %% load image
    thisImage = imread([experimentFolder filesep fileList{iFile}(1:(end-4)) '.png'], 'png');
    
    for iBoundary = 1:length(boundaryData)
        disp(['### Particle: ' fileList{iFile} ' number ' num2str(iBoundary)]);
        % calculate
        tic
        thisData.bnd = obtainBnd([...
            boundaryData(iBoundary).cartX; boundaryData(iBoundary).cartY], ...
            'cart', thisImage);
        thisData.fit = fitShapeBnd(thisData.bnd, ...
                'minConvexity', primaryLimitStruct.convexity, ...
                'minDiameter', primaryLimitStruct.diameter, ...
                'maxSphericity', primaryLimitStruct.sphericity, ...
                'maxSphericity2', primaryLimitStruct.sphericity2, ...
                'maxAngleGap', primaryLimitStruct.angleGap/360*2*pi);
        thisTime = toc
               
        % evaluation data
        thisEvalData.sphericity = thisData.fit.sphericity;
        thisEvalData.sphericity2 = thisData.fit.sphericity2;
        thisEvalData.sphericity3 = thisData.fit.sphericity3;
        thisEvalData.diameterEstimate = thisData.bnd.meanWidthHull;
        thisEvalData.angleGapDegree = thisData.bnd.largestAngleGap * 360/2/pi;
        thisEvalData.convexity = thisData.bnd.convexity;
        thisEvalData.meanWidth = thisData.fit.meanWidth;
        thisEvalData.lambda = thisData.fit.lambda;
%         thisEvalData.roundness = (2*thisEvalData.lambda)/thisEvalData.meanWidth;
        thisEvalData.roughness = thisData.fit.roughness;
        thisEvalData.squareRoughness = thisData.fit.squareRoughness;
        thisEvalData.relRoughness = thisEvalData.roughness / thisEvalData.meanWidth;
        thisEvalData.roundness = 2*thisEvalData.lambda / thisEvalData.meanWidth;
        thisEvalData.brightFraction = thisData.fit.brightFraction;
            
        % evaluate (yes/no) primary
        if thisData.fit.flagIgnore
            thisEvalData.ignorePrimary = 1;
            disp('Particle is ignored PRIMARILY')
        else
            thisEvalData.ignorePrimary = 0;
        end
        
        % get image
        thisData.image = thisImage(...
            max(boundaryData(iBoundary).cartY):-1:min(boundaryData(iBoundary).cartY), ...
            min(boundaryData(iBoundary).cartX):max(boundaryData(iBoundary).cartX));
        
        % evaluate (yes/no) secondary
        thisEvalData.ignoreSecondary = 0;
        if thisEvalData.meanWidth < secondaryLimitStruct.minSize
            disp('Particle would be ignored according to MIN SIZE')
            thisEvalData.ignoreSecondary = 1;
        end
%         if thisEvalData.meanWidth > secondaryLimitStruct.maxSize
%             disp('Particle would be ignored according to MAX SIZE')
%             thisEvalData.ignoreSecondary = 1;
%         end
%         if thisEvalData.convexity < secondaryLimitStruct.convexity
%             disp('Particle would be ignored according to CONVEXITY')
%             thisEvalData.ignoreSecondary = 1;
%         end
%         if thisEvalData.sphericity < secondaryLimitStruct.sphericity
%             disp('Particle would be ignored according to SPHERICITY')
%             thisEvalData.ignoreSecondary = 1;
%         end
        if thisEvalData.relRoughness > secondaryLimitStruct.relRoughness
            disp('Particle would be ignored according to RELATIVE ROUGHNESS')
            thisEvalData.ignoreSecondary = 1;
        end
        if thisEvalData.roughness > secondaryLimitStruct.roughness
            disp('Particle would be ignored according to ROUGHNESS')
            thisEvalData.ignoreSecondary = 1;
        end
%         if thisEvalData.roundness > secondaryLimitStruct.maxRoundness
%             disp('Particle would be ignored according to MAX ROUNDNESS')
%             thisEvalData.ignoreSecondary = 1;
%         end
%         if thisEvalData.sphericity < secondaryLimitStruct.sphericity
%             disp('Particle would be ignored according to SPHERICITY-1')
% %             thisEvalData.ignoreSecondary = 1;
%         end
%         if thisEvalData.sphericity2 < secondaryLimitStruct.sphericity2
%             disp('Particle would be ignored according to SPHERICITY-2')
% %             thisEvalData.ignoreSecondary = 1;
%         end
%         if thisEvalData.sphericity3 > secondaryLimitStruct.sphericity3
%             disp('Particle would be ignored according to SPHERICITY-3')
%             thisEvalData.ignoreSecondary = 1;
%         end
        if thisEvalData.brightFraction < secondaryLimitStruct.minBrightFraction
            disp('Particle would be ignored according to BRIGHT-FRACTION')
            thisEvalData.ignoreSecondary = 1;
        end
                  
        %% plot
        if (~thisEvalData.ignorePrimary && ~thisEvalData.ignoreSecondary) || ...
                (thisEvalData.ignorePrimary && flagViewPrimaryIgnored) || ...
                (thisEvalData.ignoreSecondary && flagViewSecondaryIgnored && ~thisEvalData.ignorePrimary)
        % view fit
        axes(hAxisFit),cla,hold on
        viewBoundary(thisData.bnd);
        viewBoundaryShape(thisData.fit);
%         title(['round=' num2str(2*thisEvalData.lambda / ...
%                                 thisEvalData.meanWidth)])
        title(['sph=[' num2str(thisEvalData.sphericity) ' ' ...
                       num2str(thisEvalData.sphericity2) ' ' ...
                       num2str(thisEvalData.sphericity3) ']'])
        title(['BRF/ PH3 = ' num2str(thisEvalData.brightFraction) '/' ...
                             num2str(thisEvalData.sphericity3)])  
        switch thisData.fit.flagIgnore
            case 6
                title('ignored by BRIGHTFRACTION')
            case 2
                title('ignored by SIZE')
        end
                   
        % view image
        axes(hAxisImage),cla,hold on
        imshow(thisData.image)
%         title(['lam=' num2str(thisEvalData.lambda)])
%         title(['frac=' num2str(thisEvalData.brightFraction)])
        if thisEvalData.ignorePrimary
            set(hFigure, 'Color', [1 0.6 0.6])
%             axes(hAxisFit),title('');
        elseif thisEvalData.ignoreSecondary
            set(hFigure, 'Color', [1 0.8 0.8])
%             axes(hAxisFit),title('');
        else
            set(hFigure, 'Color', [0.8 1 0.8])
            axes(hAxisFit)
%             title(['w=' num2str(thisEvalData.meanWidth) ', ' ...
%                    'lam=' num2str(thisEvalData.lambda) ', ' ...
%                    'r=' num2str(thisEvalData.roundness)]);
        end
        end
        
        %% save
        if (~thisEvalData.ignorePrimary && ~thisEvalData.ignoreSecondary) || ...
                (thisEvalData.ignorePrimary && flagSavePrimaryIgnored) || ...
                (thisEvalData.ignoreSecondary && flagSaveSecondaryIgnored && ~thisEvalData.ignorePrimary)
            
        print('-dpng', [analyseFolder filesep ...
            fileList{iFile}(1:(end-4)) '--' num2str(iBoundary)]);
        
        fileID = fopen([analyseFolder filesep ...
            fileList{iFile}(1:(end-4)) '--' num2str(iBoundary) '.txt'], 'w');
        fwrite(fileID, evalc('disp(thisEvalData)'));
        fwrite(fileID, ['roundness: ' num2str(...
            2*thisEvalData.lambda / thisEvalData.meanWidth) 10 10] );
        fwrite(fileID, evalc('disp(primaryLimitStruct)'));
        fwrite(fileID, evalc('disp(secondaryLimitStruct)'));
        fclose(fileID);
%         pause
        end
%         pause
    end
    
end