diff --git a/Matlab_3D/0-gt.vol b/Matlab_3D/0-gt.vol deleted file mode 100644 index a323ed5..0000000 Binary files a/Matlab_3D/0-gt.vol and /dev/null differ diff --git a/Matlab_3D/ivote3.m b/Matlab_3D/ivote3.m index bbdecf6..a6c76d0 100644 --- a/Matlab_3D/ivote3.m +++ b/Matlab_3D/ivote3.m @@ -4,11 +4,11 @@ disp('***************** NEW RUN *********************'); total = tic; % ******* Initialize voting parameters ************************************** rmax = [10 10 10]; %maximum radius of the cell -rmin = [1 1 1]; +%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 = 1; +t0 = 0; sigma = [5, 5, 5]; % t = 0.1; d_ang= ang / (iter+2); @@ -23,7 +23,7 @@ d_ang= ang / (iter+2); % X = S(1); % Y = S(2); % Z = S(3); -filename = '128-128-128/nissl-float-128.128.128.vol'; %'nissl-float-128.128.128.vol'; +filename = 'nissl-float-128.128.128.vol'; %'nissl-float-128.128.128.vol'; X = 128; Y = 128; Z = 128; @@ -50,16 +50,16 @@ Isize = size(I); %set a threshold for the gradient magnitude It = Imag > t0; -fidt = fopen('128-128-128/It.vol', 'w'); -fwrite(fidt, It, 'single'); -fclose(fidt); +% 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; +% 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) @@ -74,9 +74,9 @@ rangez = -rmax(3):rmax(3); 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_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) -M_dist2 = (mx.^2/rmin(1)^2 + my.^2/rmin(2)^2 + mz.^2/rmin(3)^2) >= 1 ; -M_dist = M_dist1 .* M_dist2; +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; @@ -91,8 +91,8 @@ g_v_prime = zeros(nV, ceil(rmax(1)*rmax(2)*rmax(3)*ang)); % vote tic; -mask = zeros(Isize); -mask1 = zeros(Isize); +% mask = zeros(Isize); +% mask1 = zeros(Isize); %for each iteration (in iterative voting) for itr = 1 :iter @@ -146,19 +146,29 @@ for itr = 1 :iter global_pz = vz + mz(pi); %convert the global 3D index of each point into a global 1D index - global_pi = sub2ind(Isize, global_px, global_py, global_pz); + 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); +% 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)]); + t_v1 = toc +% disp(['voting done. time =',num2str(t_v1)]); % update the voting direction if ang>=d_ang @@ -189,22 +199,23 @@ for itr = 1 :iter end - tdir1 = toc; - display (['updating dir done. time = ', num2str(tdir1)]); + 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); +% hv = reshape(Ivote, [X*Y*Z, 1]); +% hist(hv, 250); %% -t = 300; +t = 2000; conn = [5 5 5]; Icenter = local_max(Ivote, conn, t); -fidc = fopen(sprintf('std3.2-r10.10-v8/out%d-t%d.vol',t,t0), 'w'); +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)); diff --git a/Matlab_3D/validation.m b/Matlab_3D/validation.m index 5b2f62d..99b6aba 100644 --- a/Matlab_3D/validation.m +++ b/Matlab_3D/validation.m @@ -1,4 +1,3 @@ - clear all; disp('***************** NEW RUN *********************'); X = 128; @@ -6,21 +5,22 @@ Y = 128; Z = 128; D = 10; t0=1; -r1=10; +r1=12; r2=10; -t=2100; +t=2300; itr=8; vote=10; std = [5 5]; -gt_filename = '0-gt.vol'; +gt_filename = 'gt2.vol'; % 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); -out_filename = sprintf('D:/build/ivote3-bld/shared2D-v8/out%d.vol',t); +out_filename = sprintf('D:/build/ivote3-bld/shared2D-v8/centers.%d.vol',t); % 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); -txt_filename = sprintf('D:/build/ivote3-bld/shared2D-v8/t%d-D%d.txt',t,D); +txt_filename = sprintf('0-t%d--1.txt',t); spec = sprintf('Nissl-std%d.%d-r%d.%d-t%d-out%d.%d',std(1), std(2),r1,r2,t,itr,vote); fid0 = fopen(gt_filename); gt = fread(fid0,[X Y*Z], 'single'); fclose(fid0); +% store x,y, and z directions of cell centers defined by the ground-trurth in an array gt = reshape(gt, [X Y Z]); [gtx, gty, gtz] = ind2sub(size(gt),find(gt)); ref = [gtx gty gtz]; @@ -28,6 +28,7 @@ ref = [gtx gty gtz]; fid1 = fopen(out_filename); out = fread(fid1,[X Y*Z], 'single'); fclose(fid1); +% store x,y, and z directions of detected cells by the code in an array out = reshape(out, [X Y Z]); [cx, cy, cz] = ind2sub(size(out),find(out)); cent = [cx cy cz]; @@ -35,67 +36,135 @@ cent = [cx cy cz]; [idx, dist] = knnsearch(cent,ref); % 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. ind_eq_idx = zeros(numel(idx), 10); -t=0; +%set number of cells with the same detected cell center to zero. +noc=0; for i=1:numel(idx) + % check the current element hasn't considered yet s_ind = sum(ind_eq_idx==i); if s_ind==0 + % check how many of the elemnts assigned to the same detected cell s = sum(idx==idx(i)); if s>1 - t = t+1; - id_ref(t,1) = i; - id_ref(t,2) = s; + noc = noc+1; + %save the index and number of gt's elements with the same assigned detected cells. + id_ref(noc,1) = i; + id_ref(noc,2) = s; ind1 = find(idx==idx(i)); ind_eq_idx(i, 1:numel(ind1)) = ind1; end end end +% determine those indices which hs assigned to only one detected cell b = idx; u_eq_idx = unique(ind_eq_idx(ind_eq_idx>0)); b(u_eq_idx)=0; +% set the number of over sefgmented cells to zero. oseg=0; - for k=1:size(id_ref,1) k1 = id_ref(k,1); k2 = id_ref(k,2); + %find the minimum distance to the detected cell happened for which of the gt's element l = ind_eq_idx(k1,1:k2); [~, local_id_min] = min(dist(l)); + % set the element with minimum distance to the corresponding detected cell b(l(local_id_min)) = idx(k1); + %remove the proper element from the list of indices with same designated detected cells u_eq_idx (u_eq_idx == l(local_id_min))=0; + % check if the indices remained in the list has distance less than the value is set for validation ly = l (l~=l(local_id_min)); distl = dist(ly); sl = sum(distl0 oseg = oseg+1; end end - -% -u_idx1 = u_eq_idx (u_eq_idx>0); +%**** +% check if the oversegmented detected cells could be assigned to the gt's element with no detected cell! +%find the indices of gt's element with no detected cell +u_idx1 = find(b==0); +%find the coordinates of gt's ref1 = ref(u_idx1,:); +% find the indices of detected cells which are assigned to a gt's element ub = unique(b(b>0)); +% set the number of detected cell which are not assigned to a gt's element to zero mm=0; -for jj=1:size(cent,1) - z = sum(ub(:)==jj); +% find the indices of detected cell which are not assigned to a gt's element +for j=1:size(cent,1) + z = sum(ub(:)==j); if z==0 mm = mm+1; - cent_id1(mm) = jj; + cent_id1(mm) = j; end end - +% find the coordinated of of detected cells cent1 = cent(cent_id1,:); +% check if there are osme low enough distances [idx1, dist1] = knnsearch(cent1,ref1); -min1= min(dist1(:)); -if min1 < D - miss1 = sum(dist1(:)1 + noc_rc = noc_rc+1; + %save the index and number of gt's elements with the same assigned detected cells. + id_ref_rc(noc_rc,1) = i_rc; + id_ref_rc(noc_rc,2) = s_rc; + ind1_rc = find(idx1==idx1(i_rc)); + ind_eq_idx_rc(i_rc, 1:numel(ind1_rc)) = ind1_rc; + end + end +end +% determine those indices which hs assigned to only one detected cell +b_rc = idx1; +u_eq_idx_rc = unique(ind_eq_idx_rc(ind_eq_idx_rc>0)); +b_rc(u_eq_idx_rc)=0; +% set the number of over sefgmented cells to zero. +oseg_rc=0; +for k_rc=1:size(id_ref_rc,1) + k1_rc = id_ref_rc(k_rc,1); + k2_rc = id_ref_rc(k_rc,2); + %find the minimum distance to the detected cell happened for which of the gt's element + l_rc = ind_eq_idx_rc(k1_rc,1:k2_rc); + [~, local_id_min_rc] = min(dist1(l_rc)); + % set the element with minimum distance to the corresponding detected cell + b_rc(l_rc(local_id_min_rc)) = idx1(k1_rc); + %remove the proper element from the list of indices with same designated detected cells + u_eq_idx_rc (u_eq_idx_rc == l_rc(local_id_min_rc))=0; + % check if the indices remained in the list has distance less than the value is set for validation + ly_rc = l_rc (l_rc~=l_rc(local_id_min_rc)); + distl_rc = dist1(ly_rc); + sl_rc = sum(distl_rc0 + oseg_rc = oseg_rc+1; + end end -% + +%****** +% b include the gt's element and its detected cells, for those element with no cell detection, b has zero value b_ind = find(b==0); b_dist = dist; +% remove the disatances for those elements with no detected cells from distance array b_dist(b_ind)=-1; -TP = sum(b_dist>=0) - sum(b_dist>=D); + +b_ind_rc = find(b_rc==0); +b_dist_rc = dist1; +b_dist_rc(b_ind_rc)=-1; +% calculate Ttrue Positive, number of detected cells that have low enough distance to one of the gt's elements. +TP = sum(b_dist>=0) - sum(b_dist>D) + (sum(b_dist_rc>=0) - sum(b_dist_rc>D)); +% calculate False Negative, number of gt's elements with no detected cells. FN = size(ref, 1) - TP; +% calculate False Positive, number of detected cells, which their distance to the gt's element is long +% or considered as over segmentation. FP = size(cent,1) - TP; +% compute the validation factors. Precision = TP/(TP + FP); Recall = TP/(TP+FN); Accuracy= TP/(TP+FP+FN); diff --git a/cpp/cudafunc.cu b/cpp/cudafunc.cu index dcca2e7..5a89eaf 100644 --- a/cpp/cudafunc.cu +++ b/cpp/cudafunc.cu @@ -1,6 +1,3 @@ -//#include "cuda_fp16.h" -//#include "float_to_half.cuh" -//#include "half_to_float.cuh" #include "gaussian_blur3.cuh" #include "gradient3.cuh" #include "mag3.cuh" @@ -9,7 +6,7 @@ #include "local_max3.cuh" -void ivote3(float* center, float* img, float sigma[], float anisotropy, float phi, float d_phi, unsigned int r[], +void ivote3(float* img, float sigma[], float anisotropy, float phi, float d_phi, unsigned int r[], int iter, float t, unsigned int conn[], unsigned int x, unsigned int y, unsigned int z){ @@ -47,8 +44,8 @@ void ivote3(float* center, float* img, float sigma[], float anisotropy, float ph gpu_vote3(gpu_vote, gpu_grad, cos_phi, r, x, y, z); cudaDeviceSynchronize(); - if (i==7) - cudaMemcpy(img, gpu_vote, bytes, cudaMemcpyDeviceToHost); + /*if (i==7) + cudaMemcpy(img, gpu_vote, bytes, cudaMemcpyDeviceToHost);*/ if (phi >= d_phi){ gpu_update_dir3(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 } cudaFree(gpu_grad); - //cudaMemcpy(center, gpu_grad, bytes, cudaMemcpyDeviceToHost); + cudaMemcpy(img, gpu_vote, bytes, cudaMemcpyDeviceToHost); //allocate space on the gpu for the final detected cells. - float* gpu_output; - cudaMalloc(&gpu_output, bytes); + //float* gpu_output; + //cudaMalloc(&gpu_output, bytes); + + ////call the local max function + //gpu_local_max3(gpu_output, gpu_vote, t, conn, x, y, z); + + ////copy the final result to the cpu. + //cudaMemcpy(center, gpu_output, bytes, cudaMemcpyDeviceToHost); + // + // + cudaFree(gpu_vote); + //cudaFree(gpu_output); + +} + +void lmax(float* out, float* in, float t, unsigned int conn[], unsigned int x, unsigned int y, unsigned int z){ + unsigned int bytes = x * y * z * sizeof(float); + + //assign memory on gpu for the input data. + float* gpuV; + cudaMalloc(&gpuV, bytes); + + //copy the image data to the GPU. + cudaMemcpy(gpuV, in, bytes, cudaMemcpyHostToDevice); + + float* gpuOut; + cudaMalloc(&gpuOut, bytes); //call the local max function - gpu_local_max3(gpu_output, gpu_vote, t, conn, x, y, z); + gpu_local_max3(gpuOut, gpuV, t, conn, x, y, z); //copy the final result to the cpu. - cudaMemcpy(center, gpu_output, bytes, cudaMemcpyDeviceToHost); - - - cudaFree(gpu_vote); - cudaFree(gpu_output); + cudaMemcpy(out, gpuOut, bytes, cudaMemcpyDeviceToHost); + cudaFree(gpuV); + cudaFree(gpuOut); } \ No newline at end of file diff --git a/cpp/gaussian_blur3.cuh b/cpp/gaussian_blur3.cuh index 142ac6f..fa083ca 100644 --- a/cpp/gaussian_blur3.cuh +++ b/cpp/gaussian_blur3.cuh @@ -4,9 +4,7 @@ #include #include #include -#include #include -#include "cuda_fp16.h" #define pi 3.14159 diff --git a/cpp/mag3.cuh b/cpp/mag3.cuh index ca0bb06..a86fb91 100644 --- a/cpp/mag3.cuh +++ b/cpp/mag3.cuh @@ -4,7 +4,7 @@ #include #include #include -#include +//#include #include template diff --git a/cpp/main.cpp b/cpp/main.cpp index 5d5e11e..47d76d8 100644 --- a/cpp/main.cpp +++ b/cpp/main.cpp @@ -1,18 +1,18 @@ #include +#include #include #include #include #include #include -#include -#include #include #include #define pi 3.14159 -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[], +void ivote3(float* img, float std[], float anisotropy, float phi, float d_phi, unsigned int r[], int iter, float t, unsigned int conn[], unsigned int x, unsigned int y, unsigned int z); +void lmax(float* center, float* vote, float t1, unsigned int conn[], unsigned int x, unsigned int y, unsigned int z); void invert_data(float* cpuI, unsigned int x, unsigned int y, unsigned int z){ for(int ix = 0; ix < x; ix++){ @@ -99,10 +99,10 @@ int main(int argc, char** argv){ //set the anisotropy float anisotropy = args["anisotropy"].as_float(); unsigned int rmax = 10 ; - unsigned int r[3] = { rmax, rmax, rmax}; + unsigned int r[3] = { 12, rmax, rmax}; float std = 5; float sigma[3] = { std, std, std}; - unsigned int nlmax = 5; + unsigned int nlmax = 1; unsigned int conn[3] = { nlmax, nlmax, nlmax}; float phi_deg = 25.0; float phi = phi_deg * pi /180; @@ -123,7 +123,7 @@ int main(int argc, char** argv){ invert_data(cpuI, x, y, z); //write a new file from the cpuI. - std::ofstream original("shared2D-v8/inv-128.vol", std::ofstream::out | std::ofstream::binary); + std::ofstream original("shared2D-v1/inv-128.vol", std::ofstream::out | std::ofstream::binary); original.write((char*)cpuI, bytes); original.close(); @@ -131,8 +131,8 @@ int main(int argc, char** argv){ float* cpu_out = (float*) malloc(bytes); // call the ivote function - ivote3(cpu_out, cpuI, sigma, anisotropy, phi, d_phi, r, iter, t, conn, x, y, z); - + //ivote3(cpu_out, cpuI, sigma, anisotropy, phi, d_phi, r, iter, t, conn, x, y, z); + ivote3(cpuI, sigma, anisotropy, phi, d_phi, r, iter, t, conn, x, y, z); //write the blurred file from the cpuI. std::ofstream fblur("shared2D-v8/vote8.vol", std::ofstream::out | std::ofstream::binary); fblur.write((char*)cpuI, bytes); @@ -145,27 +145,38 @@ int main(int argc, char** argv){ fgx.close(); */ //write the output file. - std::ofstream fo("shared2D-v8/" + OutName.str(), std::ofstream::out | std::ofstream::binary); - fo.write((char*)cpu_out, bytes); - fo.close(); + //for (int t0=2000; t0<=2500; t0+=100){ + int t0 = t; + lmax(cpu_out, cpuI, t, conn, x, y, z); + //std::ofstream fo("shared2D-v8/" + OutName.str(), std::ofstream::out | std::ofstream::binary); + std::ofstream fo("shared2D-v8/" + OutName.str()+std::to_string(t0)+".vol", std::ofstream::out | std::ofstream::binary); + fo.write((char*)cpu_out, bytes); + fo.close(); // creat a file for saving the list centers - - std::ofstream list("shared2D-v8/center.txt"); + + std::ofstream list("shared2D-v8/" + OutName.str()+std::to_string(t0)+".obj"); + // set the number of detected cells to zero. + int nod = 0; if (list.is_open()){ - for (int ix=0; ix # include #include -#include -#include #include "cpyToshare.cuh" // 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. diff --git a/cpp/vote3.cuh b/cpp/vote3.cuh index 2e6b9ce..5af3cef 100644 --- a/cpp/vote3.cuh +++ b/cpp/vote3.cuh @@ -4,7 +4,6 @@ #include #include #include -#include #include #include "cpyToshare.cuh" -- libgit2 0.21.4