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 = 0; 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 = '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_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; %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_pi0 = zeros(npts); idxx=0; for id_gl=1:npts if(global_px(id_gl)>0 && global_px(id_gl)<=X && global_py(id_gl)>0 && global_py(id_gl)<=Y && global_pz(id_gl)>0 && global_pz(id_gl)<=Z) idxx = idxx + 1; global_pi0(idxx) = sub2ind(Isize, global_px(id_gl), global_py(id_gl), global_pz(id_gl)); end end global_pi = global_pi0(global_pi0~=0); npts = numel(global_pi); % 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 = 2000; conn = [5 5 5]; Icenter = local_max(Ivote, conn, t); fidc = fopen(sprintf('out%d.vol',t), 'w'); fwrite(fidc, single(Icenter), '*single'); fclose(fidc); tot = toc(total) nnz(Icenter) % [cxx1, cyy1, czz1] = ind2sub(size(Icenter),find(Icenter)); % % center = Ivote1; % % center(center