clc; clear; disp('***************** NEW RUN *********************'); total = tic; % ******* Initialize voting parameters ************************************** rmax = [15 15 15]; %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); % Iblur = I; % %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); 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_dist = (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; update_rx = zeros(Isize); update_rx(:) = rmax(1); update_ry = zeros(Isize); update_ry(:) = rmax(2); update_rz = zeros(Isize); update_rz(:) = rmax(3); %% %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 :2 %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; %%%%%%% rx = update_rx(vi); ry = update_ry(vi); rz = update_rz(vi); mx_v = mx(rmax(1)-rx+1:rmax(1)-rx+2*rx+1, rmax(2)-ry+1:rmax(2)-ry+2*ry+1, rmax(3)-rz+1:rmax(3)-rz+2*rz+1); my_v = my(rmax(1)-rx+1:rmax(1)-rx+2*rx+1, rmax(2)-ry+1:rmax(2)-ry+2*ry+1, rmax(3)-rz+1:rmax(3)-rz+2*rz+1); mz_v = mz(rmax(1)-rx+1:rmax(1)-rx+2*rx+1, rmax(2)-ry+1:rmax(2)-ry+2*ry+1, rmax(3)-rz+1:rmax(3)-rz+2*rz+1); M_dist = (mx_v.^2/rx^2 + my_v.^2/ry^2 + mz_v.^2/rz^2) <= 1 ; LV_x_v = LV_x(rmax(1)-rx+1:rmax(1)-rx+2*rx+1, rmax(2)-ry+1:rmax(2)-ry+2*ry+1, rmax(3)-rz+1:rmax(3)-rz+2*rz+1); LV_y_v = LV_y(rmax(1)-rx+1:rmax(1)-rx+2*rx+1, rmax(2)-ry+1:rmax(2)-ry+2*ry+1, rmax(3)-rz+1:rmax(3)-rz+2*rz+1); LV_z_v = LV_z(rmax(1)-rx+1:rmax(1)-rx+2*rx+1, rmax(2)-ry+1:rmax(2)-ry+2*ry+1, rmax(3)-rz+1:rmax(3)-rz+2*rz+1); %%%%%%% %calculate the angle between the voter direction and the pixel direction cos_diff = (LV_x_v .* dx + LV_y_v .* dy + LV_z_v .* 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_v(pi); global_py = vy + my_v(pi); global_pz = vz + mz_v(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; % if itr ==3 % if mod(v, 5000)==0 % mask(global_pi)= mask(global_pi) + 0.1; % mask (vi) = mask(vi) + 0.5; % end % end % if itr ==6 % if mod(v, 5000)==0 % mask1(global_pi)= mask1(global_pi) + 0.1; % mask1 (vi) = mask1(vi) + 0.5; % end % end % if itr==1 % for ix = -12:12 % for iy = -12:12 % for iz = -12:12 % mask(vx+ix, vy+iy, vz+iz) = M(ix+13,iy+13,iz+13)+ mask(vx+ix, vy+iy, vz+iz); % end % end % end % end % if itr==2 % for ix = -12:12 % for iy = -12:12 % for iz = -12:12 % mask1(vx+ix, vy+iy, vz+iz) = M(ix+13,iy+13,iz+13)+ mask1(vx+ix, vy+iy, vz+iz); % end % end % end % end end fid = fopen(sprintf('std5.5-r10.10-v8-update-rmax/vote%d',itr), 'w'); fwrite(fid, single(Ivote), '*single'); % if itr ==1 % fid = fopen(sprintf('check vote/00-nissl-vote%d.vol',itr), 'w'); % fwrite(fid, single(Ivote), '*single'); % fclose(fid); % end % if itr ==3 % fid = fopen(sprintf('128-128-128/00-nissl-vote%d.vol',itr), 'w'); % fwrite(fid, single(Ivote), '*single'); % fclose(fid); % Ivote8 = (Ivote); % end % if itr ==9 % fid = fopen(sprintf('128-128-128/00-nissl-vote%d.vol',itr), 'w'); % fwrite(fid, single(Ivote), '*single'); % fclose(fid); % end % if itr ==10 % fid = fopen(sprintf('128-128-128/00-nissl-vote%d.vol',itr), 'w'); % fwrite(fid, single(Ivote), '*single'); % fclose(fid); % end % Ivote1 = single(Ivote); % fwrite(fid, Ivote1, '*single'); % % % elseif itr ==2 % Ivote2 = single(Ivote); % fwrite(fid, Ivote2, '*single'); % % % elseif itr ==3 % Ivote3 = single(Ivote); % fwrite(fid, Ivote3, '*single'); % % % elseif itr ==4 % Ivote4 = single(Ivote); % fwrite(fid, Ivote4, '*single'); % % % elseif itr == 5 % Ivote5 = single(Ivote); % fwrite(fid, Ivote5, '*single'); % % % elseif itr == 6 % fwrite(fid, single(Ivote), '*single'); % elseif itr == 7 % fwrite(fid, single(Ivote), '*single'); % elseif itr == 8 % fwrite(fid, single(Ivote), '*single'); % elseif itr == 9 % fwrite(fid, single(Ivote), '*single'); % end fclose(fid); t_v1 = toc; disp(['voting done. time =',num2str(t_v1)]); sum(validPixels(:)==0) % update the voting direction if ang>=d_ang tic; Igrad_x = zeros(Isize); Igrad_y = zeros(Isize); Igrad_z = zeros(Isize); vv=0; nrv = zeros; for v = 1: nV % coordinates of the current voter vx = Itx(v); vy = Ity(v); vz = Itz(v); if (validPixels(v)>0) %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; update_rx(vx,vy,vz) = 2 * abs(g_px - vx); update_ry(vx,vy,vz) = 2 * abs(g_py - vy); update_rz(vx,vy,vz) = 2 * abs(g_pz - vz); else vv=vv+1; nrv(vv) = v; end end ni = nnz(nrv); if ni>0 for j=1:ni It(Itx(nrv(j)), Ity(nrv(j)), Itz(nrv(j)))=0; end clear Itx; clear Ity; clear Itz; [Itx,Ity,Itz] = ind2sub(size(It),find(It)); Vi =(find(It)); nV = nnz(It) end update_rx(update_rx>rmax(1)) = rmax(1); update_ry(update_ry>rmax(2)) = rmax(2); update_rz(update_rz>rmax(3)) = rmax(3); 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 = 0; conn = [5 5 5]; Icenter = local_max(Ivote, conn, t); fidc = fopen(sprintf('std5.5-r10.10-v8-update-rmax/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