clear;close all;clc

n = './samplemovie/';
refimgMreB = [n 'AVG_Pos1.tif']; % location of the bright field image

%%
refimg = refimgMreB;
IMG = double(imread(refimg));
IMG_vec = reshape(IMG,[],1);

figure(100),hold on
subplot(1,2,2),imagesc(IMG);

objS = 10;
I = zeros(size(IMG)+2*objS);
I = I +median(IMG_vec);
I(objS+1:size(IMG,1)+objS, objS+1:size(IMG,2)+objS) = IMG;
I = bpass(I,.1,objS,-Inf);
I_vec = reshape(I,[],1);
clear L
L = imgaussfilt(I,2);
L = del2(L, 10);

L = abs(L);
L = mat2gray(L);
L = imgaussfilt(L,1);

level = graythresh(L);
BW = im2bw(L, level);
BW = BW(objS+1:size(IMG,1)+objS, objS+1:size(IMG,2)+objS);
BW = imdilate(BW,strel('disk',10));
BW = imerode(BW,strel('disk',10));
BW = imerode(BW,strel('disk',5));
BW = imdilate(BW,strel('disk',5));
BW = bwlabel(BW);

test = regionprops(BW, 'Perimeter', 'Area');
for i = 1: max(max(BW))
    ratio(i) =  test(i).Perimeter / test(i).Area;
    if ratio(i) > 0.35
        BW(BW == i) = 0;
    end
    if max(max(I(BW == i))) < 0.0*max(max(I))
        BW(BW == i) = 0;
    end
end
toosmall = find([test.Area]<250);
toosmall = ismember(BW,toosmall);
BW = BW.*~toosmall;
subplot(1,2,1),
imagesc(BW.*IMG)
BW = bwlabel(BW);
nbCells = max(max(BW));
cells = regionprops(BW, 'Orientation', 'MajorAxisLength', 'MinorAxisLength', ...
    'Eccentricity', 'Centroid', 'Area');

alpha= cat(1, cells.Orientation);
alpha = alpha.*(pi/180);
AxMin = cat(1, cells.MinorAxisLength);
AxMaj = cat(1, cells.MajorAxisLength);
Centr = cat(1, cells.Centroid);
hold on
scatter(Centr(:,1), Centr(:,2), 'xk')

for i = 1:nbCells
    Alpha = alpha(i);
    ry = 10;
    rx = 30;
    line([Centr(i,1), Centr(i,1) + cos(Alpha)*rx], ...
        [Centr(i,2), Centr(i,2) - sin(Alpha)*rx],...
        'Color', 'k', 'LineWidth', 1)
    line([Centr(i,1), Centr(i,1) + sin(Alpha)*ry], ...
        [Centr(i,2), Centr(i,2) + cos(Alpha)*ry],...
        'Color', 'k', 'LineWidth', 1)
    hold on
end
save('mask.mat', 'BW', 'cells', 'refimg', 'refimgMreB') 
% save the resulting mask as "BW"

%%