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

BW_vec = reshape(BW,[],1);
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 cell 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));%mean(bg_vec);
mx = pkfnd(imgfilter,3.5*mean(sub_img_vec),objS2/2);
imgfilter = imgfilter2;
mxAll = mx;
mx = [];
for imx = 1:size(mxAll, 1)
    if sum(sum(BW(mxAll(imx ,2)-2:mxAll(imx ,2)+2, mxAll(imx ,1)-2:mxAll(imx ,1)+2))) > 0
        if (imgfilter2(mxAll(imx ,2), mxAll(imx ,1)) > 5*mean(sub_img_vec2(sub_img_vec2>0))) && ...
                (img(mxAll(imx ,2), mxAll(imx ,1)) > 1.3*mean(sub_img_raw_vec(sub_img_raw_vec>0)))
            mx = [mx; mxAll(imx ,:)];
        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
        indices = find(xout(:,1)< 0.06 & xout(:,2) < 8);
        
        xout = xout(indices, :);
        out = out(indices, :);
        out(:, 4) = ones(size(out, 1),1)*bg;
    end
else out = [];
    xout = [];
end
XOUT = [XOUT; xout];
OUT = [OUT; out];

xout = XOUT;
out = OUT;