diff options
Diffstat (limited to 'Wrappers/Python')
| -rw-r--r-- | Wrappers/Python/ccpi/filters/regularizers.py | 44 | ||||
| -rw-r--r-- | Wrappers/Python/demo/test_cpu_regularizers.py | 55 | ||||
| -rw-r--r-- | Wrappers/Python/setup-regularizers.py.in | 6 | ||||
| -rw-r--r-- | Wrappers/Python/src/cpu_regularizers.cpp | 291 | ||||
| -rw-r--r-- | Wrappers/Python/src/cpu_regularizers.pyx | 106 | ||||
| -rw-r--r-- | Wrappers/Python/src/gpu_regularizers.pyx | 151 | ||||
| -rw-r--r-- | Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py | 219 | ||||
| -rw-r--r-- | Wrappers/Python/test/test_gpu_regularizers.py | 113 |
8 files changed, 588 insertions, 397 deletions
diff --git a/Wrappers/Python/ccpi/filters/regularizers.py b/Wrappers/Python/ccpi/filters/regularizers.py new file mode 100644 index 0000000..d6dfa8c --- /dev/null +++ b/Wrappers/Python/ccpi/filters/regularizers.py @@ -0,0 +1,44 @@ +""" +script which assigns a proper device core function based on a flag ('cpu' or 'gpu') +""" + +from ccpi.filters.cpu_regularizers_cython import TV_ROF_CPU, TV_FGP_CPU +from ccpi.filters.gpu_regularizers import TV_ROF_GPU, TV_FGP_GPU + +def ROF_TV(inputData, regularization_parameter, iterations, + time_marching_parameter,device='cpu'): + if device == 'cpu': + return TV_ROF_CPU(inputData, + regularization_parameter, + iterations, + time_marching_parameter) + elif device == 'gpu': + return TV_ROF_GPU(inputData, + regularization_parameter, + iterations, + time_marching_parameter) + else: + raise ValueError('Unknown device {0}. Expecting gpu or cpu'\ + .format(device)) + +def FGP_TV(inputData, regularization_parameter,iterations, + tolerance_param, methodTV, nonneg, printM, device='cpu'): + if device == 'cpu': + return TV_FGP_CPU(inputData, + regularization_parameter, + iterations, + tolerance_param, + methodTV, + nonneg, + printM) + elif device == 'gpu': + return TV_FGP_GPU(inputData, + regularization_parameter, + iterations, + tolerance_param, + methodTV, + nonneg, + printM) + else: + raise ValueError('Unknown device {0}. Expecting gpu or cpu'\ + .format(device))
\ No newline at end of file diff --git a/Wrappers/Python/demo/test_cpu_regularizers.py b/Wrappers/Python/demo/test_cpu_regularizers.py index 5908c3c..76b9de7 100644 --- a/Wrappers/Python/demo/test_cpu_regularizers.py +++ b/Wrappers/Python/demo/test_cpu_regularizers.py @@ -11,11 +11,8 @@ import numpy as np import os from enum import Enum import timeit -from ccpi.filters.cpu_regularizers_boost import SplitBregman_TV , FGP_TV ,\ - LLT_model, PatchBased_Regul ,\ - TGV_PD -from ccpi.filters.cpu_regularizers_cython import ROF_TV - +from ccpi.filters.cpu_regularizers_boost import SplitBregman_TV, LLT_model, PatchBased_Regul, TGV_PD +from ccpi.filters.regularizers import ROF_TV, FGP_TV ############################################################################### #https://stackoverflow.com/questions/13875989/comparing-image-in-url-to-image-in-filesystem-in-python/13884956#13884956 #NRMSE a normalization of the root of the mean squared error @@ -112,12 +109,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 @@ -128,23 +123,26 @@ imgplot = plt.imshow(splitbregman,\ ) ###################### FGP_TV ######################################### -# u = FGP_TV(single(u0), 0.05, 100, 1e-04); + start_time = timeit.default_timer() pars = {'algorithm' : FGP_TV , \ - 'input' : u0, - 'regularization_parameter':0.05, \ - 'number_of_iterations' :200 ,\ - 'tolerance_constant':1e-4,\ - 'TV_penalty': 0 -} + 'input' : u0,\ + 'regularization_parameter':0.07, \ + 'number_of_iterations' :300 ,\ + 'tolerance_constant':0.00001,\ + 'methodTV': 0 ,\ + 'nonneg': 0 ,\ + 'printingOut': 0 + } -out = FGP_TV (pars['input'], +fgp = FGP_TV(pars['input'], pars['regularization_parameter'], pars['number_of_iterations'], pars['tolerance_constant'], - pars['TV_penalty']) + pars['methodTV'], + pars['nonneg'], + pars['printingOut'], 'cpu') -fgp = out[0] rms = rmse(Im, fgp) pars['rmse'] = rms @@ -165,8 +163,8 @@ imgplot = plt.imshow(fgp, \ a.text(0.05, 0.95, txtstr, transform=a.transAxes, fontsize=14, verticalalignment='top', bbox=props) -###################### LLT_model ######################################### +###################### LLT_model ######################################### start_time = timeit.default_timer() pars = {'algorithm': LLT_model , \ @@ -201,8 +199,6 @@ a.text(0.05, 0.95, txtstr, transform=a.transAxes, fontsize=14, imgplot = plt.imshow(llt,\ cmap="gray" ) - - # ###################### PatchBased_Regul ######################################### # # Quick 2D denoising example in Matlab: # # Im = double(imread('lena_gray_256.tif'))/255; % loading image @@ -284,16 +280,15 @@ start_time = timeit.default_timer() pars = {'algorithm': ROF_TV , \ 'input' : u0,\ - 'regularization_parameter':1,\ - 'marching_step': 0.003,\ + 'regularization_parameter':0.07,\ + 'marching_step': 0.0025,\ 'number_of_iterations': 300 } rof = ROF_TV(pars['input'], - pars['number_of_iterations'], - pars['regularization_parameter'], - pars['marching_step'] - ) -#tgv = out + pars['regularization_parameter'], + pars['number_of_iterations'], + pars['marching_step'], 'cpu') + rms = rmse(Im, rof) pars['rmse'] = rms @@ -307,9 +302,7 @@ props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) # place a text box in upper left in axes coords a.text(0.05, 0.95, txtstr, transform=a.transAxes, fontsize=14, verticalalignment='top', bbox=props) -imgplot = plt.imshow(tgv, cmap="gray") - - +imgplot = plt.imshow(rof, cmap="gray") plt.show() diff --git a/Wrappers/Python/setup-regularizers.py.in b/Wrappers/Python/setup-regularizers.py.in index a125261..8655a2e 100644 --- a/Wrappers/Python/setup-regularizers.py.in +++ b/Wrappers/Python/setup-regularizers.py.in @@ -34,7 +34,9 @@ extra_libraries = ['cilreg'] extra_include_dirs += [os.path.join(".." , ".." , "Core"), os.path.join(".." , ".." , "Core", "regularizers_CPU"), - os.path.join(".." , ".." , "Core", "regularizers_GPU") , + os.path.join(".." , ".." , "Core", "regularizers_GPU" , "Diffus_HO" ) , + os.path.join(".." , ".." , "Core", "regularizers_GPU" , "NL_Regul" ) , + os.path.join(".." , ".." , "Core", "regularizers_GPU" , "TV_ROF" ) , "."] if platform.system() == 'Windows': @@ -81,4 +83,4 @@ setup( ) -@SETUP_GPU_WRAPPERS@
\ No newline at end of file +@SETUP_GPU_WRAPPERS@ diff --git a/Wrappers/Python/src/cpu_regularizers.cpp b/Wrappers/Python/src/cpu_regularizers.cpp index e311570..3529ebd 100644 --- a/Wrappers/Python/src/cpu_regularizers.cpp +++ b/Wrappers/Python/src/cpu_regularizers.cpp @@ -27,7 +27,6 @@ limitations under the License. #include "boost/tuple/tuple.hpp" #include "SplitBregman_TV_core.h" -#include "FGP_TV_core.h" #include "LLT_model_core.h" #include "PatchBased_Regul_core.h" #include "TGV_PD_core.h" @@ -303,292 +302,6 @@ bp::list SplitBregman_TV(np::ndarray input, double d_mu, int iter, double d_epsi } - - -bp::list FGP_TV(np::ndarray input, double d_mu, int iter, double d_epsil, int methTV) { - - // the result is in the following list - bp::list result; - - int number_of_dims, dimX, dimY, dimZ, ll, j, count; - float *A, *D = NULL, *D_old = NULL, *P1 = NULL, *P2 = NULL, *P3 = NULL, *P1_old = NULL, *P2_old = NULL, *P3_old = NULL, *R1 = NULL, *R2 = NULL, *R3 = NULL; - float lambda, tk, tkp1, re, re1, re_old, epsil, funcval; - - //number_of_dims = mxGetNumberOfDimensions(prhs[0]); - //dim_array = mxGetDimensions(prhs[0]); - - number_of_dims = input.get_nd(); - int dim_array[3]; - - dim_array[0] = input.shape(0); - dim_array[1] = input.shape(1); - if (number_of_dims == 2) { - dim_array[2] = -1; - } - else { - dim_array[2] = input.shape(2); - } - // Parameter handling is be done in Python - ///*Handling Matlab input data*/ - //if ((nrhs < 2) || (nrhs > 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')"); - - ///*Handling Matlab input data*/ - //A = (float *)mxGetData(prhs[0]); /*noisy image (2D/3D) */ - A = reinterpret_cast<float *>(input.get_data()); - - //mu = (float)mxGetScalar(prhs[1]); /* regularization parameter */ - lambda = (float)d_mu; - - //iter = 35; /* default iterations number */ - - //epsil = 0.0001; /* default tolerance constant */ - epsil = (float)d_epsil; - //methTV = 0; /* default isotropic TV penalty */ - //if ((nrhs == 3) || (nrhs == 4) || (nrhs == 5)) iter = (int)mxGetScalar(prhs[2]); /* iterations number */ - //if ((nrhs == 4) || (nrhs == 5)) epsil = (float)mxGetScalar(prhs[3]); /* tolerance constant */ - //if (nrhs == 5) { - // char *penalty_type; - // penalty_type = mxArrayToString(prhs[4]); /* choosing TV penalty: 'iso' or 'l1', 'iso' is the default */ - // if ((strcmp(penalty_type, "l1") != 0) && (strcmp(penalty_type, "iso") != 0)) mexErrMsgTxt("Choose TV type: 'iso' or 'l1',"); - // if (strcmp(penalty_type, "l1") == 0) methTV = 1; /* enable 'l1' penalty */ - // mxFree(penalty_type); - //} - //if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) { mexErrMsgTxt("The input image must be in a single precision"); } - - //plhs[1] = mxCreateNumericMatrix(1, 1, mxSINGLE_CLASS, mxREAL); - bp::tuple shape1 = bp::make_tuple(dim_array[0], dim_array[1]); - np::dtype dtype = np::dtype::get_builtin<float>(); - np::ndarray out1 = np::zeros(shape1, dtype); - - //float *funcvalA = (float *)mxGetData(plhs[1]); - float * funcvalA = reinterpret_cast<float *>(out1.get_data()); - //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 = 1; - re_old = 0.0f; - - 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));*/ - - bp::tuple shape = bp::make_tuple(dim_array[0], dim_array[1]); - np::dtype dtype = np::dtype::get_builtin<float>(); - - - np::ndarray npD = np::zeros(shape, dtype); - np::ndarray npD_old = np::zeros(shape, dtype); - np::ndarray npP1 = np::zeros(shape, dtype); - np::ndarray npP2 = np::zeros(shape, dtype); - np::ndarray npP1_old = np::zeros(shape, dtype); - np::ndarray npP2_old = np::zeros(shape, dtype); - np::ndarray npR1 = np::zeros(shape, dtype); - np::ndarray npR2 = np::zeros(shape, dtype); - - D = reinterpret_cast<float *>(npD.get_data()); - D_old = reinterpret_cast<float *>(npD_old.get_data()); - P1 = reinterpret_cast<float *>(npP1.get_data()); - P2 = reinterpret_cast<float *>(npP2.get_data()); - P1_old = reinterpret_cast<float *>(npP1_old.get_data()); - P2_old = reinterpret_cast<float *>(npP2_old.get_data()); - R1 = reinterpret_cast<float *>(npR1.get_data()); - R2 = reinterpret_cast<float *>(npR2.get_data()); - - /* begin iterations */ - for (ll = 0; ll<iter; ll++) { - /* computing the gradient of the objective function */ - Obj_func2D(A, D, R1, R2, lambda, dimX, dimY); - - /*Taking a step towards minus of the gradient*/ - Grad_func2D(P1, P2, D, R1, R2, lambda, dimX, dimY); - - /* projection step */ - Proj_func2D(P1, P2, methTV, dimX, dimY); - - /*updating R and t*/ - tkp1 = (1.0f + sqrt(1.0f + 4.0f*tk*tk))*0.5f; - Rupd_func2D(P1, P1_old, P2, P2_old, R1, R2, tkp1, tk, dimX, dimY); - - /* calculate norm */ - re = 0.0f; re1 = 0.0f; - for (j = 0; j<dimX*dimY*dimZ; j++) - { - re += pow(D[j] - D_old[j], 2); - re1 += pow(D[j], 2); - } - re = sqrt(re) / sqrt(re1); - if (re < epsil) count++; - if (count > 4) { - Obj_func2D(A, D, P1, P2, lambda, dimX, dimY); - funcval = 0.0f; - for (j = 0; j<dimX*dimY*dimZ; j++) funcval += pow(D[j], 2); - //funcvalA[0] = sqrt(funcval); - float fv = sqrt(funcval); - std::memcpy(funcvalA, &fv, sizeof(float)); - break; - } - - /* check that the residual norm is decreasing */ - if (ll > 2) { - if (re > re_old) { - Obj_func2D(A, D, P1, P2, lambda, dimX, dimY); - funcval = 0.0f; - for (j = 0; j<dimX*dimY*dimZ; j++) funcval += pow(D[j], 2); - //funcvalA[0] = sqrt(funcval); - float fv = sqrt(funcval); - std::memcpy(funcvalA, &fv, sizeof(float)); - 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_func2D(A, D, P1, P2, lambda, dimX, dimY); - funcval = 0.0f; - for (j = 0; j<dimX*dimY*dimZ; j++) funcval += pow(D[j], 2); - //funcvalA[0] = sqrt(funcval); - float fv = sqrt(funcval); - std::memcpy(funcvalA, &fv, sizeof(float)); - } - } - //printf("FGP-TV iterations stopped at iteration %i with the function value %f \n", ll, funcvalA[0]); - result.append<np::ndarray>(npD); - result.append<np::ndarray>(out1); - result.append<int>(ll); - } - 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));*/ - bp::tuple shape = bp::make_tuple(dim_array[0], dim_array[1], dim_array[2]); - np::dtype dtype = np::dtype::get_builtin<float>(); - - np::ndarray npD = np::zeros(shape, dtype); - np::ndarray npD_old = np::zeros(shape, dtype); - np::ndarray npP1 = np::zeros(shape, dtype); - np::ndarray npP2 = np::zeros(shape, dtype); - np::ndarray npP3 = np::zeros(shape, dtype); - np::ndarray npP1_old = np::zeros(shape, dtype); - np::ndarray npP2_old = np::zeros(shape, dtype); - np::ndarray npP3_old = np::zeros(shape, dtype); - np::ndarray npR1 = np::zeros(shape, dtype); - np::ndarray npR2 = np::zeros(shape, dtype); - np::ndarray npR3 = np::zeros(shape, dtype); - - D = reinterpret_cast<float *>(npD.get_data()); - D_old = reinterpret_cast<float *>(npD_old.get_data()); - P1 = reinterpret_cast<float *>(npP1.get_data()); - P2 = reinterpret_cast<float *>(npP2.get_data()); - P3 = reinterpret_cast<float *>(npP3.get_data()); - P1_old = reinterpret_cast<float *>(npP1_old.get_data()); - P2_old = reinterpret_cast<float *>(npP2_old.get_data()); - P3_old = reinterpret_cast<float *>(npP3_old.get_data()); - R1 = reinterpret_cast<float *>(npR1.get_data()); - R2 = reinterpret_cast<float *>(npR2.get_data()); - R3 = reinterpret_cast<float *>(npR3.get_data()); - /* begin iterations */ - for (ll = 0; ll<iter; ll++) { - /* computing the gradient of the objective function */ - Obj_func3D(A, D, R1, R2, R3, lambda, dimX, dimY, dimZ); - /*Taking a step towards minus of the gradient*/ - Grad_func3D(P1, P2, P3, D, R1, R2, R3, lambda, dimX, dimY, dimZ); - - /* projection step */ - Proj_func3D(P1, P2, P3, dimX, dimY, dimZ); - - /*updating R and t*/ - tkp1 = (1.0f + sqrt(1.0f + 4.0f*tk*tk))*0.5f; - Rupd_func3D(P1, P1_old, P2, P2_old, P3, P3_old, R1, R2, R3, tkp1, tk, dimX, dimY, dimZ); - - /* calculate norm - stopping rules*/ - re = 0.0f; re1 = 0.0f; - for (j = 0; j<dimX*dimY*dimZ; j++) - { - re += pow(D[j] - D_old[j], 2); - re1 += pow(D[j], 2); - } - re = sqrt(re) / sqrt(re1); - /* stop if the norm residual is less than the tolerance EPS */ - if (re < epsil) count++; - if (count > 3) { - Obj_func3D(A, D, P1, P2, P3, lambda, dimX, dimY, dimZ); - funcval = 0.0f; - for (j = 0; j<dimX*dimY*dimZ; j++) funcval += pow(D[j], 2); - //funcvalA[0] = sqrt(funcval); - float fv = sqrt(funcval); - std::memcpy(funcvalA, &fv, sizeof(float)); - break; - } - - /* check that the residual norm is decreasing */ - if (ll > 2) { - if (re > re_old) { - Obj_func3D(A, D, P1, P2, P3, lambda, dimX, dimY, dimZ); - funcval = 0.0f; - for (j = 0; j<dimX*dimY*dimZ; j++) funcval += pow(D[j], 2); - //funcvalA[0] = sqrt(funcval); - float fv = sqrt(funcval); - std::memcpy(funcvalA, &fv, sizeof(float)); - 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); - copyIm(P3, P3_old, dimX, dimY, dimZ); - tk = tkp1; - - if (ll == (iter - 1)) { - Obj_func3D(A, D, P1, P2, P3, lambda, dimX, dimY, dimZ); - funcval = 0.0f; - for (j = 0; j<dimX*dimY*dimZ; j++) funcval += pow(D[j], 2); - //funcvalA[0] = sqrt(funcval); - float fv = sqrt(funcval); - std::memcpy(funcvalA, &fv, sizeof(float)); - } - - } - //printf("FGP-TV iterations stopped at iteration %i with the function value %f \n", ll, funcvalA[0]); - result.append<np::ndarray>(npD); - result.append<np::ndarray>(out1); - result.append<int>(ll); - } - - return result; -} - bp::list LLT_model(np::ndarray input, double d_lambda, double d_tau, int iter, double d_epsil, int switcher) { // the result is in the following list bp::list result; @@ -1022,9 +735,6 @@ bp::list TGV_PD(np::ndarray input, double d_lambda, double d_alpha1, double d_al result.append<np::ndarray>(npU); } - - - return result; } @@ -1040,7 +750,6 @@ BOOST_PYTHON_MODULE(cpu_regularizers_boost) np::dtype dt2 = np::dtype::get_builtin<uint16_t>(); def("SplitBregman_TV", SplitBregman_TV); - def("FGP_TV", FGP_TV); def("LLT_model", LLT_model); def("PatchBased_Regul", PatchBased_Regul); def("TGV_PD", TGV_PD); diff --git a/Wrappers/Python/src/cpu_regularizers.pyx b/Wrappers/Python/src/cpu_regularizers.pyx index e7ff78c..60e8627 100644 --- a/Wrappers/Python/src/cpu_regularizers.pyx +++ b/Wrappers/Python/src/cpu_regularizers.pyx @@ -11,62 +11,108 @@ 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. -Author: Edoardo Pasca +Author: Edoardo Pasca, Daniil Kazantsev """ import cython import numpy as np cimport numpy as np -cdef extern float TV_ROF(float *A, float *B, int dimX, int dimY, int dimZ, - int iterationsNumb, float tau, float flambda); +cdef extern float TV_ROF_CPU_main(float *Input, float *Output, float lambdaPar, int iterationsNumb, float tau, int dimX, int dimY, int dimZ); +cdef extern float TV_FGP_CPU_main(float *Input, float *Output, float lambdaPar, int iterationsNumb, float epsil, int methodTV, int nonneg, int printM, int dimX, int dimY, int dimZ); +def TV_ROF_CPU(inputData, regularization_parameter, iterationsNumb, marching_step_parameter): + if inputData.ndim == 2: + return TV_ROF_2D(inputData, regularization_parameter, iterationsNumb, marching_step_parameter) + elif inputData.ndim == 3: + return TV_ROF_3D(inputData, regularization_parameter, iterationsNumb, marching_step_parameter) + +def TV_ROF_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, + float regularization_parameter, + int iterationsNumb, + float marching_step_parameter): + cdef long dims[2] + dims[0] = inputData.shape[0] + dims[1] = inputData.shape[1] + + cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \ + np.zeros([dims[0],dims[1]], dtype='float32') + + # Run ROF iterations for 2D data + TV_ROF_CPU_main(&inputData[0,0], &outputData[0,0], regularization_parameter, iterationsNumb, marching_step_parameter, dims[0], dims[1], 1) + + return outputData + +def TV_ROF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, + int iterationsNumb, + float regularization_parameter, + float marching_step_parameter): + cdef long dims[3] + dims[0] = inputData.shape[0] + dims[1] = inputData.shape[1] + dims[2] = inputData.shape[2] + + cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \ + np.zeros([dims[0],dims[1],dims[2]], dtype='float32') + + # Run ROF iterations for 3D data + TV_ROF_CPU_main(&inputData[0,0,0], &outputData[0,0,0], regularization_parameter, iterationsNumb, marching_step_parameter, dims[0], dims[1], dims[2]) -def ROF_TV(inputData, iterations, regularization_parameter, - marching_step_parameter): + return outputData + +def TV_FGP_CPU(inputData, regularization_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, printM): if inputData.ndim == 2: - return ROF_TV_2D(inputData, iterations, regularization_parameter, - marching_step_parameter) + return TV_FGP_2D(inputData, regularization_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, printM) elif inputData.ndim == 3: - return ROF_TV_3D(inputData, iterations, regularization_parameter, - marching_step_parameter) + return TV_FGP_3D(inputData, regularization_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, printM) -def ROF_TV_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, - int iterations, +def TV_FGP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, float regularization_parameter, - float marching_step_parameter - ): + int iterationsNumb, + float tolerance_param, + int methodTV, + int nonneg, + int printM): cdef long dims[2] dims[0] = inputData.shape[0] dims[1] = inputData.shape[1] - cdef np.ndarray[np.float32_t, ndim=2, mode="c"] B = \ + cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \ np.zeros([dims[0],dims[1]], dtype='float32') #/* Run ROF iterations for 2D data */ - TV_ROF(&inputData[0,0], &B[0,0], dims[0], dims[1], 1, iterations, - marching_step_parameter, regularization_parameter) + TV_FGP_CPU_main(&inputData[0,0], &outputData[0,0], regularization_parameter, + iterationsNumb, + tolerance_param, + methodTV, + nonneg, + printM, + dims[0], dims[1], 1) - return B - + return outputData -def ROF_TV_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, - int iterations, +def TV_FGP_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, float regularization_parameter, - float marching_step_parameter - ): - cdef long dims[2] + int iterationsNumb, + float tolerance_param, + int methodTV, + int nonneg, + int printM): + cdef long dims[3] dims[0] = inputData.shape[0] dims[1] = inputData.shape[1] dims[2] = inputData.shape[2] - cdef np.ndarray[np.float32_t, ndim=3, mode="c"] B = \ - np.zeros([dims[0],dims[1],dims[2]], dtype='float32') + cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \ + np.zeros([dims[0], dims[1], dims[2]], dtype='float32') #/* Run ROF iterations for 3D data */ - TV_ROF(&inputData[0,0, 0], &B[0,0, 0], dims[0], dims[1], dims[2], iterations, - marching_step_parameter, regularization_parameter) - - return B - + TV_FGP_CPU_main(&inputData[0,0,0], &outputData[0,0,0], regularization_parameter, + iterationsNumb, + tolerance_param, + methodTV, + nonneg, + printM, + dims[0], dims[1], dims[2]) + return outputData diff --git a/Wrappers/Python/src/gpu_regularizers.pyx b/Wrappers/Python/src/gpu_regularizers.pyx index 5a5d274..f96156a 100644 --- a/Wrappers/Python/src/gpu_regularizers.pyx +++ b/Wrappers/Python/src/gpu_regularizers.pyx @@ -11,7 +11,7 @@ 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. -Author: Edoardo Pasca +Author: Edoardo Pasca, Daniil Kazantsev """ import cython @@ -25,12 +25,16 @@ cdef extern void NLM_GPU_kernel(float *A, float* B, float *Eucl_Vec, int N, int M, int Z, int dimension, int SearchW, int SimilW, int SearchW_real, float denh2, float lambdaf); -cdef extern void TV_ROF_GPU(float* A, float* B, int N, int M, int Z, int iter, float tau, float lambdaf); +cdef extern void TV_ROF_GPU_main(float* Input, float* Output, float lambdaPar, int iter, float tau, int N, int M, int Z); +cdef extern void TV_FGP_GPU_main(float *Input, float *Output, float lambdaPar, int iter, float epsil, int methodTV, int nonneg, int printM, int N, int M, int Z); + +# padding function cdef extern float pad_crop(float *A, float *Ap, int OldSizeX, int OldSizeY, int OldSizeZ, int NewSizeX, int NewSizeY, int NewSizeZ, int padXY, int switchpad_crop); - + +#Diffusion 4th order regularizer def Diff4thHajiaboli(inputData, edge_preserv_parameter, iterations, @@ -48,7 +52,7 @@ def Diff4thHajiaboli(inputData, iterations, time_marching_parameter, regularization_parameter) - +# patch-based nonlocal regularization def NML(inputData, SearchW_real, SimilW, @@ -66,23 +70,48 @@ def NML(inputData, SimilW, h, lambdaf) - -def ROF_TV_GPU(inputData, + +# Total-variation Rudin-Osher-Fatemi (ROF) +def TV_ROF_GPU(inputData, + regularization_parameter, iterations, - time_marching_parameter, - regularization_parameter): + time_marching_parameter): if inputData.ndim == 2: return ROFTV2D(inputData, - iterations, - time_marching_parameter, - regularization_parameter) + regularization_parameter, + iterations, + time_marching_parameter) elif inputData.ndim == 3: return ROFTV3D(inputData, + regularization_parameter, iterations, - time_marching_parameter, - regularization_parameter) - - + time_marching_parameter) + +# Total-variation Fast-Gradient-Projection (FGP) +def TV_FGP_GPU(inputData, + regularization_parameter, + iterations, + tolerance_param, + methodTV, + nonneg, + printM): + if inputData.ndim == 2: + return FGPTV2D(inputData, + regularization_parameter, + iterations, + tolerance_param, + methodTV, + nonneg, + printM) + elif inputData.ndim == 3: + return FGPTV3D(inputData, + regularization_parameter, + iterations, + tolerance_param, + methodTV, + nonneg, + printM) +#****************************************************************# def Diff4thHajiaboli2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, float edge_preserv_parameter, int iterations, @@ -329,48 +358,106 @@ def NML3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, switchpad_crop) return B - + +# Total-variation ROF def ROFTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, + float regularization_parameter, int iterations, - float time_marching_parameter, - float regularization_parameter): + float time_marching_parameter): cdef long dims[2] dims[0] = inputData.shape[0] dims[1] = inputData.shape[1] - cdef np.ndarray[np.float32_t, ndim=2, mode="c"] B = \ + cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \ np.zeros([dims[0],dims[1]], dtype='float32') # Running CUDA code here - TV_ROF_GPU( - &inputData[0,0], &B[0,0], - dims[0], dims[1], 0, + TV_ROF_GPU_main( + &inputData[0,0], &outputData[0,0], + regularization_parameter, iterations , time_marching_parameter, - regularization_parameter); + dims[0], dims[1], 1); - return B + return outputData def ROFTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, + float regularization_parameter, int iterations, - float time_marching_parameter, - float regularization_parameter): + float time_marching_parameter): cdef long dims[3] dims[0] = inputData.shape[0] dims[1] = inputData.shape[1] dims[2] = inputData.shape[2] - cdef np.ndarray[np.float32_t, ndim=3, mode="c"] B = \ + cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \ np.zeros([dims[0],dims[1],dims[2]], dtype='float32') # Running CUDA code here - TV_ROF_GPU( - &inputData[0,0,0], &B[0,0,0], - dims[0], dims[1], dims[2], + TV_ROF_GPU_main( + &inputData[0,0,0], &outputData[0,0,0], + regularization_parameter, iterations , time_marching_parameter, - regularization_parameter); + dims[0], dims[1], dims[2]); - return B + return outputData + +# Total-variation Fast-Gradient-Projection (FGP) +def FGPTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, + float regularization_parameter, + int iterations, + float tolerance_param, + int methodTV, + int nonneg, + int printM): + + cdef long dims[2] + dims[0] = inputData.shape[0] + dims[1] = inputData.shape[1] + + cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \ + np.zeros([dims[0],dims[1]], dtype='float32') + + # Running CUDA code here + TV_FGP_GPU_main(&inputData[0,0], &outputData[0,0], + regularization_parameter, + iterations, + tolerance_param, + methodTV, + nonneg, + printM, + dims[0], dims[1], 1); + + return outputData + +def FGPTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, + float regularization_parameter, + int iterations, + float tolerance_param, + int methodTV, + int nonneg, + int printM): + + cdef long dims[3] + dims[0] = inputData.shape[0] + dims[1] = inputData.shape[1] + dims[2] = inputData.shape[2] + + cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \ + np.zeros([dims[0],dims[1],dims[2]], dtype='float32') + + # Running CUDA code here + TV_FGP_GPU_main( + &inputData[0,0,0], &outputData[0,0,0], + regularization_parameter , + iterations, + tolerance_param, + methodTV, + nonneg, + printM, + dims[0], dims[1], dims[2]); + + return outputData diff --git a/Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py b/Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py new file mode 100644 index 0000000..63be1a0 --- /dev/null +++ b/Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py @@ -0,0 +1,219 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Thu Feb 22 11:39:43 2018 + +Testing CPU implementation against the GPU one + +@author: Daniil Kazantsev +""" + +import matplotlib.pyplot as plt +import numpy as np +import os +import timeit +from ccpi.filters.regularizers import ROF_TV, FGP_TV + +############################################################################### +def printParametersToString(pars): + txt = r'' + for key, value in pars.items(): + if key== 'algorithm' : + txt += "{0} = {1}".format(key, value.__name__) + elif key == 'input': + txt += "{0} = {1}".format(key, np.shape(value)) + else: + txt += "{0} = {1}".format(key, value) + txt += '\n' + return txt +############################################################################### +def rmse(im1, im2): + a, b = im1.shape + rmse = np.sqrt(np.sum((im1 - im2) ** 2) / float(a * b)) + return rmse + +filename = os.path.join(".." , ".." , ".." , "data" ,"lena_gray_512.tif") + +# read image +Im = plt.imread(filename) +Im = np.asarray(Im, dtype='float32') + +Im = Im/255 +perc = 0.075 +u0 = Im + np.random.normal(loc = Im , + scale = perc * Im , + size = np.shape(Im)) +# map the u0 u0->u0>0 +f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1) +u0 = f(u0).astype('float32') + + +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +print ("____________ROF-TV bench___________________") +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") + +## plot +fig = plt.figure(1) +plt.suptitle('Comparison of ROF-TV regularizer using CPU and GPU implementations') +a=fig.add_subplot(1,4,1) +a.set_title('Noisy Image') +imgplot = plt.imshow(u0,cmap="gray") + +# set parameters +pars = {'algorithm': ROF_TV, \ + 'input' : u0,\ + 'regularization_parameter':0.04,\ + 'number_of_iterations': 1200,\ + 'time_marching_parameter': 0.0025 + } +print ("#############ROF TV CPU####################") +start_time = timeit.default_timer() +rof_cpu = ROF_TV(pars['input'], + pars['regularization_parameter'], + pars['number_of_iterations'], + pars['time_marching_parameter'],'cpu') +rms = rmse(Im, rof_cpu) +pars['rmse'] = rms + +txtstr = printParametersToString(pars) +txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) +print (txtstr) +a=fig.add_subplot(1,4,2) + +# these are matplotlib.patch.Patch properties +props = dict(boxstyle='round', facecolor='wheat', alpha=0.75) +# place a text box in upper left in axes coords +a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, + verticalalignment='top', bbox=props) +imgplot = plt.imshow(rof_cpu, cmap="gray") +plt.title('{}'.format('CPU results')) + + +print ("##############ROF TV GPU##################") +start_time = timeit.default_timer() +rof_gpu = ROF_TV(pars['input'], + pars['regularization_parameter'], + pars['number_of_iterations'], + pars['time_marching_parameter'],'gpu') + +rms = rmse(Im, rof_gpu) +pars['rmse'] = rms +pars['algorithm'] = ROF_TV +txtstr = printParametersToString(pars) +txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) +print (txtstr) +a=fig.add_subplot(1,4,3) + +# these are matplotlib.patch.Patch properties +props = dict(boxstyle='round', facecolor='wheat', alpha=0.75) +# place a text box in upper left in axes coords +a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, + verticalalignment='top', bbox=props) +imgplot = plt.imshow(rof_gpu, cmap="gray") +plt.title('{}'.format('GPU results')) + + +print ("--------Compare the results--------") +tolerance = 1e-05 +diff_im = np.zeros(np.shape(rof_cpu)) +diff_im = abs(rof_cpu - rof_gpu) +diff_im[diff_im > tolerance] = 1 +a=fig.add_subplot(1,4,4) +imgplot = plt.imshow(diff_im, vmin=0, vmax=1, cmap="gray") +plt.title('{}'.format('Pixels larger threshold difference')) +if (diff_im.sum() > 1): + print ("Arrays do not match!") +else: + print ("Arrays match") + +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +print ("____________FGP-TV bench___________________") +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") + +## plot +fig = plt.figure(2) +plt.suptitle('Comparison of FGP-TV regularizer using CPU and GPU implementations') +a=fig.add_subplot(1,4,1) +a.set_title('Noisy Image') +imgplot = plt.imshow(u0,cmap="gray") + +# set parameters +pars = {'algorithm' : FGP_TV, \ + 'input' : u0,\ + 'regularization_parameter':0.04, \ + 'number_of_iterations' :1200 ,\ + 'tolerance_constant':0.00001,\ + 'methodTV': 0 ,\ + 'nonneg': 0 ,\ + 'printingOut': 0 + } + +print ("#############FGP TV CPU####################") +start_time = timeit.default_timer() +fgp_cpu = FGP_TV(pars['input'], + pars['regularization_parameter'], + pars['number_of_iterations'], + pars['tolerance_constant'], + pars['methodTV'], + pars['nonneg'], + pars['printingOut'],'cpu') + + +rms = rmse(Im, fgp_cpu) +pars['rmse'] = rms + +txtstr = printParametersToString(pars) +txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) +print (txtstr) +a=fig.add_subplot(1,4,2) + +# these are matplotlib.patch.Patch properties +props = dict(boxstyle='round', facecolor='wheat', alpha=0.75) +# place a text box in upper left in axes coords +a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, + verticalalignment='top', bbox=props) +imgplot = plt.imshow(fgp_cpu, cmap="gray") +plt.title('{}'.format('CPU results')) + + +print ("##############FGP TV GPU##################") +start_time = timeit.default_timer() +fgp_gpu = FGP_TV(pars['input'], + pars['regularization_parameter'], + pars['number_of_iterations'], + pars['tolerance_constant'], + pars['methodTV'], + pars['nonneg'], + pars['printingOut'],'gpu') + +rms = rmse(Im, fgp_gpu) +pars['rmse'] = rms +pars['algorithm'] = FGP_TV +txtstr = printParametersToString(pars) +txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) +print (txtstr) +a=fig.add_subplot(1,4,3) + +# these are matplotlib.patch.Patch properties +props = dict(boxstyle='round', facecolor='wheat', alpha=0.75) +# place a text box in upper left in axes coords +a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, + verticalalignment='top', bbox=props) +imgplot = plt.imshow(fgp_gpu, cmap="gray") +plt.title('{}'.format('GPU results')) + + +print ("--------Compare the results--------") +tolerance = 1e-05 +diff_im = np.zeros(np.shape(rof_cpu)) +diff_im = abs(fgp_cpu - fgp_gpu) +diff_im[diff_im > tolerance] = 1 +a=fig.add_subplot(1,4,4) +imgplot = plt.imshow(diff_im, vmin=0, vmax=1, cmap="gray") +plt.title('{}'.format('Pixels larger threshold difference')) +if (diff_im.sum() > 1): + print ("Arrays do not match!") +else: + print ("Arrays match") + + diff --git a/Wrappers/Python/test/test_gpu_regularizers.py b/Wrappers/Python/test/test_gpu_regularizers.py index 735a25d..640b3f9 100644 --- a/Wrappers/Python/test/test_gpu_regularizers.py +++ b/Wrappers/Python/test/test_gpu_regularizers.py @@ -6,14 +6,13 @@ Created on Tue Jan 30 10:24:26 2018 @author: ofn77899 """ - - import matplotlib.pyplot as plt import numpy as np import os from enum import Enum import timeit from ccpi.filters.gpu_regularizers import Diff4thHajiaboli, NML +from ccpi.filters.regularizers import ROF_TV, FGP_TV ############################################################################### def printParametersToString(pars): txt = r'' @@ -52,14 +51,15 @@ u0 = f(u0).astype('float32') ## plot fig = plt.figure() -a=fig.add_subplot(2,3,1) +a=fig.add_subplot(2,4,1) a.set_title('noise') imgplot = plt.imshow(u0,cmap="gray") ## Diff4thHajiaboli start_time = timeit.default_timer() -pars = {'algorithm' : Diff4thHajiaboli , \ +pars = { +'algorithm' : Diff4thHajiaboli , \ 'input' : u0, 'edge_preserv_parameter':0.1 , \ 'number_of_iterations' :250 ,\ @@ -78,7 +78,7 @@ pars['rmse'] = rms txtstr = printParametersToString(pars) txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) print (txtstr) -a=fig.add_subplot(2,3,2) +a=fig.add_subplot(2,4,2) # these are matplotlib.patch.Patch properties props = dict(boxstyle='round', facecolor='wheat', alpha=0.75) @@ -87,14 +87,15 @@ a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=12, verticalalignment='top', bbox=props) imgplot = plt.imshow(d4h, cmap="gray") -a=fig.add_subplot(2,3,5) +a=fig.add_subplot(2,4,6) # 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 a.text(0.05, 0.95, 'd4h - u0', transform=a.transAxes, fontsize=12, verticalalignment='top', bbox=props) -imgplot = plt.imshow((d4h - u0)**2, cmap="gray") +imgplot = plt.imshow((d4h - u0)**2, vmin=0, vmax=0.03, cmap="gray") +plt.colorbar(ticks=[0, 0.03], orientation='vertical') ## Patch Based Regul NML @@ -109,6 +110,7 @@ pars = {'algorithm' : NML , \ } """ pars = { +'algorithm' : NML , \ 'input' : u0, 'regularization_parameter': 0.01,\ 'searching_window_ratio':3, \ @@ -126,7 +128,7 @@ pars['rmse'] = rms txtstr = printParametersToString(pars) txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) print (txtstr) -a=fig.add_subplot(2,3,3) +a=fig.add_subplot(2,4,3) # these are matplotlib.patch.Patch properties props = dict(boxstyle='round', facecolor='wheat', alpha=0.75) @@ -135,14 +137,103 @@ a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=12, verticalalignment='top', bbox=props) imgplot = plt.imshow(nml, cmap="gray") -a=fig.add_subplot(2,3,6) +a=fig.add_subplot(2,4,7) # 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 a.text(0.05, 0.95, 'nml - u0', transform=a.transAxes, fontsize=14, verticalalignment='top', bbox=props) -imgplot = plt.imshow((nml - u0)**2, cmap="gray") +imgplot = plt.imshow((nml - u0)**2, vmin=0, vmax=0.03, cmap="gray") +plt.colorbar(ticks=[0, 0.03], orientation='vertical') -plt.show() + +## Rudin-Osher-Fatemi (ROF) TV regularization +start_time = timeit.default_timer() +pars = { +'algorithm' : ROF_TV , \ + 'input' : u0, + 'regularization_parameter': 0.04,\ + 'number_of_iterations':300,\ + 'time_marching_parameter': 0.0025 + + } + +rof_tv = TV_ROF_GPU(pars['input'], + pars['regularization_parameter'], + pars['number_of_iterations'], + pars['time_marching_parameter'],'gpu') + +rms = rmse(Im, rof_tv) +pars['rmse'] = rms +txtstr = printParametersToString(pars) +txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) +print (txtstr) +a=fig.add_subplot(2,4,4) + +# these are matplotlib.patch.Patch properties +props = dict(boxstyle='round', facecolor='wheat', alpha=0.75) +# place a text box in upper left in axes coords +a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=12, + verticalalignment='top', bbox=props) +imgplot = plt.imshow(rof_tv, cmap="gray") + +a=fig.add_subplot(2,4,8) + +# 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 +a.text(0.05, 0.95, 'rof_tv - u0', transform=a.transAxes, fontsize=14, + verticalalignment='top', bbox=props) +imgplot = plt.imshow((rof_tv - u0)**2, vmin=0, vmax=0.03, cmap="gray") +plt.colorbar(ticks=[0, 0.03], orientation='vertical') +plt.show() + +## Fast-Gradient Projection TV regularization +""" +start_time = timeit.default_timer() + +pars = {'algorithm' : FGP_TV, \ + 'input' : u0,\ + 'regularization_parameter':0.04, \ + 'number_of_iterations' :1200 ,\ + 'tolerance_constant':0.00001,\ + 'methodTV': 0 ,\ + 'nonneg': 0 ,\ + 'printingOut': 0 + } + +fgp_gpu = FGP_TV(pars['input'], + pars['regularization_parameter'], + pars['number_of_iterations'], + pars['tolerance_constant'], + pars['methodTV'], + pars['nonneg'], + pars['printingOut'],'gpu') + +rms = rmse(Im, fgp_gpu) +pars['rmse'] = rms +txtstr = printParametersToString(pars) +txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) +print (txtstr) +a=fig.add_subplot(2,4,4) + +# these are matplotlib.patch.Patch properties +props = dict(boxstyle='round', facecolor='wheat', alpha=0.75) +# place a text box in upper left in axes coords +a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=12, + verticalalignment='top', bbox=props) +imgplot = plt.imshow(fgp_gpu, cmap="gray") + +a=fig.add_subplot(2,4,8) + +# 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 +a.text(0.05, 0.95, 'fgp_gpu - u0', transform=a.transAxes, fontsize=14, + verticalalignment='top', bbox=props) +imgplot = plt.imshow((fgp_gpu - u0)**2, vmin=0, vmax=0.03, cmap="gray") +plt.colorbar(ticks=[0, 0.03], orientation='vertical') +plt.show() +""" |