Commit 46820b25c7c04467a21376027e76ceb732285a30

Authored by Laila Saadatifard
1 parent a60dc7ad

create an obj file as the output of Ivote3

Matlab_3D/0-gt.vol deleted
No preview for this file type
Matlab_3D/ivote3.m
@@ -4,11 +4,11 @@ disp('***************** NEW RUN *********************'); @@ -4,11 +4,11 @@ disp('***************** NEW RUN *********************');
4 total = tic; 4 total = tic;
5 % ******* Initialize voting parameters ************************************** 5 % ******* Initialize voting parameters **************************************
6 rmax = [10 10 10]; %maximum radius of the cell 6 rmax = [10 10 10]; %maximum radius of the cell
7 -rmin = [1 1 1]; 7 +%rmin = [1 1 1];
8 ang_deg = 25.1; %half the angular range of the voting area 8 ang_deg = 25.1; %half the angular range of the voting area
9 ang = ang_deg * pi / 180; 9 ang = ang_deg * pi / 180;
10 iter = 8 ; %number of voting iterations 10 iter = 8 ; %number of voting iterations
11 -t0 = 1; 11 +t0 = 0;
12 sigma = [5, 5, 5]; 12 sigma = [5, 5, 5];
13 % t = 0.1; 13 % t = 0.1;
14 d_ang= ang / (iter+2); 14 d_ang= ang / (iter+2);
@@ -23,7 +23,7 @@ d_ang= ang / (iter+2); @@ -23,7 +23,7 @@ d_ang= ang / (iter+2);
23 % X = S(1); 23 % X = S(1);
24 % Y = S(2); 24 % Y = S(2);
25 % Z = S(3); 25 % Z = S(3);
26 -filename = '128-128-128/nissl-float-128.128.128.vol'; %'nissl-float-128.128.128.vol'; 26 +filename = 'nissl-float-128.128.128.vol'; %'nissl-float-128.128.128.vol';
27 X = 128; 27 X = 128;
28 Y = 128; 28 Y = 128;
29 Z = 128; 29 Z = 128;
@@ -50,16 +50,16 @@ Isize = size(I); @@ -50,16 +50,16 @@ Isize = size(I);
50 50
51 %set a threshold for the gradient magnitude 51 %set a threshold for the gradient magnitude
52 It = Imag > t0; 52 It = Imag > t0;
53 -fidt = fopen('128-128-128/It.vol', 'w');  
54 -fwrite(fidt, It, 'single');  
55 -fclose(fidt); 53 +% fidt = fopen('128-128-128/It.vol', 'w');
  54 +% fwrite(fidt, It, 'single');
  55 +% fclose(fidt);
56 %Set the boundaries of the threshold image to zero 56 %Set the boundaries of the threshold image to zero
57 -It(1:rmax(1), :, :) = 0;  
58 -It(X - rmax(1):X, :,:) = 0;  
59 -It(:, 1:rmax(2), :) = 0;  
60 -It(:, Y - rmax(2):Y,:) = 0;  
61 -It(:, :, 1:rmax(3)) = 0;  
62 -It(:,:, Z - rmax(3):Z) = 0; 57 +% It(1:rmax(1), :, :) = 0;
  58 +% It(X - rmax(1):X, :,:) = 0;
  59 +% It(:, 1:rmax(2), :) = 0;
  60 +% It(:, Y - rmax(2):Y,:) = 0;
  61 +% It(:, :, 1:rmax(3)) = 0;
  62 +% It(:,:, Z - rmax(3):Z) = 0;
63 63
64 %get the indices of all of the nonzero values in the threshold image 64 %get the indices of all of the nonzero values in the threshold image
65 % (voter positions) 65 % (voter positions)
@@ -74,9 +74,9 @@ rangez = -rmax(3):rmax(3); @@ -74,9 +74,9 @@ rangez = -rmax(3):rmax(3);
74 m_mag = (sqrt(mx.^2 + my.^2 + mz.^2)); %create a template describing the distance from the center of a small cube 74 m_mag = (sqrt(mx.^2 + my.^2 + mz.^2)); %create a template describing the distance from the center of a small cube
75 75
76 % create a mask for the voting area 76 % create a mask for the voting area
77 -M_dist1 = (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)  
78 -M_dist2 = (mx.^2/rmin(1)^2 + my.^2/rmin(2)^2 + mz.^2/rmin(3)^2) >= 1 ;  
79 -M_dist = M_dist1 .* M_dist2; 77 +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)
  78 +%M_dist2 = (mx.^2/rmin(1)^2 + my.^2/rmin(2)^2 + mz.^2/rmin(3)^2) >= 1 ;
  79 +%M_dist = M_dist1 .* M_dist2;
80 % calculate the direction vector between a pixel and voter 80 % calculate the direction vector between a pixel and voter
81 LV_x = mx./m_mag; 81 LV_x = mx./m_mag;
82 LV_y = my./m_mag; 82 LV_y = my./m_mag;
@@ -91,8 +91,8 @@ g_v_prime = zeros(nV, ceil(rmax(1)*rmax(2)*rmax(3)*ang)); @@ -91,8 +91,8 @@ g_v_prime = zeros(nV, ceil(rmax(1)*rmax(2)*rmax(3)*ang));
91 91
92 % vote 92 % vote
93 tic; 93 tic;
94 -mask = zeros(Isize);  
95 -mask1 = zeros(Isize); 94 +% mask = zeros(Isize);
  95 +% mask1 = zeros(Isize);
96 96
97 %for each iteration (in iterative voting) 97 %for each iteration (in iterative voting)
98 for itr = 1 :iter 98 for itr = 1 :iter
@@ -146,19 +146,29 @@ for itr = 1 :iter @@ -146,19 +146,29 @@ for itr = 1 :iter
146 global_pz = vz + mz(pi); 146 global_pz = vz + mz(pi);
147 147
148 %convert the global 3D index of each point into a global 1D index 148 %convert the global 3D index of each point into a global 1D index
149 - global_pi = sub2ind(Isize, global_px, global_py, global_pz); 149 + global_pi0 = zeros(npts);
  150 + idxx=0;
  151 + for id_gl=1:npts
  152 + 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)
  153 + idxx = idxx + 1;
  154 + global_pi0(idxx) = sub2ind(Isize, global_px(id_gl), global_py(id_gl), global_pz(id_gl));
  155 + end
  156 + end
  157 + global_pi = global_pi0(global_pi0~=0);
  158 + npts = numel(global_pi);
  159 +% global_pi = sub2ind(Isize, global_px, global_py, global_pz);
150 160
151 g_v_prime (v, 1:npts) = global_pi; 161 g_v_prime (v, 1:npts) = global_pi;
152 162
153 163
154 Ivote( global_pi ) = Ivote( global_pi ) + vmag; 164 Ivote( global_pi ) = Ivote( global_pi ) + vmag;
155 end 165 end
156 - fid = fopen(sprintf('128-128-128/nissl-vote%d',itr), 'w');  
157 - fwrite(fid, single(Ivote), '*single');  
158 - fclose(fid); 166 +% fid = fopen(sprintf('128-128-128/nissl-vote%d',itr), 'w');
  167 +% fwrite(fid, single(Ivote), '*single');
  168 +% fclose(fid);
159 169
160 - t_v1 = toc;  
161 - disp(['voting done. time =',num2str(t_v1)]); 170 + t_v1 = toc
  171 +% disp(['voting done. time =',num2str(t_v1)]);
162 172
163 % update the voting direction 173 % update the voting direction
164 if ang>=d_ang 174 if ang>=d_ang
@@ -189,22 +199,23 @@ for itr = 1 :iter @@ -189,22 +199,23 @@ for itr = 1 :iter
189 end 199 end
190 200
191 201
192 - tdir1 = toc;  
193 - display (['updating dir done. time = ', num2str(tdir1)]); 202 + tdir1 = toc
  203 +% display (['updating dir done. time = ', num2str(tdir1)]);
194 ang = ang - d_ang; 204 ang = ang - d_ang;
195 end 205 end
196 206
197 end 207 end
198 208
199 -hv = reshape(Ivote, [X*Y*Z, 1]);  
200 -hist(hv, 250); 209 +% hv = reshape(Ivote, [X*Y*Z, 1]);
  210 +% hist(hv, 250);
201 %% 211 %%
202 -t = 300; 212 +t = 2000;
203 conn = [5 5 5]; 213 conn = [5 5 5];
204 Icenter = local_max(Ivote, conn, t); 214 Icenter = local_max(Ivote, conn, t);
205 -fidc = fopen(sprintf('std3.2-r10.10-v8/out%d-t%d.vol',t,t0), 'w'); 215 +fidc = fopen(sprintf('out%d.vol',t), 'w');
206 fwrite(fidc, single(Icenter), '*single'); 216 fwrite(fidc, single(Icenter), '*single');
207 fclose(fidc); 217 fclose(fidc);
  218 +tot = toc(total)
208 nnz(Icenter) 219 nnz(Icenter)
209 % [cxx1, cyy1, czz1] = ind2sub(size(Icenter),find(Icenter)); 220 % [cxx1, cyy1, czz1] = ind2sub(size(Icenter),find(Icenter));
210 221
Matlab_3D/validation.m
1 -  
2 clear all; 1 clear all;
3 disp('***************** NEW RUN *********************'); 2 disp('***************** NEW RUN *********************');
4 X = 128; 3 X = 128;
@@ -6,21 +5,22 @@ Y = 128; @@ -6,21 +5,22 @@ Y = 128;
6 Z = 128; 5 Z = 128;
7 D = 10; 6 D = 10;
8 t0=1; 7 t0=1;
9 -r1=10; 8 +r1=12;
10 r2=10; 9 r2=10;
11 -t=2100; 10 +t=2300;
12 itr=8; 11 itr=8;
13 vote=10; 12 vote=10;
14 std = [5 5]; 13 std = [5 5];
15 -gt_filename = '0-gt.vol'; 14 +gt_filename = 'gt2.vol';
16 % out_filename = sprintf('128-128-128/0-nissl-std%d.%d-t0%d-r%d.%d-t%d-out%d.%d.vol',std(1), std(2),t0,r1,r2,t,itr,vote); 15 % out_filename = sprintf('128-128-128/0-nissl-std%d.%d-t0%d-r%d.%d-t%d-out%d.%d.vol',std(1), std(2),t0,r1,r2,t,itr,vote);
17 -out_filename = sprintf('D:/build/ivote3-bld/shared2D-v8/out%d.vol',t); 16 +out_filename = sprintf('D:/build/ivote3-bld/shared2D-v8/centers.%d.vol',t);
18 % txt_filename = sprintf('128-128-128/0-validation-nissl-std%d.%d-r%d.%d-t%d-out%d.%d-D%d.txt',std(1), std(2),r1,r2,t,itr,vote,D); 17 % txt_filename = sprintf('128-128-128/0-validation-nissl-std%d.%d-r%d.%d-t%d-out%d.%d-D%d.txt',std(1), std(2),r1,r2,t,itr,vote,D);
19 -txt_filename = sprintf('D:/build/ivote3-bld/shared2D-v8/t%d-D%d.txt',t,D); 18 +txt_filename = sprintf('0-t%d--1.txt',t);
20 spec = sprintf('Nissl-std%d.%d-r%d.%d-t%d-out%d.%d',std(1), std(2),r1,r2,t,itr,vote); 19 spec = sprintf('Nissl-std%d.%d-r%d.%d-t%d-out%d.%d',std(1), std(2),r1,r2,t,itr,vote);
21 fid0 = fopen(gt_filename); 20 fid0 = fopen(gt_filename);
22 gt = fread(fid0,[X Y*Z], 'single'); 21 gt = fread(fid0,[X Y*Z], 'single');
23 fclose(fid0); 22 fclose(fid0);
  23 +% store x,y, and z directions of cell centers defined by the ground-trurth in an array
24 gt = reshape(gt, [X Y Z]); 24 gt = reshape(gt, [X Y Z]);
25 [gtx, gty, gtz] = ind2sub(size(gt),find(gt)); 25 [gtx, gty, gtz] = ind2sub(size(gt),find(gt));
26 ref = [gtx gty gtz]; 26 ref = [gtx gty gtz];
@@ -28,6 +28,7 @@ ref = [gtx gty gtz]; @@ -28,6 +28,7 @@ ref = [gtx gty gtz];
28 fid1 = fopen(out_filename); 28 fid1 = fopen(out_filename);
29 out = fread(fid1,[X Y*Z], 'single'); 29 out = fread(fid1,[X Y*Z], 'single');
30 fclose(fid1); 30 fclose(fid1);
  31 +% store x,y, and z directions of detected cells by the code in an array
31 out = reshape(out, [X Y Z]); 32 out = reshape(out, [X Y Z]);
32 [cx, cy, cz] = ind2sub(size(out),find(out)); 33 [cx, cy, cz] = ind2sub(size(out),find(out));
33 cent = [cx cy cz]; 34 cent = [cx cy cz];
@@ -35,67 +36,135 @@ cent = [cx cy cz]; @@ -35,67 +36,135 @@ cent = [cx cy cz];
35 [idx, dist] = knnsearch(cent,ref); 36 [idx, dist] = knnsearch(cent,ref);
36 % 2- find the same detected cells for different elements in the ground truth and assign that to which gt's element with the lowest distance. 37 % 2- find the same detected cells for different elements in the ground truth and assign that to which gt's element with the lowest distance.
37 ind_eq_idx = zeros(numel(idx), 10); 38 ind_eq_idx = zeros(numel(idx), 10);
38 -t=0; 39 +%set number of cells with the same detected cell center to zero.
  40 +noc=0;
39 for i=1:numel(idx) 41 for i=1:numel(idx)
  42 + % check the current element hasn't considered yet
40 s_ind = sum(ind_eq_idx==i); 43 s_ind = sum(ind_eq_idx==i);
41 if s_ind==0 44 if s_ind==0
  45 + % check how many of the elemnts assigned to the same detected cell
42 s = sum(idx==idx(i)); 46 s = sum(idx==idx(i));
43 if s>1 47 if s>1
44 - t = t+1;  
45 - id_ref(t,1) = i;  
46 - id_ref(t,2) = s; 48 + noc = noc+1;
  49 + %save the index and number of gt's elements with the same assigned detected cells.
  50 + id_ref(noc,1) = i;
  51 + id_ref(noc,2) = s;
47 ind1 = find(idx==idx(i)); 52 ind1 = find(idx==idx(i));
48 ind_eq_idx(i, 1:numel(ind1)) = ind1; 53 ind_eq_idx(i, 1:numel(ind1)) = ind1;
49 end 54 end
50 end 55 end
51 end 56 end
  57 +% determine those indices which hs assigned to only one detected cell
52 b = idx; 58 b = idx;
53 u_eq_idx = unique(ind_eq_idx(ind_eq_idx>0)); 59 u_eq_idx = unique(ind_eq_idx(ind_eq_idx>0));
54 b(u_eq_idx)=0; 60 b(u_eq_idx)=0;
  61 +% set the number of over sefgmented cells to zero.
55 oseg=0; 62 oseg=0;
56 -  
57 for k=1:size(id_ref,1) 63 for k=1:size(id_ref,1)
58 k1 = id_ref(k,1); 64 k1 = id_ref(k,1);
59 k2 = id_ref(k,2); 65 k2 = id_ref(k,2);
  66 + %find the minimum distance to the detected cell happened for which of the gt's element
60 l = ind_eq_idx(k1,1:k2); 67 l = ind_eq_idx(k1,1:k2);
61 [~, local_id_min] = min(dist(l)); 68 [~, local_id_min] = min(dist(l));
  69 + % set the element with minimum distance to the corresponding detected cell
62 b(l(local_id_min)) = idx(k1); 70 b(l(local_id_min)) = idx(k1);
  71 + %remove the proper element from the list of indices with same designated detected cells
63 u_eq_idx (u_eq_idx == l(local_id_min))=0; 72 u_eq_idx (u_eq_idx == l(local_id_min))=0;
  73 + % check if the indices remained in the list has distance less than the value is set for validation
64 ly = l (l~=l(local_id_min)); 74 ly = l (l~=l(local_id_min));
65 distl = dist(ly); 75 distl = dist(ly);
66 sl = sum(distl<D); 76 sl = sum(distl<D);
  77 + % if the distance is low enought, consider the corresponding cell as oversegmented
67 if sl>0 78 if sl>0
68 oseg = oseg+1; 79 oseg = oseg+1;
69 end 80 end
70 end 81 end
71 -  
72 -%  
73 -u_idx1 = u_eq_idx (u_eq_idx>0); 82 +%****
  83 +% check if the oversegmented detected cells could be assigned to the gt's element with no detected cell!
  84 +%find the indices of gt's element with no detected cell
  85 +u_idx1 = find(b==0);
  86 +%find the coordinates of gt's
74 ref1 = ref(u_idx1,:); 87 ref1 = ref(u_idx1,:);
  88 +% find the indices of detected cells which are assigned to a gt's element
75 ub = unique(b(b>0)); 89 ub = unique(b(b>0));
  90 +% set the number of detected cell which are not assigned to a gt's element to zero
76 mm=0; 91 mm=0;
77 -for jj=1:size(cent,1)  
78 - z = sum(ub(:)==jj); 92 +% find the indices of detected cell which are not assigned to a gt's element
  93 +for j=1:size(cent,1)
  94 + z = sum(ub(:)==j);
79 if z==0 95 if z==0
80 mm = mm+1; 96 mm = mm+1;
81 - cent_id1(mm) = jj; 97 + cent_id1(mm) = j;
82 end 98 end
83 end 99 end
84 - 100 +% find the coordinated of of detected cells
85 cent1 = cent(cent_id1,:); 101 cent1 = cent(cent_id1,:);
86 102
  103 +% check if there are osme low enough distances
87 [idx1, dist1] = knnsearch(cent1,ref1); 104 [idx1, dist1] = knnsearch(cent1,ref1);
88 -min1= min(dist1(:));  
89 -if min1 < D  
90 - miss1 = sum(dist1(:)<D); 105 +ind_eq_idx_rc = zeros(numel(idx1), 10);
  106 +%set number of cells with the same detected cell center to zero.
  107 +noc_rc=0;
  108 +for i_rc=1:numel(idx1)
  109 + % check the current element hasn't considered yet
  110 + s_ind_rc = sum(ind_eq_idx_rc==i_rc);
  111 + if s_ind_rc==0
  112 + % check how many of the elemnts assigned to the same detected cell
  113 + s_rc = sum(idx1==idx1(i_rc));
  114 + if s_rc>1
  115 + noc_rc = noc_rc+1;
  116 + %save the index and number of gt's elements with the same assigned detected cells.
  117 + id_ref_rc(noc_rc,1) = i_rc;
  118 + id_ref_rc(noc_rc,2) = s_rc;
  119 + ind1_rc = find(idx1==idx1(i_rc));
  120 + ind_eq_idx_rc(i_rc, 1:numel(ind1_rc)) = ind1_rc;
  121 + end
  122 + end
  123 +end
  124 +% determine those indices which hs assigned to only one detected cell
  125 +b_rc = idx1;
  126 +u_eq_idx_rc = unique(ind_eq_idx_rc(ind_eq_idx_rc>0));
  127 +b_rc(u_eq_idx_rc)=0;
  128 +% set the number of over sefgmented cells to zero.
  129 +oseg_rc=0;
  130 +for k_rc=1:size(id_ref_rc,1)
  131 + k1_rc = id_ref_rc(k_rc,1);
  132 + k2_rc = id_ref_rc(k_rc,2);
  133 + %find the minimum distance to the detected cell happened for which of the gt's element
  134 + l_rc = ind_eq_idx_rc(k1_rc,1:k2_rc);
  135 + [~, local_id_min_rc] = min(dist1(l_rc));
  136 + % set the element with minimum distance to the corresponding detected cell
  137 + b_rc(l_rc(local_id_min_rc)) = idx1(k1_rc);
  138 + %remove the proper element from the list of indices with same designated detected cells
  139 + u_eq_idx_rc (u_eq_idx_rc == l_rc(local_id_min_rc))=0;
  140 + % check if the indices remained in the list has distance less than the value is set for validation
  141 + ly_rc = l_rc (l_rc~=l_rc(local_id_min_rc));
  142 + distl_rc = dist1(ly_rc);
  143 + sl_rc = sum(distl_rc<D);
  144 + % if the distance is low enought, consider the corresponding cell as oversegmented
  145 + if sl_rc>0
  146 + oseg_rc = oseg_rc+1;
  147 + end
91 end 148 end
92 -% 149 +
  150 +%******
  151 +% b include the gt's element and its detected cells, for those element with no cell detection, b has zero value
93 b_ind = find(b==0); 152 b_ind = find(b==0);
94 b_dist = dist; 153 b_dist = dist;
  154 +% remove the disatances for those elements with no detected cells from distance array
95 b_dist(b_ind)=-1; 155 b_dist(b_ind)=-1;
96 -TP = sum(b_dist>=0) - sum(b_dist>=D); 156 +
  157 +b_ind_rc = find(b_rc==0);
  158 +b_dist_rc = dist1;
  159 +b_dist_rc(b_ind_rc)=-1;
  160 +% calculate Ttrue Positive, number of detected cells that have low enough distance to one of the gt's elements.
  161 +TP = sum(b_dist>=0) - sum(b_dist>D) + (sum(b_dist_rc>=0) - sum(b_dist_rc>D));
  162 +% calculate False Negative, number of gt's elements with no detected cells.
97 FN = size(ref, 1) - TP; 163 FN = size(ref, 1) - TP;
  164 +% calculate False Positive, number of detected cells, which their distance to the gt's element is long
  165 +% or considered as over segmentation.
98 FP = size(cent,1) - TP; 166 FP = size(cent,1) - TP;
  167 +% compute the validation factors.
99 Precision = TP/(TP + FP); 168 Precision = TP/(TP + FP);
100 Recall = TP/(TP+FN); 169 Recall = TP/(TP+FN);
101 Accuracy= TP/(TP+FP+FN); 170 Accuracy= TP/(TP+FP+FN);
1 -//#include "cuda_fp16.h"  
2 -//#include "float_to_half.cuh"  
3 -//#include "half_to_float.cuh"  
4 #include "gaussian_blur3.cuh" 1 #include "gaussian_blur3.cuh"
5 #include "gradient3.cuh" 2 #include "gradient3.cuh"
6 #include "mag3.cuh" 3 #include "mag3.cuh"
@@ -9,7 +6,7 @@ @@ -9,7 +6,7 @@
9 #include "local_max3.cuh" 6 #include "local_max3.cuh"
10 7
11 8
12 -void ivote3(float* center, float* img, float sigma[], float anisotropy, float phi, float d_phi, unsigned int r[], 9 +void ivote3(float* img, float sigma[], float anisotropy, float phi, float d_phi, unsigned int r[],
13 int iter, float t, unsigned int conn[], unsigned int x, unsigned int y, unsigned int z){ 10 int iter, float t, unsigned int conn[], unsigned int x, unsigned int y, unsigned int z){
14 11
15 12
@@ -47,8 +44,8 @@ void ivote3(float* center, float* img, float sigma[], float anisotropy, float ph @@ -47,8 +44,8 @@ void ivote3(float* center, float* img, float sigma[], float anisotropy, float ph
47 44
48 gpu_vote3<float>(gpu_vote, gpu_grad, cos_phi, r, x, y, z); 45 gpu_vote3<float>(gpu_vote, gpu_grad, cos_phi, r, x, y, z);
49 cudaDeviceSynchronize(); 46 cudaDeviceSynchronize();
50 - if (i==7)  
51 - cudaMemcpy(img, gpu_vote, bytes, cudaMemcpyDeviceToHost); 47 + /*if (i==7)
  48 + cudaMemcpy(img, gpu_vote, bytes, cudaMemcpyDeviceToHost);*/
52 49
53 if (phi >= d_phi){ 50 if (phi >= d_phi){
54 gpu_update_dir3<float>(gpu_grad, gpu_vote, cos_phi, r, x, y, z); 51 gpu_update_dir3<float>(gpu_grad, gpu_vote, cos_phi, r, x, y, z);
@@ -60,20 +57,43 @@ void ivote3(float* center, float* img, float sigma[], float anisotropy, float ph @@ -60,20 +57,43 @@ void ivote3(float* center, float* img, float sigma[], float anisotropy, float ph
60 } 57 }
61 58
62 cudaFree(gpu_grad); 59 cudaFree(gpu_grad);
63 - //cudaMemcpy(center, gpu_grad, bytes, cudaMemcpyDeviceToHost); 60 + cudaMemcpy(img, gpu_vote, bytes, cudaMemcpyDeviceToHost);
64 61
65 //allocate space on the gpu for the final detected cells. 62 //allocate space on the gpu for the final detected cells.
66 - float* gpu_output;  
67 - cudaMalloc(&gpu_output, bytes); 63 + //float* gpu_output;
  64 + //cudaMalloc(&gpu_output, bytes);
  65 +
  66 + ////call the local max function
  67 + //gpu_local_max3<float>(gpu_output, gpu_vote, t, conn, x, y, z);
  68 +
  69 + ////copy the final result to the cpu.
  70 + //cudaMemcpy(center, gpu_output, bytes, cudaMemcpyDeviceToHost);
  71 + //
  72 + //
  73 + cudaFree(gpu_vote);
  74 + //cudaFree(gpu_output);
  75 +
  76 +}
  77 +
  78 +void lmax(float* out, float* in, float t, unsigned int conn[], unsigned int x, unsigned int y, unsigned int z){
  79 + unsigned int bytes = x * y * z * sizeof(float);
  80 +
  81 + //assign memory on gpu for the input data.
  82 + float* gpuV;
  83 + cudaMalloc(&gpuV, bytes);
  84 +
  85 + //copy the image data to the GPU.
  86 + cudaMemcpy(gpuV, in, bytes, cudaMemcpyHostToDevice);
  87 +
  88 + float* gpuOut;
  89 + cudaMalloc(&gpuOut, bytes);
68 90
69 //call the local max function 91 //call the local max function
70 - gpu_local_max3<float>(gpu_output, gpu_vote, t, conn, x, y, z); 92 + gpu_local_max3<float>(gpuOut, gpuV, t, conn, x, y, z);
71 93
72 //copy the final result to the cpu. 94 //copy the final result to the cpu.
73 - cudaMemcpy(center, gpu_output, bytes, cudaMemcpyDeviceToHost);  
74 -  
75 -  
76 - cudaFree(gpu_vote);  
77 - cudaFree(gpu_output); 95 + cudaMemcpy(out, gpuOut, bytes, cudaMemcpyDeviceToHost);
78 96
  97 + cudaFree(gpuV);
  98 + cudaFree(gpuOut);
79 } 99 }
80 \ No newline at end of file 100 \ No newline at end of file
cpp/gaussian_blur3.cuh
@@ -4,9 +4,7 @@ @@ -4,9 +4,7 @@
4 #include <iostream> 4 #include <iostream>
5 #include <cuda.h> 5 #include <cuda.h>
6 #include <stim/cuda/cudatools.h> 6 #include <stim/cuda/cudatools.h>
7 -#include <stim/cuda/sharedmem.cuh>  
8 #include <stim/cuda/cudatools/error.h> 7 #include <stim/cuda/cudatools/error.h>
9 -#include "cuda_fp16.h"  
10 8
11 #define pi 3.14159 9 #define pi 3.14159
12 10
@@ -4,7 +4,7 @@ @@ -4,7 +4,7 @@
4 #include <iostream> 4 #include <iostream>
5 #include <cuda.h> 5 #include <cuda.h>
6 #include <stim/cuda/cudatools.h> 6 #include <stim/cuda/cudatools.h>
7 -#include <stim/cuda/sharedmem.cuh> 7 +//#include <stim/cuda/sharedmem.cuh>
8 #include <stim/cuda/cudatools/error.h> 8 #include <stim/cuda/cudatools/error.h>
9 9
10 template<typename T> 10 template<typename T>
1 #include <iostream> 1 #include <iostream>
  2 +#include <string>
2 #include <fstream> 3 #include <fstream>
3 #include <cuda_runtime.h> 4 #include <cuda_runtime.h>
4 #include <stim/math/vector.h> 5 #include <stim/math/vector.h>
5 #include <stim/parser/arguments.h> 6 #include <stim/parser/arguments.h>
6 #include <stim/parser/filename.h> 7 #include <stim/parser/filename.h>
7 -#include <stim/grids/image_stack.h>  
8 -#include <stim/grids/grid.h>  
9 #include <stim/visualization/colormap.h> 8 #include <stim/visualization/colormap.h>
10 #include <stim/image/image.h> 9 #include <stim/image/image.h>
11 #define pi 3.14159 10 #define pi 3.14159
12 11
13 12
14 -void ivote3(float* center, float* img, float std[], float anisotropy, float phi, float d_phi, unsigned int r[], int iter, float t, unsigned int conn[], 13 +void ivote3(float* img, float std[], float anisotropy, float phi, float d_phi, unsigned int r[], int iter, float t, unsigned int conn[],
15 unsigned int x, unsigned int y, unsigned int z); 14 unsigned int x, unsigned int y, unsigned int z);
  15 +void lmax(float* center, float* vote, float t1, unsigned int conn[], unsigned int x, unsigned int y, unsigned int z);
16 16
17 void invert_data(float* cpuI, unsigned int x, unsigned int y, unsigned int z){ 17 void invert_data(float* cpuI, unsigned int x, unsigned int y, unsigned int z){
18 for(int ix = 0; ix < x; ix++){ 18 for(int ix = 0; ix < x; ix++){
@@ -99,10 +99,10 @@ int main(int argc, char** argv){ @@ -99,10 +99,10 @@ int main(int argc, char** argv){
99 //set the anisotropy 99 //set the anisotropy
100 float anisotropy = args["anisotropy"].as_float(); 100 float anisotropy = args["anisotropy"].as_float();
101 unsigned int rmax = 10 ; 101 unsigned int rmax = 10 ;
102 - unsigned int r[3] = { rmax, rmax, rmax}; 102 + unsigned int r[3] = { 12, rmax, rmax};
103 float std = 5; 103 float std = 5;
104 float sigma[3] = { std, std, std}; 104 float sigma[3] = { std, std, std};
105 - unsigned int nlmax = 5; 105 + unsigned int nlmax = 1;
106 unsigned int conn[3] = { nlmax, nlmax, nlmax}; 106 unsigned int conn[3] = { nlmax, nlmax, nlmax};
107 float phi_deg = 25.0; 107 float phi_deg = 25.0;
108 float phi = phi_deg * pi /180; 108 float phi = phi_deg * pi /180;
@@ -123,7 +123,7 @@ int main(int argc, char** argv){ @@ -123,7 +123,7 @@ int main(int argc, char** argv){
123 invert_data(cpuI, x, y, z); 123 invert_data(cpuI, x, y, z);
124 124
125 //write a new file from the cpuI. 125 //write a new file from the cpuI.
126 - std::ofstream original("shared2D-v8/inv-128.vol", std::ofstream::out | std::ofstream::binary); 126 + std::ofstream original("shared2D-v1/inv-128.vol", std::ofstream::out | std::ofstream::binary);
127 original.write((char*)cpuI, bytes); 127 original.write((char*)cpuI, bytes);
128 original.close(); 128 original.close();
129 129
@@ -131,8 +131,8 @@ int main(int argc, char** argv){ @@ -131,8 +131,8 @@ int main(int argc, char** argv){
131 float* cpu_out = (float*) malloc(bytes); 131 float* cpu_out = (float*) malloc(bytes);
132 132
133 // call the ivote function 133 // call the ivote function
134 - ivote3(cpu_out, cpuI, sigma, anisotropy, phi, d_phi, r, iter, t, conn, x, y, z);  
135 - 134 + //ivote3(cpu_out, cpuI, sigma, anisotropy, phi, d_phi, r, iter, t, conn, x, y, z);
  135 + ivote3(cpuI, sigma, anisotropy, phi, d_phi, r, iter, t, conn, x, y, z);
136 //write the blurred file from the cpuI. 136 //write the blurred file from the cpuI.
137 std::ofstream fblur("shared2D-v8/vote8.vol", std::ofstream::out | std::ofstream::binary); 137 std::ofstream fblur("shared2D-v8/vote8.vol", std::ofstream::out | std::ofstream::binary);
138 fblur.write((char*)cpuI, bytes); 138 fblur.write((char*)cpuI, bytes);
@@ -145,27 +145,38 @@ int main(int argc, char** argv){ @@ -145,27 +145,38 @@ int main(int argc, char** argv){
145 fgx.close(); 145 fgx.close();
146 */ 146 */
147 //write the output file. 147 //write the output file.
148 - std::ofstream fo("shared2D-v8/" + OutName.str(), std::ofstream::out | std::ofstream::binary);  
149 - fo.write((char*)cpu_out, bytes);  
150 - fo.close(); 148 + //for (int t0=2000; t0<=2500; t0+=100){
  149 + int t0 = t;
  150 + lmax(cpu_out, cpuI, t, conn, x, y, z);
  151 + //std::ofstream fo("shared2D-v8/" + OutName.str(), std::ofstream::out | std::ofstream::binary);
  152 + std::ofstream fo("shared2D-v8/" + OutName.str()+std::to_string(t0)+".vol", std::ofstream::out | std::ofstream::binary);
  153 + fo.write((char*)cpu_out, bytes);
  154 + fo.close();
151 155
152 // creat a file for saving the list centers 156 // creat a file for saving the list centers
153 -  
154 - std::ofstream list("shared2D-v8/center.txt"); 157 +
  158 + std::ofstream list("shared2D-v8/" + OutName.str()+std::to_string(t0)+".obj");
  159 + // set the number of detected cells to zero.
  160 + int nod = 0;
155 if (list.is_open()){ 161 if (list.is_open()){
156 162
157 - for (int ix=0; ix<x; ix++){ 163 + for (int iz=0; iz<z; iz++){
158 for (int iy=0; iy<y; iy++){ 164 for (int iy=0; iy<y; iy++){
159 - for (int iz=0; iz<z; iz++){ 165 + for (int ix=0; ix<x; ix++){
160 166
161 int idx = iz * x * y + iy * x + ix; 167 int idx = iz * x * y + iy * x + ix;
162 if (cpu_out[idx]==1){ 168 if (cpu_out[idx]==1){
163 - list << ix << "\t" << iy << "\t"<< iz << '\n' ; 169 + nod++;
  170 + list << "v" << "\t" << ix << "\t" << iy << "\t"<< iz << '\n' ;
164 171
165 } 172 }
166 } 173 }
167 } 174 }
168 } 175 }
  176 + list << "p" << "\t";
  177 + for (unsigned int i_nod =1 ; i_nod <=nod; i_nod++){
  178 + list << i_nod << "\t";
  179 + }
169 180
170 list.close(); 181 list.close();
171 } 182 }
cpp/update_dir3.cuh
@@ -4,8 +4,6 @@ @@ -4,8 +4,6 @@
4 # include <iostream> 4 # include <iostream>
5 # include <cuda.h> 5 # include <cuda.h>
6 #include <stim/cuda/cudatools.h> 6 #include <stim/cuda/cudatools.h>
7 -#include <stim/cuda/sharedmem.cuh>  
8 -#include <cuda_fp16.h>  
9 #include "cpyToshare.cuh" 7 #include "cpyToshare.cuh"
10 8
11 // this kernel calculates the voting direction for the next iteration based on the angle between the location of this voter and the maximum vote value in its voting area. 9 // this kernel calculates the voting direction for the next iteration based on the angle between the location of this voter and the maximum vote value in its voting area.
@@ -4,7 +4,6 @@ @@ -4,7 +4,6 @@
4 #include <iostream> 4 #include <iostream>
5 #include <cuda.h> 5 #include <cuda.h>
6 #include <stim/cuda/cudatools.h> 6 #include <stim/cuda/cudatools.h>
7 -#include <stim/cuda/sharedmem.cuh>  
8 #include <stim/cuda/cudatools/error.h> 7 #include <stim/cuda/cudatools/error.h>
9 #include "cpyToshare.cuh" 8 #include "cpyToshare.cuh"
10 9