clear all;
disp('***************** NEW RUN *********************');
X = 128;
Y = 128;
Z = 128;
D = 10;
t0=1;
r1=12;
r2=10;
t=2300;
itr=8;
vote=10;
std = [5 5];
gt_filename = 'gt2.vol';
% out_filename = sprintf('128-128-128/0-nissl-std%d.%d-t0%d-r%d.%d-t%d-out%d.%d.vol',std(1), std(2),t0,r1,r2,t,itr,vote);
out_filename = sprintf('D:/build/ivote3-bld/shared2D-v8/centers.%d.vol',t);
% txt_filename = sprintf('128-128-128/0-validation-nissl-std%d.%d-r%d.%d-t%d-out%d.%d-D%d.txt',std(1), std(2),r1,r2,t,itr,vote,D);
txt_filename = sprintf('0-t%d--1.txt',t);
spec = sprintf('Nissl-std%d.%d-r%d.%d-t%d-out%d.%d',std(1), std(2),r1,r2,t,itr,vote);
fid0 = fopen(gt_filename);
gt = fread(fid0,[X Y*Z], 'single');
fclose(fid0);
% store x,y, and z directions of cell centers defined by the ground-trurth in an array
gt = reshape(gt, [X Y Z]);
[gtx, gty, gtz] = ind2sub(size(gt),find(gt));
ref = [gtx gty gtz];

fid1 = fopen(out_filename);
out = fread(fid1,[X Y*Z], 'single');
fclose(fid1);
% store x,y, and z directions of detected cells by the code in an array
out = reshape(out, [X Y Z]);
[cx, cy, cz] = ind2sub(size(out),find(out));
cent = [cx cy cz];
% 1-find the nearest detected cell for each element in the ground truth.
[idx, dist] = knnsearch(cent,ref);
% 2- find the same detected cells for different elements in the ground truth and assign that to which gt's element with the lowest distance.
ind_eq_idx = zeros(numel(idx), 10);
%set number of cells with the same detected cell center to zero.
noc=0;
for i=1:numel(idx)	
	% check the current element hasn't considered yet
	s_ind = sum(ind_eq_idx==i);
	if s_ind==0
		% check how many of the elemnts assigned to the same detected cell
		s = sum(idx==idx(i));		
		if s>1 
			noc = noc+1;
			%save the index and number of gt's elements with the same assigned detected cells.
			id_ref(noc,1) = i;
			id_ref(noc,2) = s;
			ind1 = find(idx==idx(i));
			ind_eq_idx(i, 1:numel(ind1)) = ind1;
		end
	end
end
% determine those indices which hs assigned to only one detected cell
b = idx;
u_eq_idx = unique(ind_eq_idx(ind_eq_idx>0));
b(u_eq_idx)=0;
% set the number of over sefgmented cells to zero.
oseg=0;
for k=1:size(id_ref,1)
	k1 = id_ref(k,1);
	k2 = id_ref(k,2);
	%find the minimum distance to the detected cell happened for which of the gt's element
	l = ind_eq_idx(k1,1:k2);
	[~, local_id_min] = min(dist(l));
	% set the element with minimum distance to the corresponding detected cell
	b(l(local_id_min)) = idx(k1);
	%remove the proper element from the list of indices with same designated detected cells
	u_eq_idx (u_eq_idx == l(local_id_min))=0;
	% check if the indices remained in the list has distance less than the value is set for validation
	ly = l (l~=l(local_id_min));
	distl = dist(ly);
	sl = sum(distl<D);
	% if the distance is low enought, consider the corresponding cell as oversegmented
	if sl>0
		oseg = oseg+1;
	end
end
%****
% check if the oversegmented detected cells could be assigned to the gt's element with no detected cell!
%find the indices of gt's element with no detected cell
u_idx1 = find(b==0);
%find the coordinates of gt's
ref1 = ref(u_idx1,:);
% find the indices of detected cells which are assigned to a gt's element
ub = unique(b(b>0));
% set the number of detected cell which are not assigned to a gt's element to zero
mm=0;
% find the indices of detected cell which are not assigned to a gt's element
for j=1:size(cent,1)
	z = sum(ub(:)==j);
	if z==0
		mm = mm+1;
		cent_id1(mm) = j;
	end
end
% find the coordinated of of detected cells
cent1 = cent(cent_id1,:);

% check if there are osme low enough distances
[idx1, dist1] = knnsearch(cent1,ref1);
ind_eq_idx_rc = zeros(numel(idx1), 10);
%set number of cells with the same detected cell center to zero.
noc_rc=0;
for i_rc=1:numel(idx1)	
	% check the current element hasn't considered yet
	s_ind_rc = sum(ind_eq_idx_rc==i_rc);
	if s_ind_rc==0
		% check how many of the elemnts assigned to the same detected cell
		s_rc = sum(idx1==idx1(i_rc));		
		if s_rc>1 
			noc_rc = noc_rc+1;
			%save the index and number of gt's elements with the same assigned detected cells.
			id_ref_rc(noc_rc,1) = i_rc;
			id_ref_rc(noc_rc,2) = s_rc;
			ind1_rc = find(idx1==idx1(i_rc));
			ind_eq_idx_rc(i_rc, 1:numel(ind1_rc)) = ind1_rc;
		end
	end
end
% determine those indices which hs assigned to only one detected cell
b_rc = idx1;
u_eq_idx_rc = unique(ind_eq_idx_rc(ind_eq_idx_rc>0));
b_rc(u_eq_idx_rc)=0;
% set the number of over sefgmented cells to zero.
oseg_rc=0;
for k_rc=1:size(id_ref_rc,1)
	k1_rc = id_ref_rc(k_rc,1);
	k2_rc = id_ref_rc(k_rc,2);
	%find the minimum distance to the detected cell happened for which of the gt's element
	l_rc = ind_eq_idx_rc(k1_rc,1:k2_rc);
	[~, local_id_min_rc] = min(dist1(l_rc));
	% set the element with minimum distance to the corresponding detected cell
	b_rc(l_rc(local_id_min_rc)) = idx1(k1_rc);
	%remove the proper element from the list of indices with same designated detected cells
	u_eq_idx_rc (u_eq_idx_rc == l_rc(local_id_min_rc))=0;
	% check if the indices remained in the list has distance less than the value is set for validation
	ly_rc = l_rc (l_rc~=l_rc(local_id_min_rc));
	distl_rc = dist1(ly_rc);
	sl_rc = sum(distl_rc<D);
	% if the distance is low enought, consider the corresponding cell as oversegmented
	if sl_rc>0
		oseg_rc = oseg_rc+1;
	end
end

%******
% b include the gt's element and its detected cells, for those element with no cell detection, b has zero value
b_ind = find(b==0);
b_dist = dist;
% remove the disatances for those elements with no detected cells from distance array
b_dist(b_ind)=-1;

b_ind_rc = find(b_rc==0);
b_dist_rc = dist1;
b_dist_rc(b_ind_rc)=-1;
% calculate Ttrue Positive, number of detected cells that have low enough distance to one of the gt's elements.
TP = sum(b_dist>=0) - sum(b_dist>D) + (sum(b_dist_rc>=0) - sum(b_dist_rc>D));
% calculate False Negative, number of gt's elements with no detected cells.
FN = size(ref, 1) - TP;
% calculate False Positive, number of detected cells, which their distance to the gt's element is long
% or considered as over segmentation.
FP = size(cent,1) - TP;
% compute the validation factors.
Precision = TP/(TP + FP);
Recall = TP/(TP+FN);
Accuracy= TP/(TP+FP+FN);

fid_text = fopen(txt_filename, 'w');
fprintf(fid_text, '********--------\tValidation Results for %s\t--------********\n\n',spec);
fprintf(fid_text, 'number of detected cells by GROUND_TRUTH = %d\n', size(ref,1));
fprintf(fid_text, 'number of detected cells by IVOTE3 = %d\n', size(cent,1));
fprintf(fid_text, 'number of OVER-SEGMENTED cells = %d\n', oseg);
fprintf(fid_text, 'TP = %f,\t\tFP = %f,\t\tFN = %f\n', TP,FP,FN);
fprintf(fid_text, 'Precision = %f\n', Precision);
fprintf(fid_text, 'Accuracy = %f\n', Accuracy);
fprintf(fid_text, 'Recall = %f\n', Recall);
fclose(fid_text);