segmentation / ivote3

Commit da849ea2cf240713b77e81b2b16fdfd30389bbf7

Authored by Laila Saadatifard
1 parent 76a0f6f9

fix the bug about the gradient anf make syntethic image for test

Showing 3 changed files with 521 additions and 302 deletions   Show diff stats
Inline Side-by-side
1 -tic;  
2 clc; 1 clc;
3 clear; 2 clear;
4 disp('***************** NEW RUN *********************'); 3 disp('***************** NEW RUN *********************');
@@ -6,193 +5,291 @@ total = tic; @@ -6,193 +5,291 @@ total = tic;
6 5
7 6
8 % ******* Initialize voting parameters ************************************** 7 % ******* Initialize voting parameters **************************************
9 -rmax = 9; %maximum radius of the cell  
10 -phi_deg = 25.1; %half the angular range of the voting area  
11 -phi = phi_deg * pi / 180;  
12 -iter = 5; %number of voting iterations 8 +rmax = 7; %maximum radius of the cell
  9 +ang_deg = 15.1; %half the angular range of the voting area
  10 +ang = ang_deg * pi / 180;
  11 +iter = 2; %number of voting iterations
13 t0 = 1.2; %threshold color 12 t0 = 1.2; %threshold color
14 sigma = [4, 4, 2]; 13 sigma = [4, 4, 2];
15 -final_t = 75;  
16 - 14 +% t = 0.1;
  15 +d_ang= ang / (iter);
17 16
18 % ******** Testing parameters ****************************************** 17 % ******** Testing parameters ******************************************
19 -p = [50, 20, 50] ;  
20 -ps = [100, 100, 50] - 1;  
21 -  
22 -volfile = 'img\nissl-rat.vol';  
23 -fid = fopen(volfile); % open the file that include the image  
24 -S = fread(fid, 3, 'int32');  
25 -X = S(1);  
26 -Y = S(2);  
27 -Z = S(3);  
28 -  
29 -% load the VOL data into a 2D matrix  
30 -I = fread(fid,[X Y*Z], 'uint8');  
31 -fclose(fid);  
32 -  
33 -%change this to a 3D matrix  
34 -I = reshape(I, [X, Y, Z]);  
35 -  
36 -%invert the intensity  
37 -I = 255 - I; 18 +% p = [100, 50, 100];
  19 +% ps = [200, 200, 100];
  20 +ps = [100, 100, 25];
  21 +r = 7;
  22 +I = syn_Img(r , ps);
  23 +% volfile = 'img\nissl-rat.vol';
  24 +% fid = fopen(volfile); % open the file that include the image
  25 +% S = fread(fid, 3, 'int32');
  26 +% X = S(1);
  27 +% Y = S(2);
  28 +% Z = S(3);
  29 +%
  30 +% % load the VOL data into a 2D matrix
  31 +% I = fread(fid,[X Y*Z], 'uint8');
  32 +% fclose(fid);
  33 +%
  34 +% %change this to a 3D matrix
  35 +% I = (reshape(I, [X, Y, Z]));
  36 +%
  37 +% % invert the intensity
  38 +% I = 255 - I;
38 39
39 %perform a gaussian blur 40 %perform a gaussian blur
40 -I_blur = gauss_blur3d(I, sigma); 41 +Iblur = gauss_blur3d(I, sigma);
41 42
42 -% compute the gradient  
43 -[Igrad_x, Igrad_y, Igrad_z] = gradient(I_blur); 43 +%%crop out a small subregion of I and Iblur
  44 +% Iblur = Iblur(p(1):p(1)+ps(1)-1, p(2):p(2)+ps(2)-1, p(3):p(3)+ps(3)-1);
  45 +% I = I(p(1):p(1)+ps(1)-1, p(2):p(2)+ps(2)-1, p(3):p(3)+ps(3)-1);
44 46
45 -%%crop out a small subregion of I  
46 -Isub_x = Igrad_x(p(1):p(1)+ps(1), p(2):p(2)+ps(2), p(3):p(3)+ps(3));  
47 -Isub_y = Igrad_y(p(1):p(1)+ps(1), p(2):p(2)+ps(2), p(3):p(3)+ps(3));  
48 -Isub_z = Igrad_z(p(1):p(1)+ps(1), p(2):p(2)+ps(2), p(3):p(3)+ps(3)); 47 +%
  48 +% compute the gradient
  49 +[Igrad_x, Igrad_y, Igrad_z] = gradient(Iblur);
49 50
50 %calculate the gradient magnitude 51 %calculate the gradient magnitude
51 -Imag = sqrt(Isub_x .^ 2 + Isub_y .^ 2 + Isub_z .^2); 52 +Imag = sqrt(Igrad_x .^ 2 + Igrad_y .^ 2 + Igrad_z .^2);
  53 +I_theta1 = acos(Igrad_x./Imag);
  54 +I_theta2 = acos(Igrad_y./Imag);
  55 +I_theta3 = acos(Igrad_z./Imag);
52 56
53 %set a threshold for the gradient magnitude 57 %set a threshold for the gradient magnitude
54 It = Imag > t0; 58 It = Imag > t0;
55 59
  60 +%Set the boundaries of the threshold image to zero
  61 +It(1:rmax, :, :) = 0;
  62 +It(ps(1) - rmax:ps(1), :,:) = 0;
  63 +It(:, 1:rmax, :) = 0;
  64 +It(:, ps(2) - rmax:ps(2),:) = 0;
  65 +It(:, :, 1:rmax) = 0;
  66 +It(:,:, ps(3) - rmax:ps(3)) = 0;
  67 +
  68 +%get the indices of all of the nonzero values in the threshold image
  69 +% (voter positions)
  70 +[Itx,Ity,Itz] = ind2sub(size(It),find(It));
  71 +Vi = find(It);
  72 +nV = nnz(It);
  73 +
  74 +% create a meshgrid describing coordinates relative to the voter position
  75 +range = -rmax:rmax; %create an array of values between -rmax and rmax
  76 +[mx, my, mz] = meshgrid(range, range, range); %create a template describing local pixel position in a small cube
  77 +m_mag = sqrt(mx.^2 + my.^2 + mz.^2); %create a template describing the distance from the center of a small cube
56 78
  79 +% create a mask for the voting area
  80 +M_dist = m_mag <= rmax; %mask for the voting area distance (all values < rmax from the center)
  81 +
  82 +%calculate the direction vector between a pixel and voter
  83 +LV_x = mx./m_mag;
  84 +LV_y = my./m_mag;
  85 +LV_z = mz./m_mag;
  86 +
  87 +%number of pixels in the voting area of each voter (initialize to zero)
  88 +validPixels = (zeros(nV,1));
  89 +
  90 +%indices of pixels in the voting area of each voter
  91 +% indices reference the 3D image
  92 +g_v_prime = (zeros(nV, (rmax^3)));
  93 +
  94 +%--------------Display current data-------------------------
57 subplot(3, 3, 1), 95 subplot(3, 3, 1),
58 -imshow( squeeze( I(:, :, ceil(size(I, 3)/2)) ) ./ 255 );  
59 -%colormap(gray); 96 +imagesc(squeeze(I(:, :, ceil(size(I, 3)/2))));
60 97
61 subplot(3, 3, 2), 98 subplot(3, 3, 2),
62 -imshow( squeeze( I(:, ceil(size(I, 2)/2), :) ) ./ 255 ); 99 +imagesc(squeeze(I(:, ceil(size(I, 2)/2), :)));
63 100
64 subplot(3, 3, 3), 101 subplot(3, 3, 3),
65 -imshow( squeeze( I(ceil(size(I, 1)/2), :, :) ) ./ 255 );  
66 - 102 +imagesc(squeeze(I( ceil(size(I, 1)/2), :, :)));
  103 +%
67 subplot(3, 3, 4), 104 subplot(3, 3, 4),
68 -imagesc( squeeze( Imag(:, :, ceil(size(Imag, 3)/2)) ) ); 105 +imagesc(squeeze(Iblur(:, :, ceil(size(Iblur, 3)/2))));
69 106
70 subplot(3, 3, 5), 107 subplot(3, 3, 5),
71 -imagesc( squeeze( Imag(:, ceil(size(Imag, 2)/2), :) ) ); 108 +imagesc(squeeze(Iblur(:, ceil(size(Iblur, 2)/2), :)));
72 109
73 subplot(3, 3, 6), 110 subplot(3, 3, 6),
74 -imagesc( squeeze( Imag(ceil(size(Imag, 1)/2), :, :) ) ); 111 +imagesc(squeeze(Iblur(ceil(size(Iblur, 1)/2), :, :)));
75 112
76 subplot(3, 3, 7), 113 subplot(3, 3, 7),
77 -imagesc( squeeze( It(:, :, ceil(size(It, 3)/2)) ) ); 114 +imagesc(squeeze(Imag(:, :, ceil(size(Imag, 3)/2))));
78 115
79 subplot(3, 3, 8), 116 subplot(3, 3, 8),
80 -imagesc( squeeze( It(:, ceil(size(It, 2)/2), :) ) ); 117 +imagesc(squeeze(Imag(:, ceil(size(Imag, 2)/2), :)));
81 118
82 subplot(3, 3, 9), 119 subplot(3, 3, 9),
83 -imagesc( squeeze( It(ceil(size(It, 1)/2), :, :) ) ); 120 +imagesc(squeeze(Imag(ceil(size(Imag, 1)/2), :, :)));
84 121
85 -%  
86 -% thresholding  
87 -It_x = size(It,1);  
88 -It_y = size(It,2);  
89 -It_z = size(It,3);  
90 -%  
91 -It(1:rmax, :, :) = 0;  
92 -It(It_x - rmax:It_x, :,:) = 0;  
93 -S(:, 1:rmax, :) = 0;  
94 -S(:, It_y - rmax:It_y,:) = 0;  
95 -S(:, :, 1:rmax) = 0;  
96 -S(:,:, It_z - rmax:It_z) = 0; 122 +colormap(gray);
97 123
98 -%  
99 -[Itx,Ity,Itz] = ind2sub(size(It),find(It)); 124 +%% vote
  125 +tic;
  126 +%for each iteration (in iterative voting)
  127 +for itr = 1 : iter
100 128
101 -nV = nnz(It); 129 + %initialize the vote image to zero
  130 + Ivote = (zeros(size(I)));
  131 +
  132 + %for each voter (nonzero pixels in the threshold image It)
  133 + for v = 1: nV
  134 +
  135 + %get the cartesian coordinates of the voter v in the main image I
  136 + vx = Itx(v);
  137 + vy = Ity(v);
  138 + vz = Itz(v);
  139 + vi = Vi(v);
  140 +
  141 + %retreive the gradient magnitude at the voter position
  142 + %vmag = Imag(vx,vy,vz);
  143 + vmag = Imag(vi);
  144 +
  145 + %retrieve the gradient
  146 + gx = Igrad_x(vi);
  147 + gy = Igrad_y(vi);
  148 + gz = Igrad_z(vi);
  149 +
  150 + %calculate the gradient magnitude
  151 + dir_mag = sqrt (gx^2 + gy^2 + gz^2);
  152 +
  153 + %calculate the normalized gradient direction
  154 + ngx = gx / dir
  155 +
  156 + %gdir_x = Igrad_x(vx, vy, vz) / dir_mag;
  157 + %gdir_y = Igrad_y(vx, vy, vz) / dir_mag;
  158 + %gdir_z = Igrad_z(vx, vy, vz) / dir_mag;
  159 +
  160 + %calculate the angle between the voter direction and the pixel direction
  161 + cos_diff = LV_x .* gdir_x + LV_y .* gdir_y + LV_z .* gdir_z;
  162 + ang_diff = acos(cos_diff);
  163 +
  164 + %create an angular mask for the voting area
  165 + M_angle = cos_diff > cos(ang);
  166 +
  167 + %combine the two masks to mask out the voting angle
  168 + M = M_angle .* M_dist;
  169 +
  170 + % get the coordinates of each pixel in the final voter mask M
  171 + pi = find(M);
  172 +
  173 + %calculate the number of pixels in the voting region
  174 + npts = nnz(M);
  175 + validPixels(v) = npts;
  176 +
  177 + %convert every index in the voting area from a local 3D index to a global 3D index (into the original image I)
  178 + global_px = vx + mx(pi);
  179 + global_py = vy + my(pi);
  180 + global_pz = vz + mz(pi);
  181 +
  182 + %convert the global 3D index of each point into a global 1D index
  183 + %global_pi = sub2ind(ps, global_px, global_py, global_pz);
  184 + %global_pi = (global_pz-1)*ps(1)*ps(2) + (global_py-1)*ps(1) + global_px;
  185 + global_pi = sub2ind(ps, global_px, global_py, global_pz);
  186 +
  187 + g_v_prime (v, 1:npts) = global_pi;
  188 +
  189 +
  190 + Ivote( global_pi ) = Ivote( global_pi ) + vmag;
  191 +
  192 + end
  193 +
  194 + if itr ==1
  195 + Ivote1 = Ivote;
  196 +
  197 + elseif itr ==2
  198 + Ivote2 = Ivote;
  199 +
  200 + elseif itr ==3
  201 + Ivote3 = Ivote;
  202 +
  203 + elseif itr ==4
  204 + Ivote4 = Ivote;
  205 +
  206 + elseif itr == 5
  207 + Ivote5 = Ivote;
  208 + end
  209 + t_v1 = toc;
  210 + disp(['voting done. time =',num2str(t_v1)]);
  211 +
  212 + figure, imagesc(squeeze(Ivote(:, :, ceil(size(Ivote, 3)/2),:))); colormap(gray);
  213 +
  214 + % update the voting direction
  215 + tic;
  216 + for v = 1: nV
  217 + % coordinates of the current voter
  218 + vx = Itx(v);
  219 + vy = Ity(v);
  220 + vz = Itz(v);
  221 +
  222 + %get the local value of the voting image
  223 + local_Ivote = Ivote(g_v_prime(v,1:validPixels(v)));
  224 +
  225 + %find the index of the maximum value
  226 + [~, local_max_idx] = max(local_Ivote);
  227 +
  228 + %convert this into a global subscript
  229 + [gx, gy, gz] = ind2sub(size(Ivote), g_v_prime(v,local_max_idx));
  230 +
  231 + %compute the vector from the voter position to this position
  232 + dx = gx - vx;
  233 + dy = gy - vy;
  234 + dz = gz - vz;
  235 + Igrad_x(vx, vy, vz) = dx;
  236 + Igrad_y(vx, vy, vz) = dy;
  237 + Igrad_z(vx, vy, vz) = dz;
  238 +
  239 + end
  240 +
  241 + tdir1 = toc;
  242 + display (['updating dir done. time = ', num2str(tdir1)]);
  243 + ang = ang - d_ang;
  244 + end
  245 +
102 246
  247 +%%
  248 +t = 10;
  249 +out = Ivote;
  250 +out(out<t) = 0;
  251 +% out(out>=t) = 255;
  252 +out = imregionalmax(out);
103 253
104 -% create a meshgrid describing coordinates relative to the voter position  
105 -range = -rmax:rmax;  
106 -[mx, my, mz] = meshgrid(range, range, range);  
107 -m_mag = sqrt(mx.^2 + my.^2 + mz.^2);  
108 -%create a mask for the voting area  
109 -M_dist = mx.^2 + my.^2 + mz.^2 < rmax^2; 254 +out1(:,:,:,1) = mat2gray(Iblur);
  255 +out1(:,:,:,2) = mat2gray(out);
  256 +out1(:,:,:,3) = mat2gray(Iblur);
110 257
111 -pxPerRow = size(I_blur,1);  
112 -pxPerCol = size(I_blur,2);  
113 -validPoints = zeros(nV,1);  
114 258
115 -g_v_prime = zeros(nV, ceil(rmax^2*phi/3));  
116 259
117 -%% 260 +figure(2);
118 261
119 -disp('done.');  
120 -%% vote  
121 -for itr = 1 : iter  
122 - Ivote = zeros(size(Iblur));  
123 - for v = 1: nV  
124 - % coordinates of the current voter  
125 - vx = Itx(v);  
126 - vy = Ity(v);  
127 - vz = Itz(v);  
128 -  
129 - vgrad = [Isub_x(vx,vy,vz), Isub_y(vx,vy,vz), Isub_z(vx,vy,vz)];  
130 - %determine the gradient magnitude at the current voter location  
131 - vmag = Imag(vx,vy,vz);  
132 - ang_diff = zeros(2*rmax+1, 2*rmax+1, 2*rmax+1);  
133 - for i = 1 : 2*rmax+1  
134 - for j = 1 : 2*rmax+1  
135 - for k = 1 : 2*rmax+1  
136 - a = [mx(i,j,k), my(i,j,k), mz(i,j,k)];  
137 - c = dot(vgrad, a);  
138 - ang_diff(i,j,k) = acos(c/(m_mag(i,j,k)* vmag ));  
139 - end  
140 - end  
141 - end  
142 -  
143 - M_diff = min(2*pi - ang_diff,ang_diff)<phi;  
144 - %compute the final mask  
145 - M = (M_dist .* M_diff);  
146 -  
147 - % get the coordinates of each pixel in the final mask  
148 - [vx_prime,vy_prime,vz_prime] = ind2sub(size(M),find(M));  
149 -  
150 - %transform the local coordinates of the pixels in the voting region  
151 - %to the global coordinates of the image  
152 -  
153 - npts = numel(vx_prime);  
154 - validPoints(v) = npts;  
155 -  
156 - g_v_prime(v,1:npts) = vx + (vx_prime - (rmax + 1)) + (vy + (vy_prime - (rmax + 1))-1).*pxPerRow + (vz + (vz_prime - (rmax + 1))-1).*pxPerRow*pxPerCol;  
157 -  
158 - for n=1: npts  
159 - %  
160 - Ivote( g_v_prime(v,n)) = Ivote( g_v_prime(v,n)) + vmag;  
161 - end  
162 -  
163 - end  
164 -  
165 -disp('voting done.');  
166 - %% update the voting direction  
167 -  
168 - for v = 1: nV  
169 - % coordinates of the current voter  
170 - vx = sx(v);  
171 - vy = sy(v);  
172 - vz = sz(v);  
173 -  
174 - %get the local value of the voting image  
175 - local_Ivote = Ivote(g_v_prime(v,1:validPoints(v)));  
176 -  
177 - %find the index of the maximum value  
178 - [~, local_max_idx] = max(local_Ivote);  
179 -  
180 - %convert this into a global subscript  
181 - gz = ceil(g_v_prime(v,local_max_idx)/(pxPerRow*pxPerCol)) ;  
182 - gy = ceil((g_v_prime(v,local_max_idx)-(gz - 1)*pxPerRow*pxPerCol)/pxPerRow);  
183 - gx = g_v_prime(v,local_max_idx)- pxPerRow*((gy - 1) + (gz - 1) * pxPerCol);  
184 -  
185 - %compute the vector from the voter position to this position  
186 - new_vx = gx - vx;  
187 - new_vy = gy - vy;  
188 - new_vz = gz - vz;  
189 - Igrad_x(vx,vy,vz) = new_vx;  
190 - Igrad_y(vx,vy,vz) = new_vy;  
191 - Igrad_z(vx,vy,vz) = new_vz;  
192 - end  
193 -end  
194 -cell_center = Ivote;  
195 -fido = fopen('out\outImg.vol', 'w+');  
196 -fwrite(fido, cell_center);  
197 -fclose(fido);  
198 -t1 = toc 262 +subplot(2, 3, 1),
  263 +imagesc(squeeze(Iblur(:, :, ceil(size(Iblur, 3)/2))));
  264 +
  265 +subplot(2, 3, 2),
  266 +imagesc(squeeze(Iblur(:, ceil(size(Iblur, 2)/2), :)));
  267 +
  268 +subplot(2, 3, 3),
  269 +imagesc(squeeze(Iblur( ceil(size(Iblur, 1)/2), :, :)));
  270 +
  271 +subplot(2, 3, 4),
  272 +imagesc( squeeze( out1(:, :, ceil(size(Ivote, 3)/2),:) ) );
  273 +
  274 +subplot(2, 3, 5),
  275 +imagesc(squeeze(out1(:, ceil(size(Ivote, 2)/2), :, :)));
  276 +
  277 +subplot(2, 3, 6),
  278 +imagesc(squeeze(out1( ceil(size(Ivote, 1)/2), :, :, :)));
  279 +
  280 +colormap(gray);
  281 +%%
  282 +% figure, imagesc(squeeze(Iblur(:, :, ceil(size(Iblur, 3)/2)))); colormap(gray);
  283 +% figure, imagesc(squeeze(Ivote(:, :, ceil(size(Ivote, 3)/2),:))); colormap(gray);
  284 +% figure, imagesc(squeeze(Iblur(:, ceil(size(Iblur, 2)/2), :))); colormap(gray);
  285 +% figure, imagesc(squeeze(I(:, ceil(size(Ivote, 2)/2), :))); colormap(gray);
  286 +% % figure, imagesc(squeeze(Iblur( ceil(size(Iblur, 1)/2), :, :))); colormap(gray);
  287 +% figure, imagesc(squeeze(Ivote( ceil(size(Ivote, 1)/2), :, :))); colormap(gray);
  288 +%
  289 +fid1 = fopen('ivote1.vol', 'w');
  290 +fwrite(fid1, Ivote1);
  291 +fclose(fid1);
  292 +
  293 +% fid0 = fopen('iblur.vol', 'w');
  294 +% fwrite(fid0, Iblur);
  295 +% fclose(fid0);
1 -tic;  
2 clc; 1 clc;
3 clear; 2 clear;
4 disp('***************** NEW RUN *********************'); 3 disp('***************** NEW RUN *********************');
@@ -6,193 +5,293 @@ total = tic; @@ -6,193 +5,293 @@ total = tic;
6 5
7 6
8 % ******* Initialize voting parameters ************************************** 7 % ******* Initialize voting parameters **************************************
9 -rmax = 9; %maximum radius of the cell  
10 -phi_deg = 25.1; %half the angular range of the voting area  
11 -phi = phi_deg * pi / 180;  
12 -iter = 5; %number of voting iterations 8 +rmax = 7; %maximum radius of the cell
  9 +ang_deg = 15.1; %half the angular range of the voting area
  10 +ang = ang_deg * pi / 180;
  11 +iter = 2; %number of voting iterations
13 t0 = 1.2; %threshold color 12 t0 = 1.2; %threshold color
14 sigma = [4, 4, 2]; 13 sigma = [4, 4, 2];
15 -final_t = 75;  
16 - 14 +% t = 0.1;
  15 +d_ang= ang / (iter);
17 16
18 % ******** Testing parameters ****************************************** 17 % ******** Testing parameters ******************************************
19 -p = [50, 20, 50] ;  
20 -ps = [100, 100, 50] - 1;  
21 -  
22 -volfile = 'img\nissl-rat.vol';  
23 -fid = fopen(volfile); % open the file that include the image  
24 -S = fread(fid, 3, 'int32');  
25 -X = S(1);  
26 -Y = S(2);  
27 -Z = S(3);  
28 -  
29 -% load the VOL data into a 2D matrix  
30 -I = fread(fid,[X Y*Z], 'uint8');  
31 -fclose(fid);  
32 -  
33 -%change this to a 3D matrix  
34 -I = reshape(I, [X, Y, Z]);  
35 -  
36 -%invert the intensity  
37 -I = 255 - I; 18 +% p = [100, 50, 100];
  19 +% ps = [200, 200, 100];
  20 +ps = [150, 100, 25];
  21 +r = 7;
  22 +I = syn_Img(r , ps);
  23 +% volfile = 'img\nissl-rat.vol';
  24 +% fid = fopen(volfile); % open the file that include the image
  25 +% S = fread(fid, 3, 'int32');
  26 +% X = S(1);
  27 +% Y = S(2);
  28 +% Z = S(3);
  29 +%
  30 +% % load the VOL data into a 2D matrix
  31 +% I = fread(fid,[X Y*Z], 'uint8');
  32 +% fclose(fid);
  33 +%
  34 +% %change this to a 3D matrix
  35 +% I = (reshape(I, [X, Y, Z]));
  36 +%
  37 +% % invert the intensity
  38 +% I = 255 - I;
38 39
39 %perform a gaussian blur 40 %perform a gaussian blur
40 -I_blur = gauss_blur3d(I, sigma); 41 +Iblur = gauss_blur3d(I, sigma);
41 42
42 -% compute the gradient  
43 -[Igrad_x, Igrad_y, Igrad_z] = gradient(I_blur); 43 +%%crop out a small subregion of I and Iblur
  44 +% Iblur = Iblur(p(1):p(1)+ps(1)-1, p(2):p(2)+ps(2)-1, p(3):p(3)+ps(3)-1);
  45 +% I = I(p(1):p(1)+ps(1)-1, p(2):p(2)+ps(2)-1, p(3):p(3)+ps(3)-1);
44 46
45 -%%crop out a small subregion of I  
46 -Isub_x = Igrad_x(p(1):p(1)+ps(1), p(2):p(2)+ps(2), p(3):p(3)+ps(3));  
47 -Isub_y = Igrad_y(p(1):p(1)+ps(1), p(2):p(2)+ps(2), p(3):p(3)+ps(3));  
48 -Isub_z = Igrad_z(p(1):p(1)+ps(1), p(2):p(2)+ps(2), p(3):p(3)+ps(3)); 47 +%
  48 +% compute the gradient
  49 +[Igrad_y, Igrad_x, Igrad_z] = gradient(Iblur);
49 50
50 %calculate the gradient magnitude 51 %calculate the gradient magnitude
51 -Imag = sqrt(Isub_x .^ 2 + Isub_y .^ 2 + Isub_z .^2); 52 +Imag = sqrt(Igrad_x .^ 2 + Igrad_y .^ 2 + Igrad_z .^2);
  53 +I_theta1 = acos(Igrad_x./Imag);
  54 +I_theta2 = acos(Igrad_y./Imag);
  55 +I_theta3 = acos(Igrad_z./Imag);
52 56
53 %set a threshold for the gradient magnitude 57 %set a threshold for the gradient magnitude
54 It = Imag > t0; 58 It = Imag > t0;
55 59
  60 +%Set the boundaries of the threshold image to zero
  61 +It(1:rmax, :, :) = 0;
  62 +It(ps(1) - rmax:ps(1), :,:) = 0;
  63 +It(:, 1:rmax, :) = 0;
  64 +It(:, ps(2) - rmax:ps(2),:) = 0;
  65 +It(:, :, 1:rmax) = 0;
  66 +It(:,:, ps(3) - rmax:ps(3)) = 0;
  67 +
  68 +%get the indices of all of the nonzero values in the threshold image
  69 +% (voter positions)
  70 +[Itx,Ity,Itz] = ind2sub(size(It),find(It));
  71 +Vi = find(It);
  72 +nV = nnz(It);
  73 +
  74 +% create a meshgrid describing coordinates relative to the voter position
  75 +range = -rmax:rmax; %create an array of values between -rmax and rmax
  76 +[mx, my, mz] = meshgrid(range, range, range); %create a template describing local pixel position in a small cube
  77 +m_mag = sqrt(mx.^2 + my.^2 + mz.^2); %create a template describing the distance from the center of a small cube
56 78
  79 +% create a mask for the voting area
  80 +M_dist = m_mag <= rmax; %mask for the voting area distance (all values < rmax from the center)
  81 +
  82 +%calculate the direction vector between a pixel and voter
  83 +LV_x = mx./m_mag;
  84 +LV_y = my./m_mag;
  85 +LV_z = mz./m_mag;
  86 +
  87 +%number of pixels in the voting area of each voter (initialize to zero)
  88 +validPixels = (zeros(nV,1));
  89 +
  90 +%indices of pixels in the voting area of each voter
  91 +% indices reference the 3D image
  92 +g_v_prime = (zeros(nV, (rmax^3)));
  93 +
  94 +%--------------Display current data-------------------------
57 subplot(3, 3, 1), 95 subplot(3, 3, 1),
58 -imshow( squeeze( I(:, :, ceil(size(I, 3)/2)) ) ./ 255 );  
59 -%colormap(gray); 96 +imagesc(squeeze(I(:, :, ceil(size(I, 3)/2))));
60 97
61 subplot(3, 3, 2), 98 subplot(3, 3, 2),
62 -imshow( squeeze( I(:, ceil(size(I, 2)/2), :) ) ./ 255 ); 99 +imagesc(squeeze(I(:, ceil(size(I, 2)/2), :)));
63 100
64 subplot(3, 3, 3), 101 subplot(3, 3, 3),
65 -imshow( squeeze( I(ceil(size(I, 1)/2), :, :) ) ./ 255 );  
66 - 102 +imagesc(squeeze(I( ceil(size(I, 1)/2), :, :)));
  103 +%
67 subplot(3, 3, 4), 104 subplot(3, 3, 4),
68 -imagesc( squeeze( Imag(:, :, ceil(size(Imag, 3)/2)) ) ); 105 +imagesc(squeeze(Iblur(:, :, ceil(size(Iblur, 3)/2))));
69 106
70 subplot(3, 3, 5), 107 subplot(3, 3, 5),
71 -imagesc( squeeze( Imag(:, ceil(size(Imag, 2)/2), :) ) ); 108 +imagesc(squeeze(Iblur(:, ceil(size(Iblur, 2)/2), :)));
72 109
73 subplot(3, 3, 6), 110 subplot(3, 3, 6),
74 -imagesc( squeeze( Imag(ceil(size(Imag, 1)/2), :, :) ) ); 111 +imagesc(squeeze(Iblur(ceil(size(Iblur, 1)/2), :, :)));
75 112
76 subplot(3, 3, 7), 113 subplot(3, 3, 7),
77 -imagesc( squeeze( It(:, :, ceil(size(It, 3)/2)) ) ); 114 +imagesc(squeeze(Imag(:, :, ceil(size(Imag, 3)/2))));
78 115
79 subplot(3, 3, 8), 116 subplot(3, 3, 8),
80 -imagesc( squeeze( It(:, ceil(size(It, 2)/2), :) ) ); 117 +imagesc(squeeze(Imag(:, ceil(size(Imag, 2)/2), :)));
81 118
82 subplot(3, 3, 9), 119 subplot(3, 3, 9),
83 -imagesc( squeeze( It(ceil(size(It, 1)/2), :, :) ) ); 120 +imagesc(squeeze(Imag(ceil(size(Imag, 1)/2), :, :)));
84 121
85 -%  
86 -% thresholding  
87 -It_x = size(It,1);  
88 -It_y = size(It,2);  
89 -It_z = size(It,3);  
90 -%  
91 -It(1:rmax, :, :) = 0;  
92 -It(It_x - rmax:It_x, :,:) = 0;  
93 -S(:, 1:rmax, :) = 0;  
94 -S(:, It_y - rmax:It_y,:) = 0;  
95 -S(:, :, 1:rmax) = 0;  
96 -S(:,:, It_z - rmax:It_z) = 0; 122 +colormap(gray);
97 123
98 -%  
99 -[Itx,Ity,Itz] = ind2sub(size(It),find(It)); 124 +%% vote
  125 +tic;
  126 +%for each iteration (in iterative voting)
  127 +for itr = 1 : iter
100 128
101 -nV = nnz(It); 129 + %initialize the vote image to zero
  130 + Ivote = (zeros(size(I)));
  131 +
  132 + %for each voter (nonzero pixels in the threshold image It)
  133 + for v = 1: nV
  134 +
  135 + %get the cartesian coordinates of the voter v in the main image I
  136 + vx = Itx(v);
  137 + vy = Ity(v);
  138 + vz = Itz(v);
  139 + vi = Vi(v);
  140 +
  141 + %retreive the gradient magnitude at the voter position
  142 + %vmag = Imag(vx,vy,vz);
  143 + vmag = Imag(vi);
  144 +
  145 + %retrieve the gradient
  146 + gx = Igrad_x(vi);
  147 + gy = Igrad_y(vi);
  148 + gz = Igrad_z(vi);
  149 +
  150 + %calculate the gradient magnitude
  151 + dmag = sqrt (gx^2 + gy^2 + gz^2);
  152 +
  153 + %calculate the normalized gradient direction
  154 + dx = gx / dmag;
  155 + dy = gy / dmag;
  156 + dz = gz / dmag;
  157 +
  158 + %gdir_x = Igrad_x(vx, vy, vz) / dir_mag;
  159 + %gdir_y = Igrad_y(vx, vy, vz) / dir_mag;
  160 + %gdir_z = Igrad_z(vx, vy, vz) / dir_mag;
  161 +
  162 + %calculate the angle between the voter direction and the pixel direction
  163 + cos_diff = LV_x .* dx + LV_y .* dy + LV_z .* dz;
  164 + ang_diff = acos(cos_diff);
  165 +
  166 + %create an angular mask for the voting area
  167 + M_angle = cos_diff > cos(ang);
  168 +
  169 + %combine the two masks to mask out the voting angle
  170 + M = M_angle .* M_dist;
  171 +
  172 + % get the coordinates of each pixel in the final voter mask M
  173 + pi = find(M);
  174 +
  175 + %calculate the number of pixels in the voting region
  176 + npts = nnz(M);
  177 + validPixels(v) = npts;
  178 +
  179 + %convert every index in the voting area from a local 3D index to a global 3D index (into the original image I)
  180 + global_px = vx + mx(pi);
  181 + global_py = vy + my(pi);
  182 + global_pz = vz + mz(pi);
  183 +
  184 + %convert the global 3D index of each point into a global 1D index
  185 + %global_pi = sub2ind(ps, global_px, global_py, global_pz);
  186 + %global_pi = (global_pz-1)*ps(1)*ps(2) + (global_py-1)*ps(1) + global_px;
  187 + global_pi = sub2ind(ps, global_px, global_py, global_pz);
  188 +
  189 + g_v_prime (v, 1:npts) = global_pi;
  190 +
  191 +
  192 + Ivote( global_pi ) = Ivote( global_pi ) + vmag;
  193 +
  194 + end
  195 +
  196 + if itr ==1
  197 + Ivote1 = Ivote;
  198 +
  199 + elseif itr ==2
  200 + Ivote2 = Ivote;
  201 +
  202 + elseif itr ==3
  203 + Ivote3 = Ivote;
  204 +
  205 + elseif itr ==4
  206 + Ivote4 = Ivote;
  207 +
  208 + elseif itr == 5
  209 + Ivote5 = Ivote;
  210 + end
  211 + t_v1 = toc;
  212 + disp(['voting done. time =',num2str(t_v1)]);
  213 +
  214 + figure, imagesc(squeeze(Ivote(:, :, ceil(size(Ivote, 3)/2),:))); colormap(gray);
  215 +
  216 + % update the voting direction
  217 + tic;
  218 + for v = 1: nV
  219 + % coordinates of the current voter
  220 + vx = Itx(v);
  221 + vy = Ity(v);
  222 + vz = Itz(v);
  223 +
  224 + %get the local value of the voting image
  225 + local_Ivote = Ivote(g_v_prime(v,1:validPixels(v)));
  226 +
  227 + %find the index of the maximum value
  228 + [~, local_max_idx] = max(local_Ivote);
  229 +
  230 + %convert this into a global subscript
  231 + [gx, gy, gz] = ind2sub(size(Ivote), g_v_prime(v,local_max_idx));
  232 +
  233 + %compute the vector from the voter position to this position
  234 + dx = gx - vx;
  235 + dy = gy - vy;
  236 + dz = gz - vz;
  237 + Igrad_x(vx, vy, vz) = dx;
  238 + Igrad_y(vx, vy, vz) = dy;
  239 + Igrad_z(vx, vy, vz) = dz;
  240 +
  241 + end
  242 +
  243 + tdir1 = toc;
  244 + display (['updating dir done. time = ', num2str(tdir1)]);
  245 + ang = ang - d_ang;
  246 + end
  247 +
102 248
  249 +%%
  250 +t = 10;
  251 +out = Ivote;
  252 +out(out<t) = 0;
  253 +% out(out>=t) = 255;
  254 +out = imregionalmax(out);
103 255
104 -% create a meshgrid describing coordinates relative to the voter position  
105 -range = -rmax:rmax;  
106 -[mx, my, mz] = meshgrid(range, range, range);  
107 -m_mag = sqrt(mx.^2 + my.^2 + mz.^2);  
108 -%create a mask for the voting area  
109 -M_dist = mx.^2 + my.^2 + mz.^2 < rmax^2; 256 +out1(:,:,:,1) = mat2gray(Iblur);
  257 +out1(:,:,:,2) = mat2gray(out);
  258 +out1(:,:,:,3) = mat2gray(Iblur);
110 259
111 -pxPerRow = size(I_blur,1);  
112 -pxPerCol = size(I_blur,2);  
113 -validPoints = zeros(nV,1);  
114 260
115 -g_v_prime = zeros(nV, ceil(rmax^2*phi/3));  
116 261
117 -%% 262 +figure(2);
118 263
119 -disp('done.');  
120 -%% vote  
121 -for itr = 1 : iter  
122 - Ivote = zeros(size(Iblur));  
123 - for v = 1: nV  
124 - % coordinates of the current voter  
125 - vx = Itx(v);  
126 - vy = Ity(v);  
127 - vz = Itz(v);  
128 -  
129 - vgrad = [Isub_x(vx,vy,vz), Isub_y(vx,vy,vz), Isub_z(vx,vy,vz)];  
130 - %determine the gradient magnitude at the current voter location  
131 - vmag = Imag(vx,vy,vz);  
132 - ang_diff = zeros(2*rmax+1, 2*rmax+1, 2*rmax+1);  
133 - for i = 1 : 2*rmax+1  
134 - for j = 1 : 2*rmax+1  
135 - for k = 1 : 2*rmax+1  
136 - a = [mx(i,j,k), my(i,j,k), mz(i,j,k)];  
137 - c = dot(vgrad, a);  
138 - ang_diff(i,j,k) = acos(c/(m_mag(i,j,k)* vmag ));  
139 - end  
140 - end  
141 - end  
142 -  
143 - M_diff = min(2*pi - ang_diff,ang_diff)<phi;  
144 - %compute the final mask  
145 - M = (M_dist .* M_diff);  
146 -  
147 - % get the coordinates of each pixel in the final mask  
148 - [vx_prime,vy_prime,vz_prime] = ind2sub(size(M),find(M));  
149 -  
150 - %transform the local coordinates of the pixels in the voting region  
151 - %to the global coordinates of the image  
152 -  
153 - npts = numel(vx_prime);  
154 - validPoints(v) = npts;  
155 -  
156 - g_v_prime(v,1:npts) = vx + (vx_prime - (rmax + 1)) + (vy + (vy_prime - (rmax + 1))-1).*pxPerRow + (vz + (vz_prime - (rmax + 1))-1).*pxPerRow*pxPerCol;  
157 -  
158 - for n=1: npts  
159 - %  
160 - Ivote( g_v_prime(v,n)) = Ivote( g_v_prime(v,n)) + vmag;  
161 - end  
162 -  
163 - end  
164 -  
165 -disp('voting done.');  
166 - %% update the voting direction  
167 -  
168 - for v = 1: nV  
169 - % coordinates of the current voter  
170 - vx = sx(v);  
171 - vy = sy(v);  
172 - vz = sz(v);  
173 -  
174 - %get the local value of the voting image  
175 - local_Ivote = Ivote(g_v_prime(v,1:validPoints(v)));  
176 -  
177 - %find the index of the maximum value  
178 - [~, local_max_idx] = max(local_Ivote);  
179 -  
180 - %convert this into a global subscript  
181 - gz = ceil(g_v_prime(v,local_max_idx)/(pxPerRow*pxPerCol)) ;  
182 - gy = ceil((g_v_prime(v,local_max_idx)-(gz - 1)*pxPerRow*pxPerCol)/pxPerRow);  
183 - gx = g_v_prime(v,local_max_idx)- pxPerRow*((gy - 1) + (gz - 1) * pxPerCol);  
184 -  
185 - %compute the vector from the voter position to this position  
186 - new_vx = gx - vx;  
187 - new_vy = gy - vy;  
188 - new_vz = gz - vz;  
189 - Igrad_x(vx,vy,vz) = new_vx;  
190 - Igrad_y(vx,vy,vz) = new_vy;  
191 - Igrad_z(vx,vy,vz) = new_vz;  
192 - end  
193 -end  
194 -cell_center = Ivote;  
195 -fido = fopen('out\outImg.vol', 'w+');  
196 -fwrite(fido, cell_center);  
197 -fclose(fido);  
198 -t1 = toc 264 +subplot(2, 3, 1),
  265 +imagesc(squeeze(Iblur(:, :, ceil(size(Iblur, 3)/2))));
  266 +
  267 +subplot(2, 3, 2),
  268 +imagesc(squeeze(Iblur(:, ceil(size(Iblur, 2)/2), :)));
  269 +
  270 +subplot(2, 3, 3),
  271 +imagesc(squeeze(Iblur( ceil(size(Iblur, 1)/2), :, :)));
  272 +
  273 +subplot(2, 3, 4),
  274 +imagesc( squeeze( out1(:, :, ceil(size(Ivote, 3)/2),:) ) );
  275 +
  276 +subplot(2, 3, 5),
  277 +imagesc(squeeze(out1(:, ceil(size(Ivote, 2)/2), :, :)));
  278 +
  279 +subplot(2, 3, 6),
  280 +imagesc(squeeze(out1( ceil(size(Ivote, 1)/2), :, :, :)));
  281 +
  282 +colormap(gray);
  283 +%%
  284 +% figure, imagesc(squeeze(Iblur(:, :, ceil(size(Iblur, 3)/2)))); colormap(gray);
  285 +% figure, imagesc(squeeze(Ivote(:, :, ceil(size(Ivote, 3)/2),:))); colormap(gray);
  286 +% figure, imagesc(squeeze(Iblur(:, ceil(size(Iblur, 2)/2), :))); colormap(gray);
  287 +% figure, imagesc(squeeze(I(:, ceil(size(Ivote, 2)/2), :))); colormap(gray);
  288 +% % figure, imagesc(squeeze(Iblur( ceil(size(Iblur, 1)/2), :, :))); colormap(gray);
  289 +% figure, imagesc(squeeze(Ivote( ceil(size(Ivote, 1)/2), :, :))); colormap(gray);
  290 +%
  291 +fid1 = fopen('ivote1.vol', 'w');
  292 +fwrite(fid1, Ivote1);
  293 +fclose(fid1);
  294 +
  295 +% fid0 = fopen('iblur.vol', 'w');
  296 +% fwrite(fid0, Iblur);
  297 +% fclose(fid0);
matlab/syn_Img.m 0 → 100644
  Show/Hide comments   View file @da849ea
  1 +
  2 +function I1 = syn_Img(r, ps)
  3 +
  4 +range = -r:r;
  5 +[cx, cy, cz] = meshgrid(range, range, range);
  6 +c1 = zeros(2*r+1,2*r+1,2*r+1);
  7 +c1(cx.^2 + cy.^2 + cz.^2 < r^2) = 255;
  8 +I1 = zeros(ps);
  9 +
  10 +sx = [8 36 91 45 50 76 85 67 51 21];
  11 +sy = [23 12 56 43 87 63 89 51 61 92];
  12 +sz = [8 10 9 14 17 10 11 15 13 10 11];
  13 +
  14 +for kk = 1:10
  15 +% set = [randi([r+1,ps(1)-r]) randi([r+1,ps(2)-r]) randi([r+1,ps(3)-r])];
  16 +% I1(set(1)-r:set(1)+r, set(2)-r:set(2)+r, set(3)-r:set(3)+r)= c1;
  17 + I1(sx(kk)-r:sx(kk)+r, sy(kk)-r:sy(kk)+r, sz(kk)-r:sz(kk)+r)= c1;
  18 +end
  19 +
  20 +% -------------------------------------------------------------------
  21 +
  22 +
  23 +