clc;
clear;
disp('***************** NEW RUN *********************');
% ******* Initialize voting parameters **************************************
rmin = 2;       %minimum radius of a cell
rmax = 30;		%maximum radius of the cell
phi = 25;		%half the angular range of the voting area
iterNum =6;	%number of voting iterations

%%%%%%%%%%
% cell_t = 4;
% sigmagk =0.8;
% sigma =0.75;
% final_t = 250;
cell_t = 1;
sigmagk = 0.8;
sigma = 2.5;
       % ******************Loading input Image**********************************
% I = imread('test/test10.png');
I = imread('newcrp.bmp');

X = size(I, 1);
Y = size(I, 2);
% **************** Gradient Magnitude and Direction for Image **************
% make a 2D gaussian kernel has its size defined by sigma (the std)
gX = sigma * 4;
gY = gX;

[gx, gy] = meshgrid(-gX:gX,-gY:gY);
a = 1/(2*pi*sigma^2);
b = -1/(2*sigma^2);
c = gx.^2 + gy.^2;
G = a*exp(b*c); % the Gaussian function

[dgx,dgy] = gradient(G); % derivatives of gaussian kernel in x and y directions.

%convolve I with the derivatives of the gaussian
Igx = conv2(double(I),dgx,'same');
Igy = conv2(double(I),dgy,'same');

M = sqrt(Igx.*Igx+Igy.*Igy); % compute the gradient magnitude.
D = atan2(Igx,Igy);

% ************************************** Thresholding **************************************
Mt = M > cell_t;

%find nonzero values in the threshold image (return 2D indices)
[sx, sy] = find(Mt);

%find the 1D index into the gradient direction D
idx = sub2ind(size(D), sx, sy);
utheta = D(idx);    %get the gradient direction for each coordinate

% compute the 2D Gaussian Kernel with size of [7 7] and yariance of 1.
Xgk = sigmagk * 4;
Ygk = Xgk;

[gkx, gky] = meshgrid(-Xgk:Xgk,-Ygk:Ygk);
agk = 1/(2*pi*sigmagk^2);
bgk = -1/(2*sigmagk^2);
cgk = gkx.^2 + gky.^2;
gaussiankernel = agk*exp(bgk*cgk); % the Gaussian function

% Voting
for i=1 : iterNum
    if phi>=0
        [VVV, utheta, phi] = voting (rmin, rmax, X, Y, phi, utheta, sx, sy, length(sx),M, gaussiankernel);
    end
end

%% ************************************** Detecting Cell Centers **************************************
final_t = 620;% 400-425 the best result with this data set
cellcenter = VVV;
cellcenter (VVV <final_t) = 0;
cellcenter = imregionalmax(cellcenter);

% ********************* Create an Overlay *************************************************************
output = cellcenter;

 

[tx, ty] = find(output);
% figure; imshow(mat2gray(output)); title('voting results');
% Ground Truth
gt = imread('GT.bmp');
gt(gt<255) = 0;
gt(gt==255)=1;
[gtx, gty] = find(gt);
% figure; imshow(mat2gray(gt)); title('Ground Truth');

% Validation
gt1 = [gtx(:), gty(:)];
out1 = [tx(:), ty(:)];
[IDX, D] = knnsearch(out1,gt1);
idx = IDX;
idx(D>5) = 0;
tp = nnz(idx);
lgt = length(gt1);
lout = length(out1);
fp = lout - tp;
fn = lgt - tp;
precision(final_t) = tp/(tp+fp);
recall(final_t) = tp/(tp+fn);
accuracy= tp/(tp+fp+fn)