function  [Objects] = f_imageAnalysis(Label, chHoechst, chGFAP, chTuj1, chS100beta, PreviewPath, MetaTable)
    %vol(chHoechst, 0, 500)% Hoechst
    %vol(chTuj1, 0, 5000)
    %vol(chGFAP, 0, 1000)

    if false % debug tool
        chHoechst = chHoechst(700:2700,3000:5000,:);
        chGFAP = chGFAP(700:2700,3000:5000,:);
        chTuj1 = chTuj1(700:2700,3000:5000,:);
        chS100beta = chS100beta(700:2700,3000:5000,:);
    end
    
    %% Max projections
    %NucMax = max(chHoechst, [],3); imtool(NucMax,[])
    
    %% Segment Organoid
    chOrganoid = imlincomb(0.25, chHoechst, 0.25, chGFAP, 0.25, chTuj1, 0.25, chS100beta);
    %chOrganoid = chHoechst; % other markers are condition dependent
    %vol(chOrganoid, 0, 1000)
    OrganoidMask = chOrganoid > 120; % vol(OrganoidMask,0,1)
    OrganoidMask = bwareaopen(OrganoidMask, 1e5);
    OrganoidMask = Iris_KeepbiggestObct3D(OrganoidMask);
    OrganoidMask = imclose(OrganoidMask, strel('sphere', 5));
    OrganoidMask = imfill(OrganoidMask, 'holes');
    
    %% Segment Dead core
    % strategy resize (shrink) DoG (expand interploate) 
    ImHoechstShrinked = imresize3(chHoechst, [size(chHoechst,1)/10, size(chHoechst,2)/10, size(chHoechst,3)]);% vol(ImHoechstShrinked, 0, 5000) % imtool(max(ch1,[],3))
    DeadCoreDoG = imfilter(ImHoechstShrinked, fspecial('gaussian', 1001, 1) - fspecial('gaussian', 1001, 201), 0); % vol(DeadCoreDoG,0,10000)sparse mosaic has zeros anyway
    DeadCoreMaskShrinked = DeadCoreDoG > 2000; % vol(DeadCoreMaskShrinked, 0,1) it(max(DeadCoreDoG,[],3))
    DeadCoreMaskShrinked = bwareaopen(DeadCoreMaskShrinked, 500);
    DeadCoreMaskShrinked = imclose(DeadCoreMaskShrinked, strel('disk', 7));
    DeadCoreMask = logical(imresize3(uint8(DeadCoreMaskShrinked), size(chHoechst))); % vol(DeadCoreMask, 0, 1) it(max(DeadCoreMask,[],3))
    % vol(DeadCoreDoG, 0, 3000)

    OrganoidMaskWOCore = OrganoidMask & ~DeadCoreMask;

    
    %% Segment nuclei
    % vol(chHoechst, 0, 5000)
    NucDoG = imfilter(double(chHoechst), fspecial('gaussian', 55, 1) - fspecial('gaussian', 55, 11), 'symmetric');%vol(NucDoG, 0, 100, 'hot')
    NucleiMask = NucDoG > 30; %vol(NucleiMask)
    NucleiMask = bwareaopen(NucleiMask, 100);%vol(NucleiMask) it(max(NucleiMask,[],3))
    NucleiMask = NucleiMask | DeadCoreMask;
    
    ch1LP = imfilter(chHoechst, fspecial('gaussian', 11, 1), 'symmetric');
    NucMaskHigh =  (ch1LP > 2800) .* NucleiMask; %vol(NucMaskHigh, 0, 1)
    NucMaskAlive = NucleiMask & ~NucMaskHigh; % vol(NucMaskAlive)

    NucleiWODeadCore = NucleiMask & ~DeadCoreMask;

    %% GFAP
    %vol(chGFAP,0,5000)
    %GFAPDoG = gather(imfilter(gpuArray(chGFAP), fspecial('gaussian', 355, 1) - fspecial('gaussian', 355, 150), 'symmetric'));
    GFAPDoG = imfilter(chGFAP, fspecial('gaussian', 355, 1) - fspecial('gaussian', 355, 150), 'symmetric');
    %vol(GFAPDoG, 0, 2000, 'hot') it(max(GFAPDoG,[],3))
    GFAPMask = GFAPDoG > 1000; % vol(GFAPMask, 0,1)
    GFAPMask = GFAPMask & chGFAP > 1000;
    GFAPMask = bwareaopen(GFAPMask, 600);
    GFAPMask = GFAPMask & OrganoidMask;
    GFAPMaskWODeadCore = GFAPMask & ~DeadCoreMask;
    
    %% Tuj1 % vol(chTuj1,0,5000)
    %Tuj1DoG = gather(imfilter(gpuArray(chTuj1), fspecial('gaussian', 355, 1) - fspecial('gaussian', 355, 150), 'symmetric'));
    Tuj1DoG = imfilter(chTuj1, fspecial('gaussian', 355, 1) - fspecial('gaussian', 355, 150), 'symmetric');
    %vol(Tuj1DoG, 0, 2000, 'hot') it(max(Tuj1DoG,[],3))
    Tuj1Mask = Tuj1DoG > 200; % vol(Tuj1Mask, 0,1)
    Tuj1Mask = bwareaopen(Tuj1Mask, 1000);
    Tuj1Mask = Tuj1Mask & OrganoidMask;
    Tuj1MaskWODeadCore = Tuj1Mask & ~DeadCoreMask;
    
    %% S100beta
    %vol(chS100beta, 0, 2000)
    %S100betaDoG = gather(imfilter(gpuArray(chS100beta), fspecial('gaussian', 355, 1) - fspecial('gaussian', 355, 150), 'symmetric'));
    S100betaDoG = imfilter(chS100beta, fspecial('gaussian', 355, 1) - fspecial('gaussian', 355, 150), 'symmetric');
    %vol(S100betaDoG, 0, 500, 'hot') it(max(S100betaDoG,[],3))
    S100betaMask = S100betaDoG > 350; % vol(S100betaMask, 0,1)
    S100betaMask = bwareaopen(S100betaMask, 300);
    S100betaMask = S100betaMask & OrganoidMask;
    S100betaMaskWODeadCore = S100betaMask & ~DeadCoreMask;
    
    
    GFAPandS100betaMask = S100betaMask & GFAPMask;
    
    %S100betaBackground = median(chS100beta(OrganoidMask & ~S100betaMask));
    %chS100betaCorr = chS100beta - S100betaBackground;
    %vol(chS100betaCorr,0,2000)
    
%     %% reiterate to refine by alleviating organoid border contrast causing FP
%     S100betaDoG = imfilter(chS100betaCorr, fspecial('gaussian', 355, 1) - fspecial('gaussian', 355, 150), 'symmetric');
%     S100betaMask = S100betaDoG > 600; % vol(S100betaMask, 0,1)
%     S100betaMask = bwareaopen(S100betaMask, 500);
%     S100betaMask = S100betaMask & OrganoidMask;
%     S100betaMaskWODeadCore = S100betaMask & ~DeadCoreMask;
%     
%     
%     GFAPandS100betaMask = S100betaMask & GFAPMask;

    
    %% Feature extraction
    Objects = table();
    Objects.OrganoidIdx = MetaTable.Idx;
    Objects.Well = MetaTable.Well;
    Objects.AreaName = MetaTable.AreaName;
    Objects.OrganoidID = MetaTable.OrganoidID;
    
    Objects.OrganoidVol = sum(OrganoidMask(:));
    Objects.DeadCoreVol = sum(DeadCoreMask(:));
    Objects.OrganoidWOCoreVol = sum(OrganoidMaskWOCore(:));
    %%% Objects.DeadCoreVolNorm = Objects.DeadCoreVol / Objects.OrganoidVol;
    Objects.NucleiWODeadCoreVol = sum(NucleiWODeadCore(:));
    
    Objects.Nuclei = sum(NucleiMask(:));
    Objects.GFAP = sum(GFAPMask(:));
    Objects.S100beta = sum(S100betaMask(:));
    Objects.Tuj1 = sum(Tuj1Mask(:));
    
    Objects.NucleiDead = sum(NucMaskHigh(:));
	Objects.GFAPByNucleiAlive = sum(GFAPMask(:)) / sum(NucMaskAlive(:));
	Objects.S100betaByNucleiAlive = sum(S100betaMask(:)) / sum(NucMaskAlive(:));
    
    %Objects.GFAPMask = sum(GFAPMask(:));
    Objects.GFAPNorm = Objects.GFAP / Objects.OrganoidVol;
    Objects.GFAPWODeadCore = sum(GFAPMaskWODeadCore(:));
    Objects.GFAPByNuclei = Objects.GFAP / Objects.Nuclei;
    %%% Objects.GFAPMaskWODeadCoreNorm = Objects.GFAPMaskWODeadCore / Objects.OrganoidWOCoreVol;
    
    %Objects.S100betaMask = sum(S100betaMask(:));
    Objects.S100betaNorm = Objects.S100beta / Objects.OrganoidVol;
    Objects.S100betaWODeadCore = sum(S100betaMaskWODeadCore(:));
    Objects.S100betaByNuclei = Objects.S100beta / Objects.Nuclei;
    %%% Objects.S100betaMaskWODeadCoreNorm = Objects.S100betaMaskWODeadCore / Objects.OrganoidWOCoreVol;
    
    %Objects.Tuj1Mask = sum(Tuj1Mask(:));
    Objects.Tuj1Norm = Objects.Tuj1 / Objects.OrganoidVol;
    Objects.Tuj1WODeadCore = sum(Tuj1MaskWODeadCore(:));
    Objects.Tuj1ByNuclei = Objects.Tuj1 / Objects.Nuclei;
    %%% Objects.Tuj1MaskWODeadCoreNorm = Objects.Tuj1MaskWODeadCore / Objects.OrganoidWOCoreVol;
    
    Objects.HoechstMeanInNuclei = mean(chHoechst(NucleiMask));
    Objects.GFAPMeanInGFAP = mean(chGFAP(GFAPMask));
    Objects.S100betaMeanInS100beta = mean(chS100beta(S100betaMask));
    Objects.Tuj1MeanInTuj1 = mean(chTuj1(Tuj1Mask));
    
    Objects.GFAPoverlapS100beta = sum(GFAPandS100betaMask(:));
    %%% Objects.GFAPoverlapS100betaNormToOrganoid = Objects.GFAPoverlapS100beta / Objects.OrganoidVol;
    Objects.GFAPoverlapS100betaByS100beta = Objects.GFAPoverlapS100beta / Objects.S100beta;
    Objects.GFAPoverlapS100betaByGFAP = Objects.GFAPoverlapS100beta / Objects.GFAP;
    Objects.GFAPoverlapS100betaByNuclei = Objects.GFAPoverlapS100beta / Objects.Nuclei;
    Objects.GFAPoverlapS100betaByNucleiAlive = Objects.GFAPoverlapS100beta / sum(NucMaskAlive(:));
    
    %% Previews 
    % Scalebar
    imSize = [size(chHoechst, 1), size(chHoechst, 2)];
    [BarMask, BarCenter] = f_barMask(200, 0.32393102760889064, imSize, imSize(1)-200, 200, 25);
    %it(BarMask)

    RGB = cat(3, imadjust(max(chS100beta, [], 3), [0 0.08], [0 1]), imadjust(max(chGFAP, [], 3), [0 0.3], [0 1]), imadjust(max(chHoechst, [], 3), [0 0.03], [0 1]));
    RGB = imoverlay2(RGB, BarMask, [1 1 1]);
    %imtool(RGB)
   
    PreviewNuc = f_imoverlayAlpha(imadjust(max(chHoechst,[],3), [0 0.03], [0 1]), max(NucleiMask,[],3), [0 0 1], 0.1);
    PreviewNuc = imoverlay2(PreviewNuc, BarMask, [1 1 1]);
    %imtool(PreviewNuc,[])
    
    PreviewNucDead = imoverlay2(imadjust(max(chHoechst,[],3),[0 0.07]), bwperim(max(NucMaskHigh,[],3)), [0 0 1]);
    PreviewNucDead = imoverlay2(PreviewNucDead, BarMask, [1 1 1]); %imtool(PreviewNucDead)
    
    PreviewGFAP = f_imoverlayAlpha(imadjust(max(chGFAP,[],3), [0 0.3], [0 1]), max(GFAPMask,[],3), [0 1 0], 0.1);
    %PreviewGFAP = imoverlay2(imadjust(max(chGFAP,[],3), [0 0.5], [0 1]), max(GFAPMask,[],3), [0 1 0]);
    PreviewGFAP = imoverlay2(PreviewGFAP, BarMask, [1 1 1]);
    %imtool(PreviewGFAP,[])

    PreviewS100beta = f_imoverlayAlpha(imadjust(max(chS100beta,[],3), [0 0.05], [0 1]), max(S100betaMask,[],3), [1 0 0], 0.1);
    PreviewS100beta = imoverlay2(PreviewS100beta, BarMask, [1 1 1]);
    %imtool(PreviewS100beta,[])
    
    PreviewTuj1 = imoverlay2(imadjust(max(chTuj1,[],3), [0 0.05], [0 1]), bwperim(max(Tuj1Mask,[],3)), [0 1 0]);
    PreviewTuj1 = imoverlay2(PreviewTuj1, BarMask, [1 1 1]);
    %imtool(PreviewTuj1,[])

    %% Write Previews 
    IdentityString = ['PreviewOrganoid', '_', Objects.OrganoidID{:}]
    imwrite(RGB, [PreviewPath, filesep, IdentityString, 'RGB.png'])
    imwrite(PreviewNuc, [PreviewPath, filesep, IdentityString, 'Nuc.png'])
    imwrite(PreviewGFAP, [PreviewPath, filesep, IdentityString, '.GFAP.png'])
    imwrite(PreviewS100beta, [PreviewPath, filesep, IdentityString, 'S100beta.png'])
    imwrite(PreviewTuj1, [PreviewPath, filesep, IdentityString, 'Tuj1.png'])
    imwrite(PreviewNucDead, [PreviewPath, filesep, IdentityString, 'PreviewNucDead.png'])

end
