function [FieldAssignmentsAll, OrganoidStencilAll, Wells] = Iris_GetOrganoidFields(InPath, MetaData, CoordinatesAll, illustrate)
%extract organoid IDs from CV8000 first pass and second pass integration
%   author: Paul Antony 20210213
%   Coordinates: example: readtable([InPath, filesep, 'Encoding.csv'], 'ReadVariableNames', true);
%   illustrate: true or false
%   use case: /mnt/IrisHCS/Library/HcsProjectCode/SemraSmajic/Organoid_Neurons_Dapi_488Map2_546Vimentin_647Olig2

FieldAssignmentsAll = {};
OrganoidStencilAll = {};
Wells = unique(MetaData.InfoTable{:}.Well);
for W = 1:numel(Wells) 
    WellsTimeLineDecoded = MetaData.TimeLines{1,1}.Wells(:,:);
    WellsTimeLineDecoded = rowfun(@(a,b) sprintf('%s%0.2d', char('A' + a - 1), b), WellsTimeLineDecoded, 'OutputFormat', 'cell', 'ExtractCellContents', true);
    TimeLineThis = find(strcmp(Wells{W}, WellsTimeLineDecoded));
    overViewPath = [InPath, filesep, Wells{W}, '.png'];
    Coordinates = CoordinatesAll(strcmp(CoordinatesAll.Well, Wells{W}),:);
    if illustrate
        overview = imread(overViewPath); % imtool(overview,[])
        for c = 1:height(Coordinates)
            if mod(c,100) == 0
                overview = insertText(overview, [Coordinates{c, 'FieldX'}, Coordinates{c, 'FieldY'}], Coordinates{c, 'FieldID'});
            end
        end
        imtool(overview)

        figure;
        scatter(MetaData.TimeLines{1,1}.Fields{TimeLineThis,1}.X, MetaData.TimeLines{1,1}.Fields{TimeLineThis,1}.Y, '.r')
    end

    %% Compare first pass proposed coordinates with second pass used coordinates
    try
        CoordinatesProposedSecondPass = readtable([InPath, filesep, 'Encoding.csv'], 'ReadVariableNames', true);
    catch
        CoordinatesProposedSecondPass = readtable(['.', filesep, 'Encoding.csv'], 'ReadVariableNames', true);
    end
    
    CoordinatesProposedSecondPass = CoordinatesProposedSecondPass(strcmp(CoordinatesProposedSecondPass.Well, Wells{W}),:);
    CoordinatesProposedSecondPass = [CoordinatesProposedSecondPass.FieldX, CoordinatesProposedSecondPass.FieldY];
    CoordinatesProposedSecondPassNorm = [mat2gray(CoordinatesProposedSecondPass(:,1)), mat2gray(-CoordinatesProposedSecondPass(:,2))];
    CoordinatesUsedSecondPass = [MetaData.TimeLines{1,1}.Fields{TimeLineThis,1}.X, MetaData.TimeLines{1,1}.Fields{TimeLineThis,1}.Y, MetaData.TimeLines{1,1}.Fields{TimeLineThis,1}.Field];
    CoordinatesUsedSecondPassNorm = [mat2gray(CoordinatesUsedSecondPass(:,2)), mat2gray(CoordinatesUsedSecondPass(:,1))];
    if illustrate
        figure
        scatter(CoordinatesProposedSecondPassNorm(:,1), CoordinatesProposedSecondPassNorm(:,2), '.', 'MarkerFaceColor', [0.5, 0.5, 0.5], 'MarkerEdgeColor', [0.5, 0.5, 0.5])
        hold on
        scatter(CoordinatesUsedSecondPassNorm(:,2), CoordinatesUsedSecondPassNorm(:,1), '.g')
    end

    %% Create Organoid stencil and fit proposed coordinates to used coordinates
    try
        SeedData = readtable([InPath, filesep, 'Encoding.csv'], 'ReadVariableNames', true);
    catch
        SeedData = readtable(['.', filesep, 'Encoding.csv'], 'ReadVariableNames', true);
    end
    
    SeedData = SeedData(strcmp(SeedData.Well, Wells{W}), :);
    colRange = max(SeedData.FieldX) - min(SeedData.FieldX);
    rowRange = max(SeedData.FieldY) - min(SeedData.FieldY);
    SeedData.colNorm = round(colRange * CoordinatesProposedSecondPassNorm(:,1));
    SeedData.rowNorm = round(rowRange * CoordinatesProposedSecondPassNorm(:,2));
    OrganoidStencil = zeros([rowRange + 1, colRange + 1], 'uint16');
    for i = 1:height(SeedData)
        SeedThis = SeedData(i,:);
        OrganoidStencil(SeedThis.rowNorm + 1, SeedThis.colNorm + 1) = SeedThis.OrganoidID;
    end

    % dilate dynamically to fit count of organoids
    OrganoidCount = max(SeedData.OrganoidID);
    while true
        OrganoidStencil = imdilate(OrganoidStencil, strel('disk', 1));
        [~, countCCs] = bwlabeln(OrganoidStencil);
        %imtool(OrganoidStencil, [])
        if countCCs <= OrganoidCount
            break
        end
    end

    OrganoidStencil = imfill(OrganoidStencil, 'holes');
    %imtool(OrganoidStencil,[])

    % Encode for gray tone dilation fitting
    OrganoidStencil(OrganoidStencil > 0) = OrganoidStencil(OrganoidStencil > 0) + 1;

    % add used second pass coordinates
    SecondPassDecoded = table();
    SecondPassDecoded.rowNorm = round(rowRange * CoordinatesUsedSecondPassNorm(:,1));
    SecondPassDecoded.colNorm = round(colRange * CoordinatesUsedSecondPassNorm(:,2));
    for i = 1:height(SecondPassDecoded)
        SeedThis = SecondPassDecoded(i,:);
        OrganoidStencil(SeedThis.rowNorm + 1, SeedThis.colNorm + 1) = 1; % lowest non zero
    end
    %imtool(OrganoidStencil,[])

    % gray tone dilate dynamically
    while true
        OrganoidStencil = imdilate(OrganoidStencil, strel('disk', 1));
        [~, countCCs] = bwlabeln(OrganoidStencil);
        for i = 1:height(SecondPassDecoded)
            SeedThis = SecondPassDecoded(i,:);
            OrganoidStencil(SeedThis.rowNorm + 1, SeedThis.colNorm + 1) = 1; % lowest non zero
        end
        %imtool(OrganoidStencil, [])
        if countCCs <= OrganoidCount
            break
        end
    end

    [OrganoidStencilGoast, countCCs] = bwlabeln(OrganoidStencil);
    OrganoidStencilRecon = OrganoidStencil;
    for i = 1:countCCs
        OrganoidGoastThisVec = OrganoidStencil(OrganoidStencilGoast == i);
        OrganoidStencilRecon(OrganoidStencilGoast == i) = max(OrganoidGoastThisVec);
    end

    %imtool(flip(OrganoidStencil,1),[]) % scatter uses y with origin bottom left
    %imtool(flip(OrganoidStencilRecon,1),[]) % scatter uses y with origin bottom left

    % Decode for gray tone dilation fitting
    OrganoidStencil = OrganoidStencilRecon;
    OrganoidStencil(OrganoidStencil > 0) = OrganoidStencil(OrganoidStencil > 0) - 1;

    % Join
    FieldAssignments = array2table(CoordinatesUsedSecondPass, 'VariableNames', {'X','Y', 'Field'});
    FieldAssignments.rowNorm = SecondPassDecoded.rowNorm;
    FieldAssignments.colNorm = SecondPassDecoded.colNorm;
    OrganoidStencilQC = OrganoidStencil;
    for i = 1:height(FieldAssignments)
        FieldAssignments(i, 'OrganoidID') = {OrganoidStencil(FieldAssignments.rowNorm(i)+1, FieldAssignments.colNorm(i)+1)};
        OrganoidStencilQC(FieldAssignments.rowNorm(i)+1, FieldAssignments.colNorm(i)+1) = OrganoidCount + 1;
    end
    % imtool(flip(OrganoidStencilQC, 1),[])
    FieldAssignmentsAll{W} =  FieldAssignments;
    OrganoidStencilAll{W} = OrganoidStencil;
end

end

