clear;clc;close all;

parentdir = './';
Folder = [parentdir 'samplemovie' filesep]; 
name01 = [Folder 'Pos0' filesep]; % location of the tif images
load('mask.mat') % load the mask

IMG = BW;
nbCells = size(cells, 1);

BW = bwmorph(BW,'thicken', 0);
BW = bwmorph(BW,'close', 1);
BW = bwlabel(BW);
BW = BW(1:512, 1:512);

figure,hold on,imagesc(BW)

fileList = dir([name01,'*.tif']); % list of tif images to perform peak detection
%%
name = {fileList.name};

Int = []; BG = [];
xxs = []; yys = [];

PosList = [];
i = 0;
penalty = 0;
somme = zeros(size(IMG));

pZ = 0;
mi = []; ma = mi; mit = ma; med = ma;
while i < length(fileList)-1
    i = i+1;
    img = zeros(size(IMG));
    clear imgfilter beadPos mask a b res img_vec imgfilter_vec ...
        xout mx diff Peaks Peaksz out
    img = double(imread([name01,name{i}]));
    img_vec = reshape(img,[],1);
    somme = somme + img;
    
    peakFinding;
     
    if (not(isempty(xout)))
        xout = xout(((xout(:, 5) < 512) & (xout(:, 6) < 512) &...
            (xout(:, 5) > 0) & (xout(:, 6) >0)), :);
        out = out(((xout(:, 5) < 512) & (xout(:, 6) < 512) & ...
            (xout(:, 5) > 0) & (xout(:, 6) >0)), :);
        pZ = pZ + 1;
        penalty = 0;
        if isempty(PosList)
            l = 0;
        else
            l = length(PosList(:,1));
        end
        k = length(xout(:,5));
        PosList(l+1:l+k,1) = xout(:,5);
        PosList(l+1:l+k,2) = xout(:,6);
        PosList(l+1:l+k,3) = str2num(name{i}(5:13));
        gaussInt(l+1:l+k,1) = xout(:,7);
        BG(l+1:l+k, 1) = out(:, 4);
        c = 1;
        mi = []; ma = mi; mit = ma; med = ma;
        for m = 1:k
            mitte = round([xout(m,5) xout(m,6)]);
            if (mitte(1)-c > 0) && (mitte(2)-c > 0) && (mitte(2)+c <= ...
                    length(img(:,1))) && (mitte(1)+c <= length(img(1,:)))
                aus = img(mitte(2)-c:mitte(2)+c, mitte(1)-c:mitte(1)+c);
                aus_vec = reshape(aus,[],1);
                med = [med median(aus_vec)];
                mit = [mit mean2(aus)];
                ma = [ma max(max(aus))];
                mi = [mi min(min(aus))];
            else
                intens = img(mitte(2), mitte(1));
                med = [med intens];
                mit = [mit intens];
                ma = [ma intens];
                mi = [mi intens];
            end
        end
        Int(l+1:l+k,1) = mit';
        Int(l+1:l+k,2) = med';
        Int(l+1:l+k,3) = ma';
        Int(l+1:l+k,4) = mi';
        
    else penalty = penalty + 1;
    end
    disp(i)
end

mkdir([parentdir 'tracks'])
save([parentdir  'tracks' filesep 'sampledata.mat'], ...
    'PosList', 'refimg', 'Int', 'cells', 'BW', 'refimgMreB', 'BG')

% saves the list of detected peaks in "PosList"
% PosList: x | y | image ID

