%% User inputs
%ChannelCount = 3;

try
    delete(gcp('nocreate'))
    myCluster = parcluster('local');
    AvailableWorkers = myCluster.NumWorkers;
    if AvailableWorkers >= 112
        pool = parpool(112)
    else
        pool = parpool(7)
    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 = '/work/projects/lcsb_hcs/Data/ClaudiaSaraiva/TH_Tuj1_Map2/CS_20230303_AxPD9-AxPD18_TH-488_MAP2-647_TUJ1-568_daneurons_20230303_110744/03-03-23_11-20-31_in';
     %'/work/projects/lcsb_hcs/Data/ClaudiaSaraiva/TH_Tuj1_Map2/CS_20210917_AxPD9_TH-488_MAP2-647_TUJ1-568_axpd9_20210917_182459_in';
%end '/mnt/IrisHCS/Data/ClaudiaSaraiva/TH_Tuj1_Map2/CS_20210917_AxPD9_TH-488_MAP2-647_TUJ1-568_axpd9_20210917_182459_in'

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/lscratch/users/csaraiva/TH_Tuj1_MAP2_AxPD18/60Days/CS_20230303_AxPD9-AxPD18_TH-488_MAP2-647_TUJ1-568_daneurons_20230303_110744/03-03-23_11-20-31_out';
%end

f_LogDependenciesLinux(mfilename, 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

[FieldAssignmentsAll, OrganoidStencilsAll, WellsAll] = Iris_GetOrganoidFields(InPath, MetaData, Coordinates, false);

ObjectsAll = {};
for organoid = 61:71 %41:50 %51:60  %1:size(OrganoidsToAnalyze,1)
%for organoid = 55 % A04 OR6

    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, 10000)
        %imtool(max(ch1, [], 3),[])
        
        % Map2 647
        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, 5000)
        %imtool(max(ch2, [], 3),[])
        
        % Tuj1 (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, 5000)
        %imtool(max(ch3, [], 3),[])
        
                
        % TH 488 (green)
        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, 5000)
        %imtool(max(ch4, [], 3),[])
        

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


        clear 'ch1Ims' 'ch2Ims' 'ch3Ims' 'ch4Ims' 'ch5Ims'
        chFindBestPlanes = uint32(ch1)+ uint32(ch2) + uint32(ch3); % 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 = bwareafilt(GoodPlanes, 1);
        clear 'chFindBestPlanes'
        ch1 = ch1(:,:,GoodPlanes); % it(max(ch1,[],3))
        ch2 = ch2(:,:,GoodPlanes); % it(max(ch2,[],3))
        ch3 = ch3(:,:,GoodPlanes); % it(max(ch3,[],3))
        ch4 = ch4(:,:,GoodPlanes); % it(max(ch4,[],3))
        
        combineIms = false
        if combineIms % semi manual helper
            ch1L = ch1(4500:6500, 500:2500, :);% it(max(ch1L,[],3)) % organoid 9
            ch2L = ch2(4500:6500, 500:2500, :);% it(max(ch2L,[],3))
            ch3L = ch3(4500:6500, 500:2500, :);% it(max(ch3L,[],3))
            ch4L = ch4(4500:6500, 500:2500, :);% it(max(ch4L,[],3))
            
            ch1R = ch1(6000:8000, 4000:6000, :);% it(max(ch1R,[],3)) % organoid 25
            ch2R = ch2(6000:8000, 4000:6000, :);% it(max(ch2R,[],3))
            ch3R = ch3(6000:8000, 4000:6000, :);% it(max(ch3R,[],3))
            ch4R = ch4(6000:8000, 4000:6000, :);% it(max(ch4R,[],3))
            
            ch1 = [ch1L, ch1R];
            ch2 = [ch2L, ch2R];
            ch3 = [ch3L, ch3R];
            ch4 = [ch4L, ch4R];
            
            %% write tiff
            
            Ims = {ch1, ch2, ch3, ch4};
            Labels = {'ch1', 'ch2', 'ch3', 'ch4'}
            
            for i = 1:4
                t = Tiff(['/', Labels{i}, '.tif'], 'w');
                %
                tagstruct.ImageLength = size(Ims{i},1); 
                tagstruct.ImageWidth = size(Ims{i},2);
                tagstruct.Photometric = Tiff.Photometric.MinIsBlack;
                tagstruct.BitsPerSample = 16;
                tagstruct.SamplesPerPixel = size(Ims{i},3);
                tagstruct.PlanarConfiguration = Tiff.PlanarConfiguration.Chunky; 
                tagstruct.Software = 'MATLAB'; 

                setTag(t,tagstruct)
                write(t, Ims{i});
                close(t)
            end

            %tReadObject = Tiff(['/mnt/irisgpfs/projects/lcsb_hcs/Data/ClaudiaSaraiva/GFAP_S100b_NeuN/CS_20210924_AxPD9_NeuN-488_GFAP-647_S100b-568_axpd9-ast_20210929_145617/09-29-21_15-26-44_in/', Labels{i}, '.tif'], 'r');
            %ImTest = read(tReadObject);
            
            save(['/mnt/irisgpfs/projects/lcsb_hcs/Data/ClaudiaSaraiva/TH_Tuj1_Map2/CS_20210917_AxPD9_TH-488_MAP2-647_TUJ1-568_axpd9_20210917_182459_in', Labels{i}, '.mat'], 'organoid', 'PreviewPath', 'OrganoidsToAnalyze')
            %/mnt/irisgpfs/projects/lcsb_hcs/Data/ClaudiaSaraiva/GFAP_S100b_NeuN/CS_20210924_AxPD9_NeuN-488_GFAP-647_S100b-568_axpd9-ast_20210929_145617/09-29-21_15-26-44_in
            
        end
        

        %ObjectsThisOrganoid = f_imageAnalysis(num2str(organoid), ch1, ch2, ch3, ch4, ch5, PreviewPath, OrganoidsToAnalyze(organoid, :));
        ObjectsThisOrganoid = f_imageAnalysis(num2str(organoid), ch1, ch2, ch3, ch4, PreviewPath, OrganoidsToAnalyze(organoid, :));
        ObjectsAll{organoid} = ObjectsThisOrganoid;
    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'])
