objS2 = 15;

w1 = fspecial('gaussian',[5 5],0.4);
img = imfilter(img,w1,'replicate');

I = zeros(size(img)+2*objS2);
I = I + median(img_vec);
I(objS2+1:size(img,1)+objS2, objS2+1:size(img,2)+objS2) = img;

w1 = -fspecial('log',[17 17],2);
w2 = -fspecial('log',[17 17],1.8);
w3 = -fspecial('log',[17 17],1.6);

imgfilter = imfilter(img,w1,'replicate');
imgfilter2 = imfilter(img,w2,'replicate');
imgfilter3 = imfilter(img,w3,'replicate');

imgfilter = max(max(imgfilter, imgfilter2), imgfilter3);

imgfilter(imgfilter<0) = 0;

pk_th = quantile(imgfilter, 0.98, 'All')
img30per = quantile(img, 0.30, 'All');

Mask = imgfilter(:,:)>pk_th;
Mask = imfilter(Mask/2,[0,1,0;1,1,1;0,1,0],'replicate')*2;
Mask = Mask >= 3;
Mask = bwareafilt(Mask, [4,40],4);

%%
imgfilter_vec = reshape(imgfilter,[],1);

% background
BW_vec = reshape(BW,[],1);
bg_vec = imgfilter(BW_vec == 0);

XOUT = [];
OUT = [];

mask = zeros(size(BW));
mask(BW > 0) = 1;
mask(BW == j) = 1;
area = bwarea(mask);
sub_img = mask.*imgfilter;
sub_img_raw = mask.*img;

sub_img_vec = reshape(sub_img,[],1);
sub_img_vec = sub_img_vec(sub_img_vec > 0);
%consider only pixels that are inside the cell to calculate mean and std
sub_img_raw_vec = sub_img_raw(sub_img_raw>0);
% Peak find
% out=pkfnd(im,th,sz)
signal = imgfilter(sub_img > 0);
bg = median(sub_img_raw_vec(sub_img_raw_vec>0));%mean(bg_vec);
median(sub_img_raw_vec(sub_img_raw_vec>0))
mx = pkfnd(imgfilter.*Mask,pk_th, objS2/2); % pk_th,objS2/2);
imgfilter = imgfilter2;

imgfilter = img; imgfilter2 = img;

raw_median = median(sub_img_raw_vec(sub_img_raw_vec>0));
filter_median = median(sub_img_vec);

raw_noise = std(sub_img_raw_vec);
filter_noise = std(sub_img_vec);

mxAll = mx;
mx = [];

filter_th = 4; raw_th = 3;
mx_raw_intensity = []; mx_filter_intensity = []; mx_center_intensity = [];
mx_local_BG = []; mx_local_BG_std = [];

bg_disk = fspecial('disk',4)>0;
bg_disk(3:7, 3:7) = bg_disk(3:7, 3:7) - (fspecial('disk',2)>0);

for imx = 1:size(mxAll, 1)
    % get the coordinates of the peak to check if it meets conditions
    % to be a good peak (high enough filtered and raw intensity)
    x_mx = mxAll(imx ,1); y_mx = mxAll(imx ,2);
    if BW(x_mx, y_mx) > 0
        
        sub_img = img(x_mx-4:x_mx+4, y_mx-4:y_mx+4);
        
        raw_intensity = mean([img(x_mx-1:x_mx+1, y_mx)', img(x_mx, y_mx-1:y_mx+1)]);
        mx_loc_BG = mean(sub_img(bg_disk == 1));
        mx_loc_BG_std = std(sub_img(bg_disk == 1));
        filter_intensity = imgfilter2(x_mx, y_mx);
        if raw_intensity - mx_loc_BG > max(raw_th*mx_loc_BG_std, raw_median/2)
            
            mx = [mx; mxAll(imx ,:)];
            mx_local_BG = [mx_local_BG, mx_loc_BG];
            mx_local_BG_std = [mx_local_BG_std, mx_loc_BG_std];
            mx_raw_intensity = [mx_raw_intensity; raw_intensity];
            mx_filter_intensity = [mx_filter_intensity, filter_intensity];
            mx_center_intensity = [mx_center_intensity, img(x_mx, y_mx)];
        end
    end
end

new_center_int = [];

if not(isempty(mx))
    opt.interactive = 0;
    opt.method = 6;
    
    [out, xout] = centfind(imgfilter2,mx,9,1,opt);
    
    for kk = 1:size(out,1)
        x = out(kk,1); y = out(kk, 2);
        new_center_int = [new_center_int, mean2(img(floor(x):ceil(x), floor(y):ceil(y)))];
    end
    if opt.method == 6
        out(:, 4) = ones(size(out, 1),1)*img30per;
        out(:, 5) = mx_local_BG;
        out(:, 6) = mx_local_BG_std;
        out(:, 7) = ones(size(out, 1),1)*filter_noise;
        out(:, 8) = new_center_int; % out(:, 8) = mx_raw_intensity;
        out(:, 9) = mx_filter_intensity;
        out(:, 10) = mx_center_intensity;
    end
else out = [];
    xout = [];
end
XOUT = [XOUT; xout];
OUT = [OUT; out];


xout = XOUT;
out = OUT;