From 309d84445b5ec2980db7c79832958a6481343f17 Mon Sep 17 00:00:00 2001 From: Daniil Kazantsev Date: Tue, 6 Mar 2018 12:46:09 +0000 Subject: FGP/ROF updates --- Core/regularizers_CPU/FGP_TV_core.c | 139 +++++++++----------- Core/regularizers_CPU/FGP_TV_core.h | 8 +- Core/regularizers_CPU/utils.c | 13 +- .../Matlab/mex_compile/regularizers_CPU/FGP_TV.c | 142 ++------------------- Wrappers/Python/demo/test_cpu_regularizers.py | 10 +- Wrappers/Python/src/cpu_regularizers.pyx | 5 +- 6 files changed, 89 insertions(+), 228 deletions(-) diff --git a/Core/regularizers_CPU/FGP_TV_core.c b/Core/regularizers_CPU/FGP_TV_core.c index 9c0fcfc..7a8ce48 100644 --- a/Core/regularizers_CPU/FGP_TV_core.c +++ b/Core/regularizers_CPU/FGP_TV_core.c @@ -46,7 +46,7 @@ float TV_FGP_CPU_main(float *Input, float *Output, float lambdaPar, int iteratio int count = 0; if (dimZ <= 1) { - /*2D case */ + /*2D case */ float *Output_prev=NULL, *P1=NULL, *P2=NULL, *P1_prev=NULL, *P2_prev=NULL, *R1=NULL, *R2=NULL; DimTotal = dimX*dimY; @@ -65,28 +65,28 @@ float TV_FGP_CPU_main(float *Input, float *Output, float lambdaPar, int iteratio Obj_func2D(Input, Output, R1, R2, lambdaPar, dimX, dimY); /* apply nonnegativity */ - if (nonneg == 1) for(j=0; j 4) break; + if (count > 4) break; /*storing old values*/ copyIm(Output, Output_prev, dimX, dimY, 1); @@ -120,21 +120,21 @@ float TV_FGP_CPU_main(float *Input, float *Output, float lambdaPar, int iteratio Obj_func3D(Input, Output, R1, R2, R3, lambdaPar, dimX, dimY, dimZ); /* apply nonnegativity */ - if (nonneg == 1) for(j=0; j 1.0f) { - sq_denom = 1.0f/sqrt(denom); - P1[index] = P1[index]*sq_denom; - P2[index] = P2[index]*sq_denom; + sq_denom = 1.0f/sqrtf(denom); + P1[i] = P1[i]*sq_denom; + P2[i] = P2[i]*sq_denom; } - }} + } } else { /* anisotropic TV*/ -#pragma omp parallel for shared(P1,P2) private(index,i,j,val1,val2) - for(i=0; i 1.0f) { - sq_denom = 1.0f/sqrt(denom); - P1[index] = P1[index]*sq_denom; - P2[index] = P2[index]*sq_denom; - P3[index] = P3[index]*sq_denom; - } - }}} + sq_denom = 1.0f/sqrtf(denom); + P1[i] = P1[i]*sq_denom; + P2[i] = P2[i]*sq_denom; + P3[i] = P3[i]*sq_denom; + } + } } else { /* anisotropic TV*/ -#pragma omp parallel for shared(P1,P2,P3) private(index,i,j,k,val1,val2,val3) - for(i=0; i /* Copy Image */ float copyIm(float *A, float *U, int dimX, int dimY, int dimZ) @@ -34,7 +35,7 @@ float copyIm(float *A, float *U, int dimX, int dimY, int dimZ) float TV_energy2D(float *U, float *U0, float *E_val, float lambda, int dimX, int dimY) { int i, j, i1, j1, index; - float NOMx_2, NOMy_2, E_Grad, E_Data; + float NOMx_2, NOMy_2, E_Grad=0.0f, E_Data=0.0f; /* first calculate \grad U_xy*/ for(j=0; j 5)) mexErrMsgTxt("At least 2 parameters is required: Image(2D/3D), Regularization parameter. The full list of parameters: Image(2D/3D), Regularization parameter, iterations number, tolerance, penalty type ('iso' or 'l1')"); - A = (float *) mxGetData(prhs[0]); /*noisy image (2D/3D) */ + Input = (float *) mxGetData(prhs[0]); /*noisy image (2D/3D) */ lambda = (float) mxGetScalar(prhs[1]); /* regularization parameter */ - iter = 50; /* default iterations number */ + iter = 300; /* default iterations number */ epsil = 0.0001; /* default tolerance constant */ methTV = 0; /* default isotropic TV penalty */ @@ -78,139 +78,17 @@ void mexFunction( if (strcmp(penalty_type, "l1") == 0) methTV = 1; /* enable 'l1' penalty */ mxFree(penalty_type); } - /*output function value (last iteration) */ - plhs[1] = mxCreateNumericMatrix(1, 1, mxSINGLE_CLASS, mxREAL); - float *funcvalA = (float *) mxGetData(plhs[1]); if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); } /*Handling Matlab output data*/ - dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2]; - - tk = 1.0f; - tkp1=1.0f; - count = 0; - re_old = 0.0f; + dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2]; if (number_of_dims == 2) { dimZ = 1; /*2D case*/ - D = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); - D_old = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); - P1 = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); - P2 = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); - P1_old = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); - P2_old = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); - R1 = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); - R2 = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); - - /* begin iterations */ - for(ll=0; ll 4) { - Obj_func_CALC2D(A, D, funcvalA, lambda, dimX, dimY); - break; } - - /* check that the residual norm is decreasing */ - if (ll > 2) { - if (re > re_old) { - Obj_func_CALC2D(A, D, funcvalA, lambda, dimX, dimY); - break; }} - re_old = re; - /*printf("%f %i %i \n", re, ll, count); */ - - /*storing old values*/ - copyIm(D, D_old, dimX, dimY, dimZ); - copyIm(P1, P1_old, dimX, dimY, dimZ); - copyIm(P2, P2_old, dimX, dimY, dimZ); - tk = tkp1; - - /* calculating the objective function value */ - if (ll == (iter-1)) Obj_func_CALC2D(A, D, funcvalA, lambda, dimX, dimY); - } - printf("FGP-TV iterations stopped at iteration %i with the function value %f \n", ll, funcvalA[0]); - } - if (number_of_dims == 3) { - D = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); - D_old = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); - P1 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); - P2 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); - P3 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); - P1_old = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); - P2_old = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); - P3_old = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); - R1 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); - R2 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); - R3 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); - - /* begin iterations */ - for(ll=0; ll 3) { - Obj_func_CALC3D(A, D, funcvalA, lambda, dimX, dimY, dimZ); - break;} - - /* check that the residual norm is decreasing */ - if (ll > 2) { - if (re > re_old) { - Obj_func_CALC3D(A, D, funcvalA, lambda, dimX, dimY, dimZ); - }} - re_old = re; - /*printf("%f %i %i \n", re, ll, count); */ - - /*storing old values*/ - copyIm(D, D_old, dimX, dimY, dimZ); - copyIm(P1, P1_old, dimX, dimY, dimZ); - copyIm(P2, P2_old, dimX, dimY, dimZ); - copyIm(P3, P3_old, dimX, dimY, dimZ); - tk = tkp1; - - if (ll == (iter-1)) Obj_func_CALC3D(A, D, funcvalA, lambda, dimX, dimY, dimZ); - } - printf("FGP-TV iterations stopped at iteration %i with the function value %f \n", ll, funcvalA[0]); - } + Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); + TV_FGP_CPU_main(Input, Output, lambda, iter, epsil, methTV, 0, 0, dimX, dimY, dimZ); + } + if (number_of_dims == 3) { + } } diff --git a/Wrappers/Python/demo/test_cpu_regularizers.py b/Wrappers/Python/demo/test_cpu_regularizers.py index f1eb3c3..7f08605 100644 --- a/Wrappers/Python/demo/test_cpu_regularizers.py +++ b/Wrappers/Python/demo/test_cpu_regularizers.py @@ -111,12 +111,10 @@ rms = rmse(Im, splitbregman) pars['rmse'] = rms txtstr = printParametersToString(pars) txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) -print (txtstr) - +print (txtstr) a=fig.add_subplot(2,4,2) - # these are matplotlib.patch.Patch properties props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) # place a text box in upper left in axes coords @@ -130,14 +128,14 @@ imgplot = plt.imshow(splitbregman,\ start_time = timeit.default_timer() pars = {'algorithm' : TV_FGP_CPU , \ - 'input' : u0, + 'input' : u0,\ 'regularization_parameter':0.07, \ 'number_of_iterations' :300 ,\ 'tolerance_constant':0.00001,\ 'methodTV': 0 ,\ 'nonneg': 0 ,\ - 'printingOut': 0 -} + 'printingOut': 0 + } out = TV_FGP_CPU (pars['input'], pars['regularization_parameter'], diff --git a/Wrappers/Python/src/cpu_regularizers.pyx b/Wrappers/Python/src/cpu_regularizers.pyx index d62ca59..60e8627 100644 --- a/Wrappers/Python/src/cpu_regularizers.pyx +++ b/Wrappers/Python/src/cpu_regularizers.pyx @@ -93,7 +93,7 @@ def TV_FGP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, return outputData def TV_FGP_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, - float regularization_parameter, + float regularization_parameter, int iterationsNumb, float tolerance_param, int methodTV, @@ -109,11 +109,10 @@ def TV_FGP_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, #/* Run ROF iterations for 3D data */ TV_FGP_CPU_main(&inputData[0,0,0], &outputData[0,0,0], regularization_parameter, - iterationsNumb, + iterationsNumb, tolerance_param, methodTV, nonneg, printM, dims[0], dims[1], dims[2]) - return outputData -- cgit v1.2.3