segmentation / ivote3

clc;
clear;
disp('***************** NEW RUN *********************');
total = tic;


% ******* Initialize voting parameters **************************************
rmax = [9 9 5];		%maximum radius of the cell
ang_deg = 20.1;		%half the angular range of the voting area
ang = ang_deg * pi / 180;
iter = 5;	%number of voting iterations
t0 = 1.0; %threshold color
sigma = [3, 3, 1.5];
% t = 0.1;

d_ang= ang / (iter);

% ******** Testing parameters ******************************************
p = [50, 50, 150];
ps = [400, 400, 200];
% ps = [100, 50, 40];
% I = syn_Img(rmax , ps);
volfile = 'nissl-rat.vol';
fid = fopen(volfile); % open the file that include the image
S = fread(fid, 3, 'int32');
X = S(1);
Y = S(2);
Z = S(3);
% load the VOL data into a 2D matrix
I = fread(fid,[X Y*Z], 'uint8');
fclose(fid);
%change this to a 3D matrix
I = (reshape(I, [X, Y, Z]));
% invert the intensity
I = (255 - I);

%perform a gaussian blur
Iblur = gauss_blur3d(I, sigma);

%crop out a small subregion of I and Iblur
Iblur = Iblur(p(1):p(1)+ps(1)-1, p(2):p(2)+ps(2)-1, p(3):p(3)+ps(3)-1);
I = I(p(1):p(1)+ps(1)-1, p(2):p(2)+ps(2)-1, p(3):p(3)+ps(3)-1);
%
% compute the gradient
[Igrad_y, Igrad_x, Igrad_z] = gradient(Iblur);

%calculate the gradient magnitude
Imag = sqrt(Igrad_x .^ 2 + Igrad_y .^ 2 + Igrad_z .^2);
Isize = size(I);
I = single(I);
Iblur = single(Iblur);

%
%set a threshold for the gradient magnitude
It = Imag > t0;

%Set the boundaries of the threshold image to zero
It(1:rmax(1), :, :) = 0;
It(ps(1) - rmax(1):ps(1), :,:) = 0;
It(:, 1:rmax(2), :) = 0;
It(:, ps(2) - rmax(2):ps(2),:) = 0;
It(:, :, 1:rmax(3)) = 0;
It(:,:, ps(3) - rmax(3):ps(3)) = 0;
%
%get the indices of all of the nonzero values in the threshold image
%   (voter positions)
[Itx,Ity,Itz] = ind2sub(size(It),find(It));
Vi =(find(It));
nV = nnz(It);
%
% create a meshgrid describing coordinates relative to the voter position
rangex = -rmax(1):rmax(1);                                 %create an array of values between -rmax and rmax
rangey = -rmax(2):rmax(2);    
rangez = -rmax(3):rmax(3);    
[mx, my, mz] = meshgrid(rangex, rangey, rangez);       %create a template describing local pixel position in a small cube
m_mag = (sqrt(mx.^2 + my.^2 + mz.^2));                %create a template describing the distance from the center of a small cube

% create a mask for the voting area
M_dist =  (mx.^2/rmax(1)^2 + my.^2/rmax(2)^2 + mz.^2/rmax(3)^2)   <= 1;                             %mask for the voting area distance (all values < rmax from the center)
 
% calculate the direction vector between a pixel and voter
LV_x = mx./m_mag;
LV_y = my./m_mag;
LV_z = mz./m_mag;

%number of pixels in the voting area of each voter (initialize to zero)
validPixels = (zeros(nV,1)); 
%%
%indices of pixels in the voting area of each voter
%   indices reference the 3D image
g_v_prime = zeros(nV, ceil(rmax(1)*rmax(2)*rmax(3)*ang));


%% vote
tic;
%for each iteration (in iterative voting)
for itr = 1 : iter+1

	%initialize the vote image to zero
	Ivote = zeros(Isize);

	%for each voter (nonzero pixels in the threshold image It)
	for v = 1: nV

		%get the cartesian coordinates of the voter v in the main image I
		vx = Itx(v);
		vy = Ity(v);
		vz = Itz(v);
		vi = Vi(v);
		
		%retreive the gradient magnitude at the voter position
		vmag = Imag(vi);
		
		%retrieve the gradient
		gx = Igrad_x(vi);
		gy = Igrad_y(vi);
		gz = Igrad_z(vi);
		
		%calculate the gradient magnitude
		dmag = sqrt (gx^2 + gy^2 + gz^2);
		
		%calculate the normalized gradient direction
		dx = gx / dmag;
		dy = gy / dmag;
		dz = gz / dmag;
		
		%calculate the angle between the voter direction and the pixel direction
		cos_diff = LV_x .* dx + LV_y .* dy + LV_z .* dz;
		ang_diff = acos(cos_diff);

		%create an angular mask for the voting area
		M_angle = cos_diff >= cos(ang);

		%combine the two masks to mask out the voting angle
		M = M_angle .* M_dist;		

		% get the coordinates of each pixel in the final voter mask M
		pi = find(M);
		
		%calculate the number of pixels in the voting region
		npts =  nnz(M);
		validPixels(v) = npts; 

		%convert every index in the voting area from a local 3D index to a global 3D index (into the original image I)
		global_px = vx + mx(pi);
		global_py = vy + my(pi);
		global_pz = vz + mz(pi);
		
		%convert the global 3D index of each point into a global 1D index
		global_pi = sub2ind(ps, global_px, global_py, global_pz);

		g_v_prime (v, 1:npts) = global_pi;


		Ivote( global_pi ) = Ivote( global_pi ) + vmag;

	end

	if itr ==1
		Ivote1 = single(Ivote);
	
	elseif itr ==2
		Ivote2 = single(Ivote);
	
	elseif itr ==3
		Ivote3 = single(Ivote);
	
	elseif itr ==4
		Ivote4 = single(Ivote);
		
	elseif itr == 5
		Ivote5 = single(Ivote);
	end
	t_v1 = toc;   
	disp(['voting done.  time =',num2str(t_v1)]);

	% update the voting direction
	if ang>0
		tic;
		for v = 1: nV
			% coordinates of the current voter
			vx = Itx(v);
			vy = Ity(v);
			vz = Itz(v);

			%get the local value of the voting image
			local_Ivote = Ivote(g_v_prime(v,1:validPixels(v)));

			%find the index of the maximum value
			[~, local_max_idx] = max(local_Ivote);

			%convert this into a global subscript
			[g_px, g_py, g_pz] = ind2sub(size(Ivote), g_v_prime(v,local_max_idx));

			%compute the vector from the voter position to this position
			Igrad_x(vx, vy, vz) = g_px - vx;
			Igrad_y(vx, vy, vz) = g_py - vy;
			Igrad_z(vx, vy, vz) = g_pz - vz;

		end
	

	tdir1 = toc;
	display (['updating dir done.  time = ', num2str(tdir1)]);
	ang = ang - d_ang;
	end
 end
 

%%
t = 350;
conn = [5 5 3];
Icenter = local_max(Ivote, conn, t);
% center = Ivote1;
% center(center<t) = 0;
% center = imregionalmax(center);
% cn = nnz(center);
% [cx, cy, cz] = ind2sub(size(center), find(center));
% Icenter = zeros(size(center));
% for cc =1:cn
% 	Icenter(cx(cc), cy(cc), cz(cc)) = 255;
% end

% fid_Ic = fopen('image_center2-300.vol', 'w');
% fwrite(fid_Ic, Icenter);
% fclose(fid_Ic);
cn = nnz(Icenter);
[cx, cy, cz] = ind2sub(size(Icenter), find(Icenter));
Ic2d = zeros(size(Icenter,1), size(Icenter,2));
for cc =1:cn
	Ic2d(cx(cc), cy(cc)) = 1;
end
I2d = max(I, [], 3);
% figure(1),imagesc(I2d); colormap(gray);
% figure(2),imagesc(Ic2d); colormap(gray);
%
out1(:,:,1) = mat2gray(I2d);
out1(:,:,2) = mat2gray(Ic2d);
out1(:,:,3) = mat2gray(I2d);
figure(1), imagesc((out1));
%%% % imwrite(mat2gray(c2d), 'vote.bmp');
%%
% figure(1); imagesc(squeeze(I(:,:,ceil(size(I,3)/2)))), colormap(gray);
%  figure(33); imagesc(squeeze(Ivote3(:,:,ceil(size(Ivote,3)/2)))), colormap(gray);