segmentation / ivote3

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


% ******* Initialize voting parameters **************************************
rmax = 10;		%maximum radius of the cell
ang_deg = 25.1;		%half the angular range of the voting area
ang = ang_deg * pi / 180;
iter = 6;	%number of voting iterations
t0 = 1.2; %threshold color
sigma = [4, 4, 2];
% t = 0.1;
d_ang= ang / (iter);

% ******** Testing parameters ******************************************
p = [100, 50, 100];
ps = [100, 100, 50];
% 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);
I_theta1 = acos(Igrad_x./Imag);
I_theta2 = acos(Igrad_y./Imag);
I_theta3 = acos(Igrad_z./Imag);

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

%Set the boundaries of the threshold image to zero
It(1:rmax, :, :) = 0;
It(ps(1) - rmax:ps(1), :,:) = 0;
It(:, 1:rmax, :) = 0;
It(:, ps(2) - rmax:ps(2),:) = 0;
It(:, :, 1:rmax) = 0;
It(:,:, ps(3) - rmax: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
range = -rmax:rmax;                                 %create an array of values between -rmax and rmax
[mx, my, mz] = meshgrid(range, range, range);       %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 =  m_mag <= rmax;                             %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, (rmax^3)));


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

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

	%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(vx,vy,vz);
		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 = Ivote;
	
	elseif itr ==2
		Ivote2 = Ivote;
	
	elseif itr ==3
		Ivote3 = Ivote;
	
	elseif itr ==4
		Ivote4 = Ivote;
		
	elseif itr == 5
		Ivote5 = Ivote;
	end
	t_v1 = toc;   
	disp(['voting done.  time =',num2str(t_v1)]);

	% update the voting direction
	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
 

%%
t = 250;
out = Ivote;
out(out<t) = 0;

fid_out = fopen('D:\ivote3_files\out.vol', 'w');
fwrite(fid_out, out);
fclose(fid_out);