matlab/main.m

clc;
clear;
disp('***************** NEW RUN *********************');
total = tic;
% ******* Initialize voting parameters **************************************
rmax = 7;		%maximum radius of the cell
ang_deg = 15.1;		%half the angular range of the voting area
ang = ang_deg * pi / 180;
iter = 2;	%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 = [200, 200, 100];
ps = [150, 100, 25];
r = 7;
I = syn_Img(r , ps);
% volfile = 'img\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)));
%--------------Display current data-------------------------
subplot(3, 3, 1),
imagesc(squeeze(I(:, :, ceil(size(I, 3)/2))));
subplot(3, 3, 2),
imagesc(squeeze(I(:, ceil(size(I, 2)/2), :)));
subplot(3, 3, 3),
imagesc(squeeze(I( ceil(size(I, 1)/2), :, :)));
%
subplot(3, 3, 4),
imagesc(squeeze(Iblur(:, :, ceil(size(Iblur, 3)/2))));
subplot(3, 3, 5),
imagesc(squeeze(Iblur(:, ceil(size(Iblur, 2)/2), :)));
subplot(3, 3, 6),
imagesc(squeeze(Iblur(ceil(size(Iblur, 1)/2), :, :)));
subplot(3, 3, 7),
imagesc(squeeze(Imag(:, :, ceil(size(Imag, 3)/2))));
subplot(3, 3, 8),
imagesc(squeeze(Imag(:, ceil(size(Imag, 2)/2), :)));
subplot(3, 3, 9),
imagesc(squeeze(Imag(ceil(size(Imag, 1)/2), :, :)));
colormap(gray);
%% 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;
		
		%gdir_x = Igrad_x(vx, vy, vz) / dir_mag;
		%gdir_y = Igrad_y(vx, vy, vz) / dir_mag;
		%gdir_z = Igrad_z(vx, vy, vz) / dir_mag;
		%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);
		%global_pi = (global_pz-1)*ps(1)*ps(2) + (global_py-1)*ps(1) + global_px;
		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)]);
	figure, imagesc(squeeze(Ivote(:, :, ceil(size(Ivote, 3)/2),:))); colormap(gray);
	% 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
		[gx, gy, gz] = ind2sub(size(Ivote), g_v_prime(v,local_max_idx));
		%compute the vector from the voter position to this position
		dx = gx - vx;
		dy = gy - vy;
		dz = gz - vz;
		Igrad_x(vx, vy, vz) = dx;
		Igrad_y(vx, vy, vz) = dy;
		Igrad_z(vx, vy, vz) = dz;
	end
	tdir1 = toc;
	display (['updating dir done.  time = ', num2str(tdir1)]);
	ang = ang - d_ang;
 end
 
%%
t = 10;
out = Ivote;
out(out<t) = 0;
% out(out>=t) = 255;
out = imregionalmax(out);
out1(:,:,:,1) = mat2gray(Iblur);
out1(:,:,:,2) = mat2gray(out);
out1(:,:,:,3) = mat2gray(Iblur);
figure(2);
subplot(2, 3, 1),
imagesc(squeeze(Iblur(:, :, ceil(size(Iblur, 3)/2))));
subplot(2, 3, 2),
imagesc(squeeze(Iblur(:, ceil(size(Iblur, 2)/2), :)));
subplot(2, 3, 3),
imagesc(squeeze(Iblur( ceil(size(Iblur, 1)/2), :, :)));
subplot(2, 3, 4),
imagesc( squeeze( out1(:, :, ceil(size(Ivote, 3)/2),:) ) );
subplot(2, 3, 5),
imagesc(squeeze(out1(:, ceil(size(Ivote, 2)/2), :, :)));
subplot(2, 3, 6),
imagesc(squeeze(out1( ceil(size(Ivote, 1)/2), :, :, :)));
colormap(gray);
%%
% figure, imagesc(squeeze(Iblur(:, :, ceil(size(Iblur, 3)/2)))); colormap(gray);
% figure, imagesc(squeeze(Ivote(:, :, ceil(size(Ivote, 3)/2),:))); colormap(gray);
%  figure, imagesc(squeeze(Iblur(:, ceil(size(Iblur, 2)/2), :))); colormap(gray);
%  figure, imagesc(squeeze(I(:, ceil(size(Ivote, 2)/2), :))); colormap(gray);
% % figure, imagesc(squeeze(Iblur( ceil(size(Iblur, 1)/2), :, :))); colormap(gray);
% figure, imagesc(squeeze(Ivote( ceil(size(Ivote, 1)/2), :, :))); colormap(gray);
% 
fid1 = fopen('ivote1.vol', 'w');
fwrite(fid1, Ivote1);
fclose(fid1);
% fid0 = fopen('iblur.vol', 'w');
% fwrite(fid0, Iblur);
% fclose(fid0);