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