From 4bcfee09d4b8fc23ea231521b4ceb7aaeecf2811 Mon Sep 17 00:00:00 2001 From: dkazanc Date: Thu, 14 Feb 2019 17:41:00 +0000 Subject: first TGVD version --- Core/regularisers_CPU/TGV_core.c | 202 ++++++++++++++------- Core/regularisers_CPU/TGV_core.h | 2 +- Wrappers/Matlab/demos/demoMatlab_3Ddenoise.m | 15 +- Wrappers/Matlab/demos/demoMatlab_denoise.m | 16 +- Wrappers/Matlab/mex_compile/regularisers_CPU/TGV.c | 29 +-- 5 files changed, 170 insertions(+), 94 deletions(-) diff --git a/Core/regularisers_CPU/TGV_core.c b/Core/regularisers_CPU/TGV_core.c index 964c8a1..01306f3 100644 --- a/Core/regularisers_CPU/TGV_core.c +++ b/Core/regularisers_CPU/TGV_core.c @@ -121,7 +121,7 @@ float TGV_main(float *U0, float *U, float lambda, float alpha1, float alpha0, in DualP_3D(U, V1, V2, V3, P1, P2, P3, (long)(dimX), (long)(dimY), (long)(dimZ), sigma); /*Projection onto convex set for P*/ - ProjP_3D(P1, P2, P3, (long)(dimX), (long)(dimY), alpha1); + ProjP_3D(P1, P2, P3, (long)(dimX), (long)(dimY), (long)(dimZ), alpha1); /* Calculate Dual Variable Q */ DualQ_3D(V1, V2, V3, Q1, Q2, Q3, Q4, Q5, Q6, (long)(dimX), (long)(dimY), (long)(dimZ), sigma); @@ -189,10 +189,9 @@ float ProjP_2D(float *P1, float *P2, long dimX, long dimY, float alpha1) #pragma omp parallel for shared(P1,P2) private(i,j,index,grad_magn) for(i=0; i 1.0) { + index = j*dimX+i; + grad_magn = (sqrtf(pow(P1[index],2) + pow(P2[index],2)))/alpha1; + if (grad_magn > 1.0f) { P1[index] /= grad_magn; P2[index] /= grad_magn; } @@ -207,21 +206,14 @@ float DualQ_2D(float *V1, float *V2, float *Q1, float *Q2, float *Q3, long dimX, #pragma omp parallel for shared(Q1,Q2,Q3,V1,V2) private(i,j,index,q1,q2,q11,q22) for(i=0; i 1.0) { + if (grad_magn > 1.0f) { Q1[index] /= grad_magn; Q2[index] /= grad_magn; Q3[index] /= grad_magn; @@ -279,37 +271,29 @@ float newU(float *U, float *U_old, long dimX, long dimY) float UpdV_2D(float *V1, float *V2, float *P1, float *P2, float *Q1, float *Q2, float *Q3, long dimX, long dimY, float tau) { long i, j, index; - float q1, q11, q2, q22, div1, div2; -#pragma omp parallel for shared(V1,V2,P1,P2,Q1,Q2,Q3) private(i, j, index, q1, q11, q2, q22, div1, div2) + float q1, q3_x, q3_y, q2, div1, div2; +#pragma omp parallel for shared(V1,V2,P1,P2,Q1,Q2,Q3) private(i, j, index, q1, q3_x, q3_y, q2, div1, div2) for(i=0; i 1.0f) { P1[index] /= grad_magn; P2[index] /= grad_magn; @@ -357,39 +338,29 @@ float ProjP_3D(float *P1, float *P2, float *P3, long dimX, long dimY, long dimZ, /*Calculating dual variable Q (using forward differences)*/ float DualQ_3D(float *V1, float *V2, float *V3, float *Q1, float *Q2, float *Q3, float *Q4, float *Q5, float *Q6, long dimX, long dimY, long dimZ, float sigma) { - long i,j,k, index; - float q1, q2, q3, q4, q5, q6, q11, q22, q33, q44, q55, q66; -#pragma omp parallel for shared(Q1,Q2,Q3,Q4,Q5,Q6,V1,V2,V3) private(i,j,k,index,q1,q2,q3,q4,q5,q6,q11,q22,q33,q44,q55,q66) + long i,j,k,index; + float q1, q2, q3, q11, q22, q33, q44, q55, q66; +#pragma omp parallel for shared(Q1,Q2,Q3,Q4,Q5,Q6,V1,V2,V3) private(i,j,k,index,q1,q2,q3,q11,q22,q33,q44,q55,q66) for(i=0; i 1.0f) { + Q1[index] /= grad_magn; + Q2[index] /= grad_magn; + Q3[index] /= grad_magn; + Q4[index] /= grad_magn; + Q5[index] /= grad_magn; + Q6[index] /= grad_magn; + } + }}} + return 1; +} +/* Divergence and projection for P*/ +float DivProjP_3D(float *U, float *U0, float *P1, float *P2, float *P3, long dimX, long dimY, long dimZ, float lambda, float tau) +{ + long i,j,k,index; + float P_v1, P_v2, P_v3, div; +#pragma omp parallel for shared(U,U0,P1,P2,P3) private(i,j,k,index,P_v1,P_v2,P_v3,div) + for(i=0; i>>>>>>>>>>>>> MULTI-CHANNEL priors <<<<<<<<<<<<<<< % fprintf('Denoise a volume using the FGP-dTV model (CPU) \n'); diff --git a/Wrappers/Matlab/demos/demoMatlab_denoise.m b/Wrappers/Matlab/demos/demoMatlab_denoise.m index 3506cca..14d3096 100644 --- a/Wrappers/Matlab/demos/demoMatlab_denoise.m +++ b/Wrappers/Matlab/demos/demoMatlab_denoise.m @@ -60,20 +60,20 @@ figure; imshow(u_sb, [0 1]); title('SB-TV denoised image (CPU)'); % figure; imshow(u_sbG, [0 1]); title('SB-TV denoised image (GPU)'); %% fprintf('Denoise using the TGV model (CPU) \n'); -lambda_TGV = 0.04; % regularisation parameter -alpha1 = 1; % parameter to control the first-order term -alpha0 = 0.7; % parameter to control the second-order term -iter_TGV = 500; % number of Primal-Dual iterations for TGV +lambda_TGV = 0.045; % regularisation parameter +alpha1 = 1.0; % parameter to control the first-order term +alpha0 = 2.0; % parameter to control the second-order term +iter_TGV = 2000; % number of Primal-Dual iterations for TGV tic; u_tgv = TGV(single(u0), lambda_TGV, alpha1, alpha0, iter_TGV); toc; rmseTGV = (RMSE(u_tgv(:),Im(:))); fprintf('%s %f \n', 'RMSE error for TGV is:', rmseTGV); figure; imshow(u_tgv, [0 1]); title('TGV denoised image (CPU)'); %% % fprintf('Denoise using the TGV model (GPU) \n'); -% lambda_TGV = 0.04; % regularisation parameter -% alpha1 = 1; % parameter to control the first-order term -% alpha0 = 0.7; % parameter to control the second-order term -% iter_TGV = 500; % number of Primal-Dual iterations for TGV +% lambda_TGV = 0.045; % regularisation parameter +% alpha1 = 1.0; % parameter to control the first-order term +% alpha0 = 2.0; % parameter to control the second-order term +% iter_TGV = 2000; % number of Primal-Dual iterations for TGV % tic; u_tgv_gpu = TGV_GPU(single(u0), lambda_TGV, alpha1, alpha0, iter_TGV); toc; % rmseTGV_gpu = (RMSE(u_tgv_gpu(:),Im(:))); % fprintf('%s %f \n', 'RMSE error for TGV is:', rmseTGV_gpu); diff --git a/Wrappers/Matlab/mex_compile/regularisers_CPU/TGV.c b/Wrappers/Matlab/mex_compile/regularisers_CPU/TGV.c index 5459bf5..aa4eed4 100644 --- a/Wrappers/Matlab/mex_compile/regularisers_CPU/TGV.c +++ b/Wrappers/Matlab/mex_compile/regularisers_CPU/TGV.c @@ -21,14 +21,14 @@ limitations under the License. #include "TGV_core.h" /* C-OMP implementation of Primal-Dual denoising method for - * Total Generilized Variation (TGV)-L2 model [1] (2D case only) + * Total Generilized Variation (TGV)-L2 model [1] (2D/3D) * * Input Parameters: - * 1. Noisy image (2D) (required) - * 2. lambda - regularisation parameter (required) - * 3. parameter to control the first-order term (alpha1) (default - 1) - * 4. parameter to control the second-order term (alpha0) (default - 0.5) - * 5. Number of Chambolle-Pock (Primal-Dual) iterations (default is 300) + * 1. Noisy image/volume (2D/3D) + * 2. lambda - regularisation parameter + * 3. parameter to control the first-order term (alpha1) + * 4. parameter to control the second-order term (alpha0) + * 5. Number of Chambolle-Pock (Primal-Dual) iterations * 6. Lipshitz constant (default is 12) * * Output: @@ -44,7 +44,7 @@ void mexFunction( { int number_of_dims, iter; - mwSize dimX, dimY; + mwSize dimX, dimY, dimZ; const mwSize *dim_array; float *Input, *Output=NULL, lambda, alpha0, alpha1, L2; @@ -55,7 +55,7 @@ void mexFunction( /*Handling Matlab input data*/ if ((nrhs < 2) || (nrhs > 6)) mexErrMsgTxt("At least 2 parameters is required, all parameters are: Image(2D), Regularisation parameter, alpha0, alpha1, iterations number, Lipshitz Constant"); - Input = (float *) mxGetData(prhs[0]); /*noisy image (2D) */ + Input = (float *) mxGetData(prhs[0]); /*noisy image/volume */ lambda = (float) mxGetScalar(prhs[1]); /* regularisation parameter */ alpha1 = 1.0f; /* parameter to control the first-order term */ alpha0 = 0.5f; /* parameter to control the second-order term */ @@ -69,12 +69,15 @@ void mexFunction( if (nrhs == 6) L2 = (float) mxGetScalar(prhs[5]); /* Lipshitz constant */ /*Handling Matlab output data*/ - dimX = dim_array[0]; dimY = dim_array[1]; + dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2]; if (number_of_dims == 2) { - Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); - /* running the function */ - TGV_main(Input, Output, lambda, alpha1, alpha0, iter, L2, dimX, dimY); + dimZ = 1; /*2D case*/ + Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); } - if (number_of_dims == 3) {mexErrMsgTxt("Only 2D images accepted");} + if (number_of_dims == 3) { + Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); + } + /* running the function */ + TGV_main(Input, Output, lambda, alpha1, alpha0, iter, L2, dimX, dimY, dimZ); } -- cgit v1.2.3