function  [Objects] = f_imageAnalysis(Label, chHoechst, chMap2, chTuj1, chTH, PreviewPath, MetaTable)
    
     if false % debug tool
        chHoechst = chHoechst(700:2700,3000:5000,:);
        chTH = chTH(700:2700,3000:5000,:);
        chTuj1 = chTuj1(700:2700,3000:5000,:);
        chMap2 = chMap2(700:2700,3000:5000,:);
    end
    
    %% Organoid segmentation
%     chCombo = imlincomb(0.25, chHoechst, 0.25, chMap2, 0.25, chTuj1, 0.25, chTH);%vol(chCombo, 0, 500)
%     OrganoidMask = chCombo > 115;% vol(OrganoidMask)
%     OrganoidMask = bwareaopen(OrganoidMask, 1E6);
%     OrganoidMask = imclose(OrganoidMask, strel('sphere', 3));

    chOrganoid = imlincomb(0.25, chHoechst, 0.25, chMap2, 0.25, chTuj1, 0.25, chTH);
    %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 > 1200; % vol(DeadCoreMaskShrinked, 0,1)
    DeadCoreMaskShrinked = bwareaopen(DeadCoreMaskShrinked, 500);
    DeadCoreMaskShrinked = imclose(DeadCoreMaskShrinked, strel('disk', 7));
    DeadCoreMask = logical(imresize3(uint8(DeadCoreMaskShrinked), size(chHoechst))); % vol(DeadCoreMask, 0, 1)
    % 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 > 35; %vol(NucleiMask)
    NucleiMask = bwareaopen(NucleiMask, 100);%vol(NucleiMask)
    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;
    

    %% Segment Map2
   % vol(chMap2, 0, 3000) 
%     %% Map2 & Tuj1
%     % detect large structures
    for p=1:size(chMap2, 3)
        Map2_FT(:,:,p) = f_LPF_by_FFT(chMap2(:,:,p), 'Butterworth', [15,1], 0);
    end % vol(chMap2_FT * 1000, 0, 100, hot)

    Map2Mask = Map2_FT > 0.009; %0.008; % vol(Map2Mask)
    Map2Mask = bwareaopen(Map2Mask, 500); %27
    
    Map2MaskWODeadCore = Map2Mask & ~DeadCoreMask;
    
    clear('Map2_FT') 

    %Map2DoG = imfilter(chMap2, fspecial('gaussian', 90, 1) - fspecial('gaussian', 90, 11), 'symmetric'));
    %Map2DoG = gather(imfilter(gpuArray(chMap2), fspecial('gaussian', 90, 1) - fspecial('gaussian', 90, 11), 'symmetric'));
%     Map2DoGMask = Map2DoG > 0; % vol(Map2DoGMask, 0,1)    
%     Map2GlobalMask = chMap2 > 1600; % vol(chMap2,0, 10000) % vol(Map2GlobalMask,0,1)
%     Map2Mask = Map2DoGMask & Map2GlobalMask;
%     Map2Mask = bwareaopen(Map2Mask, 1000);
   

    %% TH
    % vol(chTH, 0, 3000)
    for p=1:size(chTH, 3)
        TH_FT(:,:,p) = f_LPF_by_FFT(chTH(:,:,p), 'Butterworth', [15,1], 0);
    end % vol(TH_FT * 1000, 0, 100, hot)

    %THMask = TH_FT > 0.005; % vol(THMask)
    THMask = TH_FT > 0.0025; % vol(THMask)
    
    %THMask = bwareaopen(THMask, 27);
    THMask = bwareaopen(THMask, 200);
    THMaskWODeadCore = THMask & ~DeadCoreMask;
    THandMap2Mask =  THMask & Map2Mask; %vol(THandMap2)
    THandMap2MaskWODeadCore = THandMap2Mask & ~DeadCoreMask;
    clear('TH_FT') 
    
    %% Tuj1
%     %Tuj1DoG = gather(imfilter(gpuArray(chTuj1), fspecial('gaussian', 90, 1) - fspecial('gaussian', 90, 11), 'symmetric'));
%     %Tuj1DoG = imfilter(chTuj1, fspecial('gaussian', 90, 1) - fspecial('gaussian', 90, 11), 'symmetric'));
%     Tuj1DoGMask = Tuj1DoG > 0; % vol(Tuj1DoGMask, 0,1)    
%     Tuj1GlobalMask = chTuj1 > 1600; % vol(chTuj1,0, 10000) % vol(Tuj1GlobalMask,0,1)
%     Tuj1Mask = Tuj1DoGMask & Tuj1GlobalMask;
%     Tuj1Mask = bwareaopen(Tuj1Mask, 1000);

% vol(chTuj1, 0, 3000)
   for p=1:size(chTuj1, 3)
        Tuj1_FT(:,:,p) = f_LPF_by_FFT(chTuj1(:,:,p), 'Butterworth', [15,1], 0);
    end % vol(chTuj1_FT * 1000, 0, 100, hot)

    Tuj1Mask = Tuj1_FT > 0.008; % vol(Tuj1Mask) 0.008;
    Tuj1Mask = bwareaopen(Tuj1Mask, 50); %27
    
    Tuj1MaskWODeadCore = Tuj1Mask & ~DeadCoreMask;
    THandTuj1Mask =  THMask & Tuj1Mask; %vol(THandTuj1Mask)
    THandTuj1MaskWODeadCore = THandTuj1Mask & ~DeadCoreMask;
    
    clear('Tuj1_FT') 

    

    %% TH skeleton3D EPFL
    disp('Start skel')
    tic
    skelTH = Skeleton3D(THMask);
    toc
    disp('Skel done')
    %vol(skelTH, 0, 1)
    [AdjacencyMatrixTH, nodeTH, linkTH] = Skel2Graph3D(skelTH,0);                       
    %imtool(AdjacencyMatrixTH, [])
    NodeTH = zeros(size(THMask), 'uint8');
    NodeIdxs = vertcat(nodeTH(:).idx);
    NodeTH(NodeIdxs) = 1;
    %vol(NodeTH)    
    if size(NodeIdxs, 1) == 0
        return
    end
    NodeTHPreview = uint8(skelTH) + NodeTH + uint8(THMask); 
    NodeTHPreview2D = max(NodeTHPreview, [], 3);
    %it(NodeTHPreview2D)
    %vol(NodeTHPreview, 0, 3, 'jet')    
    NodeDegreeVectorTH = sum(AdjacencyMatrixTH, 1);

    ZeroNodeExplanationNeeded = 0;
    if ZeroNodeExplanationNeeded
        ZeroNodes = find(NodeDegreeVectorTH == 0);
        ZeroNodesLinIdx = vertcat(nodeTH(ZeroNodes).idx);
        ZeroNodeMask = zeros(size(THMaskClipped), 'uint8');
        ZeroNodeMask(ZeroNodesLinIdx) = 1; %vol(ZeroNodeMask)
        NodePreviewZeroCase = uint8(skelTH) + NodeMaskTH + 10*uint8(ZeroNodeMask) + uint8(THMask);
    end  

    %% TH Fragmentation

    % Define structuring element for surface detection
    Conn6 = strel('sphere', 1); % 6 connectivity
    % Detect surface
    SurfaceTH = THMask & ~(imerode(THMask, Conn6));
    %vol(SurfaceTH)
    
    
    %% Feture extraction
    %%ObjectsThisOrganoid = table();
    
    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.NucleiDead = sum(NucMaskHigh(:));
    Objects.NucleiNorm = Objects.Nuclei / Objects.OrganoidVol;
    
    Objects.Map2 = sum(Map2Mask(:));
    Objects.Map2ByNuclei = Objects.Map2 / Objects.Nuclei;
    Objects.Map2ByNucleiAlive = Objects.Map2 / sum(NucMaskAlive(:));
    Objects.Map2Norm = Objects.Map2 / Objects.OrganoidVol;
    Objects.Map2WODeadCore = sum(Map2MaskWODeadCore(:));
    %Objects.Map2MaskWODeadCoreNorm = Objects.Map2MaskWODeadCore / Objects.OrganoidWOCoreVol;

    Objects.TH = sum(THMask(:));
    Objects.THByNuclei = Objects.TH / Objects.Nuclei;
    Objects.THByNucleiAlive = Objects.TH / sum(NucMaskAlive(:));
    Objects.THNorm = Objects.TH / Objects.OrganoidVol;
    Objects.THMaskWODeadCore = sum(THMaskWODeadCore(:));
    %Objects.THMaskWODeadCoreNorm = Objects.THMaskWODeadCore / Objects.OrganoidWOCoreVol;
    
    Objects.Tuj1 = sum(Tuj1Mask(:));
    Objects.Tuj1ByNuclei = Objects.Tuj1 / Objects.Nuclei;
    Objects.Tuj1ByNucleiAlive = Objects.Tuj1 / sum(NucMaskAlive(:));
    Objects.Tuj1Norm = Objects.Tuj1 / Objects.OrganoidVol;
    Objects.Tuj1WODeadCore = sum(Tuj1MaskWODeadCore(:));
    %Objects.Tuj1MaskWODeadCoreNorm = Objects.Tuj1MaskWODeadCore / Objects.OrganoidWOCoreVol;
    
    %Objects.THMask = sum(THMask(:));
    %Objects.TotalNucMask = sum(NucleiMask(:));
    %Objects.Tuj1Area = sum(Tuj1Mask(:));
    
    Objects.THoverlapMap2 = sum(THandMap2Mask(:));
    Objects.THoverlapMap2Norm = Objects.THoverlapMap2 / Objects.OrganoidVol;
    Objects.THandMap2WODeadCore = sum(THandMap2MaskWODeadCore(:));
    %Objects.THandMap2MaskWODeadCoreNorm = Objects.THoverlapMap2 / Objects.OrganoidWOCoreVol;
    Objects.THoverlapMap2ByNuclei = sum(THandMap2Mask(:)) / Objects.Nuclei;
    Objects.THoverlapMap2ByNucleiAlive = sum(THandMap2Mask(:)) / sum(NucMaskAlive(:));
    Objects.THoverlapMap2ByMap2 = sum(THandMap2Mask(:)) / Objects.Map2;
    
    Objects.THoverlapTuj1 = sum(THandTuj1Mask(:));
    Objects.THoverlapTuj1Norm = Objects.THoverlapTuj1 / Objects.OrganoidVol;
    Objects.THandTuj1WODeadCore = sum(THandTuj1MaskWODeadCore(:));
    Objects.THoverlapTuj1ByNuclei = sum(THandTuj1Mask(:)) / Objects.Nuclei;
    Objects.THoverlapTuj1ByNucleiAlive = sum(THandTuj1Mask(:)) / sum(NucMaskAlive(:));
    Objects.THoverlapTuj1ByTuj1 = sum(THandTuj1Mask(:)) / Objects.Tuj1;
      
    
    Objects.HoechstMeanInNuclei = mean(chHoechst(NucleiMask));
    Objects.Map2MeanInMap2 = mean(chMap2(Map2Mask));
    Objects.THMeanInTH = mean(chTH(THMask));
    Objects.Tuj1MeanInTH = mean(chTuj1(Tuj1Mask));
    
%     Objects.MeanOrganoidNuclei = mean(chHoechst(OrganoidMask));
%     Objects.MeanOrganoidMap2 = mean(chMap2(OrganoidMask));
%     Objects.MeanOrganoidTuj1 = mean(chTuj1(OrganoidMask));
%     Objects.MeanOrganoidTH = mean(chTH(OrganoidMask));
    

    Objects.SkelTH = sum(skelTH(:));
    Objects.SkelTHNorm = Objects.SkelTH / Objects.OrganoidVol; 
    Objects.Nodes = size(nodeTH, 2);
    Objects.Links = size(linkTH, 2);
    Objects.NodesNorm = Objects.Nodes / Objects.SkelTH;
    Objects.LinksNorm = Objects.Links / Objects.SkelTH;
    Objects.THFragmentation = sum(SurfaceTH(:)) / sum(THMask(:));
    
    
    
    %% Previews 
    % Scalebar
    imSize = [size(chHoechst, 1), size(chHoechst, 2)];
    [BarMask, BarCenter] = f_barMask(200, 0.32393102760889064, imSize, imSize(1)-200, 200, 25);
    %it(BarMask)
    %chEmpty = zeros(size(chHoechst),'uint16');

    PreviewHoechst = imoverlay2(imadjust(max(chHoechst,[],3),[0 0.07]), bwperim(max(NucleiMask,[],3)), [0 0 1]);
    PreviewHoechst = imoverlay2(PreviewHoechst, BarMask, [1 1 1]); %imtool(PreviewHoechst)
    
    PreviewHoechstNoMask = imoverlay2(imadjust(max(chHoechst,[],3),[0.002 0.07]), [0 0 1]); %, bwperim(max(NucleiMask,[],3)), [0 0 1]);
    PreviewHoechstNoMask = imoverlay2(PreviewHoechstNoMask, BarMask, [1 1 1]); %imtool(PreviewHoechstNoMask)
    
    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)
    
    %PreviewMap2 = imoverlay2(imadjust(max(chMap2,[],3),[0 0.2]), bwperim(max(Map2Mask,[],3)), [1 0 1]);
    PreviewMap2 = imoverlay2(imadjust(max(chMap2,[],3),[0 0.15]), bwperim(max(Map2Mask,[],3)), [1 0 1]);
    %PreviewMap2 = f_imoverlayAlpha(imadjust(max(chMap2,[],3),[0 0.15]), max(Map2Mask,[],3), [1 0 1], 0.02);
    PreviewMap2 = imoverlay2(PreviewMap2, BarMask, [1 1 1]); %imtool(PreviewMap2)
    
    PreviewMap2NoMask = imoverlay2(imadjust(max(chMap2,[],3),[0.002 0.15]), [1 0 1]);
    PreviewMap2NoMask = imoverlay2(PreviewMap2NoMask, BarMask, [1 1 1]); %imtool(PreviewMap2NoMask)
    
    %PreviewTuj1 = imoverlay2(imadjust(max(chTuj1,[],3),[0 0.02]), bwperim(max(Tuj1Mask,[],3)), [1 0 0]);
    PreviewTuj1 = imoverlay2(imadjust(max(chTuj1,[],3),[0 0.1]), bwperim(max(Tuj1Mask,[],3)), [1 0 0]);
    PreviewTuj1 = imoverlay2(PreviewTuj1, BarMask, [1 1 1]); %imtool(PreviewTuj1)
    
    PreviewTuj1NoMask = imoverlay2(imadjust(max(chTuj1,[],3),[0.002 0.09]), [1 0 0]);
    PreviewTuj1NoMask = imoverlay2(PreviewTuj1NoMask, BarMask, [1 1 1]); %imtool(PreviewTuj1NoMask)
    
    %PreviewTH = imoverlay2(imadjust(max(chTH,[],3),[0 0.02]), bwperim(max(THMask,[],3)), [0 1 0]);
    PreviewTH = imoverlay2(imadjust(max(chTH,[],3),[0 0.05]), bwperim(max(THMask,[],3)), [0 1 0]);
    PreviewTH = imoverlay2(PreviewTH, BarMask, [1 1 1]); %imtool(PreviewTH)
    
    PreviewTHNoMask = imoverlay2(imadjust(max(chTH,[],3),[0.002 0.02]), [0 1 0]);
    PreviewTHNoMask = imoverlay2(PreviewTHNoMask, BarMask, [1 1 1]); %imtool(PreviewTHNoMask)
    
    RGB = cat(3, imadjust(max(chMap2, [], 3), [0 0.2], [0 1]), imadjust(max(chTH, [], 3), [0 0.02], [0 1]), imadjust(max(chHoechst, [], 3), [0 0.07], [0 1]));
    RGB = imoverlay2(RGB, BarMask, [1 1 1]);
    %imtool(RGB)

    RGBpaper = cat(3, imadjust(max(chTuj1, [], 3), [0.002 0.09], [0 1]), imadjust(max(chTH, [], 3), [0.002 0.02], [0 1]), imadjust(max(chHoechst, [], 3), [0.002 0.07], [0 1]));
    RGBpaper = imoverlay2(RGBpaper, BarMask, [1 1 1]);
    %imtool(RGB)
    %% Write Previews 
    IdentityString = ['Organoid', '_', Objects.OrganoidID{:}{:}];
    imwrite(PreviewHoechst, [PreviewPath, filesep, IdentityString, '_Hoechst.png'])
    imwrite(PreviewHoechstNoMask, [PreviewPath, filesep, IdentityString, '_HoechstNoMask.png'])
    imwrite(PreviewMap2, [PreviewPath, filesep, IdentityString, '_Map2.png'])
    imwrite(PreviewMap2NoMask, [PreviewPath, filesep, IdentityString, '_Map2NoMask.png'])
    imwrite(PreviewTuj1, [PreviewPath, filesep, IdentityString, '_Tuj1.png'])
    imwrite(PreviewTuj1NoMask, [PreviewPath, filesep, IdentityString, '_Tuj1NoMask.png'])
    imwrite(PreviewTH, [PreviewPath, filesep, IdentityString, '_TH.png'])
    imwrite(PreviewTHNoMask, [PreviewPath, filesep, IdentityString, '_THNoMask.png'])
    imwrite(RGB, [PreviewPath, filesep, IdentityString, 'RGB.png'])
    imwrite(RGBpaper, [PreviewPath, filesep, IdentityString, 'RGBpaper.png'])
    imwrite(PreviewNucDead, [PreviewPath, filesep, IdentityString, 'PreviewNucDead.png'])
    
end

    
