segmentation / ivote3

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


% ******* Initialize voting parameters **************************************
rmax = 9;		%maximum radius of the cell
ang_deg = 15.1;		%half the angular range of the voting area
ang = ang_deg * pi / 180;
iter = 6;	%number of voting iterations
t0 = 1.0; %threshold color
sigma = [2, 2, 1];
% t = 0.1;
iter = iter-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);

%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+1

	%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
	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 = 50;
% center = Ivote;
% % center(center<t) = 0;
% % center = imregionalmax(center);
% cn = nnz(center);
% [cx, cy, cz] = ind2sub(size(center), find(center));
% c2d = zeros(size(center,1), size(center,2));
% for cc =1:cn
% 	c2d(cx(cc), cy(cc)) = 1;
% end
% I2d = max(Iblur, [], 3);
% % figure(1),imagesc(I2d); colormap(gray);
% figure(2),imagesc(c2d); colormap(gray);
% 
% % out(:,:,1) = mat2gray(I2d);
% % out(:,:,2) = mat2gray(c2d);
% % out(:,:,3) = mat2gray(I2d);
% % figure (5), imagesc(out); 
% imwrite(mat2gray(c2d), 'vote.bmp');
figure(6); imagesc(squeeze(Ivote1(:,ceil(size(Ivote1,2)/2),:))); colormap(gray);