function [] = f_imageAnalysis(ch1, ch2, WellThis, FieldThis, MesFile, PreviewPath)
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here

    % vol(ch1,100,10000)

    %% segment nuclei
    NucMask = zeros(size(ch1), 'logical');
    NucDoG = imfilter(ch1, fspecial('gaussian', 101, 3) - fspecial('gaussian', 101, 21), 'symmetric');
    % vol(NucDoG, 0, 200)
    parfor t = 1:size(ch1, 3)
        t
        %NucMask(:,:,t) = medfilt2(ch1(:,:,t)) > 430; % imtool(NucMask(:,:,t)) % imtool(ch1(:,:,t),[])% imtool(medfilt2(ch1(:,:,t)),[])
        NucMask(:,:,t) = medfilt2(NucDoG(:,:,t)) > 20; % imtool(NucMask(:,:,t)) % imtool(ch1(:,:,t),[])% imtool(medfilt2(ch1(:,:,t)),[])
        NucMask(:,:,t) = bwareaopen(NucMask(:,:,t), 50);
        %% Split Nuclei
        D = bwdist(~NucMask(:,:,t));
        %imtool(D,[]);
        WI = imcomplement(D);
        W = watershed(imhmin(WI, 0.5));%it(W)
        NucStencil = uint16(NucMask(:,:,t)) .* uint16(W);
        NucMask(:,:,t) = NucStencil > 0;
        NucMask(:,:,t) = bwareaopen(NucMask(:,:,t), 50);
        % imtool(f_imoverlayIris(imadjust(ch1(:,:,1), [0 0.05], [0 1]), bwperim(NucMask(:,:,1)), [0 0 1]))
        % imtool(f_imoverlayIris(imadjust(max(ch1,[],3), [0 0.05], [0 1]), bwperim(NucMask(:,:,1)), [0 0 1]))
        % vol(NucMask, 0, 1)
    end
    
    
%     parfor t = 1:size(ch1, 3)
%         t
%         %NucMask(:,:,t) = medfilt2(ch1(:,:,t)) > 430; % imtool(NucMask(:,:,t)) % imtool(ch1(:,:,t),[])% imtool(medfilt2(ch1(:,:,t)),[])
%         NucMask(:,:,t) = medfilt2(ch1(:,:,t)) > 420; % imtool(NucMask(:,:,t)) % imtool(ch1(:,:,t),[])% imtool(medfilt2(ch1(:,:,t)),[])
%         NucMask(:,:,t) = bwareaopen(NucMask(:,:,t), 50);
%         %% Split Nuclei
%         D = bwdist(~NucMask(:,:,t));
%         %imtool(D,[]);
%         WI = imcomplement(D);
%         W = watershed(imhmin(WI, 0.5));%it(W)
%         NucStencil = uint16(NucMask(:,:,t)) .* uint16(W);
%         NucMask(:,:,t) = NucStencil > 0;
%         NucMask(:,:,t) = bwareaopen(NucMask(:,:,t), 50);
%     end
    %vol(NucMask, 0,1)
   
    Tracks = zeros(size(ch1), 'uint16'); % imtool(ch1(:,:,1))
    GoodTracks = [];
    [LM, countNucT1] = bwlabeln(bwareafilt(imclearborder(NucMask(:,:,1)), [100, 300])); % imtool(LM,[])
    for nucT1 = 1:countNucT1
        disp(nucT1)
        [r,c] = find(LM == nucT1);
        idx = sub2ind(size(Tracks), r, c, repmat(1, size(r)));
        Tracks(idx) = nucT1;
        ThisNucCandidateLargestReconstructed = LM == nucT1; % imtool(ThisNucCandidateLargestReconstructed)
        for t = 2:300
            NucTLastThis = ThisNucCandidateLargestReconstructed;
            OverlapInTime = NucTLastThis & NucMask(:,:,t);
            ThisNucCandidates = OverlapInTime;
            ThisNucCandidateLargest = bwareafilt(ThisNucCandidates, 1); % imtool(ThisNucCandidateLargest)
            ThisNucCandidateLargestReconstructed = imreconstruct(ThisNucCandidateLargest, NucMask(:,:,t));%imtool(OverlapInTime)
            if sum(OverlapInTime(:)) < 0.8 * sum(ThisNucCandidateLargestReconstructed(:))
                break
            end
            if sum(OverlapInTime(:)) > 1.2 * sum(ThisNucCandidateLargestReconstructed(:))
                break
            end

            [r,c] = find(ThisNucCandidateLargestReconstructed);
            idx = sub2ind(size(Tracks), r, c, repmat(t, size(r)));
            Tracks(idx) = nucT1;
            if t == 300
                GoodTracks = [GoodTracks, nucT1];
            end
        end
    end
    
    %vol(Tracks)
    %vi(Tracks)
    
    TracksSelect = uint16(ismember(Tracks, GoodTracks)) .* Tracks;
    %vol(TracksSelect)
    %vol(TracksSelect(:,:,end))
    
    
    for track = 1:numel(GoodTracks)
        track
        trackDilatedThis = imdilate(TracksSelect == GoodTracks(track), strel('disk', 3)); % vol(trackDilatedThis)
        PathThis = [PreviewPath(1:end-9), filesep, WellThis, '_', num2str(track), '.csv']
        DataThis = {};
        for t = 1:300
            FrameFluo4 = ch2(:,:,t);
            DataThis{t} = mean(FrameFluo4(trackDilatedThis(:,:,t)));
        end
        DataThis = cell2table(DataThis', 'VariableNames', {'Fluo4MeanIntensity'});
        writetable(DataThis, PathThis)
    end


    %% Previews

    imSize = size(NucMask);
    [BarMask, BarCenter] = f_barMask(20, 2 * 0.32393102760889064, imSize, imSize(1)-50, 50, 20); % 20x binning 2
    %it(BarMask)

    NucPreview = f_imoverlayIris(imadjust(ch1(:,:,1), [0 0.03], [0 1]), bwperim(TracksSelect(:,:,1)), [0 0 1]);
    NucPreview = f_imoverlayIris(NucPreview, BarMask, [1 1 1]);
    %imtool(NucPreview)
    NucPreviewPath = [PreviewPath, filesep, WellThis, '_', FieldThis, '_Nuc.png'];
    imwrite(NucPreview, NucPreviewPath)
    
    TrackPreview = f_imoverlayIris(imadjust(max(ch1,[],3), [0 0.03], [0 1]), bwperim(max(TracksSelect,[],3)), [0 0 1]);
    TrackPreview = f_imoverlayIris(TrackPreview, BarMask, [1 1 1]);
    %imtool(TrackPreview,[])
    

end

