From d4d12945f5396b10259fdeeb2b09f99b0e2c6afd Mon Sep 17 00:00:00 2001 From: algol Date: Thu, 22 Feb 2018 12:34:39 +0000 Subject: GPU test fixed, cpu vs gpu test added for ROF-TV --- Wrappers/Python/demo/test_cpu_regularizers.py | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) (limited to 'Wrappers/Python/demo') diff --git a/Wrappers/Python/demo/test_cpu_regularizers.py b/Wrappers/Python/demo/test_cpu_regularizers.py index 5908c3c..d147b85 100644 --- a/Wrappers/Python/demo/test_cpu_regularizers.py +++ b/Wrappers/Python/demo/test_cpu_regularizers.py @@ -284,9 +284,9 @@ start_time = timeit.default_timer() pars = {'algorithm': ROF_TV , \ 'input' : u0,\ - 'regularization_parameter':1,\ - 'marching_step': 0.003,\ - 'number_of_iterations': 300 + 'regularization_parameter':25,\ + 'marching_step': 0.001,\ + 'number_of_iterations': 350 } rof = ROF_TV(pars['input'], pars['number_of_iterations'], @@ -307,9 +307,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() -- cgit v1.2.3 From 0ebf1f949b7150692e356627c455e905b642af97 Mon Sep 17 00:00:00 2001 From: Daniil Kazantsev Date: Sun, 25 Feb 2018 17:54:03 +0000 Subject: updates to CPU/GPU core for ROF and demos --- Core/regularizers_CPU/ROF_TV_core.c | 9 ++++----- Core/regularizers_GPU/TV_ROF/TV_ROF_GPU.cu | 19 +++++++++---------- Wrappers/Python/demo/test_cpu_regularizers.py | 6 +++--- Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py | 6 +++--- Wrappers/Python/test/test_gpu_regularizers.py | 6 +++--- 5 files changed, 22 insertions(+), 24 deletions(-) (limited to 'Wrappers/Python/demo') diff --git a/Core/regularizers_CPU/ROF_TV_core.c b/Core/regularizers_CPU/ROF_TV_core.c index ba7fe48..fd47c3f 100644 --- a/Core/regularizers_CPU/ROF_TV_core.c +++ b/Core/regularizers_CPU/ROF_TV_core.c @@ -268,7 +268,7 @@ float TV_kernel(float *D1, float *D2, float *D3, float *B, float *A, float lambd dv2 = D2[index] - D2[(dimX*dimY)*k + j*dimX+i2]; dv3 = D3[index] - D3[(dimX*dimY)*k2 + j*dimX+i]; - B[index] = B[index] + tau*lambda*(dv1 + dv2 + dv3) + tau*(A[index] - B[index]); + B[index] = B[index] + tau*(2.0f*lambda*(dv1 + dv2 + dv3) - (B[index] - A[index])); }}} } else { @@ -284,10 +284,9 @@ float TV_kernel(float *D1, float *D2, float *D3, float *B, float *A, float lambd /* divergence components */ dv1 = D1[index] - D1[j2*dimX + i]; - dv2 = D2[index] - D2[j*dimX + i2]; - - //B[index] = B[index] + tau*lambda*(dv1 + dv2) + tau*(A[index] - B[index]); - B[index] = B[index] + tau*((dv1 + dv2) - lambda*(B[index] - A[index])); + dv2 = D2[index] - D2[j*dimX + i2]; + + B[index] = B[index] + tau*(2.0f*lambda*(dv1 + dv2) - (B[index] - A[index])); }} } return *B; diff --git a/Core/regularizers_GPU/TV_ROF/TV_ROF_GPU.cu b/Core/regularizers_GPU/TV_ROF/TV_ROF_GPU.cu index 897b5d0..b67b53b 100755 --- a/Core/regularizers_GPU/TV_ROF/TV_ROF_GPU.cu +++ b/Core/regularizers_GPU/TV_ROF/TV_ROF_GPU.cu @@ -54,7 +54,7 @@ limitations under the License. #define BLKXSIZE2D 16 #define BLKYSIZE2D 16 -#define EPS 1.0e-4 +#define EPS 1.0e-5 #define idivup(a, b) ( ((a)%(b) != 0) ? (a)/(b)+1 : (a)/(b) ) @@ -71,7 +71,7 @@ __host__ __device__ int sign (float x) /* differences 1 */ __global__ void D1_func2D(float* Input, float* D1, int N, int M) { - int i1, j1, i2, j2; + int i1, j1, i2; float NOMx_1,NOMy_1,NOMy_0,denom1,denom2,T1; int i = blockDim.x * blockIdx.x + threadIdx.x; int j = blockDim.y * blockIdx.y + threadIdx.y; @@ -84,7 +84,6 @@ __host__ __device__ int sign (float x) i1 = i + 1; if (i1 >= N) i1 = i-1; i2 = i - 1; if (i2 < 0) i2 = i+1; j1 = j + 1; if (j1 >= M) j1 = j-1; - j2 = j - 1; if (j2 < 0) j2 = j+1; /* Forward-backward differences */ NOMx_1 = Input[j1*N + i] - Input[index]; /* x+ */ @@ -102,7 +101,7 @@ __host__ __device__ int sign (float x) /* differences 2 */ __global__ void D2_func2D(float* Input, float* D2, int N, int M) { - int i1, j1, i2, j2; + int i1, j1, j2; float NOMx_1,NOMy_1,NOMx_0,denom1,denom2,T2; int i = blockDim.x * blockIdx.x + threadIdx.x; int j = blockDim.y * blockIdx.y + threadIdx.y; @@ -113,7 +112,6 @@ __host__ __device__ int sign (float x) /* boundary conditions (Neumann reflections) */ i1 = i + 1; if (i1 >= N) i1 = i-1; - i2 = i - 1; if (i2 < 0) i2 = i+1; j1 = j + 1; if (j1 >= M) j1 = j-1; j2 = j - 1; if (j2 < 0) j2 = j+1; @@ -149,7 +147,7 @@ __host__ __device__ int sign (float x) dv1 = D1[index] - D1[j2*N + i]; dv2 = D2[index] - D2[j*N + i2]; - Update[index] = Update[index] + tau*((dv1 + dv2) - lambda*(Update[index] - Input[index])); + Update[index] = Update[index] + tau*(2.0f*lambda*(dv1 + dv2) - (Update[index] - Input[index])); } } @@ -299,21 +297,22 @@ __host__ __device__ int sign (float x) dv2 = D2[index] - D2[(dimX*dimY)*k + j*dimX+i2]; dv3 = D3[index] - D3[(dimX*dimY)*k2 + j*dimX+i]; - Update[index] = Update[index] + tau*lambda*(dv1 + dv2 + dv3) + tau*(Update[index] - Input[index]); + Update[index] = Update[index] + tau*(2.0f*lambda*(dv1 + dv2 + dv3) - (Update[index] - Input[index])); } } ///////////////////////////////////////////////// // HOST FUNCTION -extern "C" void TV_ROF_GPU_kernel(float* Input, float* Output, int N, int M, int Z, int iter, float tau, float lambda) +extern "C" void TV_ROF_GPU(float* Input, float* Output, int N, int M, int Z, int iter, float tau, float lambda) { // set up device int dev = 0; CHECK(cudaSetDevice(dev)); - float *d_input, *d_update, *d_D1, *d_D2; + float *d_input, *d_update, *d_D1, *d_D2; + if (Z == 0) Z = 1; CHECK(cudaMalloc((void**)&d_input,N*M*Z*sizeof(float))); CHECK(cudaMalloc((void**)&d_update,N*M*Z*sizeof(float))); CHECK(cudaMalloc((void**)&d_D1,N*M*Z*sizeof(float))); @@ -346,7 +345,7 @@ extern "C" void TV_ROF_GPU_kernel(float* Input, float* Output, int N, int M, int CHECK(cudaFree(d_D3)); } else { - // TV - 2D case + // TV - 2D case dim3 dimBlock(BLKXSIZE2D,BLKYSIZE2D); dim3 dimGrid(idivup(N,BLKXSIZE2D), idivup(M,BLKYSIZE2D)); diff --git a/Wrappers/Python/demo/test_cpu_regularizers.py b/Wrappers/Python/demo/test_cpu_regularizers.py index d147b85..b08c418 100644 --- a/Wrappers/Python/demo/test_cpu_regularizers.py +++ b/Wrappers/Python/demo/test_cpu_regularizers.py @@ -284,9 +284,9 @@ start_time = timeit.default_timer() pars = {'algorithm': ROF_TV , \ 'input' : u0,\ - 'regularization_parameter':25,\ - 'marching_step': 0.001,\ - 'number_of_iterations': 350 + 'regularization_parameter':0.04,\ + 'marching_step': 0.0025,\ + 'number_of_iterations': 300 } rof = ROF_TV(pars['input'], pars['number_of_iterations'], diff --git a/Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py b/Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py index 8c91c73..d742a7f 100644 --- a/Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py +++ b/Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py @@ -58,8 +58,8 @@ imgplot = plt.imshow(u0,cmap="gray") # set parameters pars = {'algorithm': ROF_TV , \ 'input' : u0,\ - 'regularization_parameter':12,\ - 'time_marching_parameter': 0.001,\ + 'regularization_parameter':0.04,\ + 'time_marching_parameter': 0.0025,\ 'number_of_iterations': 600 } print ("#################ROF TV CPU#####################") @@ -120,4 +120,4 @@ plt.title('{}'.format('Pixels larger threshold difference')) if (diff_im.sum() > 1): print ("Arrays do not match!") else: - print ("Arrays match") \ No newline at end of file + print ("Arrays match") diff --git a/Wrappers/Python/test/test_gpu_regularizers.py b/Wrappers/Python/test/test_gpu_regularizers.py index 29f5bad..c473270 100644 --- a/Wrappers/Python/test/test_gpu_regularizers.py +++ b/Wrappers/Python/test/test_gpu_regularizers.py @@ -154,9 +154,9 @@ start_time = timeit.default_timer() pars = { 'algorithm' : GPU_ROF_TV , \ 'input' : u0, - 'regularization_parameter': 25,\ - 'time_marching_parameter': 0.001, \ - 'number_of_iterations':350 + 'regularization_parameter': 0.04,\ + 'time_marching_parameter': 0.0025, \ + 'number_of_iterations':300 } rof_tv = GPU_ROF_TV(pars['input'], -- cgit v1.2.3 From 74ff5b5f319b2f7ea3e078c62041ec8a1bb28335 Mon Sep 17 00:00:00 2001 From: algol Date: Mon, 5 Mar 2018 18:12:01 +0000 Subject: Cmake/Cython fixes to compile ROF-FGP --- Core/CMakeLists.txt | 3 ++- Core/regularizers_CPU/FGP_TV_core.c | 12 +++++------ Core/regularizers_CPU/FGP_TV_core.h | 2 +- Core/regularizers_CPU/ROF_TV_core.c | 4 ++-- Core/regularizers_CPU/ROF_TV_core.h | 2 +- Core/regularizers_GPU/TV_FGP/TV_FGP_GPU_core.cu | 12 +++++------ Core/regularizers_GPU/TV_FGP/TV_FGP_GPU_core.h | 2 +- Wrappers/Python/demo/test_cpu_regularizers.py | 25 ++++++++++++---------- Wrappers/Python/src/cpu_regularizers.cpp | 1 - Wrappers/Python/src/cpu_regularizers.pyx | 25 +++++++++++----------- Wrappers/Python/src/gpu_regularizers.pyx | 2 +- .../Python/test/test_cpu_vs_gpu_regularizers.py | 9 ++++---- 12 files changed, 51 insertions(+), 48 deletions(-) (limited to 'Wrappers/Python/demo') diff --git a/Core/CMakeLists.txt b/Core/CMakeLists.txt index 3e3f89e..68b7111 100644 --- a/Core/CMakeLists.txt +++ b/Core/CMakeLists.txt @@ -129,7 +129,8 @@ if (CUDA_FOUND) CUDA_ADD_LIBRARY(cilregcuda SHARED ${CMAKE_CURRENT_SOURCE_DIR}/regularizers_GPU/Diffus_HO/Diff4th_GPU_kernel.cu ${CMAKE_CURRENT_SOURCE_DIR}/regularizers_GPU/NL_Regul/NLM_GPU_kernel.cu - ${CMAKE_CURRENT_SOURCE_DIR}/regularizers_GPU/TV_ROF/TV_ROF_GPU.cu + ${CMAKE_CURRENT_SOURCE_DIR}/regularizers_GPU/TV_ROF/TV_ROF_GPU_core.cu + ${CMAKE_CURRENT_SOURCE_DIR}/regularizers_GPU/TV_FGP/TV_FGP_GPU_core.cu ) if (UNIX) message ("I'd install into ${CMAKE_INSTALL_PREFIX}/lib") diff --git a/Core/regularizers_CPU/FGP_TV_core.c b/Core/regularizers_CPU/FGP_TV_core.c index 2f1439d..5365631 100644 --- a/Core/regularizers_CPU/FGP_TV_core.c +++ b/Core/regularizers_CPU/FGP_TV_core.c @@ -23,7 +23,7 @@ limitations under the License. * * Input Parameters: * 1. Noisy image/volume - * 2. lambda - regularization parameter + * 2. lambdaPar - regularization parameter * 3. Number of iterations * 4. eplsilon: tolerance constant * 5. TV-type: methodTV - 'iso' (0) or 'l1' (1) @@ -37,7 +37,7 @@ limitations under the License. * [1] Amir Beck and Marc Teboulle, "Fast Gradient-Based Algorithms for Constrained Total Variation Image Denoising and Deblurring Problems" */ -float TV_FGP_CPU(float *Input, float *Output, float lambda, int iter, float epsil, int methodTV, int nonneg, int printM, int dimX, int dimY, int dimZ) +float TV_FGP_CPU_main(float *Input, float *Output, float lambdaPar, int iter, float epsil, int methodTV, int nonneg, int printM, int dimX, int dimY, int dimZ) { int ll, j, DimTotal; float re, re1; @@ -62,13 +62,13 @@ float TV_FGP_CPU(float *Input, float *Output, float lambda, int iter, float epsi for(ll=0; ll (y)) ? (x) : (y)) #define MIN(x, y) (((x) < (y)) ? (x) : (y)) @@ -46,7 +46,7 @@ int sign(float x) { */ /* Running iterations of TV-ROF function */ -float TV_ROF_CPU(float *Input, float *Output, int dimX, int dimY, int dimZ, int iterationsNumb, float tau, float lambda) +float TV_ROF_CPU_main(float *Input, float *Output, int dimX, int dimY, int dimZ, int iterationsNumb, float tau, float lambda) { float *D1, *D2, *D3; int i, DimTotal; diff --git a/Core/regularizers_CPU/ROF_TV_core.h b/Core/regularizers_CPU/ROF_TV_core.h index b32d0d5..63c0419 100644 --- a/Core/regularizers_CPU/ROF_TV_core.h +++ b/Core/regularizers_CPU/ROF_TV_core.h @@ -47,7 +47,7 @@ limitations under the License. extern "C" { #endif CCPI_EXPORT float TV_kernel(float *D1, float *D2, float *D3, float *B, float *A, float lambda, float tau, int dimY, int dimX, int dimZ); -CCPI_EXPORT float TV_ROF_CPU(float *Input, float *Output, int dimX, int dimY, int dimZ, int iterationsNumb, float tau, float lambda); +CCPI_EXPORT float TV_ROF_CPU_main(float *Input, float *Output, int dimX, int dimY, int dimZ, int iterationsNumb, float tau, float lambda); CCPI_EXPORT float D1_func(float *A, float *D1, int dimY, int dimX, int dimZ); CCPI_EXPORT float D2_func(float *A, float *D2, int dimY, int dimX, int dimZ); CCPI_EXPORT float D3_func(float *A, float *D3, int dimY, int dimX, int dimZ); diff --git a/Core/regularizers_GPU/TV_FGP/TV_FGP_GPU_core.cu b/Core/regularizers_GPU/TV_FGP/TV_FGP_GPU_core.cu index 0533a85..cbe9b5e 100755 --- a/Core/regularizers_GPU/TV_FGP/TV_FGP_GPU_core.cu +++ b/Core/regularizers_GPU/TV_FGP/TV_FGP_GPU_core.cu @@ -25,7 +25,7 @@ limitations under the License. * * Input Parameters: * 1. Noisy image/volume - * 2. lambda - regularization parameter + * 2. lambdaPar - regularization parameter * 3. Number of iterations * 4. eplsilon: tolerance constant * 5. TV-type: methodTV - 'iso' (0) or 'l1' (1) @@ -341,7 +341,7 @@ __global__ void copy_kernel3D(float *Input, float* Output, int N, int M, int Z, /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/ ////////////MAIN HOST FUNCTION /////////////// -extern "C" void TV_FGP_GPU(float *Input, float *Output, float lambda, int iter, float epsil, int methodTV, int nonneg, int printM, int dimX, int dimY, int dimZ) +extern "C" void TV_FGP_GPU(float *Input, float *Output, float lambdaPar, int iter, float epsil, int methodTV, int nonneg, int printM, int dimX, int dimY, int dimZ) { int deviceCount = -1; // number of devices cudaGetDeviceCount(&deviceCount); @@ -383,13 +383,13 @@ extern "C" void TV_FGP_GPU(float *Input, float *Output, float lambda, int iter, cudaMemset(R2, 0, ImSize*sizeof(float)); /********************** Run CUDA 2D kernel here ********************/ - multip = (1.0f/(8.0f*lambda)); + multip = (1.0f/(8.0f*lambdaPar)); /* The main kernel */ for (i = 0; i < iter; i++) { /* computing the gradient of the objective function */ - Obj_func2D_kernel<<>>(d_input, d_update, R1, R2, dimX, dimY, ImSize, lambda); + Obj_func2D_kernel<<>>(d_input, d_update, R1, R2, dimX, dimY, ImSize, lambdaPar); checkCudaErrors( cudaDeviceSynchronize() ); checkCudaErrors(cudaPeekAtLastError() ); @@ -493,13 +493,13 @@ extern "C" void TV_FGP_GPU(float *Input, float *Output, float lambda, int iter, cudaMemset(R2, 0, ImSize*sizeof(float)); cudaMemset(R3, 0, ImSize*sizeof(float)); /********************** Run CUDA 3D kernel here ********************/ - multip = (1.0f/(8.0f*lambda)); + multip = (1.0f/(8.0f*lambdaPar)); /* The main kernel */ for (i = 0; i < iter; i++) { /* computing the gradient of the objective function */ - Obj_func3D_kernel<<>>(d_input, d_update, R1, R2, R3, dimX, dimY, dimZ, ImSize, lambda); + Obj_func3D_kernel<<>>(d_input, d_update, R1, R2, R3, dimX, dimY, dimZ, ImSize, lambdaPar); checkCudaErrors( cudaDeviceSynchronize() ); checkCudaErrors(cudaPeekAtLastError() ); diff --git a/Core/regularizers_GPU/TV_FGP/TV_FGP_GPU_core.h b/Core/regularizers_GPU/TV_FGP/TV_FGP_GPU_core.h index 15c7120..2b8fdcf 100755 --- a/Core/regularizers_GPU/TV_FGP/TV_FGP_GPU_core.h +++ b/Core/regularizers_GPU/TV_FGP/TV_FGP_GPU_core.h @@ -5,6 +5,6 @@ #ifndef _TV_FGP_GPU_ #define _TV_FGP_GPU_ -extern "C" void TV_FGP_GPU(float *Input, float *Output, float lambda, int iter, float epsil, int methodTV, int nonneg, int printM, int dimX, int dimY, int dimZ); +extern "C" void TV_FGP_GPU(float *Input, float *Output, float lambdaPar, int iter, float epsil, int methodTV, int nonneg, int printM, int dimX, int dimY, int dimZ); #endif diff --git a/Wrappers/Python/demo/test_cpu_regularizers.py b/Wrappers/Python/demo/test_cpu_regularizers.py index b08c418..53b8538 100644 --- a/Wrappers/Python/demo/test_cpu_regularizers.py +++ b/Wrappers/Python/demo/test_cpu_regularizers.py @@ -11,10 +11,9 @@ 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 ,\ +from ccpi.filters.cpu_regularizers_boost import SplitBregman_TV, LLT_model, PatchBased_Regul ,\ TGV_PD -from ccpi.filters.cpu_regularizers_cython import ROF_TV +from ccpi.filters.cpu_regularizers_cython import TV_ROF_CPU, TV_FGP_CPU ############################################################################### #https://stackoverflow.com/questions/13875989/comparing-image-in-url-to-image-in-filesystem-in-python/13884956#13884956 @@ -128,21 +127,25 @@ 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 , \ +pars = {'algorithm' : TV_FGP_CPU , \ 'input' : u0, 'regularization_parameter':0.05, \ 'number_of_iterations' :200 ,\ - 'tolerance_constant':1e-4,\ - 'TV_penalty': 0 + 'tolerance_constant':1e-5,\ + 'methodTV': 0 ,\ + 'nonneg': 0 ,\ + 'printingOut': 0 } -out = FGP_TV (pars['input'], +out = TV_FGP_CPU (pars['input'], pars['regularization_parameter'], pars['number_of_iterations'], pars['tolerance_constant'], - pars['TV_penalty']) + pars['methodTV'], + pars['nonneg'], + pars['printingOut']) fgp = out[0] rms = rmse(Im, fgp) @@ -282,13 +285,13 @@ imgplot = plt.imshow(tgv, cmap="gray") start_time = timeit.default_timer() -pars = {'algorithm': ROF_TV , \ +pars = {'algorithm': TV_ROF_CPU , \ 'input' : u0,\ 'regularization_parameter':0.04,\ 'marching_step': 0.0025,\ 'number_of_iterations': 300 } -rof = ROF_TV(pars['input'], +rof = TV_ROF_CPU(pars['input'], pars['number_of_iterations'], pars['regularization_parameter'], pars['marching_step'] diff --git a/Wrappers/Python/src/cpu_regularizers.cpp b/Wrappers/Python/src/cpu_regularizers.cpp index 43d5d11..b8156d5 100644 --- a/Wrappers/Python/src/cpu_regularizers.cpp +++ b/Wrappers/Python/src/cpu_regularizers.cpp @@ -1040,7 +1040,6 @@ BOOST_PYTHON_MODULE(cpu_regularizers_boost) np::dtype dt2 = np::dtype::get_builtin(); 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 b8089a8..448da31 100644 --- a/Wrappers/Python/src/cpu_regularizers.pyx +++ b/Wrappers/Python/src/cpu_regularizers.pyx @@ -18,8 +18,8 @@ import cython import numpy as np cimport numpy as np -cdef extern float TV_ROF_CPU(float *Input, float *Output, int dimX, int dimY, int dimZ, int iterationsNumb, float tau, float flambda); -cdef extern float TV_FGP_CPU(float *Input, float *Output, float lambda, int iter, float epsil, int methodTV, int nonneg, int printM, int dimX, int dimY, int dimZ); +cdef extern float TV_ROF_CPU_main(float *Input, float *Output, int dimX, int dimY, int dimZ, int iterationsNumb, float tau, float flambda); +cdef extern float TV_FGP_CPU_main(float *Input, float *Output, float lambdaPar, int iter, float epsil, int methodTV, int nonneg, int printM, int dimX, int dimY, int dimZ); # Can we use the same name here in "def" as the C function? def TV_ROF_CPU(inputData, iterations, regularization_parameter, @@ -45,11 +45,10 @@ def TV_ROF_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, np.zeros([dims[0],dims[1]], dtype='float32') #/* Run ROF iterations for 2D data */ - TV_ROF_CPU(&inputData[0,0], &B[0,0], dims[0], dims[1], 1, iterations, + TV_ROF_CPU_main(&inputData[0,0], &B[0,0], dims[0], dims[1], 1, iterations, marching_step_parameter, regularization_parameter) - return B - + return B def TV_ROF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, int iterations, @@ -65,7 +64,7 @@ def TV_ROF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, np.zeros([dims[0],dims[1],dims[2]], dtype='float32') #/* Run ROF iterations for 3D data */ - TV_FGP_CPU(&inputData[0,0,0], &B[0,0,0], dims[0], dims[1], dims[2], iterations, + TV_ROF_CPU_main(&inputData[0,0,0], &B[0,0,0], dims[0], dims[1], dims[2], iterations, marching_step_parameter, regularization_parameter) return B @@ -88,11 +87,11 @@ def TV_FGP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 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_FGP_CPU(&inputData[0,0], &B[0,0], regularization_parameter, + TV_FGP_CPU_main(&inputData[0,0], &outputData[0,0], regularization_parameter, iterations, tolerance_param, methodTV, @@ -100,8 +99,7 @@ def TV_FGP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, printM, dims[0], dims[1], 1) - return B - + return outputData def TV_FGP_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, float regularization_parameter, @@ -115,15 +113,16 @@ def TV_FGP_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 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') #/* Run ROF iterations for 3D data */ - TV_FGP_CPU(&inputData[0,0, 0], &B[0,0, 0], iterations, + TV_FGP_CPU_main(&inputData[0,0,0], &outputData[0,0, 0], regularization_parameter, + iterations, tolerance_param, methodTV, nonneg, printM, dims[0], dims[1], dims[2]) - return B + return outputData diff --git a/Wrappers/Python/src/gpu_regularizers.pyx b/Wrappers/Python/src/gpu_regularizers.pyx index a14a20d..e99bfa7 100644 --- a/Wrappers/Python/src/gpu_regularizers.pyx +++ b/Wrappers/Python/src/gpu_regularizers.pyx @@ -26,7 +26,7 @@ cdef extern void NLM_GPU_kernel(float *A, float* B, float *Eucl_Vec, int SearchW, int SimilW, int SearchW_real, float denh2, float lambdaf); cdef extern void TV_ROF_GPU(float* Input, float* Output, int N, int M, int Z, int iter, float tau, float lambdaf); -cdef extern void TV_FGP_GPU(float *Input, float *Output, float lambda, int iter, float epsil, int methodTV, int nonneg, int printM, int N, int M, int Z); +cdef extern void TV_FGP_GPU(float *Input, float *Output, float lambdaPar, int iter, float epsil, int methodTV, int nonneg, int printM, int N, int M, int Z); cdef extern float pad_crop(float *A, float *Ap, int OldSizeX, int OldSizeY, int OldSizeZ, diff --git a/Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py b/Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py index 6344021..e162afa 100644 --- a/Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py +++ b/Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py @@ -12,7 +12,7 @@ import matplotlib.pyplot as plt import numpy as np import os import timeit -from ccpi.filters.gpu_regularizers import Diff4thHajiaboli, NML, TV_ROF_GPU +from ccpi.filters.gpu_regularizers import Diff4thHajiaboli, NML, GPU_ROF_TV from ccpi.filters.cpu_regularizers_cython import TV_ROF_CPU ############################################################################### def printParametersToString(pars): @@ -56,11 +56,11 @@ imgplot = plt.imshow(u0,cmap="gray") # set parameters -pars = {'algorithm': ROF_TV , \ +pars = {'algorithm': TV_ROF_CPU , \ 'input' : u0,\ 'regularization_parameter':0.04,\ 'time_marching_parameter': 0.0025,\ - 'number_of_iterations': 600 + 'number_of_iterations': 1200 } print ("#################ROF TV CPU#####################") start_time = timeit.default_timer() @@ -89,13 +89,14 @@ plt.title('{}'.format('CPU results')) print ("#################ROF TV GPU#####################") start_time = timeit.default_timer() -rof_gpu = TV_ROF_GPU(pars['input'], +rof_gpu = GPU_ROF_TV(pars['input'], pars['number_of_iterations'], pars['time_marching_parameter'], pars['regularization_parameter']) rms = rmse(Im, rof_gpu) pars['rmse'] = rms +pars['algorithm'] = GPU_ROF_TV txtstr = printParametersToString(pars) txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) print (txtstr) -- cgit v1.2.3 From 8d310478254f3cda63f3663729b416f425ad70b6 Mon Sep 17 00:00:00 2001 From: algol Date: Tue, 6 Mar 2018 11:45:53 +0000 Subject: work on FGP intergration --- Core/CMakeLists.txt | 2 +- Core/regularizers_CPU/ROF_TV_core.c | 22 +++++++++++----------- Core/regularizers_CPU/utils.c | 8 ++++---- Wrappers/Python/demo/test_cpu_regularizers.py | 17 +++++++++-------- Wrappers/Python/src/cpu_regularizers.pyx | 14 +++++--------- 5 files changed, 30 insertions(+), 33 deletions(-) (limited to 'Wrappers/Python/demo') diff --git a/Core/CMakeLists.txt b/Core/CMakeLists.txt index 68b7111..4e85002 100644 --- a/Core/CMakeLists.txt +++ b/Core/CMakeLists.txt @@ -88,7 +88,7 @@ add_library(cilreg SHARED ${CMAKE_CURRENT_SOURCE_DIR}/regularizers_CPU/PatchBased_Regul_core.c ${CMAKE_CURRENT_SOURCE_DIR}/regularizers_CPU/SplitBregman_TV_core.c ${CMAKE_CURRENT_SOURCE_DIR}/regularizers_CPU/TGV_PD_core.c - ${CMAKE_CURRENT_SOURCE_DIR}/regularizers_CPU/ROF_TV_core.c + ${CMAKE_CURRENT_SOURCE_DIR}/regularizers_CPU/ROF_TV_core.c ${CMAKE_CURRENT_SOURCE_DIR}/regularizers_CPU/utils.c ) target_link_libraries(cilreg ${EXTRA_LIBRARIES} ) diff --git a/Core/regularizers_CPU/ROF_TV_core.c b/Core/regularizers_CPU/ROF_TV_core.c index a59a3c6..a4f82a6 100644 --- a/Core/regularizers_CPU/ROF_TV_core.c +++ b/Core/regularizers_CPU/ROF_TV_core.c @@ -46,7 +46,7 @@ int sign(float x) { */ /* Running iterations of TV-ROF function */ -float TV_ROF_CPU_mainfloat TV_ROF_CPU_main(float *Input, float *Output, float lambdaPar, int iterationsNumb, float tau, int dimX, int dimY, int dimZ) +float TV_ROF_CPU_main(float *Input, float *Output, float lambdaPar, int iterationsNumb, float tau, int dimX, int dimY, int dimZ) { float *D1, *D2, *D3; int i, DimTotal; @@ -132,9 +132,9 @@ float D1_func(float *A, float *D1, int dimX, int dimY, int dimZ) NOMy_0 = A[index] - A[(j)*dimX + i2]; /* y- */ denom1 = NOMx_1*NOMx_1; - denom2 = 0.5*(sign(NOMy_1) + sign(NOMy_0))*(MIN(fabs(NOMy_1),fabs(NOMy_0))); + denom2 = 0.5f*(sign(NOMy_1) + sign(NOMy_0))*(MIN(fabs(NOMy_1),fabs(NOMy_0))); denom2 = denom2*denom2; - T1 = sqrt(denom1 + denom2 + EPS); + T1 = sqrtf(denom1 + denom2 + EPS); D1[index] = NOMx_1/T1; }} } @@ -170,11 +170,11 @@ float D2_func(float *A, float *D2, int dimX, int dimY, int dimZ) denom1 = NOMy_1*NOMy_1; - denom2 = 0.5*(sign(NOMx_1) + sign(NOMx_0))*(MIN(fabs(NOMx_1),fabs(NOMx_0))); + denom2 = 0.5f*(sign(NOMx_1) + sign(NOMx_0))*(MIN(fabs(NOMx_1),fabs(NOMx_0))); denom2 = denom2*denom2; - denom3 = 0.5*(sign(NOMz_1) + sign(NOMz_0))*(MIN(fabs(NOMz_1),fabs(NOMz_0))); + denom3 = 0.5f*(sign(NOMz_1) + sign(NOMz_0))*(MIN(fabs(NOMz_1),fabs(NOMz_0))); denom3 = denom3*denom3; - T2 = sqrt(denom1 + denom2 + denom3 + EPS); + T2 = sqrtf(denom1 + denom2 + denom3 + EPS); D2[index] = NOMy_1/T2; }}} } @@ -196,9 +196,9 @@ float D2_func(float *A, float *D2, int dimX, int dimY, int dimZ) /*NOMy_0 = A[(i)*dimY + j] - A[(i)*dimY + j2]; */ /* y- */ denom1 = NOMy_1*NOMy_1; - denom2 = 0.5*(sign(NOMx_1) + sign(NOMx_0))*(MIN(fabs(NOMx_1),fabs(NOMx_0))); + denom2 = 0.5f*(sign(NOMx_1) + sign(NOMx_0))*(MIN(fabs(NOMx_1),fabs(NOMx_0))); denom2 = denom2*denom2; - T2 = sqrt(denom1 + denom2 + EPS); + T2 = sqrtf(denom1 + denom2 + EPS); D2[index] = NOMy_1/T2; }} } @@ -233,11 +233,11 @@ float D3_func(float *A, float *D3, int dimY, int dimX, int dimZ) /*NOMz_0 = A[(dimX*dimY)*k + (i)*dimY + j] - A[(dimX*dimY)*k2 + (i)*dimY + j]; */ /* z- */ denom1 = NOMz_1*NOMz_1; - denom2 = 0.5*(sign(NOMx_1) + sign(NOMx_0))*(MIN(fabs(NOMx_1),fabs(NOMx_0))); + denom2 = 0.5f*(sign(NOMx_1) + sign(NOMx_0))*(MIN(fabs(NOMx_1),fabs(NOMx_0))); denom2 = denom2*denom2; - denom3 = 0.5*(sign(NOMy_1) + sign(NOMy_0))*(MIN(fabs(NOMy_1),fabs(NOMy_0))); + denom3 = 0.5f*(sign(NOMy_1) + sign(NOMy_0))*(MIN(fabs(NOMy_1),fabs(NOMy_0))); denom3 = denom3*denom3; - T3 = sqrt(denom1 + denom2 + denom3 + EPS); + T3 = sqrtf(denom1 + denom2 + denom3 + EPS); D3[index] = NOMz_1/T3; }}} return *D3; diff --git a/Core/regularizers_CPU/utils.c b/Core/regularizers_CPU/utils.c index 951fb91..cdf3d0e 100644 --- a/Core/regularizers_CPU/utils.c +++ b/Core/regularizers_CPU/utils.c @@ -45,10 +45,10 @@ float TV_energy2D(float *U, float *U0, float *E_val, float lambda, int dimX, int j1 = j + 1; if (j == dimY-1) j1 = j; /* Forward differences */ - NOMx_2 = pow(U[j1*dimX + i] - U[index],2); /* x+ */ - NOMy_2 = pow(U[j*dimX + i1] - U[index],2); /* y+ */ - E_Grad += sqrt(NOMx_2 + NOMy_2); /* gradient term energy */ - E_Data += 0.5f * lambda*(pow((U[index]-U0[index]),2)); /* fidelity term energy */ + NOMx_2 = powf(U[j1*dimX + i] - U[index],2); /* x+ */ + NOMy_2 = powf(U[j*dimX + i1] - U[index],2); /* y+ */ + E_Grad += sqrtf(NOMx_2 + NOMy_2); /* gradient term energy */ + E_Data += 0.5f * lambda*(powf((U[index]-U0[index]),2)); /* fidelity term energy */ } } E_val[0] = E_Grad + E_Data; diff --git a/Wrappers/Python/demo/test_cpu_regularizers.py b/Wrappers/Python/demo/test_cpu_regularizers.py index 53b8538..f1eb3c3 100644 --- a/Wrappers/Python/demo/test_cpu_regularizers.py +++ b/Wrappers/Python/demo/test_cpu_regularizers.py @@ -131,9 +131,9 @@ imgplot = plt.imshow(splitbregman,\ start_time = timeit.default_timer() pars = {'algorithm' : TV_FGP_CPU , \ 'input' : u0, - 'regularization_parameter':0.05, \ - 'number_of_iterations' :200 ,\ - 'tolerance_constant':1e-5,\ + 'regularization_parameter':0.07, \ + 'number_of_iterations' :300 ,\ + 'tolerance_constant':0.00001,\ 'methodTV': 0 ,\ 'nonneg': 0 ,\ 'printingOut': 0 @@ -156,7 +156,7 @@ txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) print (txtstr) -a=fig.add_subplot(2,4,3) +a=fig.add_subplot(2,4,4) # these are matplotlib.patch.Patch properties props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) @@ -168,8 +168,9 @@ 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 , \ @@ -204,7 +205,7 @@ 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: @@ -292,8 +293,8 @@ pars = {'algorithm': TV_ROF_CPU , \ 'number_of_iterations': 300 } rof = TV_ROF_CPU(pars['input'], - pars['number_of_iterations'], - pars['regularization_parameter'], + pars['regularization_parameter'], + pars['number_of_iterations'], pars['marching_step'] ) #tgv = out diff --git a/Wrappers/Python/src/cpu_regularizers.pyx b/Wrappers/Python/src/cpu_regularizers.pyx index 2654831..d62ca59 100644 --- a/Wrappers/Python/src/cpu_regularizers.pyx +++ b/Wrappers/Python/src/cpu_regularizers.pyx @@ -21,14 +21,11 @@ cimport numpy as np 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): +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) + 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) + 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, @@ -47,10 +44,9 @@ def TV_ROF_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, return outputData def TV_ROF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, - int iterations, + int iterationsNumb, float regularization_parameter, - float marching_step_parameter - ): + float marching_step_parameter): cdef long dims[3] dims[0] = inputData.shape[0] dims[1] = inputData.shape[1] -- cgit v1.2.3 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(-) (limited to 'Wrappers/Python/demo') 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 From 69ecdd57434d591eb3fa4afefb72174d3e025fb9 Mon Sep 17 00:00:00 2001 From: algol Date: Tue, 6 Mar 2018 12:48:21 +0000 Subject: FGP_CPU (Cythonized) now works in demo --- Wrappers/Python/demo/test_cpu_regularizers.py | 10 +++------- 1 file changed, 3 insertions(+), 7 deletions(-) (limited to 'Wrappers/Python/demo') diff --git a/Wrappers/Python/demo/test_cpu_regularizers.py b/Wrappers/Python/demo/test_cpu_regularizers.py index 7f08605..1d97857 100644 --- a/Wrappers/Python/demo/test_cpu_regularizers.py +++ b/Wrappers/Python/demo/test_cpu_regularizers.py @@ -137,7 +137,7 @@ pars = {'algorithm' : TV_FGP_CPU , \ 'printingOut': 0 } -out = TV_FGP_CPU (pars['input'], +fgp = TV_FGP_CPU(pars['input'], pars['regularization_parameter'], pars['number_of_iterations'], pars['tolerance_constant'], @@ -145,7 +145,6 @@ out = TV_FGP_CPU (pars['input'], pars['nonneg'], pars['printingOut']) -fgp = out[0] rms = rmse(Im, fgp) pars['rmse'] = rms @@ -154,7 +153,7 @@ txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) print (txtstr) -a=fig.add_subplot(2,4,4) +a=fig.add_subplot(2,4,3) # these are matplotlib.patch.Patch properties props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) @@ -168,7 +167,6 @@ a.text(0.05, 0.95, txtstr, transform=a.transAxes, fontsize=14, ###################### LLT_model ######################################### -""" start_time = timeit.default_timer() pars = {'algorithm': LLT_model , \ @@ -203,8 +201,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 @@ -286,7 +282,7 @@ start_time = timeit.default_timer() pars = {'algorithm': TV_ROF_CPU , \ 'input' : u0,\ - 'regularization_parameter':0.04,\ + 'regularization_parameter':0.07,\ 'marching_step': 0.0025,\ 'number_of_iterations': 300 } -- cgit v1.2.3 From 4593040d228fec5ef1f1a301e511708d47f0d55e Mon Sep 17 00:00:00 2001 From: algol Date: Tue, 6 Mar 2018 16:51:31 +0000 Subject: demos fixed #40 --- Wrappers/Python/demo/test_cpu_regularizers.py | 21 +- Wrappers/Python/src/cpu_regularizers.cpp | 289 --------------------- .../Python/test/test_cpu_vs_gpu_regularizers.py | 4 +- Wrappers/Python/test/test_gpu_regularizers.py | 59 ++++- 4 files changed, 63 insertions(+), 310 deletions(-) (limited to 'Wrappers/Python/demo') diff --git a/Wrappers/Python/demo/test_cpu_regularizers.py b/Wrappers/Python/demo/test_cpu_regularizers.py index 1d97857..76b9de7 100644 --- a/Wrappers/Python/demo/test_cpu_regularizers.py +++ b/Wrappers/Python/demo/test_cpu_regularizers.py @@ -11,10 +11,8 @@ import numpy as np import os from enum import Enum import timeit -from ccpi.filters.cpu_regularizers_boost import SplitBregman_TV, LLT_model, PatchBased_Regul ,\ - TGV_PD -from ccpi.filters.cpu_regularizers_cython import TV_ROF_CPU, TV_FGP_CPU - +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 @@ -127,7 +125,7 @@ imgplot = plt.imshow(splitbregman,\ ###################### FGP_TV ######################################### start_time = timeit.default_timer() -pars = {'algorithm' : TV_FGP_CPU , \ +pars = {'algorithm' : FGP_TV , \ 'input' : u0,\ 'regularization_parameter':0.07, \ 'number_of_iterations' :300 ,\ @@ -137,13 +135,13 @@ pars = {'algorithm' : TV_FGP_CPU , \ 'printingOut': 0 } -fgp = TV_FGP_CPU(pars['input'], +fgp = FGP_TV(pars['input'], pars['regularization_parameter'], pars['number_of_iterations'], pars['tolerance_constant'], pars['methodTV'], pars['nonneg'], - pars['printingOut']) + pars['printingOut'], 'cpu') rms = rmse(Im, fgp) pars['rmse'] = rms @@ -280,18 +278,17 @@ imgplot = plt.imshow(tgv, cmap="gray") start_time = timeit.default_timer() -pars = {'algorithm': TV_ROF_CPU , \ +pars = {'algorithm': ROF_TV , \ 'input' : u0,\ 'regularization_parameter':0.07,\ 'marching_step': 0.0025,\ 'number_of_iterations': 300 } -rof = TV_ROF_CPU(pars['input'], +rof = ROF_TV(pars['input'], pars['regularization_parameter'], pars['number_of_iterations'], - pars['marching_step'] - ) -#tgv = out + pars['marching_step'], 'cpu') + rms = rmse(Im, rof) pars['rmse'] = rms diff --git a/Wrappers/Python/src/cpu_regularizers.cpp b/Wrappers/Python/src/cpu_regularizers.cpp index b8156d5..88b7c3c 100644 --- a/Wrappers/Python/src/cpu_regularizers.cpp +++ b/Wrappers/Python/src/cpu_regularizers.cpp @@ -303,292 +303,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(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(); - //np::ndarray out1 = np::zeros(shape1, dtype); - - ////float *funcvalA = (float *)mxGetData(plhs[1]); - //float * funcvalA = reinterpret_cast(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(); - - - //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(npD.get_data()); - //D_old = reinterpret_cast(npD_old.get_data()); - //P1 = reinterpret_cast(npP1.get_data()); - //P2 = reinterpret_cast(npP2.get_data()); - //P1_old = reinterpret_cast(npP1_old.get_data()); - //P2_old = reinterpret_cast(npP2_old.get_data()); - //R1 = reinterpret_cast(npR1.get_data()); - //R2 = reinterpret_cast(npR2.get_data()); - - ///* begin iterations */ - //for (ll = 0; ll 4) { - //Obj_func2D(A, D, P1, P2, lambda, dimX, dimY); - //funcval = 0.0f; - //for (j = 0; j 2) { - //if (re > re_old) { - //Obj_func2D(A, D, P1, P2, lambda, dimX, dimY); - //funcval = 0.0f; - //for (j = 0; j(npD); - //result.append(out1); - //result.append(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(); - - //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(npD.get_data()); - //D_old = reinterpret_cast(npD_old.get_data()); - //P1 = reinterpret_cast(npP1.get_data()); - //P2 = reinterpret_cast(npP2.get_data()); - //P3 = reinterpret_cast(npP3.get_data()); - //P1_old = reinterpret_cast(npP1_old.get_data()); - //P2_old = reinterpret_cast(npP2_old.get_data()); - //P3_old = reinterpret_cast(npP3_old.get_data()); - //R1 = reinterpret_cast(npR1.get_data()); - //R2 = reinterpret_cast(npR2.get_data()); - //R3 = reinterpret_cast(npR3.get_data()); - ///* begin iterations */ - //for (ll = 0; ll 3) { - //Obj_func3D(A, D, P1, P2, P3, lambda, dimX, dimY, dimZ); - //funcval = 0.0f; - //for (j = 0; j 2) { - //if (re > re_old) { - //Obj_func3D(A, D, P1, P2, P3, lambda, dimX, dimY, dimZ); - //funcval = 0.0f; - //for (j = 0; j(npD); - //result.append(out1); - //result.append(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 +736,6 @@ bp::list TGV_PD(np::ndarray input, double d_lambda, double d_alpha1, double d_al result.append(npU); } - - - return result; } diff --git a/Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py b/Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py index a9d0f31..63be1a0 100644 --- a/Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py +++ b/Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py @@ -3,7 +3,7 @@ """ Created on Thu Feb 22 11:39:43 2018 -Testing CPU implementation against GPU one +Testing CPU implementation against the GPU one @author: Daniil Kazantsev """ @@ -12,8 +12,6 @@ import matplotlib.pyplot as plt import numpy as np import os import timeit -#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 from ccpi.filters.regularizers import ROF_TV, FGP_TV ############################################################################### diff --git a/Wrappers/Python/test/test_gpu_regularizers.py b/Wrappers/Python/test/test_gpu_regularizers.py index 04aeeb4..640b3f9 100644 --- a/Wrappers/Python/test/test_gpu_regularizers.py +++ b/Wrappers/Python/test/test_gpu_regularizers.py @@ -11,8 +11,8 @@ import numpy as np import os from enum import Enum import timeit -from ccpi.filters.gpu_regularizers import Diff4thHajiaboli, NML, TV_ROF_GPU - +from ccpi.filters.gpu_regularizers import Diff4thHajiaboli, NML +from ccpi.filters.regularizers import ROF_TV, FGP_TV ############################################################################### def printParametersToString(pars): txt = r'' @@ -151,20 +151,19 @@ plt.colorbar(ticks=[0, 0.03], orientation='vertical') ## Rudin-Osher-Fatemi (ROF) TV regularization start_time = timeit.default_timer() - pars = { -'algorithm' : TV_ROF_GPU , \ +'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']) + pars['time_marching_parameter'],'gpu') rms = rmse(Im, rof_tv) pars['rmse'] = rms @@ -190,3 +189,51 @@ a.text(0.05, 0.95, 'rof_tv - u0', transform=a.transAxes, fontsize=14, 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() +""" -- cgit v1.2.3