segmentation / ivote3

Commit c81508738f93816ed43adc8422c0ff05ab81147b

Authored by Laila Saadatifard
1 parent e4c75a97

fix the matlab ivote3 code to works on the synthetic image.

Showing 2 changed files with 38 additions and 265 deletions   Show diff stats
Inline Side-by-side
Matlab_3D/main.m
  Show/Hide comments   View file @c815087
... ... @@ -6,11 +6,11 @@ total = tic;
6 6  
7 7 % ******* Initialize voting parameters **************************************
8 8 rmax = 7; %maximum radius of the cell
9   -ang_deg = 15.1; %half the angular range of the voting area
  9 +ang_deg = 25.1; %half the angular range of the voting area
10 10 ang = ang_deg * pi / 180;
11   -iter = 3; %number of voting iterations
  11 +iter = 5; %number of voting iterations
12 12 t0 = 1.2; %threshold color
13   -sigma = [4, 4, 2];
  13 +sigma = [1, 1, 1];
14 14 % t = 0.1;
15 15 d_ang= ang / (iter);
16 16  
... ... @@ -18,8 +18,7 @@ d_ang= ang / (iter);
18 18 % p = [100, 50, 100];
19 19 % ps = [200, 200, 100];
20 20 ps = [100, 50, 40];
21   -r = 7;
22   -I = syn_Img(r , ps);
  21 +I = syn_Img(rmax , ps);
23 22 % volfile = 'img\nissl-rat.vol';
24 23 % fid = fopen(volfile); % open the file that include the image
25 24 % S = fread(fid, 3, 'int32');
... ... @@ -205,8 +204,6 @@ for itr = 1 : iter
205 204 t_v1 = toc;
206 205 disp(['voting done. time =',num2str(t_v1)]);
207 206  
208   - figure, imagesc(squeeze(Ivote(:, :, ceil(size(Ivote, 3)/2),:))); colormap(gray);
209   -
210 207 % update the voting direction
211 208 tic;
212 209 for v = 1: nV
... ... @@ -222,15 +219,12 @@ for itr = 1 : iter
222 219 [~, local_max_idx] = max(local_Ivote);
223 220  
224 221 %convert this into a global subscript
225   - [gx, gy, gz] = ind2sub(size(Ivote), g_v_prime(v,local_max_idx));
  222 + [g_px, g_py, g_pz] = ind2sub(size(Ivote), g_v_prime(v,local_max_idx));
226 223  
227 224 %compute the vector from the voter position to this position
228   - dx = gx - vx;
229   - dy = gy - vy;
230   - dz = gz - vz;
231   - Igrad_x(vx, vy, vz) = dx;
232   - Igrad_y(vx, vy, vz) = dy;
233   - Igrad_z(vx, vy, vz) = dz;
  225 + Igrad_x(vx, vy, vz) = g_px - vx;
  226 + Igrad_y(vx, vy, vz) = g_py - vy;
  227 + Igrad_z(vx, vy, vz) = g_pz - vz;
234 228  
235 229 end
236 230  
... ... @@ -241,7 +235,7 @@ for itr = 1 : iter
241 235  
242 236  
243 237 %%
244   -t = 10;
  238 +t = 2500;
245 239 out = Ivote;
246 240 out(out<t) = 0;
247 241 % out(out>=t) = 255;
... ... @@ -250,10 +244,10 @@ out = imregionalmax(out);
250 244 out1(:,:,:,1) = mat2gray(Iblur);
251 245 out1(:,:,:,2) = mat2gray(out);
252 246 out1(:,:,:,3) = mat2gray(Iblur);
  247 +% out = out1;
253 248  
254 249  
255   -
256   -figure(2);
  250 +figure(4);
257 251  
258 252 subplot(2, 3, 1),
259 253 imagesc(squeeze(Iblur(:, :, ceil(size(Iblur, 3)/2))));
... ... @@ -282,19 +276,30 @@ colormap(gray);
282 276 % % figure, imagesc(squeeze(Iblur( ceil(size(Iblur, 1)/2), :, :))); colormap(gray);
283 277 % figure, imagesc(squeeze(Ivote( ceil(size(Ivote, 1)/2), :, :))); colormap(gray);
284 278 %
285   -fid0 = fopen('iblur.vol', 'w');
286   -fwrite(fid0, Iblur);
287   -fclose(fid0);
288   -
289   -fid1 = fopen('ivote1.vol', 'w');
290   -fwrite(fid1, Ivote1);
291   -fclose(fid1);
292   -
293   -fid2 = fopen('ivote2.vol', 'w');
294   -fwrite(fid2, Ivote2);
295   -fclose(fid2);
296   -
297   -fid3 = fopen('ivote3.vol', 'w');
298   -fwrite(fid3, Ivote3);
299   -fclose(fid3);
300   -
  279 +% fid0 = fopen('D:\ivote3_files\iblur.vol', 'w');
  280 +% fwrite(fid0, Iblur);
  281 +% fclose(fid0);
  282 +%
  283 +% fid1 = fopen('D:\ivote3_files\ivote1.vol', 'w');
  284 +% fwrite(fid1, Ivote1);
  285 +% fclose(fid1);
  286 +%
  287 +% fid2 = fopen('D:\ivote3_files\ivote2.vol', 'w');
  288 +% fwrite(fid2, Ivote2);
  289 +% fclose(fid2);
  290 +%
  291 +% fid3 = fopen('D:\ivote3_files\ivote3.vol', 'w');
  292 +% fwrite(fid3, Ivote3);
  293 +% fclose(fid3);
  294 +%
  295 +% fid4 = fopen('D:\ivote3_files\ivote4.vol', 'w');
  296 +% fwrite(fid4, Ivote4);
  297 +% fclose(fid4);
  298 +%
  299 +% fid5 = fopen('D:\ivote3_files\ivote5.vol', 'w');
  300 +% fwrite(fid5, Ivote5);
  301 +% fclose(fid5);
  302 +%
  303 +% fid10 = fopen('D:\ivote3_files\ivote10.vol', 'w');
  304 +% fwrite(fid10, Ivote);
  305 +% fclose(fid10);
... ...
Matlab_3D/old-matlab.m deleted
View file @e4c75a9
1   -clc;
2   -clear;
3   -disp('***************** NEW RUN *********************');
4   -total = tic;
5   -
6   -
7   -% ******* Initialize voting parameters **************************************
8   -rmax = 9; %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 = 1; %number of voting iterations
12   -t0 = 1.2; %threshold color
13   -sigma = [4, 4, 2];
14   -% t = 10;
15   -d_ang= ang / (iter);
16   -
17   -% ******** Testing parameters ******************************************
18   -p = [100, 50, 100] ;
19   -ps = [50, 50, 25] - 1;
20   -% p = [100, 100, 100] ;
21   -% ps = [300, 300, 150] - 1;
22   -
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 = single(reshape(I, [X, Y, Z]));
36   -
37   -% invert the intensity
38   -I = 255 - I;
39   -
40   -%perform a gaussian blur
41   -I_blur = gauss_blur3d(I, sigma);
42   -
43   -% compute the gradient
44   -[Igrad_x, Igrad_y, Igrad_z] = gradient(I_blur);
45   -
46   -%%crop out a small subregion of I
47   -Isub_x = Igrad_x(p(1):p(1)+ps(1), p(2):p(2)+ps(2), p(3):p(3)+ps(3));
48   -Isub_y = Igrad_y(p(1):p(1)+ps(1), p(2):p(2)+ps(2), p(3):p(3)+ps(3));
49   -Isub_z = Igrad_z(p(1):p(1)+ps(1), p(2):p(2)+ps(2), p(3):p(3)+ps(3));
50   -
51   -%calculate the gradient magnitude
52   -Imag = sqrt(Isub_x .^ 2 + Isub_y .^ 2 + Isub_z .^2);
53   -I_theta = acos(Isub_z./Imag);
54   -I_phi = atan(Isub_y./Isub_x);
55   -
56   -
57   -%set a threshold for the gradient magnitude
58   -It = Imag > t0;
59   -
60   -%
61   -% thresholding
62   -It_x = size(It,1);
63   -It_y = size(It,2);
64   -It_z = size(It,3);
65   -%
66   -It(1:rmax, :, :) = 0;
67   -It(It_x - rmax:It_x, :,:) = 0;
68   -It(:, 1:rmax, :) = 0;
69   -It(:, It_y - rmax:It_y,:) = 0;
70   -It(:, :, 1:rmax) = 0;
71   -It(:,:, It_z - rmax:It_z) = 0;
72   -
73   -%
74   -[Itx,Ity,Itz] = (ind2sub(size(It),find(It)));
75   -Itx = single(Itx);
76   -Ity = single(Ity);
77   -Itz = single(Itz);
78   -nV = nnz(It);
79   -
80   -
81   -% create a meshgrid describing coordinates relative to the voter position
82   -range = single(-rmax):single(rmax);
83   -[mx, my, mz] = (meshgrid(range, range, range));
84   -m_mag = sqrt(mx.^2 + my.^2 + mz.^2);
85   -
86   -% create a mask for the voting area
87   -M_dist = mx.^2 + my.^2 + mz.^2 < rmax^2;
88   -M_theta = acos(mz./m_mag);
89   -M_phi = atan(my./mx);
90   -M_phi (my ==0 & mx ==0) = 0;
91   -pxPerRow = single(size(It,1));
92   -pxPerCol = single(size(It,2));
93   -validPoints = single(zeros(nV,1));
94   -
95   -% g_v_prime = zeros(nV, ceil(rmax^2*ang/3));
96   -g_v_prime = single(zeros(nV, (rmax^3)));
97   -
98   -
99   -
100   -disp('first part done.');
101   -%% vote
102   -tic;
103   - for itr = 1 : iter
104   -
105   - Ivote = single(zeros(size(It)));
106   - for v = 1: 1
107   -
108   -%
109   - % coordinates of the current voter
110   - vx = Itx(v);
111   - vy = Ity(v);
112   - vz = Itz(v);
113   -%
114   - vtheta= I_theta(vx,vy,vz);
115   - vphi = I_phi(vx,vy,vz);
116   - vmag = Imag(vx,vy,vz);
117   -% ang_diff = single(zeros([2*rmax+1 2*rmax+1 2*rmax+1]));
118   -% for i=1: 2*rmax +1
119   -% for j=1: 2*rmax+1
120   -% for k = 1:2*rmax+1
121   -% ang_diff(i,j,k) = acos(dot([Isub_x(vx,vy,vz) Isub_y(vx,vy,vz) Isub_z(vx,vy,vz)], [mx(i,j,k) my(i,j,k) mz(i,j,k)])/(vmag*norm([mx(i,j,k) my(i,j,k) mz(i,j,k)])));
122   -% end
123   -% end
124   -% end
125   - ang_diff = acos(((sin(M_theta).*sin(vtheta)).*cos(vphi - M_phi)) + cos(M_theta).*cos(vtheta));
126   - M_diff = min(2*pi - ang_diff,ang_diff)<ang;
127   -% M_diff = abs(ang_diff)<ang;
128   - %compute the final mask
129   - M = single(M_dist .* M_diff);
130   -%
131   - % get the coordinates of each pixel in the final mask
132   - [vx_prime,vy_prime,vz_prime] = ind2sub(size(M),find(M));
133   -
134   - %transform the local coordinates of the pixels in the voting region
135   - %to the global coordinates of the image
136   -
137   - npts = numel(vx_prime);
138   - validPoints(v) = npts;
139   -
140   - g_v_prime(v,1:npts) = vx + (single(vx_prime) - single(rmax + 1)) + (vy + (single(vy_prime) - single(rmax + 1))-1).*single(pxPerRow) + (vz + (single(vz_prime) - single(rmax + 1))-1).*single(pxPerRow)*single(pxPerCol);
141   -
142   - for n=1: npts
143   - %
144   - Ivote( g_v_prime(v,n)) = Ivote( single(g_v_prime(v,n))) + vmag;
145   - end
146   -
147   - end
148   -
149   -
150   - t_v1 = toc;
151   - disp(['voting done. time =',num2str(t_v1)]);
152   -
153   -
154   - % update the voting direction
155   - tic;
156   - for v = 1: nV
157   - % coordinates of the current voter
158   - vx = single(Itx(v));
159   - vy = single(Ity(v));
160   - vz = single(Itz(v));
161   -
162   - %get the local value of the voting image
163   - local_Ivote = Ivote(g_v_prime(v,1:validPoints(v)));
164   -
165   - %find the index of the maximum value
166   - [~, local_max_idx] = max(local_Ivote);
167   -
168   - %convert this into a global subscript
169   - gz = ceil(g_v_prime(v,local_max_idx)/(pxPerRow*pxPerCol)) ;
170   - gy = ceil((g_v_prime(v,local_max_idx)-(gz - 1)*pxPerRow*pxPerCol)/pxPerRow);
171   - gx = g_v_prime(v,local_max_idx)- pxPerRow*((gy - 1) + (gz - 1) * pxPerCol);
172   -
173   - %compute the vector from the voter position to this position
174   - up_vx = gx - vx;
175   - up_vy = gy - vy;
176   - up_vz = gz - vz;
177   - I_theta(vx,vy,vz) = acos(up_vz/sqrt(up_vx^2 + up_vy^2 + up_vz^2));
178   - I_phi(vx,vy,vz) = atan(up_vy/up_vx);
179   - end
180   -
181   - tdir1 = toc;
182   - display (['updating dir done. time = ', num2str(tdir1)]);
183   - ang = ang - d_ang;
184   - end
185   -
186   -
187   -%%
188   -out1(:,:,:,1) = mat2gray( I(p(1):p(1)+ps(1), p(2):p(2)+ps(2), p(3):p(3)+ps(3)));
189   -out1(:,:,:,2) = mat2gray(Ivote);
190   -out1(:,:,:,3) = mat2gray(I(p(1):p(1)+ps(1), p(2):p(2)+ps(2), p(3):p(3)+ps(3)));
191   -
192   -subplot(3, 3, 1),
193   -imagesc( squeeze( out1(:, :, ceil(size(out1, 3)/2),:) ) );
194   -
195   -subplot(3, 3, 2),
196   -imagesc( squeeze( out1(:, ceil(size(out1, 2)/2), :,:) ) );
197   -
198   -subplot(3, 3, 3),
199   -imagesc( squeeze( out1(ceil(size(out1, 1)/2), :, :,:) ) );
200   -
201   -t = 10;
202   -cell_center = Ivote;
203   -cell_center(cell_center<t) =0;
204   -cell_center(cell_center>=t) = 255;
205   -% cell_center = imregionalmax(cell_center);
206   -
207   -out2(:,:,:,1) = mat2gray( I(p(1):p(1)+ps(1), p(2):p(2)+ps(2), p(3):p(3)+ps(3)));
208   -out2(:,:,:,2) = mat2gray(cell_center);
209   -out2(:,:,:,3) = mat2gray(I(p(1):p(1)+ps(1), p(2):p(2)+ps(2), p(3):p(3)+ps(3)));
210   -
211   -
212   -subplot(3, 3, 4),
213   -imagesc( squeeze( out2(:, :, ceil(size(out2, 3)/2),:) ) );
214   -
215   -subplot(3, 3, 5),
216   -imagesc( squeeze( out2(:, ceil(size(out2, 2)/2), :,:) ) );
217   -
218   -subplot(3, 3, 6),
219   -imagesc( squeeze( out2(ceil(size(out2, 1)/2), :, :,:) ) );
220   -
221   -colormap(gray);
222   -subplot(3, 3, 7),
223   -imagesc( squeeze( Imag(:, :, ceil(size(Imag, 3)/2),:) ) );
224   -
225   -subplot(3, 3, 8),
226   -imagesc( squeeze( Imag(:, ceil(size(Imag, 2)/2), :,:) ) );
227   -
228   -subplot(3, 3, 9),
229   -imagesc( squeeze( Imag(ceil(size(Imag, 1)/2), :, :,:) ) );
230   -
231   -
232   -