function  [Objects] = f_imageAnalysis(Label, chHoechst, chMAP2, chFOXA2, chTH, PreviewPath, MetaTable)
    
    %vol(chHoechst, 0, 1000)
    %vol(chMAP2, 0, 2000)
    %vol(chTH, 0, 2000)
    %vol(chFOXA2, 0, 1000)
    
    %% Organoid segmentation
    chOrganoid = imlincomb(0.25, chHoechst, 0.25, chMAP2, 0.25, chFOXA2, 0.25, chTH);
    %vol(chOrganoid, 0, 1000)
    OrganoidMask = medfilt3(chOrganoid, [5,5,3]) > 300; % vol(OrganoidMask,0,1)
    OrganoidMask = Iris_KeepbiggestObct3D(OrganoidMask);
    OrganoidMask = imclose(OrganoidMask, strel('cuboid', [11,11,5]));
    OrganoidMask = imfill(OrganoidMask, 'holes');
    
    
    %% Segment Nuclei
    %ch1DoG = imfilter(double(ch, chMAP2, chGlut4,Hoechst), fspecial('gaussian', 55, 1) - fspecial('gaussian', 55, 11), 'symmetric');%vol(ch1DoG, 0, 100, 'hot')
    ch1DoG = imfilter(double(chHoechst), fspecial('gaussian', 55, 1) - fspecial('gaussian', 55, 11), 'symmetric');%vol(ch1DoG, 0, 100, 'hot')
    NucleiMask = ch1DoG > 120; %vol(NucleiMask)
    NucleiMask = bwareaopen(NucleiMask, 100);%vol(NucleiMask)


     %% Segment MAP2
    %MAP2Mask = chMAP2 > 1200; %vol(MAP2Mask, 0,1)
    %MAP2Mask = bwareaopen(MAP2Mask, 250);
    
    for p=1:size(chMAP2, 3)
        Map2_FT(:,:,p) = f_LPF_by_FFT(chMAP2(:,:,p), 'Butterworth', [15,1], 0);
    end % vol(Map2_FT * 1000, 0, 100, hot)

    MAP2Mask = Map2_FT > 0.015; %vol(MAP2Mask)
    MAP2Mask = bwareaopen(MAP2Mask, 27);
    
    clear('Map2_FT') 
    
  %% FOXA2
    FOXA2DoG = imfilter(chFOXA2, fspecial('gaussian', 101, 1.5) - fspecial('gaussian', 101, 33), 'symmetric');%vol(FOXA2DoG, 0, 500, 'hot')
    FOXA2Mask = FOXA2DoG > 200; %vol(NucleiMask)
    FOXA2Mask= FOXA2Mask & NucleiMask;
    FOXA2Mask = bwareaopen(FOXA2Mask, 100);%vol(FOXA2Mask)
    
%takes away the grid from the slide
    SparseMosaicEmptyRegion = chFOXA2 == 0;
    SparseMosaicEmptyRegionDilated = imdilate(SparseMosaicEmptyRegion, strel('disk', 33));
    FOXA2Mask= FOXA2Mask & OrganoidMask & ~SparseMosaicEmptyRegionDilated;
    %% TH
    %vol(chTH, 0, 1000, hot)
    for p=1:size(chTH, 3)
        TH_FT(:,:,p) = f_LPF_by_FFT(chTH(:,:,p), 'Butterworth', [15,1], 0);
    end % vol(chTH_FT * 1000, 0, 100, hot)

    THMask = TH_FT > 0.008; % vol(THMask)
    %THMask = TH_FT > 0.002; % vol(THMask)
    THMask = bwareaopen(THMask, 27);
    THandMAP2 =  THMask & MAP2Mask; %vol(THandMAP2)
    
    clear('TH_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)
    
    
    %% Feature extraction
    Objects = table();
    Objects.OrganoidIdx = {MetaTable.Idx};
    Objects.Well = {MetaTable.Well};
    Objects.AreaName = {MetaTable.AreaName};
    Objects.OrganoidID = {MetaTable.OrganoidID};
    Objects.TotalNucMask = sum(NucleiMask(:));
    Objects.MAP2Mask = sum(MAP2Mask(:));
    Objects.THMask = sum(THMask(:));
    Objects.THMaskByNuc = sum(THMask(:))/ sum(NucleiMask(:));
    Objects.THMaskByMAP2 = sum(THMask(:))/ sum(MAP2Mask(:));
    Objects.MAP2MaskByNuc = sum(MAP2Mask(:))/ sum(NucleiMask(:));
    Objects.THoverlapMAP2 = sum(THandMAP2(:));
    Objects.FOXA2Area = sum(FOXA2Mask(:));
    Objects.FOXA2MaskByNuc = sum(FOXA2Mask(:))/ sum(NucleiMask(:));
    %Objects.MeanOrganoidHoechst = mean(chHoechst(OrganoidMask));
    %Objects.MeanOrganoidMAP2 = mean(chMAP2(OrganoidMask));
    %Objects.MeanOrganoidTH= mean(chTH(OrganoidMask));
    %Objects.MeanOrganoidFOXA2 = mean(chFOXA2(OrganoidMask));
    
    Objects.SkelTH = sum(skelTH(:));
    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');

    RGB = cat(3, imadjust(max(chMAP2, [], 3), [0 0.08], [0 1]), imadjust(max(chTH, [], 3), [0 0.03], [0 1]), imadjust(max(chFOXA2, [], 3), [0 0.07], [0 1]));
    RGB = imoverlay2(RGB, BarMask, [1 1 1]);
    %imtool(RGB)
    
    PreviewNuc = imoverlay2(imadjust(max(chHoechst,[],3),[0 0.06]), BarMask, [1 1 1]);
    %imtool(PreviewNuc,[])
    
    PreviewNucMask = imoverlay2(imadjust(max(chHoechst,[],3), [0 0.06], [0 1]), bwperim(max(NucleiMask,[],3)), [0 0 1]);
    PreviewNucMask = imoverlay2(PreviewNucMask, BarMask, [1 1 1]);
    %imtool(PreviewNucMask,[])
    
    PreviewTH = imoverlay2(imadjust(max(chTH,[],3),[0 0.03]), BarMask, [1 1 1]);
    %imtool(PreviewTH,[])
    
    PreviewTHMask = imoverlay2(imadjust(max(chTH,[],3), [0 0.025], [0 1]), bwperim(max(THMask,[],3)), [0 1 0]);
    PreviewTHMask = imoverlay2(PreviewTHMask, BarMask, [1 1 1]);
    %imtool(PreviewTHMask,[])
    
    PreviewFOXA2 = imoverlay2(imadjust(max(chFOXA2,[],3),[0 0.04]), BarMask, [1 1 1]);
    %imtool(PreviewFOXA2,[])
    
    PreviewFOXA2Mask = imoverlay2(imadjust(max(chFOXA2,[],3), [0 0.04], [0 1]), bwperim(max(FOXA2Mask,[],3)), [1 0 0]);
    PreviewFOXA2Mask = imoverlay2(PreviewFOXA2Mask, BarMask, [1 1 1]);
    %imtool(PreviewFOXA2Mask,[])
    
    PreviewMAP2 = imoverlay2(imadjust(max(chMAP2,[],3),[0 0.08]), BarMask, [1 1 1]);
    %imtool(PreviewMAP2,[])
    
    PreviewMAP2Mask = imoverlay2(imadjust(max(chMAP2,[],3), [0 0.08], [0 1]), bwperim(max(MAP2Mask,[],3)), [1 0 1]);
    PreviewMAP2Mask = imoverlay2(PreviewMAP2Mask, BarMask, [1 1 1]);
    %imtool(PreviewMAP2Mask,[])
   

    %% Write Previews 
    IdentityString = ['Organoid', '_', Objects.AreaName{:}{:}, '__',Objects.OrganoidID{:}{:}];
    imwrite(RGB, [PreviewPath, filesep, IdentityString, 'RGB.png'])
    imwrite(PreviewNuc, [PreviewPath, filesep, IdentityString, 'Nuc.png'])
    imwrite(PreviewNucMask, [PreviewPath, filesep, IdentityString, 'NucMask.png'])
    imwrite(PreviewTH, [PreviewPath, filesep, IdentityString, 'TH.png'])
    imwrite(PreviewTHMask, [PreviewPath, filesep, IdentityString, 'THMask.png'])
    imwrite(PreviewFOXA2, [PreviewPath, filesep, IdentityString, 'FOXA2.png'])
    imwrite(PreviewFOXA2Mask, [PreviewPath, filesep, IdentityString, 'FOXA2Mask.png'])
    imwrite(PreviewMAP2, [PreviewPath, filesep, IdentityString, 'MAP2.png'])
    imwrite(PreviewMAP2Mask, [PreviewPath, filesep, IdentityString, 'MAP2Mask.png'])

    
end

    
