%% bandpas   s: put frame around picture and remove it afterwards to avoid cuts by bandpassfilter
objS2 = 15; 
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;

imgfilter = bpass(I,1.6,objS2, -Inf);
imgfilter2 = bpass(I,1.0,objS2, -Inf);

filter_th = 4;
% peak intensity of the bandpass image must be > filter_th*median_bandpass_intensity. 2.5
raw_th = 1.6; 
% peak intensity of the raw image must be > raw_th*average_raw_intensity. 2

%%
imgfilter = imgfilter(objS2+1:size(img,1)+objS2, objS2+1:size(img,2)+objS2);
imgfilter2 = imgfilter2(objS2+1:size(img,1)+objS2, objS2+1:size(img,2)+objS2);
imgfilter_vec = reshape(imgfilter,[],1);

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

XOUT = [];
OUT = [];
mask = zeros(size(BW));
mask(BW > 0) = 1;
sub_img = mask.*imgfilter;
sub_img2 = mask.*imgfilter2;
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 cells and positive (~ abs) to calculate mean and std
sub_img_vec2 = reshape(sub_img2,[],1);
sub_img_vec2 = sub_img_vec2(sub_img_vec2 > 0);
sub_img_raw_vec = reshape(sub_img_raw,[],1);

signal = imgfilter(sub_img > 0);
bg = median(sub_img_raw_vec(sub_img_raw_vec>0));
filter_noise_1 = median(sub_img_vec);
mx = pkfnd(imgfilter,filter_th*filter_noise_1,objS2/2);
imgfilter = imgfilter2;

filter_noise_2 = median(sub_img_vec2(sub_img_vec2>0));
raw_noise = median(sub_img_raw_vec(sub_img_raw_vec>0));

mxAll = mx; 
mx = []; mx_raw_intensity = []; mx_filter_intensity = []; mx_center_intensity = [];

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
            
            raw_intensity = mean2(img(x_mx-1:x_mx+1, y_mx-1:y_mx+1));
            filter_intensity = imgfilter2(x_mx, y_mx);
            if (filter_intensity > filter_th*filter_noise_2) && (raw_intensity > raw_th*raw_noise)
                
                mx = [mx; mxAll(imx ,:)];
                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

if not(isempty(mx))
    opt.interactive = 0;
    opt.method = 6;
    
    [out, xout] = centfind(imgfilter2,mx,9,1,opt);
    if opt.method == 6
        out(:, 4) = ones(size(out, 1),1)*bg;
        out(:, 5) = ones(size(out, 1),1)*raw_noise;
        out(:, 6) = ones(size(out, 1),1)*filter_noise_1;
        out(:, 7) = ones(size(out, 1),1)*filter_noise_2;
        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;