%% User inputs

% /mnt/irisgpfs/projects/lcsb_hcs/Library/HcsProjectCode/GemmaGomez/Organoids_SNCA_TH
%/mnt/irisgpfs/projects/lcsb_hcs/Library/HcsProjectCode/DanielaVega/Organoids_SNCA_TH_GPU_Hindbrain/

try
    delete(gcp('nocreate'))
    myCluster = parcluster('local');
    AvailableWorkers = myCluster.NumWorkers;
    if AvailableWorkers >= 7
        pool = parpool(7)
    else
        pool = parpool(1)
    end
catch ParError
    disp('ParError')
end

addpath(genpath('/work/projects/lcsb_hcs/Library/hcsforge'))
addpath(genpath('/work/projects/lcsb_hcs/Library/hcsIris'))


%if ~exist('InPath') % otherwise use the one defined via the database
    %InPath = '/mnt/AtlasHCS/YokogawaCV8000Standalone/BTSData/MeasurementData/danielavega/DV_20211122_HBO_EXP8_9_60DOD_EXP9_30DOD_3SlidesRS_map2psynsyn6030_20211122_170007/11-22-21_17-06-50'%to take data from Atlas
    InPath = '/mnt/irisgpfs/projects/lcsb_hcs/Data/GemmaGomez/GG_20241116_Rescan_MO_S1_A34_mos1a34_20241116_122057/11-16-24_12-26-12_in';    
    %InPath = '/work/projects/lcsb_hcs/Data/DanielaVega/Organoids_SNCA_TH_GPU_Hindbrain/DV_20211122_HBO_EXP8_9_60DOD_EXP9_30DOD_3SlidesRS_map2psynsyn6030_20211122_170007_in';
    %InPath = '/mnt/irisgpfs/projects/lcsb_hcs/Data/GemmaGomez/GG_20230112_Rescan_A4_MOs3060_ASM_e1819_SNCA_PSNCA_mo3060asme18e19p_20230112_150606/01-12-23_15-09-34_in';
    %InPath = '/mnt/lscratch/users/ggomez/HCSdata/01-12-23_15-09-34_in';
    
%end

MesPath = ls([InPath, '/*.mes']); MesPath = MesPath(1:end-1); % remove line break
MetaData = f_CV8000_getChannelInfo(InPath, MesPath);

%if ~exist('OutPath') % otherwise use the one defined via the database
    %OutPath = '/mnt/irisgpfs/projects/lcsb_hcs/Data/GemmaGomez/GG_20240426_Rescan_A34_MOs60_ASM_e1821_SNCA_PSNCA_moasmsncapsnca_20230426_111454/04-26-23_11-43-36_out';
    %OutPath = '/work/projects/lcsb_hcs/Data/DanielaVega/Organoids_SNCA_TH_GPU_Hindbrain/DV_20211122_HBO_EXP8_9_60DOD_EXP9_30DOD_3SlidesRS_map2psynsyn6030_20211122_170007_out'
    %OutPath = '/mnt/lscracth/projects/lcsb_hcs/Data/GemmaGomez/GG_20230112_Rescan_A4_MOs3060_ASM_e1819_SNCA_PSNCA_mo3060asme18e19p_20230112_150606/01-12-23_15-09-34_out';
    OutPath = '/mnt/lscratch/users/ggomez/CONNECT_revisions/HBOs_characterization/GG_20241116_Rescan_MO_S1_A34_mos1a34_20241116_122057/11-16-24_12-26-12_out';
%end

f_LogDependenciesLinux(mfilename, OutPath)
%FileNameShort = mfilename;
%f_LogDependenciesLinux(FileNameShort, OutPath)


%% Prepare folders
mkdir(OutPath)
ThumbnailPath = [OutPath, filesep, 'Thumbnails'];
mkdir(ThumbnailPath)
PreviewPath = [OutPath, filesep, 'Previews'];
mkdir(PreviewPath)

%% Load Metadata
MetaData = f_CV8000_getChannelInfo(InPath, MesPath);

%% Analyze selected organoids
InfoTable = MetaData.InfoTable{:};
PlaneCount = max(cell2mat(cellfun(@(x) str2double(x), InfoTable.Plane, 'UniformOutput', false)));
ChannelCount = max(cell2mat(cellfun(@(x) str2double(x), InfoTable.Channel, 'UniformOutput', false)));
% Read SlideLayout
try
    OrganoidsToAnalyze = readtable([InPath, filesep, 'SlideLayout.csv'], 'ReadVariableNames', false);
    %OrganoidsToAnalyze = readtable([InPath, filesep, 'SlideLayout.csv'], 'ReadVariableNames', false, 'delimiter', ';');
catch
    %OrganoidsToAnalyze = readtable(['.', filesep, 'SlideLayout.csv'], 'ReadVariableNames', false);
    OrganoidsToAnalyze = readtable(['.', filesep, 'SlideLayout.csv'], 'ReadVariableNames', false, 'delimiter', ';');
end
OrganoidsToAnalyze.Properties.VariableNames = {'Idx', 'Well', 'AreaName'};
OrganoidsToAnalyze.OrganoidID = rowfun(@(a,b) sprintf('%s_%s', a, b), OrganoidsToAnalyze, 'InputVariables', {'Well', 'Idx'}, 'ExtractCellContents', true, 'OutputFormat', 'cell');

try
    Coordinates = readtable([InPath, filesep, 'Encoding.csv'], 'ReadVariableNames', true);
catch
    Coordinates = readtable(['.', filesep, 'Encoding.csv'], 'ReadVariableNames', true);
end

    
%[FieldAssignments, OrganoidStencil] = Iris_GetOrganoidFields(InPath, MetaData, overviewSlide, Coordinates, false);
[FieldAssignmentsAll, OrganoidStencilsAll, WellsAll] = Iris_GetOrganoidFields(InPath, MetaData, Coordinates, false);
%[FieldAssignments, OrganoidStencils, Wells] = Iris_GetOrganoidFields(InPath, MetaData, Coordinates, true);
% imtool(flip(OrganoidStencil,1),[])
% imtool(flip(OrganoidStencils{1},1),[])
% imtool(flip(OrganoidStencils{2},1),[])

ObjectsAll = {};
for organoid = 1:size(OrganoidsToAnalyze,1)
%for organoid = 34  
    try
        tic
        disp(['Well ', OrganoidsToAnalyze{organoid, 'Well'}{:}, '__Organoid #' num2str(organoid)])
        OrganoidThis = OrganoidsToAnalyze(organoid, :);
        SlideThis = OrganoidThis.Well{:};
        Sample = OrganoidThis.AreaName{:};
        ID = OrganoidThis.Idx;

        %% coordinate slide specific data
        ThisWellString = OrganoidsToAnalyze{organoid, 'Well'}{:};
        SlideIdx = find(strcmp(WellsAll, ThisWellString));
        FieldAssignments = FieldAssignmentsAll{SlideIdx};
        OrganoidStencils = OrganoidStencilsAll{SlideIdx};

        %% Load organoid rescan image mosaic (optimized for speed)
        ImSize = size(imread(InfoTable{1,'file'}{:}));
        %FieldsThisOrganoid = FieldAssignments(FieldAssignments.OrganoidID == organoid, :);
        OrganoidIDThisSlide = OrganoidsToAnalyze{organoid,'Idx'}{:};
        OrganoidIDThisSlide = regexp(OrganoidIDThisSlide, 'OR(.*)', 'tokens')
        OrganoidIDThisSlide =  str2double(OrganoidIDThisSlide{:}{:});
        FieldsThisOrganoid = FieldAssignments(FieldAssignments.OrganoidID == OrganoidIDThisSlide, :);
        XVecThisOrganoid = sort(unique(FieldsThisOrganoid.X));
        YVecThisOrganoid = sort(unique(FieldsThisOrganoid.Y));
        ch1Meta=table();
        ch2Meta=table();
        ch3Meta=table();
        ch4Meta=table();
        ch5Meta=table();
        progressCh1 = 0;
        progressCh2 = 0;
        progressCh3 = 0;
        progressCh4 = 0;
        progressCh5 = 0;

        for ch = 1:ChannelCount
            for field = 1:height(FieldsThisOrganoid)
                for plane = 1:PlaneCount
                    FieldThis = FieldsThisOrganoid.Field(field);
                    ThisTilePath = sortrows(InfoTable(strcmp(InfoTable.Well, ThisWellString) & strcmp(InfoTable.Field, sprintf('%03d', FieldThis)) & strcmp(InfoTable.Plane, sprintf('%02d', plane)), :), 'Channel');
                    StartRowThisTile =    ((length(YVecThisOrganoid) - (find(YVecThisOrganoid == FieldsThisOrganoid.Y(field)))) * ImSize(1)) + 1;
                    StartColumnThisTile = (((find(XVecThisOrganoid == FieldsThisOrganoid.X(field))) -1) * ImSize(2)) + 1;

                    switch ch
                        case 1
                            progressCh1 = progressCh1 + 1;
                            ch1Meta(progressCh1, 'StartRowThisTile') = {StartRowThisTile};
                            ch1Meta(progressCh1, 'StartColumnThisTile') = {StartColumnThisTile};
                            ch1Meta(progressCh1, 'Plane') = {plane};
                            ch1Meta(progressCh1, 'Field') = {field};
                            ch1Meta(progressCh1, 'ThisTilePath') = ThisTilePath.file(1);
                        case 2
                            progressCh2 = progressCh2 + 1;
                            ch2Meta(progressCh2, 'StartRowThisTile') = {StartRowThisTile};
                            ch2Meta(progressCh2, 'StartColumnThisTile') = {StartColumnThisTile};
                            ch2Meta(progressCh2, 'Plane') = {plane};
                            ch2Meta(progressCh2, 'Field') = {field};
                            ch2Meta(progressCh2, 'ThisTilePath') = ThisTilePath.file(2);
                        case 3
                            progressCh3 = progressCh3 + 1;
                            ch3Meta(progressCh3, 'StartRowThisTile') = {StartRowThisTile};
                            ch3Meta(progressCh3, 'StartColumnThisTile') = {StartColumnThisTile};
                            ch3Meta(progressCh3, 'Plane') = {plane};
                            ch3Meta(progressCh3, 'Field') = {field};
                            ch3Meta(progressCh3, 'ThisTilePath') = ThisTilePath.file(3);
                        case 4
                            progressCh4 = progressCh4 + 1;
                            ch4Meta(progressCh4, 'StartRowThisTile') = {StartRowThisTile};
                            ch4Meta(progressCh4, 'StartColumnThisTile') = {StartColumnThisTile};
                            ch4Meta(progressCh4, 'Plane') = {plane};
                            ch4Meta(progressCh4, 'Field') = {field};
                            ch4Meta(progressCh4, 'ThisTilePath') = ThisTilePath.file(4);
                        case 5
                            progressCh5 = progressCh5 + 1;
                            ch5Meta(progressCh5, 'StartRowThisTile') = {StartRowThisTile};
                            ch5Meta(progressCh5, 'StartColumnThisTile') = {StartColumnThisTile};
                            ch5Meta(progressCh5, 'Plane') = {plane};
                            ch5Meta(progressCh5, 'Field') = {field};
                            ch5Meta(progressCh5, 'ThisTilePath') = ThisTilePath.file(5);
                    end

                end
            end
        end
        ch1Ims = {};
        ch2Ims = {};
        ch3Ims = {};
        ch4Ims = {};
        ch5Ims = {};

        parfor chi = 1:(height(ch1Meta))
            ch1Ims{chi} = imread(ch1Meta.ThisTilePath{chi});
            ch2Ims{chi} = imread(ch2Meta.ThisTilePath{chi});
            ch3Ims{chi} = imread(ch3Meta.ThisTilePath{chi});
            ch4Ims{chi} = imread(ch4Meta.ThisTilePath{chi});
            %ch5Ims{chi} = imread(ch5Meta.ThisTilePath{chi});
        end

        Planes = numel(unique(ch1Meta.Plane));
        % Nuclei (blue) Hoechst
        ch1 = zeros(length(YVecThisOrganoid) * ImSize(1), length(XVecThisOrganoid) * ImSize(2), Planes, 'uint16');
        for i = 1:height(ch1Meta)
            ch1(ch1Meta{i, 'StartRowThisTile'}:ch1Meta{i, 'StartRowThisTile'} + ImSize(1)-1, ch1Meta{i, 'StartColumnThisTile'}:ch1Meta{i, 'StartColumnThisTile'} + ImSize(2)-1, ch1Meta{i, 'Plane'}) = ch1Ims{i};
        end
        %vol(ch1, 0, 1000)
        %imtool(max(ch1, [], 3),[])
        
        % Map2 (green)
        ch2 = zeros(length(YVecThisOrganoid) * ImSize(1), length(XVecThisOrganoid) * ImSize(2), Planes, 'uint16');
        for i = 1:height(ch1Meta)
            ch2(ch2Meta{i, 'StartRowThisTile'}:ch2Meta{i, 'StartRowThisTile'} + ImSize(1)-1, ch2Meta{i, 'StartColumnThisTile'}:ch2Meta{i, 'StartColumnThisTile'} + ImSize(2)-1, ch2Meta{i, 'Plane'}) = ch2Ims{i};
        end
        %vol(ch2, 0, 500)
        %imtool(max(ch2, [], 3),[])
        
        % pSynuclein (red)
        ch3 = zeros(length(YVecThisOrganoid) * ImSize(1), length(XVecThisOrganoid) * ImSize(2), Planes, 'uint16');
        for i = 1:height(ch1Meta)
            ch3(ch3Meta{i, 'StartRowThisTile'}:ch3Meta{i, 'StartRowThisTile'} + ImSize(1)-1, ch3Meta{i, 'StartColumnThisTile'}:ch3Meta{i, 'StartColumnThisTile'} + ImSize(2)-1, ch3Meta{i, 'Plane'}) = ch3Ims{i};
        end
        %vol(ch3, 0, 50000)
        %imtool(max(ch3, [], 3),[])
        
        % Synuclein (deep red)
        ch4 = zeros(length(YVecThisOrganoid) * ImSize(1), length(XVecThisOrganoid) * ImSize(2), Planes, 'uint16');
        for i = 1:height(ch1Meta)
            ch4(ch4Meta{i, 'StartRowThisTile'}:ch4Meta{i, 'StartRowThisTile'} + ImSize(1)-1, ch4Meta{i, 'StartColumnThisTile'}:ch4Meta{i, 'StartColumnThisTile'} + ImSize(2)-1, ch4Meta{i, 'Plane'}) = ch4Ims{i};
        end
        %vol(ch4, 0, 500)
        %imtool(max(ch4, [], 3),[])
        
        if false
            ch1 = ch1(5000:7000, 3000:5000,:);% vol(ch1f, 0, 2000) it(max(ch1f,[],3))
            ch2 = ch2(5000:7000, 3000:5000,:);% vol(ch1f, 0, 2000)
            ch3 = ch3(5000:7000, 3000:5000,:);% vol(ch1f, 0, 2000)
            ch4 = ch4(5000:7000, 3000:5000,:);% vol(ch1f, 0, 2000)
        end

        OrganoidLoadTime = toc;
        disp(['Needed ', num2str(OrganoidLoadTime), ' seconds to load organoid images'])

        HandleOverfitting = false;
        %HandleOverfitting = true;
        if HandleOverfitting
            if organoid == 1
                ch1_i_1 = ch1;
                ch2_i_1 = ch2;
                ch3_i_1 = ch3;
                ch4_i_1 = ch4;
                
                chFindBestPlanes_i_1 = uint32(ch1_i_1)+ uint32(ch2_i_1) + uint32(ch3_i_1) + uint32(ch4_i_1); % no BF as too bright
                chFindBestPlanesSummary_i_1 = squeeze(sum(chFindBestPlanes_i_1,[1,2])); % plane sum intensity vector
                %figure; plot(chFindBestPlanesSummary)
                GoodPlanes_i_1 = chFindBestPlanesSummary_i_1 > (max(chFindBestPlanesSummary_i_1)/3); % Above 33% of maximum
                %GoodPlanes = chFindBestPlanesSummary > (max(chFindBestPlanesSummary)/10); % Above 33% of maximum
                GoodPlanes_i_1 = bwareafilt(GoodPlanes_i_1, 1);
                clear 'chFindBestPlanes'
                ch1_i_1 = ch1_i_1(:,:,GoodPlanes_i_1); % Hoechst % vol(ch1, 0, 10000) % imtool(max())
                ch2_i_1 = ch2_i_1(:,:,GoodPlanes_i_1); % TH
                ch3_i_1 = ch3_i_1(:,:,GoodPlanes_i_1); % Map2 % imtool(max(ch3_i_1,[],3))
                ch4_i_1 = ch4_i_1(:,:,GoodPlanes_i_1); % Synuclein
                
                ch1_i_1c = ch1_i_1(3000:7000, 500:3500, 1:4); %imtool(max(ch1_i_1c,[],3), [])
                ch2_i_1c = ch2_i_1(3000:7000, 500:3500, 1:4); %imtool(max(ch2_i_1c,[],3), [])
                ch3_i_1c = ch3_i_1(3000:7000, 500:3500, 1:4); %imtool(max(ch3_i_1c,[],3), [])
                ch4_i_1c = ch4_i_1(3000:7000, 500:3500, 1:4); %imtool(max(ch4_i_1c,[],3), [])
                
            end
            
            if organoid == 4
                ch1_i_4 = ch1;
                ch2_i_4 = ch2;
                ch3_i_4 = ch3;
                ch4_i_4 = ch4;
                
                chFindBestPlanes_i_4 = uint32(ch1_i_4)+ uint32(ch2_i_4) + uint32(ch3_i_4) + uint32(ch4_i_4); % no BF as too bright
                chFindBestPlanesSummary_i_4 = squeeze(sum(chFindBestPlanes_i_4,[1,2])); % plane sum intensity vector
                %figure; plot(chFindBestPlanesSummary)
                GoodPlanes_i_4 = chFindBestPlanesSummary_i_4 > (max(chFindBestPlanesSummary_i_4)/3); % Above 33% of maximum
                %GoodPlanes = chFindBestPlanesSummary > (max(chFindBestPlanesSummary)/10); % Above 33% of maximum
                GoodPlanes_i_4 = bwareafilt(GoodPlanes_i_4, 1);
                clear 'chFindBestPlanes'
                ch1_i_4 = ch1_i_4(:,:,GoodPlanes_i_4); % Hoechst % vol(ch1, 0, 10000)
                ch2_i_4 = ch2_i_4(:,:,GoodPlanes_i_4); % TH
                ch3_i_4 = ch3_i_4(:,:,GoodPlanes_i_4); % Map2 % imtool(max(ch3_i_4,[],3))
                ch4_i_4 = ch4_i_4(:,:,GoodPlanes_i_4); % Synuclein
                
                ch1_i_4c = ch1_i_4(2500:6500, 600:3600, 1:4); %imtool(max(ch1_i_4c,[],3), [])
                ch2_i_4c = ch2_i_4(2500:6500, 600:3600, 1:4); %imtool(max(ch2_i_4c,[],3), [])
                ch3_i_4c = ch3_i_4(2500:6500, 600:3600, 1:4); %imtool(max(ch3_i_4c,[],3), [])
                ch4_i_4c = ch4_i_4(2500:6500, 600:3600, 1:4); %imtool(max(ch4_i_4c,[],3), [])

            end
            
            if organoid == 5
                ch1_i_5 = ch1;
                ch2_i_5 = ch2;
                ch3_i_5 = ch3;
                ch4_i_5 = ch4;
                
                chFindBestPlanes_i_5 = uint32(ch1_i_5)+ uint32(ch2_i_5) + uint32(ch3_i_5) + uint32(ch4_i_5); % no BF as too bright
                chFindBestPlanesSummary_i_5 = squeeze(sum(chFindBestPlanes_i_5,[1,2])); % plane sum intensity vector
                %figure; plot(chFindBestPlanesSummary)
                GoodPlanes_i_5 = chFindBestPlanesSummary_i_5 > (max(chFindBestPlanesSummary_i_5)/3); % Above 33% of maximum
                %GoodPlanes = chFindBestPlanesSummary > (max(chFindBestPlanesSummary)/10); % Above 33% of maximum
                GoodPlanes_i_5 = bwareafilt(GoodPlanes_i_5, 1);
                clear 'chFindBestPlanes'
                ch1_i_5 = ch1_i_5(:,:,GoodPlanes_i_5); % Hoechst % vol(ch1, 0, 10000)
                ch2_i_5 = ch2_i_5(:,:,GoodPlanes_i_5); % TH
                ch3_i_5 = ch3_i_5(:,:,GoodPlanes_i_5); % Map2% imtool(max(ch3_i_5,[],3))
                ch4_i_5 = ch4_i_5(:,:,GoodPlanes_i_5); % Synuclein
                
                ch1_i_5c = ch1_i_5(3000:7000, 900:3900, 2:5); %imtool(max(ch1_i_5c,[],3), [])
                ch2_i_5c = ch2_i_5(3000:7000, 900:3900, 2:5); %imtool(max(ch2_i_5c,[],3), [])
                ch3_i_5c = ch3_i_5(3000:7000, 900:3900, 2:5); %imtool(max(ch3_i_5c,[],3), [])
                ch4_i_5c = ch4_i_5(3000:7000, 900:3900, 2:5); %imtool(max(ch4_i_5c,[],3), [])
            end
            ch1 = cat(2, ch1_i_1c, ch1_i_4c, ch1_i_5c); % Hoechst % vol(ch1, 0, 10000)
            ch2 = cat(2, ch2_i_1c, ch2_i_4c, ch2_i_5c); % TH
            ch3 = cat(2, ch3_i_1c, ch3_i_4c, ch3_i_5c); % Map2
            ch4 = cat(2, ch4_i_1c, ch4_i_4c, ch4_i_5c); % Synuclein
            
        end

        clear 'ch1Ims' 'ch2Ims' 'ch3Ims' 'ch4Ims' 'ch5Ims'
        chFindBestPlanes = uint32(ch1)+ uint32(ch2) + uint32(ch3) + uint32(ch4); % no BF as too bright
        chFindBestPlanesSummary = squeeze(sum(chFindBestPlanes,[1,2])); % plane sum intensity vector
        %figure; plot(chFindBestPlanesSummary)
        GoodPlanes = chFindBestPlanesSummary > (max(chFindBestPlanesSummary)/3); % Above 33% of maximum
        %GoodPlanes = chFindBestPlanesSummary > (max(chFindBestPlanesSummary)/10); % Above 33% of maximum
        GoodPlanes = bwareafilt(GoodPlanes, 1);
        clear 'chFindBestPlanes'
        ch1 = ch1(:,:,GoodPlanes); % Hoechst % vol(ch1, 0, 10000)
        ch2 = ch2(:,:,GoodPlanes); % TH
        ch3 = ch3(:,:,GoodPlanes); % Map2
        ch4 = ch4(:,:,GoodPlanes); % Synuclein

        if ~HandleOverfitting
            Objects = f_imageAnalysis(num2str(organoid), ch1, ch2, ch3, ch4, PreviewPath, OrganoidsToAnalyze(organoid, :));
            ObjectsAll{organoid} = Objects;
        end
    catch ME
        disp(['Error for organoid ', num2str(organoid)])
        errorMessage = sprintf('Error in function %s() at line %d.\n\nError Message:\n%s', ...
            ME.stack(1).name, ME.stack(1).line, ME.message);
        fprintf(1, '%s\n', errorMessage);
        continue
    end
end

data = vertcat(ObjectsAll{:})
save([OutPath, filesep, 'data.mat'], 'data')
writetable(data, [OutPath, filesep, 'data.csv'])
%final_file = [OutPath, filesep, 'data.csv']
%writematrix(data, final_file)
