From d000db76c60654cdb0b07ea7f7967ceeebfbd73a Mon Sep 17 00:00:00 2001 From: Daniil Kazantsev Date: Tue, 14 May 2019 16:13:39 +0100 Subject: fixes all matlab issues --- demos/Matlab_demos/demoMatlab_3Ddenoise.m | 200 ++++++++++++++++ demos/Matlab_demos/demoMatlab_denoise.m | 188 +++++++++++++++ demos/Matlab_demos/demoMatlab_inpaint.m | 40 ++++ demos/demoMatlab_3Ddenoise.m | 198 ---------------- demos/demoMatlab_denoise.m | 188 --------------- demos/demoMatlab_inpaint.m | 35 --- src/Core/regularisers_CPU/MASK_merge_core.c | 262 --------------------- src/Core/regularisers_CPU/MASK_merge_core.h | 56 ----- src/Matlab/mex_compile/compileCPU_mex_Linux.m | 2 +- src/Matlab/mex_compile/compileCPU_mex_WINDOWS.m | 6 +- src/Matlab/mex_compile/compileGPU_mex.m | 8 +- src/Matlab/mex_compile/regularisers_CPU/ROF_TV.c | 30 ++- .../regularisers_GPU/PatchSelect_GPU.cpp | 94 ++++++++ .../mex_compile/regularisers_GPU/ROF_TV_GPU.cpp | 27 ++- 14 files changed, 578 insertions(+), 756 deletions(-) create mode 100644 demos/Matlab_demos/demoMatlab_3Ddenoise.m create mode 100644 demos/Matlab_demos/demoMatlab_denoise.m create mode 100644 demos/Matlab_demos/demoMatlab_inpaint.m delete mode 100644 demos/demoMatlab_3Ddenoise.m delete mode 100644 demos/demoMatlab_denoise.m delete mode 100644 demos/demoMatlab_inpaint.m delete mode 100644 src/Core/regularisers_CPU/MASK_merge_core.c delete mode 100644 src/Core/regularisers_CPU/MASK_merge_core.h create mode 100644 src/Matlab/mex_compile/regularisers_GPU/PatchSelect_GPU.cpp diff --git a/demos/Matlab_demos/demoMatlab_3Ddenoise.m b/demos/Matlab_demos/demoMatlab_3Ddenoise.m new file mode 100644 index 0000000..d7ff60c --- /dev/null +++ b/demos/Matlab_demos/demoMatlab_3Ddenoise.m @@ -0,0 +1,200 @@ +% Volume (3D) denoising demo using CCPi-RGL +clear; close all +fsep = '/'; + + +Path1 = sprintf(['..' fsep '..' fsep 'src' fsep 'Matlab' fsep 'mex_compile' fsep 'installed'], 1i); +Path2 = sprintf(['..' fsep 'data' fsep], 1i); +Path3 = sprintf(['..' fsep '..' fsep 'src' fsep 'Matlab' fsep 'supp'], 1i); +addpath(Path1); +addpath(Path2); +addpath(Path3); + +N = 512; +slices = 15; +vol3D = zeros(N,N,slices, 'single'); +Ideal3D = zeros(N,N,slices, 'single'); +Im = double(imread('lena_gray_512.tif'))/255; % loading image +for i = 1:slices +vol3D(:,:,i) = Im + .05*randn(size(Im)); +Ideal3D(:,:,i) = Im; +end +vol3D(vol3D < 0) = 0; +figure; imshow(vol3D(:,:,7), [0 1]); title('Noisy image'); +%% +fprintf('Denoise a volume using the ROF-TV model (CPU) \n'); +lambda_reg = 0.03; % regularsation parameter for all methods +tau_rof = 0.0025; % time-marching constant +iter_rof = 300; % number of ROF iterations +epsil_tol = 0.0; % tolerance +tic; [u_rof,infovec] = ROF_TV(single(vol3D), lambda_reg, iter_rof, tau_rof, epsil_tol); toc; +energyfunc_val_rof = TV_energy(single(u_rof),single(vol3D),lambda_reg, 1); % get energy function value +rmse_rof = (RMSE(Ideal3D(:),u_rof(:))); +fprintf('%s %f \n', 'RMSE error for ROF is:', rmse_rof); +figure; imshow(u_rof(:,:,7), [0 1]); title('ROF-TV denoised volume (CPU)'); +%% +% fprintf('Denoise a volume using the ROF-TV model (GPU) \n'); +% lambda_reg = 0.03; % regularsation parameter for all methods +% tau_rof = 0.0025; % time-marching constant +% iter_rof = 300; % number of ROF iterations +% epsil_tol = 0.0; % tolerance +% tic; u_rofG = ROF_TV_GPU(single(vol3D), lambda_reg, iter_rof, tau_rof, epsil_tol); toc; +% rmse_rofG = (RMSE(Ideal3D(:),u_rofG(:))); +% fprintf('%s %f \n', 'RMSE error for ROF is:', rmse_rofG); +% figure; imshow(u_rofG(:,:,7), [0 1]); title('ROF-TV denoised volume (GPU)'); +%% +fprintf('Denoise a volume using the FGP-TV model (CPU) \n'); +lambda_reg = 0.03; % regularsation parameter for all methods +iter_fgp = 300; % number of FGP iterations +epsil_tol = 0.0; % tolerance +tic; [u_fgp,infovec] = FGP_TV(single(vol3D), lambda_reg, iter_fgp, epsil_tol); toc; +energyfunc_val_fgp = TV_energy(single(u_fgp),single(vol3D),lambda_reg, 1); % get energy function value +rmse_fgp = (RMSE(Ideal3D(:),u_fgp(:))); +fprintf('%s %f \n', 'RMSE error for FGP-TV is:', rmse_fgp); +figure; imshow(u_fgp(:,:,7), [0 1]); title('FGP-TV denoised volume (CPU)'); +%% +fprintf('Denoise a volume using the FGP-TV model (GPU) \n'); +% lambda_reg = 0.03; % regularsation parameter for all methods +% iter_fgp = 300; % number of FGP iterations +% epsil_tol = 0.0; % tolerance +% tic; u_fgpG = FGP_TV_GPU(single(vol3D), lambda_reg, iter_fgp, epsil_tol); toc; +% rmse_fgpG = (RMSE(Ideal3D(:),u_fgpG(:))); +% fprintf('%s %f \n', 'RMSE error for FGP-TV is:', rmse_fgpG); +% figure; imshow(u_fgpG(:,:,7), [0 1]); title('FGP-TV denoised volume (GPU)'); +%% +fprintf('Denoise a volume using the SB-TV model (CPU) \n'); +iter_sb = 150; % number of SB iterations +epsil_tol = 0.0; % tolerance +tic; [u_sb,infovec] = SB_TV(single(vol3D), lambda_reg, iter_sb, epsil_tol); toc; +energyfunc_val_sb = TV_energy(single(u_sb),single(vol3D),lambda_reg, 1); % get energy function value +rmse_sb = (RMSE(Ideal3D(:),u_sb(:))); +fprintf('%s %f \n', 'RMSE error for SB-TV is:', rmse_sb); +figure; imshow(u_sb(:,:,7), [0 1]); title('SB-TV denoised volume (CPU)'); +%% +% fprintf('Denoise a volume using the SB-TV model (GPU) \n'); +% iter_sb = 150; % number of SB iterations +% epsil_tol = 0.0; % tolerance +% tic; u_sbG = SB_TV_GPU(single(vol3D), lambda_reg, iter_sb, epsil_tol); toc; +% rmse_sbG = (RMSE(Ideal3D(:),u_sbG(:))); +% fprintf('%s %f \n', 'RMSE error for SB-TV is:', rmse_sbG); +% figure; imshow(u_sbG(:,:,7), [0 1]); title('SB-TV denoised volume (GPU)'); +%% +fprintf('Denoise a volume using the ROF-LLT model (CPU) \n'); +lambda_ROF = lambda_reg; % ROF regularisation parameter +lambda_LLT = lambda_reg*0.35; % LLT regularisation parameter +iter_LLT = 300; % iterations +tau_rof_llt = 0.0025; % time-marching constant +epsil_tol = 0.0; % tolerance +tic; [u_rof_llt, infovec] = LLT_ROF(single(vol3D), lambda_ROF, lambda_LLT, iter_LLT, tau_rof_llt, epsil_tol); toc; +rmse_rof_llt = (RMSE(Ideal3D(:),u_rof_llt(:))); +fprintf('%s %f \n', 'RMSE error for ROF-LLT is:', rmse_rof_llt); +figure; imshow(u_rof_llt(:,:,7), [0 1]); title('ROF-LLT denoised volume (CPU)'); +%% +% fprintf('Denoise a volume using the ROF-LLT model (GPU) \n'); +% lambda_ROF = lambda_reg; % ROF regularisation parameter +% lambda_LLT = lambda_reg*0.35; % LLT regularisation parameter +% iter_LLT = 300; % iterations +% tau_rof_llt = 0.0025; % time-marching constant +% epsil_tol = 0.0; % tolerance +% tic; u_rof_llt_g = LLT_ROF_GPU(single(vol3D), lambda_ROF, lambda_LLT, iter_LLT, tau_rof_llt, epsil_tol); toc; +% rmse_rof_llt = (RMSE(Ideal3D(:),u_rof_llt_g(:))); +% fprintf('%s %f \n', 'RMSE error for ROF-LLT is:', rmse_rof_llt); +% figure; imshow(u_rof_llt_g(:,:,7), [0 1]); title('ROF-LLT denoised volume (GPU)'); +%% +fprintf('Denoise a volume using Nonlinear-Diffusion model (CPU) \n'); +iter_diff = 300; % number of diffusion iterations +lambda_regDiff = 0.025; % regularisation for the diffusivity +sigmaPar = 0.015; % edge-preserving parameter +tau_param = 0.025; % time-marching constant +epsil_tol = 0.0; % tolerance +tic; [u_diff, infovec] = NonlDiff(single(vol3D), lambda_regDiff, sigmaPar, iter_diff, tau_param, 'Huber', epsil_tol); toc; +rmse_diff = (RMSE(Ideal3D(:),u_diff(:))); +fprintf('%s %f \n', 'RMSE error for Diffusion is:', rmse_diff); +figure; imshow(u_diff(:,:,7), [0 1]); title('Diffusion denoised volume (CPU)'); +%% +% fprintf('Denoise a volume using Nonlinear-Diffusion model (GPU) \n'); +% iter_diff = 300; % number of diffusion iterations +% lambda_regDiff = 0.025; % regularisation for the diffusivity +% sigmaPar = 0.015; % edge-preserving parameter +% tau_param = 0.025; % time-marching constant +% tic; u_diff_g = NonlDiff_GPU(single(vol3D), lambda_regDiff, sigmaPar, iter_diff, tau_param, 'Huber', epsil_tol); toc; +% rmse_diff = (RMSE(Ideal3D(:),u_diff_g(:))); +% fprintf('%s %f \n', 'RMSE error for Diffusion is:', rmse_diff); +% figure; imshow(u_diff_g(:,:,7), [0 1]); title('Diffusion denoised volume (GPU)'); +%% +fprintf('Denoise using Fourth-order anisotropic diffusion model (CPU) \n'); +iter_diff = 300; % number of diffusion iterations +lambda_regDiff = 3.5; % regularisation for the diffusivity +sigmaPar = 0.02; % edge-preserving parameter +tau_param = 0.0015; % time-marching constant +epsil_tol = 0.0; % tolerance +tic; u_diff4 = Diffusion_4thO(single(vol3D), lambda_regDiff, sigmaPar, iter_diff, tau_param, epsil_tol); toc; +rmse_diff4 = (RMSE(Ideal3D(:),u_diff4(:))); +fprintf('%s %f \n', 'RMSE error for Anis.Diff of 4th order is:', rmse_diff4); +figure; imshow(u_diff4(:,:,7), [0 1]); title('Diffusion 4thO denoised volume (CPU)'); +%% +% fprintf('Denoise using Fourth-order anisotropic diffusion model (GPU) \n'); +% iter_diff = 300; % number of diffusion iterations +% lambda_regDiff = 3.5; % regularisation for the diffusivity +% sigmaPar = 0.02; % edge-preserving parameter +% tau_param = 0.0015; % time-marching constant +% tic; u_diff4_g = Diffusion_4thO_GPU(single(vol3D), lambda_regDiff, sigmaPar, iter_diff, tau_param, epsil_tol); toc; +% rmse_diff4 = (RMSE(Ideal3D(:),u_diff4_g(:))); +% fprintf('%s %f \n', 'RMSE error for Anis.Diff of 4th order is:', rmse_diff4); +% figure; imshow(u_diff4_g(:,:,7), [0 1]); title('Diffusion 4thO denoised volume (GPU)'); +%% +fprintf('Denoise using the TGV model (CPU) \n'); +lambda_TGV = 0.03; % regularisation parameter +alpha1 = 1.0; % parameter to control the first-order term +alpha0 = 2.0; % parameter to control the second-order term +L2 = 12.0; % convergence parameter +iter_TGV = 500; % number of Primal-Dual iterations for TGV +epsil_tol = 0.0; % tolerance +tic; u_tgv = TGV(single(vol3D), lambda_TGV, alpha1, alpha0, iter_TGV, L2, epsil_tol); toc; +rmseTGV = RMSE(Ideal3D(:),u_tgv(:)); +fprintf('%s %f \n', 'RMSE error for TGV is:', rmseTGV); +figure; imshow(u_tgv(:,:,3), [0 1]); title('TGV denoised volume (CPU)'); +%% +% fprintf('Denoise using the TGV model (GPU) \n'); +% lambda_TGV = 0.03; % regularisation parameter +% alpha1 = 1.0; % parameter to control the first-order term +% alpha0 = 2.0; % parameter to control the second-order term +% iter_TGV = 500; % number of Primal-Dual iterations for TGV +% tic; u_tgv_gpu = TGV_GPU(single(vol3D), lambda_TGV, alpha1, alpha0, iter_TGV, L2, epsil_tol); toc; +% rmseTGV = RMSE(Ideal3D(:),u_tgv_gpu(:)); +% fprintf('%s %f \n', 'RMSE error for TGV is:', rmseTGV); +% figure; imshow(u_tgv_gpu(:,:,3), [0 1]); title('TGV denoised volume (GPU)'); +%% +%>>>>>>>>>>>>>> MULTI-CHANNEL priors <<<<<<<<<<<<<<< % +fprintf('Denoise a volume using the FGP-dTV model (CPU) \n'); + +% create another volume (reference) with slightly less amount of noise +vol3D_ref = zeros(N,N,slices, 'single'); +for i = 1:slices +vol3D_ref(:,:,i) = Im + .01*randn(size(Im)); +end +vol3D_ref(vol3D_ref < 0) = 0; +% vol3D_ref = zeros(size(Im),'single'); % pass zero reference (dTV -> TV) + +iter_fgp = 300; % number of FGP iterations +epsil_tol = 0.0; % tolerance +eta = 0.2; % Reference image gradient smoothing constant +tic; u_fgp_dtv = FGP_dTV(single(vol3D), single(vol3D_ref), lambda_reg, iter_fgp, epsil_tol, eta); toc; +figure; imshow(u_fgp_dtv(:,:,7), [0 1]); title('FGP-dTV denoised volume (CPU)'); +%% +fprintf('Denoise a volume using the FGP-dTV model (GPU) \n'); + +% create another volume (reference) with slightly less amount of noise +vol3D_ref = zeros(N,N,slices, 'single'); +for i = 1:slices +vol3D_ref(:,:,i) = Im + .01*randn(size(Im)); +end +vol3D_ref(vol3D_ref < 0) = 0; +% vol3D_ref = zeros(size(Im),'single'); % pass zero reference (dTV -> TV) + +iter_fgp = 300; % number of FGP iterations +epsil_tol = 0.0; % tolerance +eta = 0.2; % Reference image gradient smoothing constant +tic; u_fgp_dtv_g = FGP_dTV_GPU(single(vol3D), single(vol3D_ref), lambda_reg, iter_fgp, epsil_tol, eta); toc; +figure; imshow(u_fgp_dtv_g(:,:,7), [0 1]); title('FGP-dTV denoised volume (GPU)'); +%% diff --git a/demos/Matlab_demos/demoMatlab_denoise.m b/demos/Matlab_demos/demoMatlab_denoise.m new file mode 100644 index 0000000..5af927f --- /dev/null +++ b/demos/Matlab_demos/demoMatlab_denoise.m @@ -0,0 +1,188 @@ +% Image (2D) denoising demo using CCPi-RGL +clear; close all +fsep = '/'; + +Path1 = sprintf(['..' fsep '..' fsep 'src' fsep 'Matlab' fsep 'mex_compile' fsep 'installed'], 1i); +Path2 = sprintf(['..' fsep 'data' fsep], 1i); +Path3 = sprintf(['..' fsep '..' fsep 'src' fsep 'Matlab' fsep 'supp'], 1i); +addpath(Path1); +addpath(Path2); +addpath(Path3); + +Im = double(imread('lena_gray_512.tif'))/255; % loading image +u0 = Im + .05*randn(size(Im)); u0(u0 < 0) = 0; +figure; imshow(u0, [0 1]); title('Noisy image'); +%% +fprintf('Denoise using the ROF-TV model (CPU) \n'); +lambda_reg = 0.03; % regularsation parameter for all methods +iter_rof = 1500; % number of ROF iterations +tau_rof = 0.003; % time-marching constant +epsil_tol = 0.0; % tolerance / 1.0e-06 +tic; [u_rof,infovec] = ROF_TV(single(u0), lambda_reg, iter_rof, tau_rof, epsil_tol); toc; +energyfunc_val_rof = TV_energy(single(u_rof),single(u0),lambda_reg, 1); % get energy function value +rmseROF = (RMSE(u_rof(:),Im(:))); +fprintf('%s %f \n', 'RMSE error for ROF-TV is:', rmseROF); +[ssimval] = ssim(u_rof*255,single(Im)*255); +fprintf('%s %f \n', 'MSSIM error for ROF-TV is:', ssimval); +figure; imshow(u_rof, [0 1]); title('ROF-TV denoised image (CPU)'); +%% +%fprintf('Denoise using the ROF-TV model (GPU) \n'); +%tic; [u_rofG,infovec] = ROF_TV_GPU(single(u0), lambda_reg, iter_rof, tau_rof, epsil_tol); toc; +%figure; imshow(u_rofG, [0 1]); title('ROF-TV denoised image (GPU)'); +%% +fprintf('Denoise using the FGP-TV model (CPU) \n'); +lambda_reg = 0.03; +iter_fgp = 500; % number of FGP iterations +epsil_tol = 0.0; % tolerance +tic; [u_fgp,infovec] = FGP_TV(single(u0), lambda_reg, iter_fgp, epsil_tol); toc; +energyfunc_val_fgp = TV_energy(single(u_fgp),single(u0),lambda_reg, 1); % get energy function value +rmseFGP = (RMSE(u_fgp(:),Im(:))); +fprintf('%s %f \n', 'RMSE error for FGP-TV is:', rmseFGP); +[ssimval] = ssim(u_fgp*255,single(Im)*255); +fprintf('%s %f \n', 'MSSIM error for FGP-TV is:', ssimval); +figure; imshow(u_fgp, [0 1]); title('FGP-TV denoised image (CPU)'); +%% +% fprintf('Denoise using the FGP-TV model (GPU) \n'); +% tic; u_fgpG = FGP_TV_GPU(single(u0), lambda_reg, iter_fgp, epsil_tol); toc; +% figure; imshow(u_fgpG, [0 1]); title('FGP-TV denoised image (GPU)'); +%% +fprintf('Denoise using the SB-TV model (CPU) \n'); +lambda_reg = 0.03; +iter_sb = 200; % number of SB iterations +epsil_tol = 0.0; % tolerance +tic; [u_sb,infovec] = SB_TV(single(u0), lambda_reg, iter_sb, epsil_tol); toc; +energyfunc_val_sb = TV_energy(single(u_sb),single(u0),lambda_reg, 1); % get energy function value +rmseSB = (RMSE(u_sb(:),Im(:))); +fprintf('%s %f \n', 'RMSE error for SB-TV is:', rmseSB); +[ssimval] = ssim(u_sb*255,single(Im)*255); +fprintf('%s %f \n', 'MSSIM error for SB-TV is:', ssimval); +figure; imshow(u_sb, [0 1]); title('SB-TV denoised image (CPU)'); +%% +% fprintf('Denoise using the SB-TV model (GPU) \n'); +% tic; u_sbG = SB_TV_GPU(single(u0), lambda_reg, iter_sb, epsil_tol); toc; +% figure; imshow(u_sbG, [0 1]); title('SB-TV denoised image (GPU)'); +%% +fprintf('Denoise using Nonlinear-Diffusion model (CPU) \n'); +iter_diff = 450; % number of diffusion iterations +lambda_regDiff = 0.025; % regularisation for the diffusivity +sigmaPar = 0.015; % edge-preserving parameter +tau_param = 0.02; % time-marching constant +epsil_tol = 0.0; % tolerance +tic; [u_diff,infovec] = NonlDiff(single(u0), lambda_regDiff, sigmaPar, iter_diff, tau_param, 'Huber', epsil_tol); toc; +rmseDiffus = (RMSE(u_diff(:),Im(:))); +fprintf('%s %f \n', 'RMSE error for Nonlinear Diffusion is:', rmseDiffus); +[ssimval] = ssim(u_diff*255,single(Im)*255); +fprintf('%s %f \n', 'MSSIM error for NDF is:', ssimval); +figure; imshow(u_diff, [0 1]); title('Diffusion denoised image (CPU)'); +%% +%fprintf('Denoise using Nonlinear-Diffusion model (GPU) \n'); +%tic; u_diff_g = NonlDiff_GPU(single(u0), lambda_regDiff, sigmaPar, iter_diff, tau_param, 'Huber', epsil_tol); toc; +%figure; imshow(u_diff_g, [0 1]); title('Diffusion denoised image (GPU)'); +%% +fprintf('Denoise using the TGV model (CPU) \n'); +lambda_TGV = 0.035; % regularisation parameter +alpha1 = 1.0; % parameter to control the first-order term +alpha0 = 2.0; % parameter to control the second-order term +L2 = 12.0; % convergence parameter +iter_TGV = 1200; % number of Primal-Dual iterations for TGV +epsil_tol = 0.0; % tolerance +tic; [u_tgv,infovec] = TGV(single(u0), lambda_TGV, alpha1, alpha0, iter_TGV, L2, epsil_tol); toc; +figure; imshow(u_tgv, [0 1]); title('TGV denoised image (CPU)'); +rmseTGV = (RMSE(u_tgv(:),Im(:))); +fprintf('%s %f \n', 'RMSE error for TGV is:', rmseTGV); +[ssimval] = ssim(u_tgv*255,single(Im)*255); +fprintf('%s %f \n', 'MSSIM error for TGV is:', ssimval); +%% +% fprintf('Denoise using the TGV model (GPU) \n'); +% tic; u_tgv_gpu = TGV_GPU(single(u0), lambda_TGV, alpha1, alpha0, iter_TGV, L2, epsil_tol); toc; +% figure; imshow(u_tgv_gpu, [0 1]); title('TGV denoised image (GPU)'); +%% +fprintf('Denoise using the ROF-LLT model (CPU) \n'); +lambda_ROF = 0.02; % ROF regularisation parameter +lambda_LLT = 0.015; % LLT regularisation parameter +iter_LLT = 2000; % iterations +tau_rof_llt = 0.01; % time-marching constant +epsil_tol = 0.0; % tolerance +tic; [u_rof_llt,infovec] = LLT_ROF(single(u0), lambda_ROF, lambda_LLT, iter_LLT, tau_rof_llt,epsil_tol); toc; +rmseROFLLT = (RMSE(u_rof_llt(:),Im(:))); +fprintf('%s %f \n', 'RMSE error for TGV is:', rmseROFLLT); +[ssimval] = ssim(u_rof_llt*255,single(Im)*255); +fprintf('%s %f \n', 'MSSIM error for ROFLLT is:', ssimval); +figure; imshow(u_rof_llt, [0 1]); title('ROF-LLT denoised image (CPU)'); +%% +% fprintf('Denoise using the ROF-LLT model (GPU) \n'); +% tic; u_rof_llt_g = LLT_ROF_GPU(single(u0), lambda_ROF, lambda_LLT, iter_LLT, tau_rof_llt, epsil_tol); toc; +% figure; imshow(u_rof_llt_g, [0 1]); title('ROF-LLT denoised image (GPU)'); +%% +fprintf('Denoise using Fourth-order anisotropic diffusion model (CPU) \n'); +iter_diff = 800; % number of diffusion iterations +lambda_regDiff = 3; % regularisation for the diffusivity +sigmaPar = 0.03; % edge-preserving parameter +tau_param = 0.0025; % time-marching constant +epsil_tol = 0.0; % tolerance +tic; [u_diff4,infovec] = Diffusion_4thO(single(u0), lambda_regDiff, sigmaPar, iter_diff, tau_param, epsil_tol); toc; +rmseDiffHO = (RMSE(u_diff4(:),Im(:))); +fprintf('%s %f \n', 'RMSE error for Fourth-order anisotropic diffusion is:', rmseDiffHO); +[ssimval] = ssim(u_diff4*255,single(Im)*255); +fprintf('%s %f \n', 'MSSIM error for DIFF4th is:', ssimval); +figure; imshow(u_diff4, [0 1]); title('Diffusion 4thO denoised image (CPU)'); +%% +%fprintf('Denoise using Fourth-order anisotropic diffusion model (GPU) \n'); +%tic; u_diff4_g = Diffusion_4thO_GPU(single(u0), lambda_regDiff, sigmaPar, iter_diff, tau_param); toc; +%figure; imshow(u_diff4_g, [0 1]); title('Diffusion 4thO denoised image (GPU)'); +%% +fprintf('Weights pre-calculation for Non-local TV (takes time on CPU) \n'); +SearchingWindow = 7; +PatchWindow = 2; +NeighboursNumber = 20; % the number of neibours to include +h = 0.23; % edge related parameter for NLM +tic; [H_i, H_j, Weights] = PatchSelect(single(u0), SearchingWindow, PatchWindow, NeighboursNumber, h); toc; +%% +fprintf('Denoise using Non-local Total Variation (CPU) \n'); +iter_nltv = 3; % number of nltv iterations +lambda_nltv = 0.055; % regularisation parameter for nltv +tic; u_nltv = Nonlocal_TV(single(u0), H_i, H_j, 0, Weights, lambda_nltv, iter_nltv); toc; +rmse_nltv = (RMSE(u_nltv(:),Im(:))); +fprintf('%s %f \n', 'RMSE error for Non-local Total Variation is:', rmse_nltv); +[ssimval] = ssim(u_nltv*255,single(Im)*255); +fprintf('%s %f \n', 'MSSIM error for NLTV is:', ssimval); +figure; imagesc(u_nltv, [0 1]); colormap(gray); daspect([1 1 1]); title('Non-local Total Variation denoised image (CPU)'); +%% +%>>>>>>>>>>>>>> MULTI-CHANNEL priors <<<<<<<<<<<<<<< % + +fprintf('Denoise using the FGP-dTV model (CPU) \n'); +% create another image (reference) with slightly less amount of noise +u_ref = Im + .01*randn(size(Im)); u_ref(u_ref < 0) = 0; +% u_ref = zeros(size(Im),'single'); % pass zero reference (dTV -> TV) + +lambda_reg = 0.04; +iter_fgp = 1000; % number of FGP iterations +epsil_tol = 0.0; % tolerance +eta = 0.2; % Reference image gradient smoothing constant +tic; [u_fgp_dtv,infovec] = FGP_dTV(single(u0), single(u_ref), lambda_reg, iter_fgp, epsil_tol, eta); toc; +rmse_dTV= (RMSE(u_fgp_dtv(:),Im(:))); +fprintf('%s %f \n', 'RMSE error for Directional Total Variation (dTV) is:', rmse_dTV); +figure; imshow(u_fgp_dtv, [0 1]); title('FGP-dTV denoised image (CPU)'); +%% +% fprintf('Denoise using the FGP-dTV model (GPU) \n'); +% % create another image (reference) with slightly less amount of noise +% u_ref = Im + .01*randn(size(Im)); u_ref(u_ref < 0) = 0; +% % u_ref = zeros(size(Im),'single'); % pass zero reference (dTV -> TV) +% +% iter_fgp = 1000; % number of FGP iterations +% epsil_tol = 1.0e-06; % tolerance +% eta = 0.2; % Reference image gradient smoothing constant +% tic; u_fgp_dtvG = FGP_dTV_GPU(single(u0), single(u_ref), lambda_reg, iter_fgp, epsil_tol, eta); toc; +% figure; imshow(u_fgp_dtvG, [0 1]); title('FGP-dTV denoised image (GPU)'); +%% +fprintf('Denoise using the TNV prior (CPU) \n'); +slices = 5; N = 512; +vol3D = zeros(N,N,slices, 'single'); +for i = 1:slices +vol3D(:,:,i) = Im + .05*randn(size(Im)); +end +vol3D(vol3D < 0) = 0; + +iter_tnv = 200; % number of TNV iterations +tic; u_tnv = TNV(single(vol3D), lambda_reg, iter_tnv); toc; +figure; imshow(u_tnv(:,:,3), [0 1]); title('TNV denoised stack of channels (CPU)'); diff --git a/demos/Matlab_demos/demoMatlab_inpaint.m b/demos/Matlab_demos/demoMatlab_inpaint.m new file mode 100644 index 0000000..67a6a23 --- /dev/null +++ b/demos/Matlab_demos/demoMatlab_inpaint.m @@ -0,0 +1,40 @@ +% Image (2D) inpainting demo using CCPi-RGL +clear; close all + +fsep = '/'; + +Path1 = sprintf(['..' fsep '..' fsep 'src' fsep 'Matlab' fsep 'mex_compile' fsep 'installed'], 1i); +Path2 = sprintf(['..' fsep 'data' fsep], 1i); +Path3 = sprintf(['..' fsep '..' fsep 'src' fsep 'Matlab' fsep 'supp'], 1i); +addpath(Path1); +addpath(Path2); +addpath(Path3); + +load('SinoInpaint.mat'); +Sinogram = Sinogram./max(Sinogram(:)); +Sino_mask = Sinogram.*(1-single(Mask)); +figure; +subplot(1,2,1); imshow(Sino_mask, [0 1]); title('Missing data sinogram'); +subplot(1,2,2); imshow(Mask, [0 1]); title('Mask'); +%% +fprintf('Inpaint using Linear-Diffusion model (CPU) \n'); +iter_diff = 5000; % number of diffusion iterations +lambda_regDiff = 6000; % regularisation for the diffusivity +sigmaPar = 0.0; % edge-preserving parameter +tau_param = 0.000075; % time-marching constant +tic; u_diff = NonlDiff_Inp(single(Sino_mask), Mask, lambda_regDiff, sigmaPar, iter_diff, tau_param); toc; +figure; imshow(u_diff, [0 1]); title('Linear-Diffusion inpainted sinogram (CPU)'); +%% +fprintf('Inpaint using Nonlinear-Diffusion model (CPU) \n'); +iter_diff = 1500; % number of diffusion iterations +lambda_regDiff = 80; % regularisation for the diffusivity +sigmaPar = 0.00009; % edge-preserving parameter +tau_param = 0.000008; % time-marching constant +tic; u_diff = NonlDiff_Inp(single(Sino_mask), Mask, lambda_regDiff, sigmaPar, iter_diff, tau_param, 'Huber'); toc; +figure; imshow(u_diff, [0 1]); title('Non-Linear Diffusion inpainted sinogram (CPU)'); +%% +fprintf('Inpaint using Nonlocal Vertical Marching model (CPU) \n'); +Increment = 1; % linear increment for the searching window +tic; [u_nom,maskupd] = NonlocalMarching_Inpaint(single(Sino_mask), Mask, Increment); toc; +figure; imshow(u_nom, [0 1]); title('NVM inpainted sinogram (CPU)'); +%% \ No newline at end of file diff --git a/demos/demoMatlab_3Ddenoise.m b/demos/demoMatlab_3Ddenoise.m deleted file mode 100644 index 3942eea..0000000 --- a/demos/demoMatlab_3Ddenoise.m +++ /dev/null @@ -1,198 +0,0 @@ -% Volume (3D) denoising demo using CCPi-RGL -clear; close all -Path1 = sprintf(['..' filesep 'src' filesep 'Matlab' filesep 'mex_compile' filesep 'installed'], 1i); -Path2 = sprintf(['data' filesep], 1i); -Path3 = sprintf(['..' filesep 'src' filesep 'Matlab' filesep 'supp'], 1i); -addpath(Path1); -addpath(Path2); -addpath(Path3); - -N = 512; -slices = 15; -vol3D = zeros(N,N,slices, 'single'); -Ideal3D = zeros(N,N,slices, 'single'); -Im = double(imread('lena_gray_512.tif'))/255; % loading image -for i = 1:slices -vol3D(:,:,i) = Im + .05*randn(size(Im)); -Ideal3D(:,:,i) = Im; -end -vol3D(vol3D < 0) = 0; -figure; imshow(vol3D(:,:,7), [0 1]); title('Noisy image'); - -%% -fprintf('Denoise a volume using the ROF-TV model (CPU) \n'); -lambda_reg = 0.03; % regularsation parameter for all methods -tau_rof = 0.0025; % time-marching constant -iter_rof = 300; % number of ROF iterations -epsil_tol = 0.0; % tolerance -tic; [u_rof,infovec] = ROF_TV(single(vol3D), lambda_reg, iter_rof, tau_rof, epsil_tol); toc; -energyfunc_val_rof = TV_energy(single(u_rof),single(vol3D),lambda_reg, 1); % get energy function value -rmse_rof = (RMSE(Ideal3D(:),u_rof(:))); -fprintf('%s %f \n', 'RMSE error for ROF is:', rmse_rof); -figure; imshow(u_rof(:,:,7), [0 1]); title('ROF-TV denoised volume (CPU)'); -%% -% fprintf('Denoise a volume using the ROF-TV model (GPU) \n'); -% lambda_reg = 0.03; % regularsation parameter for all methods -% tau_rof = 0.0025; % time-marching constant -% iter_rof = 300; % number of ROF iterations -% epsil_tol = 0.0; % tolerance -% tic; u_rofG = ROF_TV_GPU(single(vol3D), lambda_reg, iter_rof, tau_rof, epsil_tol); toc; -% rmse_rofG = (RMSE(Ideal3D(:),u_rofG(:))); -% fprintf('%s %f \n', 'RMSE error for ROF is:', rmse_rofG); -% figure; imshow(u_rofG(:,:,7), [0 1]); title('ROF-TV denoised volume (GPU)'); -%% -fprintf('Denoise a volume using the FGP-TV model (CPU) \n'); -lambda_reg = 0.03; % regularsation parameter for all methods -iter_fgp = 300; % number of FGP iterations -epsil_tol = 0.0; % tolerance -tic; [u_fgp,infovec] = FGP_TV(single(vol3D), lambda_reg, iter_fgp, epsil_tol); toc; -energyfunc_val_fgp = TV_energy(single(u_fgp),single(vol3D),lambda_reg, 1); % get energy function value -rmse_fgp = (RMSE(Ideal3D(:),u_fgp(:))); -fprintf('%s %f \n', 'RMSE error for FGP-TV is:', rmse_fgp); -figure; imshow(u_fgp(:,:,7), [0 1]); title('FGP-TV denoised volume (CPU)'); -%% -fprintf('Denoise a volume using the FGP-TV model (GPU) \n'); -% lambda_reg = 0.03; % regularsation parameter for all methods -% iter_fgp = 300; % number of FGP iterations -% epsil_tol = 0.0; % tolerance -% tic; u_fgpG = FGP_TV_GPU(single(vol3D), lambda_reg, iter_fgp, epsil_tol); toc; -% rmse_fgpG = (RMSE(Ideal3D(:),u_fgpG(:))); -% fprintf('%s %f \n', 'RMSE error for FGP-TV is:', rmse_fgpG); -% figure; imshow(u_fgpG(:,:,7), [0 1]); title('FGP-TV denoised volume (GPU)'); -%% -fprintf('Denoise a volume using the SB-TV model (CPU) \n'); -iter_sb = 150; % number of SB iterations -epsil_tol = 0.0; % tolerance -tic; [u_sb,infovec] = SB_TV(single(vol3D), lambda_reg, iter_sb, epsil_tol); toc; -energyfunc_val_sb = TV_energy(single(u_sb),single(vol3D),lambda_reg, 1); % get energy function value -rmse_sb = (RMSE(Ideal3D(:),u_sb(:))); -fprintf('%s %f \n', 'RMSE error for SB-TV is:', rmse_sb); -figure; imshow(u_sb(:,:,7), [0 1]); title('SB-TV denoised volume (CPU)'); -%% -% fprintf('Denoise a volume using the SB-TV model (GPU) \n'); -% iter_sb = 150; % number of SB iterations -% epsil_tol = 0.0; % tolerance -% tic; u_sbG = SB_TV_GPU(single(vol3D), lambda_reg, iter_sb, epsil_tol); toc; -% rmse_sbG = (RMSE(Ideal3D(:),u_sbG(:))); -% fprintf('%s %f \n', 'RMSE error for SB-TV is:', rmse_sbG); -% figure; imshow(u_sbG(:,:,7), [0 1]); title('SB-TV denoised volume (GPU)'); -%% -fprintf('Denoise a volume using the ROF-LLT model (CPU) \n'); -lambda_ROF = lambda_reg; % ROF regularisation parameter -lambda_LLT = lambda_reg*0.35; % LLT regularisation parameter -iter_LLT = 300; % iterations -tau_rof_llt = 0.0025; % time-marching constant -epsil_tol = 0.0; % tolerance -tic; [u_rof_llt, infovec] = LLT_ROF(single(vol3D), lambda_ROF, lambda_LLT, iter_LLT, tau_rof_llt, epsil_tol); toc; -rmse_rof_llt = (RMSE(Ideal3D(:),u_rof_llt(:))); -fprintf('%s %f \n', 'RMSE error for ROF-LLT is:', rmse_rof_llt); -figure; imshow(u_rof_llt(:,:,7), [0 1]); title('ROF-LLT denoised volume (CPU)'); -%% -% fprintf('Denoise a volume using the ROF-LLT model (GPU) \n'); -% lambda_ROF = lambda_reg; % ROF regularisation parameter -% lambda_LLT = lambda_reg*0.35; % LLT regularisation parameter -% iter_LLT = 300; % iterations -% tau_rof_llt = 0.0025; % time-marching constant -% epsil_tol = 0.0; % tolerance -% tic; u_rof_llt_g = LLT_ROF_GPU(single(vol3D), lambda_ROF, lambda_LLT, iter_LLT, tau_rof_llt, epsil_tol); toc; -% rmse_rof_llt = (RMSE(Ideal3D(:),u_rof_llt_g(:))); -% fprintf('%s %f \n', 'RMSE error for ROF-LLT is:', rmse_rof_llt); -% figure; imshow(u_rof_llt_g(:,:,7), [0 1]); title('ROF-LLT denoised volume (GPU)'); -%% -fprintf('Denoise a volume using Nonlinear-Diffusion model (CPU) \n'); -iter_diff = 300; % number of diffusion iterations -lambda_regDiff = 0.025; % regularisation for the diffusivity -sigmaPar = 0.015; % edge-preserving parameter -tau_param = 0.025; % time-marching constant -epsil_tol = 0.0; % tolerance -tic; [u_diff, infovec] = NonlDiff(single(vol3D), lambda_regDiff, sigmaPar, iter_diff, tau_param, 'Huber', epsil_tol); toc; -rmse_diff = (RMSE(Ideal3D(:),u_diff(:))); -fprintf('%s %f \n', 'RMSE error for Diffusion is:', rmse_diff); -figure; imshow(u_diff(:,:,7), [0 1]); title('Diffusion denoised volume (CPU)'); -%% -% fprintf('Denoise a volume using Nonlinear-Diffusion model (GPU) \n'); -% iter_diff = 300; % number of diffusion iterations -% lambda_regDiff = 0.025; % regularisation for the diffusivity -% sigmaPar = 0.015; % edge-preserving parameter -% tau_param = 0.025; % time-marching constant -% tic; u_diff_g = NonlDiff_GPU(single(vol3D), lambda_regDiff, sigmaPar, iter_diff, tau_param, 'Huber', epsil_tol); toc; -% rmse_diff = (RMSE(Ideal3D(:),u_diff_g(:))); -% fprintf('%s %f \n', 'RMSE error for Diffusion is:', rmse_diff); -% figure; imshow(u_diff_g(:,:,7), [0 1]); title('Diffusion denoised volume (GPU)'); -%% -fprintf('Denoise using Fourth-order anisotropic diffusion model (CPU) \n'); -iter_diff = 300; % number of diffusion iterations -lambda_regDiff = 3.5; % regularisation for the diffusivity -sigmaPar = 0.02; % edge-preserving parameter -tau_param = 0.0015; % time-marching constant -epsil_tol = 0.0; % tolerance -tic; u_diff4 = Diffusion_4thO(single(vol3D), lambda_regDiff, sigmaPar, iter_diff, tau_param, epsil_tol); toc; -rmse_diff4 = (RMSE(Ideal3D(:),u_diff4(:))); -fprintf('%s %f \n', 'RMSE error for Anis.Diff of 4th order is:', rmse_diff4); -figure; imshow(u_diff4(:,:,7), [0 1]); title('Diffusion 4thO denoised volume (CPU)'); -%% -% fprintf('Denoise using Fourth-order anisotropic diffusion model (GPU) \n'); -% iter_diff = 300; % number of diffusion iterations -% lambda_regDiff = 3.5; % regularisation for the diffusivity -% sigmaPar = 0.02; % edge-preserving parameter -% tau_param = 0.0015; % time-marching constant -% tic; u_diff4_g = Diffusion_4thO_GPU(single(vol3D), lambda_regDiff, sigmaPar, iter_diff, tau_param, epsil_tol); toc; -% rmse_diff4 = (RMSE(Ideal3D(:),u_diff4_g(:))); -% fprintf('%s %f \n', 'RMSE error for Anis.Diff of 4th order is:', rmse_diff4); -% figure; imshow(u_diff4_g(:,:,7), [0 1]); title('Diffusion 4thO denoised volume (GPU)'); -%% -fprintf('Denoise using the TGV model (CPU) \n'); -lambda_TGV = 0.03; % regularisation parameter -alpha1 = 1.0; % parameter to control the first-order term -alpha0 = 2.0; % parameter to control the second-order term -L2 = 12.0; % convergence parameter -iter_TGV = 500; % number of Primal-Dual iterations for TGV -epsil_tol = 0.0; % tolerance -tic; u_tgv = TGV(single(vol3D), lambda_TGV, alpha1, alpha0, iter_TGV, L2, epsil_tol); toc; -rmseTGV = RMSE(Ideal3D(:),u_tgv(:)); -fprintf('%s %f \n', 'RMSE error for TGV is:', rmseTGV); -figure; imshow(u_tgv(:,:,3), [0 1]); title('TGV denoised volume (CPU)'); -%% -% fprintf('Denoise using the TGV model (GPU) \n'); -% lambda_TGV = 0.03; % regularisation parameter -% alpha1 = 1.0; % parameter to control the first-order term -% alpha0 = 2.0; % parameter to control the second-order term -% iter_TGV = 500; % number of Primal-Dual iterations for TGV -% tic; u_tgv_gpu = TGV_GPU(single(vol3D), lambda_TGV, alpha1, alpha0, iter_TGV, L2, epsil_tol); toc; -% rmseTGV = RMSE(Ideal3D(:),u_tgv_gpu(:)); -% fprintf('%s %f \n', 'RMSE error for TGV is:', rmseTGV); -% figure; imshow(u_tgv_gpu(:,:,3), [0 1]); title('TGV denoised volume (GPU)'); -%% -%>>>>>>>>>>>>>> MULTI-CHANNEL priors <<<<<<<<<<<<<<< % -fprintf('Denoise a volume using the FGP-dTV model (CPU) \n'); - -% create another volume (reference) with slightly less amount of noise -vol3D_ref = zeros(N,N,slices, 'single'); -for i = 1:slices -vol3D_ref(:,:,i) = Im + .01*randn(size(Im)); -end -vol3D_ref(vol3D_ref < 0) = 0; -% vol3D_ref = zeros(size(Im),'single'); % pass zero reference (dTV -> TV) - -iter_fgp = 300; % number of FGP iterations -epsil_tol = 0.0; % tolerance -eta = 0.2; % Reference image gradient smoothing constant -tic; u_fgp_dtv = FGP_dTV(single(vol3D), single(vol3D_ref), lambda_reg, iter_fgp, epsil_tol, eta); toc; -figure; imshow(u_fgp_dtv(:,:,7), [0 1]); title('FGP-dTV denoised volume (CPU)'); -%% -fprintf('Denoise a volume using the FGP-dTV model (GPU) \n'); - -% create another volume (reference) with slightly less amount of noise -vol3D_ref = zeros(N,N,slices, 'single'); -for i = 1:slices -vol3D_ref(:,:,i) = Im + .01*randn(size(Im)); -end -vol3D_ref(vol3D_ref < 0) = 0; -% vol3D_ref = zeros(size(Im),'single'); % pass zero reference (dTV -> TV) - -iter_fgp = 300; % number of FGP iterations -epsil_tol = 0.0; % tolerance -eta = 0.2; % Reference image gradient smoothing constant -tic; u_fgp_dtv_g = FGP_dTV_GPU(single(vol3D), single(vol3D_ref), lambda_reg, iter_fgp, epsil_tol, eta); toc; -figure; imshow(u_fgp_dtv_g(:,:,7), [0 1]); title('FGP-dTV denoised volume (GPU)'); -%% diff --git a/demos/demoMatlab_denoise.m b/demos/demoMatlab_denoise.m deleted file mode 100644 index 9d89138..0000000 --- a/demos/demoMatlab_denoise.m +++ /dev/null @@ -1,188 +0,0 @@ -% Image (2D) denoising demo using CCPi-RGL -clear; close all -fsep = '/'; - -Path1 = sprintf(['..' fsep 'src' fsep 'Matlab' fsep 'mex_compile' fsep 'installed'], 1i); -Path2 = sprintf(['data' fsep], 1i); -Path3 = sprintf(['..' fsep 'src' fsep 'Matlab' fsep 'supp'], 1i); -addpath(Path1); -addpath(Path2); -addpath(Path3); - -Im = double(imread('lena_gray_512.tif'))/255; % loading image -u0 = Im + .05*randn(size(Im)); u0(u0 < 0) = 0; -figure; imshow(u0, [0 1]); title('Noisy image'); -%% -fprintf('Denoise using the ROF-TV model (CPU) \n'); -lambda_reg = 0.03; % regularsation parameter for all methods -iter_rof = 2000; % number of ROF iterations -tau_rof = 0.01; % time-marching constant -epsil_tol = 0.0; % tolerance / 1.0e-06 -tic; [u_rof,infovec] = ROF_TV(single(u0), lambda_reg, iter_rof, tau_rof, epsil_tol); toc; -energyfunc_val_rof = TV_energy(single(u_rof),single(u0),lambda_reg, 1); % get energy function value -rmseROF = (RMSE(u_rof(:),Im(:))); -fprintf('%s %f \n', 'RMSE error for ROF-TV is:', rmseROF); -[ssimval] = ssim(u_rof*255,single(Im)*255); -fprintf('%s %f \n', 'MSSIM error for ROF-TV is:', ssimval); -figure; imshow(u_rof, [0 1]); title('ROF-TV denoised image (CPU)'); -%% -%fprintf('Denoise using the ROF-TV model (GPU) \n'); -%tic; [u_rofG,infovec] = ROF_TV_GPU(single(u0), lambda_reg, iter_rof, tau_rof, epsil_tol); toc; -%figure; imshow(u_rofG, [0 1]); title('ROF-TV denoised image (GPU)'); -%% -fprintf('Denoise using the FGP-TV model (CPU) \n'); -lambda_reg = 0.03; -iter_fgp = 500; % number of FGP iterations -epsil_tol = 0.0; % tolerance -tic; [u_fgp,infovec] = FGP_TV(single(u0), lambda_reg, iter_fgp, epsil_tol); toc; -energyfunc_val_fgp = TV_energy(single(u_fgp),single(u0),lambda_reg, 1); % get energy function value -rmseFGP = (RMSE(u_fgp(:),Im(:))); -fprintf('%s %f \n', 'RMSE error for FGP-TV is:', rmseFGP); -[ssimval] = ssim(u_fgp*255,single(Im)*255); -fprintf('%s %f \n', 'MSSIM error for FGP-TV is:', ssimval); -figure; imshow(u_fgp, [0 1]); title('FGP-TV denoised image (CPU)'); -%% -% fprintf('Denoise using the FGP-TV model (GPU) \n'); -% tic; u_fgpG = FGP_TV_GPU(single(u0), lambda_reg, iter_fgp, epsil_tol); toc; -% figure; imshow(u_fgpG, [0 1]); title('FGP-TV denoised image (GPU)'); -%% -fprintf('Denoise using the SB-TV model (CPU) \n'); -lambda_reg = 0.03; -iter_sb = 200; % number of SB iterations -epsil_tol = 0.0; % tolerance -tic; [u_sb,infovec] = SB_TV(single(u0), lambda_reg, iter_sb, epsil_tol); toc; -energyfunc_val_sb = TV_energy(single(u_sb),single(u0),lambda_reg, 1); % get energy function value -rmseSB = (RMSE(u_sb(:),Im(:))); -fprintf('%s %f \n', 'RMSE error for SB-TV is:', rmseSB); -[ssimval] = ssim(u_sb*255,single(Im)*255); -fprintf('%s %f \n', 'MSSIM error for SB-TV is:', ssimval); -figure; imshow(u_sb, [0 1]); title('SB-TV denoised image (CPU)'); -%% -% fprintf('Denoise using the SB-TV model (GPU) \n'); -% tic; u_sbG = SB_TV_GPU(single(u0), lambda_reg, iter_sb, epsil_tol); toc; -% figure; imshow(u_sbG, [0 1]); title('SB-TV denoised image (GPU)'); -%% -fprintf('Denoise using Nonlinear-Diffusion model (CPU) \n'); -iter_diff = 450; % number of diffusion iterations -lambda_regDiff = 0.025; % regularisation for the diffusivity -sigmaPar = 0.015; % edge-preserving parameter -tau_param = 0.02; % time-marching constant -epsil_tol = 0.0; % tolerance -tic; [u_diff,infovec] = NonlDiff(single(u0), lambda_regDiff, sigmaPar, iter_diff, tau_param, 'Huber', epsil_tol); toc; -rmseDiffus = (RMSE(u_diff(:),Im(:))); -fprintf('%s %f \n', 'RMSE error for Nonlinear Diffusion is:', rmseDiffus); -[ssimval] = ssim(u_diff*255,single(Im)*255); -fprintf('%s %f \n', 'MSSIM error for NDF is:', ssimval); -figure; imshow(u_diff, [0 1]); title('Diffusion denoised image (CPU)'); -%% -%fprintf('Denoise using Nonlinear-Diffusion model (GPU) \n'); -%tic; u_diff_g = NonlDiff_GPU(single(u0), lambda_regDiff, sigmaPar, iter_diff, tau_param, 'Huber', epsil_tol); toc; -%figure; imshow(u_diff_g, [0 1]); title('Diffusion denoised image (GPU)'); -%% -fprintf('Denoise using the TGV model (CPU) \n'); -lambda_TGV = 0.035; % regularisation parameter -alpha1 = 1.0; % parameter to control the first-order term -alpha0 = 2.0; % parameter to control the second-order term -L2 = 12.0; % convergence parameter -iter_TGV = 1200; % number of Primal-Dual iterations for TGV -epsil_tol = 0.0; % tolerance -tic; [u_tgv,infovec] = TGV(single(u0), lambda_TGV, alpha1, alpha0, iter_TGV, L2, epsil_tol); toc; -figure; imshow(u_tgv, [0 1]); title('TGV denoised image (CPU)'); -rmseTGV = (RMSE(u_tgv(:),Im(:))); -fprintf('%s %f \n', 'RMSE error for TGV is:', rmseTGV); -[ssimval] = ssim(u_tgv*255,single(Im)*255); -fprintf('%s %f \n', 'MSSIM error for TGV is:', ssimval); -%% -% fprintf('Denoise using the TGV model (GPU) \n'); -% tic; u_tgv_gpu = TGV_GPU(single(u0), lambda_TGV, alpha1, alpha0, iter_TGV, L2, epsil_tol); toc; -% figure; imshow(u_tgv_gpu, [0 1]); title('TGV denoised image (GPU)'); -%% -fprintf('Denoise using the ROF-LLT model (CPU) \n'); -lambda_ROF = 0.02; % ROF regularisation parameter -lambda_LLT = 0.015; % LLT regularisation parameter -iter_LLT = 2000; % iterations -tau_rof_llt = 0.01; % time-marching constant -epsil_tol = 0.0; % tolerance -tic; [u_rof_llt,infovec] = LLT_ROF(single(u0), lambda_ROF, lambda_LLT, iter_LLT, tau_rof_llt,epsil_tol); toc; -rmseROFLLT = (RMSE(u_rof_llt(:),Im(:))); -fprintf('%s %f \n', 'RMSE error for TGV is:', rmseROFLLT); -[ssimval] = ssim(u_rof_llt*255,single(Im)*255); -fprintf('%s %f \n', 'MSSIM error for ROFLLT is:', ssimval); -figure; imshow(u_rof_llt, [0 1]); title('ROF-LLT denoised image (CPU)'); -%% -% fprintf('Denoise using the ROF-LLT model (GPU) \n'); -% tic; u_rof_llt_g = LLT_ROF_GPU(single(u0), lambda_ROF, lambda_LLT, iter_LLT, tau_rof_llt, epsil_tol); toc; -% figure; imshow(u_rof_llt_g, [0 1]); title('ROF-LLT denoised image (GPU)'); -%% -fprintf('Denoise using Fourth-order anisotropic diffusion model (CPU) \n'); -iter_diff = 800; % number of diffusion iterations -lambda_regDiff = 3; % regularisation for the diffusivity -sigmaPar = 0.03; % edge-preserving parameter -tau_param = 0.0025; % time-marching constant -epsil_tol = 0.0; % tolerance -tic; [u_diff4,infovec] = Diffusion_4thO(single(u0), lambda_regDiff, sigmaPar, iter_diff, tau_param, epsil_tol); toc; -rmseDiffHO = (RMSE(u_diff4(:),Im(:))); -fprintf('%s %f \n', 'RMSE error for Fourth-order anisotropic diffusion is:', rmseDiffHO); -[ssimval] = ssim(u_diff4*255,single(Im)*255); -fprintf('%s %f \n', 'MSSIM error for DIFF4th is:', ssimval); -figure; imshow(u_diff4, [0 1]); title('Diffusion 4thO denoised image (CPU)'); -%% -%fprintf('Denoise using Fourth-order anisotropic diffusion model (GPU) \n'); -%tic; u_diff4_g = Diffusion_4thO_GPU(single(u0), lambda_regDiff, sigmaPar, iter_diff, tau_param); toc; -%figure; imshow(u_diff4_g, [0 1]); title('Diffusion 4thO denoised image (GPU)'); -%% -fprintf('Weights pre-calculation for Non-local TV (takes time on CPU) \n'); -SearchingWindow = 7; -PatchWindow = 2; -NeighboursNumber = 20; % the number of neibours to include -h = 0.23; % edge related parameter for NLM -tic; [H_i, H_j, Weights] = PatchSelect(single(u0), SearchingWindow, PatchWindow, NeighboursNumber, h); toc; -%% -fprintf('Denoise using Non-local Total Variation (CPU) \n'); -iter_nltv = 3; % number of nltv iterations -lambda_nltv = 0.055; % regularisation parameter for nltv -tic; u_nltv = Nonlocal_TV(single(u0), H_i, H_j, 0, Weights, lambda_nltv, iter_nltv); toc; -rmse_nltv = (RMSE(u_nltv(:),Im(:))); -fprintf('%s %f \n', 'RMSE error for Non-local Total Variation is:', rmse_nltv); -[ssimval] = ssim(u_nltv*255,single(Im)*255); -fprintf('%s %f \n', 'MSSIM error for NLTV is:', ssimval); -figure; imagesc(u_nltv, [0 1]); colormap(gray); daspect([1 1 1]); title('Non-local Total Variation denoised image (CPU)'); -%% -%>>>>>>>>>>>>>> MULTI-CHANNEL priors <<<<<<<<<<<<<<< % - -fprintf('Denoise using the FGP-dTV model (CPU) \n'); -% create another image (reference) with slightly less amount of noise -u_ref = Im + .01*randn(size(Im)); u_ref(u_ref < 0) = 0; -% u_ref = zeros(size(Im),'single'); % pass zero reference (dTV -> TV) - -lambda_reg = 0.04; -iter_fgp = 1000; % number of FGP iterations -epsil_tol = 0.0; % tolerance -eta = 0.2; % Reference image gradient smoothing constant -tic; [u_fgp_dtv,infovec] = FGP_dTV(single(u0), single(u_ref), lambda_reg, iter_fgp, epsil_tol, eta); toc; -rmse_dTV= (RMSE(u_fgp_dtv(:),Im(:))); -fprintf('%s %f \n', 'RMSE error for Directional Total Variation (dTV) is:', rmse_dTV); -figure; imshow(u_fgp_dtv, [0 1]); title('FGP-dTV denoised image (CPU)'); -%% -% fprintf('Denoise using the FGP-dTV model (GPU) \n'); -% % create another image (reference) with slightly less amount of noise -% u_ref = Im + .01*randn(size(Im)); u_ref(u_ref < 0) = 0; -% % u_ref = zeros(size(Im),'single'); % pass zero reference (dTV -> TV) -% -% iter_fgp = 1000; % number of FGP iterations -% epsil_tol = 1.0e-06; % tolerance -% eta = 0.2; % Reference image gradient smoothing constant -% tic; u_fgp_dtvG = FGP_dTV_GPU(single(u0), single(u_ref), lambda_reg, iter_fgp, epsil_tol, eta); toc; -% figure; imshow(u_fgp_dtvG, [0 1]); title('FGP-dTV denoised image (GPU)'); -%% -fprintf('Denoise using the TNV prior (CPU) \n'); -slices = 5; N = 512; -vol3D = zeros(N,N,slices, 'single'); -for i = 1:slices -vol3D(:,:,i) = Im + .05*randn(size(Im)); -end -vol3D(vol3D < 0) = 0; - -iter_tnv = 200; % number of TNV iterations -tic; u_tnv = TNV(single(vol3D), lambda_reg, iter_tnv); toc; -figure; imshow(u_tnv(:,:,3), [0 1]); title('TNV denoised stack of channels (CPU)'); diff --git a/demos/demoMatlab_inpaint.m b/demos/demoMatlab_inpaint.m deleted file mode 100644 index a85f2b9..0000000 --- a/demos/demoMatlab_inpaint.m +++ /dev/null @@ -1,35 +0,0 @@ -% Image (2D) inpainting demo using CCPi-RGL -clear; close all -Path1 = sprintf(['..' filesep 'src' filesep 'Matlab' filesep 'mex_compile' filesep 'installed'], 1i); -Path2 = sprintf(['data' filesep], 1i); -addpath(Path1); -addpath(Path2); - -load('SinoInpaint.mat'); -Sinogram = Sinogram./max(Sinogram(:)); -Sino_mask = Sinogram.*(1-single(Mask)); -figure; -subplot(1,2,1); imshow(Sino_mask, [0 1]); title('Missing data sinogram'); -subplot(1,2,2); imshow(Mask, [0 1]); title('Mask'); -%% -fprintf('Inpaint using Linear-Diffusion model (CPU) \n'); -iter_diff = 5000; % number of diffusion iterations -lambda_regDiff = 6000; % regularisation for the diffusivity -sigmaPar = 0.0; % edge-preserving parameter -tau_param = 0.000075; % time-marching constant -tic; u_diff = NonlDiff_Inp(single(Sino_mask), Mask, lambda_regDiff, sigmaPar, iter_diff, tau_param); toc; -figure; imshow(u_diff, [0 1]); title('Linear-Diffusion inpainted sinogram (CPU)'); -%% -fprintf('Inpaint using Nonlinear-Diffusion model (CPU) \n'); -iter_diff = 1500; % number of diffusion iterations -lambda_regDiff = 80; % regularisation for the diffusivity -sigmaPar = 0.00009; % edge-preserving parameter -tau_param = 0.000008; % time-marching constant -tic; u_diff = NonlDiff_Inp(single(Sino_mask), Mask, lambda_regDiff, sigmaPar, iter_diff, tau_param, 'Huber'); toc; -figure; imshow(u_diff, [0 1]); title('Non-Linear Diffusion inpainted sinogram (CPU)'); -%% -fprintf('Inpaint using Nonlocal Vertical Marching model (CPU) \n'); -Increment = 1; % linear increment for the searching window -tic; [u_nom,maskupd] = NonlocalMarching_Inpaint(single(Sino_mask), Mask, Increment); toc; -figure; imshow(u_nom, [0 1]); title('NVM inpainted sinogram (CPU)'); -%% \ No newline at end of file diff --git a/src/Core/regularisers_CPU/MASK_merge_core.c b/src/Core/regularisers_CPU/MASK_merge_core.c deleted file mode 100644 index d77c812..0000000 --- a/src/Core/regularisers_CPU/MASK_merge_core.c +++ /dev/null @@ -1,262 +0,0 @@ -/* - * This work is part of the Core Imaging Library developed by - * Visual Analytics and Imaging System Group of the Science Technology - * Facilities Council, STFC - * - * Copyright 2019 Daniil Kazantsev - * Copyright 2019 Srikanth Nagella, Edoardo Pasca - * - * Licensed under the Apache License, Version 2.0 (the "License"); - * you may not use this file except in compliance with the License. - * You may obtain a copy of the License at - * http://www.apache.org/licenses/LICENSE-2.0 - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ - -#include "MASK_merge_core.h" -#include "utils.h" - -/* A method to ensure connectivity within regions of the segmented image/volume. Here we assume - * that the MASK has been obtained using some classification/segmentation method such as k-means or gaussian - * mixture. Some pixels/voxels have been misclassified and we check the spatial dependences - * and correct the mask. We check the connectivity using the bresenham line algorithm within the non-local window - * surrounding the pixel of interest. - * - * Input Parameters: - * 1. MASK [0:255], the result of some classification algorithm (information-based preferably) - * 2. The list of classes (e.g. [3,4]) to apply the method. The given order matters. - * 3. The total number of classes in the MASK. - * 4. The size of the Correction Window inside which the method works. - - * Output: - * 1. MASK_upd - the UPDATED MASK where some regions have been corrected (merged) or removed - * 2. CORRECTEDRegions - The array of the same size as MASK where all regions which were - * changed are highlighted and the changes have been counted - */ - -float Mask_merge_main(unsigned char *MASK, unsigned char *MASK_upd, unsigned char *CORRECTEDRegions, unsigned char *SelClassesList, int SelClassesList_length, int classesNumb, int CorrectionWindow, int dimX, int dimY, int dimZ) -{ - long i,j,l; - int counterG, switcher; - long DimTotal; - unsigned char *MASK_temp, *ClassesList, CurrClass, temp; - DimTotal = (long)(dimX*dimY*dimZ); - - /* defines the list for all classes in the mask */ - ClassesList = (unsigned char*) calloc (classesNumb,sizeof(unsigned char)); - - /* find which classes (values) are present in the segmented data */ - CurrClass = MASK[0]; ClassesList[0]= MASK[0]; counterG = 1; - for(i=0; iHIGH the obtained values (classes) */ - for(i=0; i ClassesList[j+1]) { - temp = ClassesList[j+1]; - ClassesList[j+1] = ClassesList[j]; - ClassesList[j] = temp; - }}} - - MASK_temp = (unsigned char*) calloc (DimTotal,sizeof(unsigned char)); - - /* copy given MASK to MASK_upd*/ - copyIm_unchar(MASK, MASK_upd, (long)(dimX), (long)(dimY), (long)(dimZ)); - - if (dimZ == 1) { - /********************** PERFORM 2D MASK PROCESSING ************************/ - #pragma omp parallel for shared(MASK,MASK_upd) private(i,j) - for(i=0; i continue */ - if (MASK_temp[j*dimX+i] == MASK[j*dimX+i]) { - /* !One needs to work with a specific class to avoid overlaps! It is - crucial to establish relevant classes first (given as an input in SelClassesList) */ - if (MASK_temp[j*dimX+i] == ClassesList[SelClassesList[l]]) { - /* i = 258; j = 165; */ - Mask_update2D(MASK_temp, MASK_upd, CORRECTEDRegions, i, j, CorrectionWindow, (long)(dimX), (long)(dimY)); - }} - }} - /* copy the updated mask */ - copyIm_unchar(MASK_upd, MASK_temp, (long)(dimX), (long)(dimY), (long)(dimZ)); - } - } - else { - /********************** PERFORM 3D MASK PROCESSING ************************/ - - - } - - free(MASK_temp); - return 0; -} - - -/********************************************************************/ -/***************************2D Functions*****************************/ -/********************************************************************/ -float OutiersRemoval2D(unsigned char *MASK, unsigned char *MASK_upd, long i, long j, long dimX, long dimY) -{ - /*if the ROI pixel does not belong to the surrondings, turn it into the surronding*/ - long i_m, j_m, i1, j1, counter; - counter = 0; - for(i_m=-1; i_m<=1; i_m++) { - for(j_m=-1; j_m<=1; j_m++) { - i1 = i+i_m; - j1 = j+j_m; - if (((i1 >= 0) && (i1 < dimX)) && ((j1 >= 0) && (j1 < dimY))) { - if (MASK[j*dimX+i] != MASK[j1*dimX+i1]) counter++; - } - }} - if (counter >= 8) MASK_upd[j*dimX+i] = MASK[j1*dimX+i1]; - return *MASK_upd; -} - -float Mask_update2D(unsigned char *MASK_temp, unsigned char *MASK_upd, unsigned char *CORRECTEDRegions, long i, long j, int CorrectionWindow, long dimX, long dimY) -{ - long i_m, j_m, i1, j1, CounterOtherClass; - - /* STEP2: in a larger neighbourhood check that the other class is present */ - CounterOtherClass = 0; - for(i_m=-CorrectionWindow; i_m<=CorrectionWindow; i_m++) { - for(j_m=-CorrectionWindow; j_m<=CorrectionWindow; j_m++) { - i1 = i+i_m; - j1 = j+j_m; - if (((i1 >= 0) && (i1 < dimX)) && ((j1 >= 0) && (j1 < dimY))) { - if (MASK_temp[j*dimX+i] != MASK_temp[j1*dimX+i1]) CounterOtherClass++; - } - }} - if (CounterOtherClass > 0) { - /* the other class is present in the vicinity of CorrectionWindow, continue to STEP 3 */ - /* - STEP 3: Loop through all neighbours in CorrectionWindow and check the spatial connection. - Meaning that we're instrested if there are any classes between points A and B that - does not belong to A and B (A,B \in C) - */ - for(i_m=-CorrectionWindow; i_m<=CorrectionWindow; i_m++) { - for(j_m=-CorrectionWindow; j_m<=CorrectionWindow; j_m++) { - i1 = i+i_m; - j1 = j+j_m; - if (((i1 >= 0) && (i1 < dimX)) && ((j1 >= 0) && (j1 < dimY))) { - if (MASK_temp[j*dimX+i] == MASK_temp[j1*dimX+i1]) { - /* A and B points belong to the same class, do STEP 4*/ - /* STEP 4: Run Bresenham line algorithm between A and B points - and convert all points on the way to the class of A/B. */ - bresenham2D(i, j, i1, j1, MASK_temp, MASK_upd, CORRECTEDRegions, (long)(dimX), (long)(dimY)); - } - } - }} - } - return 0; -} -int bresenham2D(int i, int j, int i1, int j1, unsigned char *MASK, unsigned char *MASK_upd, unsigned char *CORRECTEDRegions, long dimX, long dimY) -{ - int n; - int x[] = {i, i1}; - int y[] = {j, j1}; - int steep = (fabs(y[1]-y[0]) > fabs(x[1]-x[0])); - int ystep = 0; - - //printf("[%i][%i][%i][%i]\n", x[1], y[1], steep, kk) ; - //if (steep == 1) {swap(x[0],y[0]); swap(x[1],y[1]);} - - if (steep == 1) { - // swaping - int a, b; - - a = x[0]; - b = y[0]; - x[0] = b; - y[0] = a; - - a = x[1]; - b = y[1]; - x[1] = b; - y[1] = a; - } - - if (x[0] > x[1]) { - int a, b; - a = x[0]; - b = x[1]; - x[0] = b; - x[1] = a; - - a = y[0]; - b = y[1]; - y[0] = b; - y[1] = a; - } //(x[0] > x[1]) - - int delx = x[1]-x[0]; - int dely = fabs(y[1]-y[0]); - int error = 0; - int x_n = x[0]; - int y_n = y[0]; - if (y[0] < y[1]) {ystep = 1;} - else {ystep = -1;} - - for(n = 0; n= delx) { - y_n = y_n + ystep; - error = error - delx; - } // (2*error >= delx) - //printf("[%i][%i][%i]\n", X_new[n], Y_new[n], n) ; - } // for(int n = 0; n -#include -#include -#include -#include "omp.h" -#include "utils.h" -#include "CCPiDefines.h" - -/* A method to ensure connectivity within regions of the segmented image/volume. Here we assume - * that the MASK has been obtained using some classification/segmentation method such as k-means or gaussian - * mixture. Some pixels/voxels have been misclassified and we check the spatial dependences - * and correct the mask. We check the connectivity using the bresenham line algorithm within the non-local window - * surrounding the pixel of interest. - * - * Input Parameters: - * 1. MASK [0:255], the result of some classification algorithm (information-based preferably) - * 2. The list of classes (e.g. [3,4]) to apply the method. The given order matters. - * 3. The total number of classes in the MASK. - * 4. The size of the Correction Window inside which the method works. - - * Output: - * 1. MASK_upd - the UPDATED MASK where some regions have been corrected (merged) or removed - * 2. CORRECTEDRegions - The array of the same size as MASK where all regions which were - * changed are highlighted and the changes have been counted - */ - - -#ifdef __cplusplus -extern "C" { -#endif -CCPI_EXPORT float Mask_merge_main(unsigned char *MASK, unsigned char *MASK_upd, unsigned char *CORRECTEDRegions, unsigned char *SelClassesList, int SelClassesList_length, int classesNumb, int CorrectionWindow, int dimX, int dimY, int dimZ); -CCPI_EXPORT float OutiersRemoval2D(unsigned char *MASK, unsigned char *MASK_upd, long i, long j, long dimX, long dimY); -CCPI_EXPORT float Mask_update2D(unsigned char *MASK_temp, unsigned char *MASK_upd, unsigned char *CORRECTEDRegions, long i, long j, int CorrectionWindow, long dimX, long dimY); -CCPI_EXPORT int bresenham2D(int i, int j, int i1, int j1, unsigned char *MASK, unsigned char *MASK_upd, unsigned char *CORRECTEDRegions, long dimX, long dimY); -#ifdef __cplusplus -} -#endif diff --git a/src/Matlab/mex_compile/compileCPU_mex_Linux.m b/src/Matlab/mex_compile/compileCPU_mex_Linux.m index 19fb1a5..2ed7ea6 100644 --- a/src/Matlab/mex_compile/compileCPU_mex_Linux.m +++ b/src/Matlab/mex_compile/compileCPU_mex_Linux.m @@ -77,5 +77,5 @@ delete PatchSelect_core* Nonlocal_TV_core* delete Diffusion_Inpaint_core* NonlocalMarching_Inpaint_core* fprintf('%s \n', '<<<<<<< Regularisers successfully compiled! >>>>>>>'); -pathA2 = sprintf(['..' fsep '..' fsep '..' fsep '..' fsep 'demos'], 1i); +pathA2 = sprintf(['..' fsep '..' fsep '..' fsep '..' fsep 'demos' fsep 'Matlab_demos'], 1i); cd(pathA2); \ No newline at end of file diff --git a/src/Matlab/mex_compile/compileCPU_mex_WINDOWS.m b/src/Matlab/mex_compile/compileCPU_mex_WINDOWS.m index 3a9e2af..7d0917f 100644 --- a/src/Matlab/mex_compile/compileCPU_mex_WINDOWS.m +++ b/src/Matlab/mex_compile/compileCPU_mex_WINDOWS.m @@ -7,7 +7,11 @@ % Here I present two ways how software can be compiled, if you have some % other suggestions/remarks please contact me at dkazanc@hotmail.com % >>>>>>>>>>>>>>>>>>>>>>>>>>>>> - +% cuda related messgaes can be solved in Linux: +%export LD_LIBRARY_PATH="$LD_LIBRARY_PATH:/usr/local/cuda/lib64" +%export CUDA_HOME=/usr/local/cuda +%export LD_PRELOAD=/home/algol/anaconda3/lib/libstdc++.so.6.0.24 +% >> fsep = '/'; pathcopyFrom = sprintf(['..' fsep '..' fsep 'Core' fsep 'regularisers_CPU'], 1i); diff --git a/src/Matlab/mex_compile/compileGPU_mex.m b/src/Matlab/mex_compile/compileGPU_mex.m index 7e15233..56fcd38 100644 --- a/src/Matlab/mex_compile/compileGPU_mex.m +++ b/src/Matlab/mex_compile/compileGPU_mex.m @@ -66,10 +66,14 @@ fprintf('%s \n', 'Compiling ROF-LLT...'); mex -g -I/usr/local/cuda-10.0/include -L/usr/local/cuda-10.0/lib64 -lcudart -lcufft -lmwgpu LLT_ROF_GPU.cpp LLT_ROF_GPU_core.o movefile('LLT_ROF_GPU.mex*',Pathmove); +fprintf('%s \n', 'Compiling PatchSelect...'); +!/usr/local/cuda/bin/nvcc -O0 -c PatchSelect_GPU_core.cu -Xcompiler -fPIC -I~/SOFT/MATLAB9/extern/include/ +mex -g -I/usr/local/cuda-10.0/include -L/usr/local/cuda-10.0/lib64 -lcudart -lcufft -lmwgpu PatchSelect_GPU.cpp PatchSelect_GPU_core.o +movefile('PatchSelect_GPU.mex*',Pathmove); delete TV_ROF_GPU_core* TV_FGP_GPU_core* TV_SB_GPU_core* dTV_FGP_GPU_core* NonlDiff_GPU_core* Diffus_4thO_GPU_core* TGV_GPU_core* LLT_ROF_GPU_core* CCPiDefines.h -delete PatchSelect_core* Nonlocal_TV_core* shared.h +delete PatchSelect_GPU_core* Nonlocal_TV_core* shared.h fprintf('%s \n', 'All successfully compiled!'); -pathA2 = sprintf(['..' fsep '..' fsep '..' fsep '..' fsep 'demos'], 1i); +pathA2 = sprintf(['..' fsep '..' fsep '..' fsep '..' fsep 'demos' fsep 'Matlab_demos'], 1i); cd(pathA2); diff --git a/src/Matlab/mex_compile/regularisers_CPU/ROF_TV.c b/src/Matlab/mex_compile/regularisers_CPU/ROF_TV.c index ffe7b91..d6bdaa6 100644 --- a/src/Matlab/mex_compile/regularisers_CPU/ROF_TV.c +++ b/src/Matlab/mex_compile/regularisers_CPU/ROF_TV.c @@ -26,7 +26,7 @@ * * Input Parameters: * 1. Noisy image/volume [REQUIRED] - * 2. lambda - regularization parameter [REQUIRED] + * 2. lambda - regularization parameter [REQUIRED] scalar or the same size as 1 * 3. Number of iterations, for explicit scheme >= 150 is recommended [REQUIRED] * 4. tau - marching step for explicit scheme, ~1 is recommended [REQUIRED] * 5. eplsilon: tolerance constant [REQUIRED] @@ -38,7 +38,7 @@ * This function is based on the paper by * [1] Rudin, Osher, Fatemi, "Nonlinear Total Variation based noise removal algorithms" * - * D. Kazantsev, 2016-18 + * D. Kazantsev, 2016-19 */ void mexFunction( @@ -49,18 +49,29 @@ void mexFunction( int number_of_dims, iter_numb; mwSize dimX, dimY, dimZ; const mwSize *dim_array_i; - float *Input, *Output=NULL, lambda, tau, epsil; + float *Input, *Output=NULL, lambda_scalar, tau, epsil; float *infovec=NULL; + float *lambda; dim_array_i = mxGetDimensions(prhs[0]); number_of_dims = mxGetNumberOfDimensions(prhs[0]); /*Handling Matlab input data*/ Input = (float *) mxGetData(prhs[0]); - lambda = (float) mxGetScalar(prhs[1]); /* regularization parameter */ + /*special check to find the input scalar or an array*/ + int mrows = mxGetM(prhs[1]); + int ncols = mxGetN(prhs[1]); + if (mrows==1 && ncols==1) { + lambda = (float*) calloc (1 ,sizeof(float)); + lambda_scalar = (float) mxGetScalar(prhs[1]); /* regularization parameter */ + lambda[0] = lambda_scalar; + } + else { + lambda = (float *) mxGetData(prhs[1]); /* regularization parameter */ + } iter_numb = (int) mxGetScalar(prhs[2]); /* iterations number */ tau = (float) mxGetScalar(prhs[3]); /* marching step parameter */ - epsil = (float) mxGetScalar(prhs[4]); /* tolerance */ + epsil = (float) mxGetScalar(prhs[4]); /* tolerance */ if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); } if(nrhs != 5) mexErrMsgTxt("Four inputs reqired: Image(2D,3D), regularization parameter, iterations number, marching step constant, tolerance"); @@ -80,6 +91,11 @@ void mexFunction( mwSize vecdim[1]; vecdim[0] = 2; infovec = (float*)mxGetPr(plhs[1] = mxCreateNumericArray(1, vecdim, mxSINGLE_CLASS, mxREAL)); - - TV_ROF_CPU_main(Input, Output, infovec, lambda, iter_numb, tau, epsil, dimX, dimY, dimZ); + + if (mrows==1 && ncols==1) { + TV_ROF_CPU_main(Input, Output, infovec, lambda, 0, iter_numb, tau, epsil, dimX, dimY, dimZ); + free(lambda); + } + else TV_ROF_CPU_main(Input, Output, infovec, lambda, 1, iter_numb, tau, epsil, dimX, dimY, dimZ); + } diff --git a/src/Matlab/mex_compile/regularisers_GPU/PatchSelect_GPU.cpp b/src/Matlab/mex_compile/regularisers_GPU/PatchSelect_GPU.cpp new file mode 100644 index 0000000..13319fa --- /dev/null +++ b/src/Matlab/mex_compile/regularisers_GPU/PatchSelect_GPU.cpp @@ -0,0 +1,94 @@ +/* + * This work is part of the Core Imaging Library developed by + * Visual Analytics and Imaging System Group of the Science Technology + * Facilities Council, STFC and Diamond Light Source Ltd. + * + * Copyright 2017 Daniil Kazantsev + * Copyright 2017 Srikanth Nagella, Edoardo Pasca + * Copyright 2018 Diamond Light Source Ltd. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * http://www.apache.org/licenses/LICENSE-2.0 + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include "matrix.h" +#include "mex.h" +#include "PatchSelect_GPU_core.h" + +/* CUDA implementation of non-local weight pre-calculation for non-local priors + * Weights and associated indices are stored into pre-allocated arrays and passed + * to the regulariser + * + * + * Input Parameters: + * 1. 2D/3D grayscale image/volume + * 2. Searching window (half-size of the main bigger searching window, e.g. 11) + * 3. Similarity window (half-size of the patch window, e.g. 2) + * 4. The number of neighbours to take (the most prominent after sorting neighbours will be taken) + * 5. noise-related parameter to calculate non-local weights + * + * Output [2D]: + * 1. AR_i - indeces of i neighbours + * 2. AR_j - indeces of j neighbours + * 3. Weights_ij - associated weights + * + * Output [3D]: + * 1. AR_i - indeces of i neighbours + * 2. AR_j - indeces of j neighbours + * 3. AR_k - indeces of j neighbours + * 4. Weights_ijk - associated weights + */ +/**************************************************/ +void mexFunction( + int nlhs, mxArray *plhs[], + int nrhs, const mxArray *prhs[]) +{ + int number_of_dims, SearchWindow, SimilarWin, NumNeighb; + mwSize dimX, dimY, dimZ; + const mwSize *dim_array; + unsigned short *H_i=NULL, *H_j=NULL, *H_k=NULL; + float *A, *Weights = NULL, h; + mwSize dim_array2[3]; /* for 2D data */ + mwSize dim_array3[4]; /* for 3D data */ + + dim_array = mxGetDimensions(prhs[0]); + number_of_dims = mxGetNumberOfDimensions(prhs[0]); + + /*Handling Matlab input data*/ + A = (float *) mxGetData(prhs[0]); /* a 2D or 3D image/volume */ + SearchWindow = (int) mxGetScalar(prhs[1]); /* Large Searching window */ + SimilarWin = (int) mxGetScalar(prhs[2]); /* Similarity window (patch-search)*/ + NumNeighb = (int) mxGetScalar(prhs[3]); /* the total number of neighbours to take */ + h = (float) mxGetScalar(prhs[4]); /* NLM parameter */ + + dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2]; + dim_array2[0] = dimX; dim_array2[1] = dimY; dim_array2[2] = NumNeighb; /* 2D case */ + dim_array3[0] = dimX; dim_array3[1] = dimY; dim_array3[2] = dimZ; dim_array3[3] = NumNeighb; /* 3D case */ + + /****************2D INPUT ***************/ + if (number_of_dims == 2) { + dimZ = 0; + H_i = (unsigned short*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array2, mxUINT16_CLASS, mxREAL)); + H_j = (unsigned short*)mxGetPr(plhs[1] = mxCreateNumericArray(3, dim_array2, mxUINT16_CLASS, mxREAL)); + Weights = (float*)mxGetPr(plhs[2] = mxCreateNumericArray(3, dim_array2, mxSINGLE_CLASS, mxREAL)); + } + /****************3D INPUT ***************/ + if (number_of_dims == 3) { + /* + H_i = (unsigned short*)mxGetPr(plhs[0] = mxCreateNumericArray(4, dim_array3, mxUINT16_CLASS, mxREAL)); + H_j = (unsigned short*)mxGetPr(plhs[1] = mxCreateNumericArray(4, dim_array3, mxUINT16_CLASS, mxREAL)); + H_k = (unsigned short*)mxGetPr(plhs[2] = mxCreateNumericArray(4, dim_array3, mxUINT16_CLASS, mxREAL)); + Weights = (float*)mxGetPr(plhs[3] = mxCreateNumericArray(4, dim_array3, mxSINGLE_CLASS, mxREAL)); + */ + mexErrMsgTxt("3D version is not yet supported"); + } + + PatchSelect_GPU_main(A, H_i, H_j, Weights, (long)(dimX), (long)(dimY), SearchWindow, SimilarWin, NumNeighb, h); + } diff --git a/src/Matlab/mex_compile/regularisers_GPU/ROF_TV_GPU.cpp b/src/Matlab/mex_compile/regularisers_GPU/ROF_TV_GPU.cpp index d9b7e83..2cd4469 100644 --- a/src/Matlab/mex_compile/regularisers_GPU/ROF_TV_GPU.cpp +++ b/src/Matlab/mex_compile/regularisers_GPU/ROF_TV_GPU.cpp @@ -25,7 +25,7 @@ * * Input Parameters: * 1. Noisy image/volume [REQUIRED] - * 2. lambda - regularization parameter [REQUIRED] + * 2. lambda - regularization parameter [REQUIRED], scalar or the same size as 1 * 3. Number of iterations, for explicit scheme >= 150 is recommended [REQUIRED] * 4. tau - marching step for explicit scheme, ~1 is recommended [REQUIRED] * 5. eplsilon: tolerance constant [REQUIRED] @@ -37,7 +37,7 @@ * This function is based on the paper by * [1] Rudin, Osher, Fatemi, "Nonlinear Total Variation based noise removal algorithms" * - * D. Kazantsev, 2016-18 + * D. Kazantsev, 2016-19 */ void mexFunction( int nlhs, mxArray *plhs[], @@ -47,15 +47,26 @@ void mexFunction( int number_of_dims, iter_numb; mwSize dimX, dimY, dimZ; const mwSize *dim_array_i; - float *Input, *Output=NULL, lambda, tau, epsil; + float *Input, *Output=NULL, lambda_scalar, tau, epsil; float *infovec=NULL; - + float *lambda; + dim_array_i = mxGetDimensions(prhs[0]); number_of_dims = mxGetNumberOfDimensions(prhs[0]); /*Handling Matlab input data*/ Input = (float *) mxGetData(prhs[0]); - lambda = (float) mxGetScalar(prhs[1]); /* regularization parameter */ + /*special check to find the input scalar or an array*/ + int mrows = mxGetM(prhs[1]); + int ncols = mxGetN(prhs[1]); + if (mrows==1 && ncols==1) { + lambda = (float*) calloc (1 ,sizeof(float)); + lambda_scalar = (float) mxGetScalar(prhs[1]); /* regularization parameter */ + lambda[0] = lambda_scalar; + } + else { + lambda = (float *) mxGetData(prhs[1]); /* regularization parameter */ + } iter_numb = (int) mxGetScalar(prhs[2]); /* iterations number */ tau = (float) mxGetScalar(prhs[3]); /* marching step parameter */ epsil = (float) mxGetScalar(prhs[4]); /* tolerance */ @@ -79,5 +90,9 @@ void mexFunction( vecdim[0] = 2; infovec = (float*)mxGetPr(plhs[1] = mxCreateNumericArray(1, vecdim, mxSINGLE_CLASS, mxREAL)); - TV_ROF_GPU_main(Input, Output, infovec, lambda, iter_numb, tau, epsil, dimX, dimY, dimZ); + if (mrows==1 && ncols==1) { + TV_ROF_GPU_main(Input, Output, infovec, lambda, 0, iter_numb, tau, epsil, dimX, dimY, dimZ); + free(lambda); + } + else TV_ROF_GPU_main(Input, Output, infovec, lambda, 1, iter_numb, tau, epsil, dimX, dimY, dimZ); } -- cgit v1.2.3