segmentation / ivote3

clear;
disp('***************** NEW RUN *********************');
total = tic;
% ******* Initialize voting parameters **************************************
rmax = [10 10 10];		%maximum radius of the cell
rmin = [1 1 1];	
ang_deg = 25.1;		%half the angular range of the voting area
ang = ang_deg * pi / 180;
iter = 8 ;	%number of voting iterations
t0 = 1;
sigma = [5, 5, 5];
% t = 0.1;
d_ang= ang / (iter+2);
% ******** 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);
filename = '128-128-128/nissl-float-128.128.128.vol'; %'nissl-float-128.128.128.vol';
X = 128;
Y = 128;
Z = 128;
fidi = fopen(filename);
% load the VOL data into a 2D matrix
I = fread(fidi,[X Y*Z], 'single');
fclose(fidi);
%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);
% 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);

% h  = reshape(Imag, [X*Y*Z, 1]);
% hist(h, 100);

%set a threshold for the gradient magnitude
It = Imag > t0;
fidt = fopen('128-128-128/It.vol', 'w');
fwrite(fidt, It, 'single');
fclose(fidt);
%Set the boundaries of the threshold image to zero
It(1:rmax(1), :, :) = 0;
It(X - rmax(1):X, :,:) = 0;
It(:, 1:rmax(2), :) = 0;
It(:, Y - rmax(2):Y,:) = 0;
It(:, :, 1:rmax(3)) = 0;
It(:,:, Z - rmax(3):Z) = 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_dist1 = (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)
M_dist2 = (mx.^2/rmin(1)^2 + my.^2/rmin(2)^2 + mz.^2/rmin(3)^2)   >= 1 ;      
M_dist = M_dist1 .* M_dist2;
% 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;
mask = zeros(Isize);
mask1 = zeros(Isize);

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

	%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;
		
		%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(Isize, global_px, global_py, global_pz);

		g_v_prime (v, 1:npts) = global_pi;


		Ivote( global_pi ) = Ivote( global_pi ) + vmag;
	end
 	fid = fopen(sprintf('128-128-128/nissl-vote%d',itr), 'w');
	fwrite(fid, single(Ivote), '*single');
	fclose(fid);
 	
	t_v1 = toc;   
	disp(['voting done.  time =',num2str(t_v1)]);
	
	% update the voting direction
	if ang>=d_ang
		tic;
		Igrad_x = zeros(Isize);
		Igrad_y = zeros(Isize);
		Igrad_z = zeros(Isize);
		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
 
hv  = reshape(Ivote, [X*Y*Z, 1]);
hist(hv, 250);
%%
t = 300;
conn = [5 5 5];
Icenter = local_max(Ivote, conn, t);
fidc = fopen(sprintf('std3.2-r10.10-v8/out%d-t%d.vol',t,t0), 'w');
fwrite(fidc, single(Icenter), '*single');
fclose(fidc);
nnz(Icenter)
% [cxx1, cyy1, czz1] = ind2sub(size(Icenter),find(Icenter));

% % 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);