From da849ea2cf240713b77e81b2b16fdfd30389bbf7 Mon Sep 17 00:00:00 2001 From: Laila Saadatifard Date: Thu, 8 Oct 2015 18:30:13 -0500 Subject: [PATCH] fix the bug about the gradient anf make syntethic image for test --- matlab/main.asv | 399 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------------------------------------------------------------------------------------------------------------------------------------------------- matlab/main.m | 401 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------------------------------------------------------------------------------------------------------------------------------------------------- matlab/syn_Img.m | 23 +++++++++++++++++++++++ 3 files changed, 521 insertions(+), 302 deletions(-) create mode 100644 matlab/syn_Img.m diff --git a/matlab/main.asv b/matlab/main.asv index bb8acd2..bc09600 100644 --- a/matlab/main.asv +++ b/matlab/main.asv @@ -1,4 +1,3 @@ -tic; clc; clear; disp('***************** NEW RUN *********************'); @@ -6,193 +5,291 @@ total = tic; % ******* Initialize voting parameters ************************************** -rmax = 9; %maximum radius of the cell -phi_deg = 25.1; %half the angular range of the voting area -phi = phi_deg * pi / 180; -iter = 5; %number of voting iterations +rmax = 7; %maximum radius of the cell +ang_deg = 15.1; %half the angular range of the voting area +ang = ang_deg * pi / 180; +iter = 2; %number of voting iterations t0 = 1.2; %threshold color sigma = [4, 4, 2]; -final_t = 75; - +% t = 0.1; +d_ang= ang / (iter); % ******** Testing parameters ****************************************** -p = [50, 20, 50] ; -ps = [100, 100, 50] - 1; - -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; +% p = [100, 50, 100]; +% ps = [200, 200, 100]; +ps = [100, 100, 25]; +r = 7; +I = syn_Img(r , 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 -I_blur = gauss_blur3d(I, sigma); +Iblur = gauss_blur3d(I, sigma); -% compute the gradient -[Igrad_x, Igrad_y, Igrad_z] = gradient(I_blur); +%%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); -%%crop out a small subregion of I -Isub_x = Igrad_x(p(1):p(1)+ps(1), p(2):p(2)+ps(2), p(3):p(3)+ps(3)); -Isub_y = Igrad_y(p(1):p(1)+ps(1), p(2):p(2)+ps(2), p(3):p(3)+ps(3)); -Isub_z = Igrad_z(p(1):p(1)+ps(1), p(2):p(2)+ps(2), p(3):p(3)+ps(3)); +% +% compute the gradient +[Igrad_x, Igrad_y, Igrad_z] = gradient(Iblur); %calculate the gradient magnitude -Imag = sqrt(Isub_x .^ 2 + Isub_y .^ 2 + Isub_z .^2); +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), -imshow( squeeze( I(:, :, ceil(size(I, 3)/2)) ) ./ 255 ); -%colormap(gray); +imagesc(squeeze(I(:, :, ceil(size(I, 3)/2)))); subplot(3, 3, 2), -imshow( squeeze( I(:, ceil(size(I, 2)/2), :) ) ./ 255 ); +imagesc(squeeze(I(:, ceil(size(I, 2)/2), :))); subplot(3, 3, 3), -imshow( squeeze( I(ceil(size(I, 1)/2), :, :) ) ./ 255 ); - +imagesc(squeeze(I( ceil(size(I, 1)/2), :, :))); +% subplot(3, 3, 4), -imagesc( squeeze( Imag(:, :, ceil(size(Imag, 3)/2)) ) ); +imagesc(squeeze(Iblur(:, :, ceil(size(Iblur, 3)/2)))); subplot(3, 3, 5), -imagesc( squeeze( Imag(:, ceil(size(Imag, 2)/2), :) ) ); +imagesc(squeeze(Iblur(:, ceil(size(Iblur, 2)/2), :))); subplot(3, 3, 6), -imagesc( squeeze( Imag(ceil(size(Imag, 1)/2), :, :) ) ); +imagesc(squeeze(Iblur(ceil(size(Iblur, 1)/2), :, :))); subplot(3, 3, 7), -imagesc( squeeze( It(:, :, ceil(size(It, 3)/2)) ) ); +imagesc(squeeze(Imag(:, :, ceil(size(Imag, 3)/2)))); subplot(3, 3, 8), -imagesc( squeeze( It(:, ceil(size(It, 2)/2), :) ) ); +imagesc(squeeze(Imag(:, ceil(size(Imag, 2)/2), :))); subplot(3, 3, 9), -imagesc( squeeze( It(ceil(size(It, 1)/2), :, :) ) ); +imagesc(squeeze(Imag(ceil(size(Imag, 1)/2), :, :))); -% -% thresholding -It_x = size(It,1); -It_y = size(It,2); -It_z = size(It,3); -% -It(1:rmax, :, :) = 0; -It(It_x - rmax:It_x, :,:) = 0; -S(:, 1:rmax, :) = 0; -S(:, It_y - rmax:It_y,:) = 0; -S(:, :, 1:rmax) = 0; -S(:,:, It_z - rmax:It_z) = 0; +colormap(gray); -% -[Itx,Ity,Itz] = ind2sub(size(It),find(It)); +%% vote +tic; +%for each iteration (in iterative voting) +for itr = 1 : iter -nV = nnz(It); + %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 + dir_mag = sqrt (gx^2 + gy^2 + gz^2); + + %calculate the normalized gradient direction + ngx = gx / dir + + %gdir_x = Igrad_x(vx, vy, vz) / dir_mag; + %gdir_y = Igrad_y(vx, vy, vz) / dir_mag; + %gdir_z = Igrad_z(vx, vy, vz) / dir_mag; + + %calculate the angle between the voter direction and the pixel direction + cos_diff = LV_x .* gdir_x + LV_y .* gdir_y + LV_z .* gdir_z; + 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); + %global_pi = (global_pz-1)*ps(1)*ps(2) + (global_py-1)*ps(1) + global_px; + 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)]); + + figure, imagesc(squeeze(Ivote(:, :, ceil(size(Ivote, 3)/2),:))); colormap(gray); + + % 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 + [gx, gy, gz] = ind2sub(size(Ivote), g_v_prime(v,local_max_idx)); + + %compute the vector from the voter position to this position + dx = gx - vx; + dy = gy - vy; + dz = gz - vz; + Igrad_x(vx, vy, vz) = dx; + Igrad_y(vx, vy, vz) = dy; + Igrad_z(vx, vy, vz) = dz; + + end + + tdir1 = toc; + display (['updating dir done. time = ', num2str(tdir1)]); + ang = ang - d_ang; + end + +%% +t = 10; +out = Ivote; +out(out=t) = 255; +out = imregionalmax(out); -% create a meshgrid describing coordinates relative to the voter position -range = -rmax:rmax; -[mx, my, mz] = meshgrid(range, range, range); -m_mag = sqrt(mx.^2 + my.^2 + mz.^2); -%create a mask for the voting area -M_dist = mx.^2 + my.^2 + mz.^2 < rmax^2; +out1(:,:,:,1) = mat2gray(Iblur); +out1(:,:,:,2) = mat2gray(out); +out1(:,:,:,3) = mat2gray(Iblur); -pxPerRow = size(I_blur,1); -pxPerCol = size(I_blur,2); -validPoints = zeros(nV,1); -g_v_prime = zeros(nV, ceil(rmax^2*phi/3)); -%% +figure(2); -disp('done.'); -%% vote -for itr = 1 : iter - Ivote = zeros(size(Iblur)); - for v = 1: nV - % coordinates of the current voter - vx = Itx(v); - vy = Ity(v); - vz = Itz(v); - - vgrad = [Isub_x(vx,vy,vz), Isub_y(vx,vy,vz), Isub_z(vx,vy,vz)]; - %determine the gradient magnitude at the current voter location - vmag = Imag(vx,vy,vz); - ang_diff = zeros(2*rmax+1, 2*rmax+1, 2*rmax+1); - for i = 1 : 2*rmax+1 - for j = 1 : 2*rmax+1 - for k = 1 : 2*rmax+1 - a = [mx(i,j,k), my(i,j,k), mz(i,j,k)]; - c = dot(vgrad, a); - ang_diff(i,j,k) = acos(c/(m_mag(i,j,k)* vmag )); - end - end - end - - M_diff = min(2*pi - ang_diff,ang_diff) 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), -imshow( squeeze( I(:, :, ceil(size(I, 3)/2)) ) ./ 255 ); -%colormap(gray); +imagesc(squeeze(I(:, :, ceil(size(I, 3)/2)))); subplot(3, 3, 2), -imshow( squeeze( I(:, ceil(size(I, 2)/2), :) ) ./ 255 ); +imagesc(squeeze(I(:, ceil(size(I, 2)/2), :))); subplot(3, 3, 3), -imshow( squeeze( I(ceil(size(I, 1)/2), :, :) ) ./ 255 ); - +imagesc(squeeze(I( ceil(size(I, 1)/2), :, :))); +% subplot(3, 3, 4), -imagesc( squeeze( Imag(:, :, ceil(size(Imag, 3)/2)) ) ); +imagesc(squeeze(Iblur(:, :, ceil(size(Iblur, 3)/2)))); subplot(3, 3, 5), -imagesc( squeeze( Imag(:, ceil(size(Imag, 2)/2), :) ) ); +imagesc(squeeze(Iblur(:, ceil(size(Iblur, 2)/2), :))); subplot(3, 3, 6), -imagesc( squeeze( Imag(ceil(size(Imag, 1)/2), :, :) ) ); +imagesc(squeeze(Iblur(ceil(size(Iblur, 1)/2), :, :))); subplot(3, 3, 7), -imagesc( squeeze( It(:, :, ceil(size(It, 3)/2)) ) ); +imagesc(squeeze(Imag(:, :, ceil(size(Imag, 3)/2)))); subplot(3, 3, 8), -imagesc( squeeze( It(:, ceil(size(It, 2)/2), :) ) ); +imagesc(squeeze(Imag(:, ceil(size(Imag, 2)/2), :))); subplot(3, 3, 9), -imagesc( squeeze( It(ceil(size(It, 1)/2), :, :) ) ); +imagesc(squeeze(Imag(ceil(size(Imag, 1)/2), :, :))); -% -% thresholding -It_x = size(It,1); -It_y = size(It,2); -It_z = size(It,3); -% -It(1:rmax, :, :) = 0; -It(It_x - rmax:It_x, :,:) = 0; -S(:, 1:rmax, :) = 0; -S(:, It_y - rmax:It_y,:) = 0; -S(:, :, 1:rmax) = 0; -S(:,:, It_z - rmax:It_z) = 0; +colormap(gray); -% -[Itx,Ity,Itz] = ind2sub(size(It),find(It)); +%% vote +tic; +%for each iteration (in iterative voting) +for itr = 1 : iter -nV = nnz(It); + %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; + + %gdir_x = Igrad_x(vx, vy, vz) / dir_mag; + %gdir_y = Igrad_y(vx, vy, vz) / dir_mag; + %gdir_z = Igrad_z(vx, vy, vz) / dir_mag; + + %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); + %global_pi = (global_pz-1)*ps(1)*ps(2) + (global_py-1)*ps(1) + global_px; + 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)]); + + figure, imagesc(squeeze(Ivote(:, :, ceil(size(Ivote, 3)/2),:))); colormap(gray); + + % 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 + [gx, gy, gz] = ind2sub(size(Ivote), g_v_prime(v,local_max_idx)); + + %compute the vector from the voter position to this position + dx = gx - vx; + dy = gy - vy; + dz = gz - vz; + Igrad_x(vx, vy, vz) = dx; + Igrad_y(vx, vy, vz) = dy; + Igrad_z(vx, vy, vz) = dz; + + end + + tdir1 = toc; + display (['updating dir done. time = ', num2str(tdir1)]); + ang = ang - d_ang; + end + +%% +t = 10; +out = Ivote; +out(out=t) = 255; +out = imregionalmax(out); -% create a meshgrid describing coordinates relative to the voter position -range = -rmax:rmax; -[mx, my, mz] = meshgrid(range, range, range); -m_mag = sqrt(mx.^2 + my.^2 + mz.^2); -%create a mask for the voting area -M_dist = mx.^2 + my.^2 + mz.^2 < rmax^2; +out1(:,:,:,1) = mat2gray(Iblur); +out1(:,:,:,2) = mat2gray(out); +out1(:,:,:,3) = mat2gray(Iblur); -pxPerRow = size(I_blur,1); -pxPerCol = size(I_blur,2); -validPoints = zeros(nV,1); -g_v_prime = zeros(nV, ceil(rmax^2*phi/3)); -%% +figure(2); -disp('done.'); -%% vote -for itr = 1 : iter - Ivote = zeros(size(Iblur)); - for v = 1: nV - % coordinates of the current voter - vx = Itx(v); - vy = Ity(v); - vz = Itz(v); - - vgrad = [Isub_x(vx,vy,vz), Isub_y(vx,vy,vz), Isub_z(vx,vy,vz)]; - %determine the gradient magnitude at the current voter location - vmag = Imag(vx,vy,vz); - ang_diff = zeros(2*rmax+1, 2*rmax+1, 2*rmax+1); - for i = 1 : 2*rmax+1 - for j = 1 : 2*rmax+1 - for k = 1 : 2*rmax+1 - a = [mx(i,j,k), my(i,j,k), mz(i,j,k)]; - c = dot(vgrad, a); - ang_diff(i,j,k) = acos(c/(m_mag(i,j,k)* vmag )); - end - end - end - - M_diff = min(2*pi - ang_diff,ang_diff)