clc; clear; disp('***************** NEW RUN *********************'); total = tic; % ******* Initialize voting parameters ************************************** rmax = 7; %maximum radius of the cell ang_deg = 25.1; %half the angular range of the voting area ang = ang_deg * pi / 180; iter = 5; %number of voting iterations t0 = 1.2; %threshold color sigma = [1, 1, 1]; % t = 0.1; d_ang= ang / (iter); % ******** Testing parameters ****************************************** % p = [100, 50, 100]; % ps = [200, 200, 100]; ps = [100, 50, 40]; I = syn_Img(rmax , 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; %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 = 2500; out = Ivote; out(out=t) = 255; out = imregionalmax(out); out1(:,:,:,1) = mat2gray(Iblur); out1(:,:,:,2) = mat2gray(out); out1(:,:,:,3) = mat2gray(Iblur); % out = out1; figure(4); 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); % % fid0 = fopen('D:\ivote3_files\iblur.vol', 'w'); % fwrite(fid0, Iblur); % fclose(fid0); % % fid1 = fopen('D:\ivote3_files\ivote1.vol', 'w'); % fwrite(fid1, Ivote1); % fclose(fid1); % % fid2 = fopen('D:\ivote3_files\ivote2.vol', 'w'); % fwrite(fid2, Ivote2); % fclose(fid2); % % fid3 = fopen('D:\ivote3_files\ivote3.vol', 'w'); % fwrite(fid3, Ivote3); % fclose(fid3); % % fid4 = fopen('D:\ivote3_files\ivote4.vol', 'w'); % fwrite(fid4, Ivote4); % fclose(fid4); % % fid5 = fopen('D:\ivote3_files\ivote5.vol', 'w'); % fwrite(fid5, Ivote5); % fclose(fid5); % % fid10 = fopen('D:\ivote3_files\ivote10.vol', 'w'); % fwrite(fid10, Ivote); % fclose(fid10);