From 5e7b28053dfe06008657bcdb68462dc3d84b8a22 Mon Sep 17 00:00:00 2001 From: Daniil Kazantsev Date: Mon, 25 Nov 2019 23:18:25 +0000 Subject: added PD_TV_2D_CPU version --- build/run.sh | 6 +- demos/demo_cpu_regularisers.py | 51 +++++++- src/Core/CMakeLists.txt | 31 ++--- src/Core/regularisers_CPU/FGP_TV_core.c | 76 ++++------- src/Core/regularisers_CPU/FGP_TV_core.h | 11 +- src/Core/regularisers_CPU/PD_TV_core.c | 166 ++++++++++++++++++++++++ src/Core/regularisers_CPU/PD_TV_core.h | 57 ++++++++ src/Core/regularisers_CPU/utils.c | 46 ++++++- src/Core/regularisers_CPU/utils.h | 1 + src/Matlab/mex_compile/compileCPU_mex_Linux.m | 7 +- src/Matlab/mex_compile/regularisers_CPU/PD_TV.c | 98 ++++++++++++++ src/Python/ccpi/filters/regularisers.py | 28 +++- src/Python/setup-regularisers.py.in | 25 ++-- src/Python/src/cpu_regularisers.pyx | 42 +++++- 14 files changed, 541 insertions(+), 104 deletions(-) create mode 100644 src/Core/regularisers_CPU/PD_TV_core.c create mode 100644 src/Core/regularisers_CPU/PD_TV_core.h create mode 100644 src/Matlab/mex_compile/regularisers_CPU/PD_TV.c diff --git a/build/run.sh b/build/run.sh index 285476b..bbbbc6a 100755 --- a/build/run.sh +++ b/build/run.sh @@ -6,11 +6,11 @@ rm -r ../build_proj mkdir ../build_proj cd ../build_proj/ #make clean -export CIL_VERSION=19.09 +export CIL_VERSION=19.10 # install Python modules without CUDA -#cmake ../ -DBUILD_PYTHON_WRAPPER=ON -DBUILD_MATLAB_WRAPPER=OFF -DBUILD_CUDA=OFF -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=./install +cmake ../ -DBUILD_PYTHON_WRAPPER=ON -DBUILD_MATLAB_WRAPPER=OFF -DBUILD_CUDA=OFF -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=./install # install Python modules with CUDA -cmake ../ -DBUILD_PYTHON_WRAPPER=ON -DBUILD_MATLAB_WRAPPER=OFF -DBUILD_CUDA=ON -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=./install +# cmake ../ -DBUILD_PYTHON_WRAPPER=ON -DBUILD_MATLAB_WRAPPER=OFF -DBUILD_CUDA=ON -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=./install # install Matlab modules without CUDA # cmake ../ -DBUILD_PYTHON_WRAPPER=OFF -DMatlab_ROOT_DIR=/dls_sw/apps/matlab/r2014a/ -DBUILD_MATLAB_WRAPPER=ON -DBUILD_CUDA=OFF -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=./install # install Matlab modules with CUDA diff --git a/demos/demo_cpu_regularisers.py b/demos/demo_cpu_regularisers.py index 8655623..d0bbe63 100644 --- a/demos/demo_cpu_regularisers.py +++ b/demos/demo_cpu_regularisers.py @@ -12,7 +12,7 @@ import matplotlib.pyplot as plt import numpy as np import os import timeit -from ccpi.filters.regularisers import ROF_TV, FGP_TV, SB_TV, TGV, LLT_ROF, FGP_dTV, TNV, NDF, Diff4th +from ccpi.filters.regularisers import ROF_TV, FGP_TV, PD_TV, SB_TV, TGV, LLT_ROF, FGP_dTV, TNV, NDF, Diff4th from ccpi.filters.regularisers import PatchSelect, NLTV from ccpi.supp.qualitymetrics import QualityTools ############################################################################### @@ -129,7 +129,7 @@ imgplot = plt.imshow(u0,cmap="gray") pars = {'algorithm' : FGP_TV, \ 'input' : u0,\ 'regularisation_parameter':0.02, \ - 'number_of_iterations' :400 ,\ + 'number_of_iterations' :1500 ,\ 'tolerance_constant':1e-06,\ 'methodTV': 0 ,\ 'nonneg': 0} @@ -159,6 +159,53 @@ a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, imgplot = plt.imshow(fgp_cpu, cmap="gray") plt.title('{}'.format('CPU results')) +#%% +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +print ("_______________PD-TV (2D)__________________") +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") + +## plot +fig = plt.figure() +plt.suptitle('Performance of PD-TV regulariser using the CPU') +a=fig.add_subplot(1,2,1) +a.set_title('Noisy Image') +imgplot = plt.imshow(u0,cmap="gray") + +# set parameters +pars = {'algorithm' : PD_TV, \ + 'input' : u0,\ + 'regularisation_parameter':0.02, \ + 'number_of_iterations' :1500 ,\ + 'tolerance_constant':1e-06,\ + 'methodTV': 0 ,\ + 'nonneg': 1, + 'lipschitz_const' : 12} + +print ("#############PD TV CPU####################") +start_time = timeit.default_timer() +(pd_cpu,info_vec_cpu) = PD_TV(pars['input'], + pars['regularisation_parameter'], + pars['number_of_iterations'], + pars['tolerance_constant'], + pars['methodTV'], + pars['nonneg'], + pars['lipschitz_const'], 'cpu') + +Qtools = QualityTools(Im, pd_cpu) +pars['rmse'] = Qtools.rmse() + +txtstr = printParametersToString(pars) +txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) +print (txtstr) +a=fig.add_subplot(1,2,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(pd_cpu, cmap="gray") +plt.title('{}'.format('CPU results')) #%% print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") print ("_______________SB-TV (2D)__________________") diff --git a/src/Core/CMakeLists.txt b/src/Core/CMakeLists.txt index eea0d63..ac7ec91 100644 --- a/src/Core/CMakeLists.txt +++ b/src/Core/CMakeLists.txt @@ -20,7 +20,7 @@ if (OPENMP_FOUND) set (CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${OpenMP_EXE_LINKER_FLAGS} ${OpenMP_CXX_FLAGS}") set (CMAKE_SHARED_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${OpenMP_SHARED_LINKER_FLAGS} ${OpenMP_CXX_FLAGS}") set (CMAKE_STATIC_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${OpenMP_STATIC_LINKER_FLAGS} ${OpenMP_CXX_FLAGS}") - + endif() ## Build the regularisers package as a library @@ -39,21 +39,21 @@ if(WIN32) set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${FLAGS}") set (CMAKE_C_FLAGS "${CMAKE_CXX_FLAGS} ${FLAGS}") set (CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} /NODEFAULTLIB:MSVCRT.lib") - + set (EXTRA_LIBRARIES) - + message("library lib: ${LIBRARY_LIB}") - + elseif(UNIX) - set (FLAGS "-O2 -funsigned-char -Wall -Wl,--no-undefined -DCCPiReconstructionIterative_EXPORTS ") + set (FLAGS "-O2 -funsigned-char -Wall -Wl,--no-undefined -DCCPiReconstructionIterative_EXPORTS ") set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${FLAGS}") set (CMAKE_C_FLAGS "${CMAKE_CXX_FLAGS} ${FLAGS}") - - set (EXTRA_LIBRARIES + + set (EXTRA_LIBRARIES "gomp" "m" ) - + endif() message("CMAKE_CXX_FLAGS ${CMAKE_CXX_FLAGS}") @@ -66,6 +66,7 @@ message("Adding regularisers as a shared library") add_library(cilreg SHARED ${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/FGP_TV_core.c ${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/SB_TV_core.c + ${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/PD_TV_core.c ${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/TGV_core.c ${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/Diffusion_core.c ${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/Diffus4th_order_core.c @@ -80,8 +81,8 @@ add_library(cilreg SHARED ${CMAKE_CURRENT_SOURCE_DIR}/inpainters_CPU/NonlocalMarching_Inpaint_core.c ) target_link_libraries(cilreg ${EXTRA_LIBRARIES} ) -include_directories(cilreg PUBLIC - ${LIBRARY_INC}/include +include_directories(cilreg PUBLIC + ${LIBRARY_INC}/include ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/ ${CMAKE_CURRENT_SOURCE_DIR}/inpainters_CPU/ ) @@ -92,14 +93,14 @@ if (UNIX) message ("I'd install into ${CMAKE_INSTALL_PREFIX}/lib") install(TARGETS cilreg LIBRARY DESTINATION lib - CONFIGURATIONS ${CMAKE_BUILD_TYPE} + CONFIGURATIONS ${CMAKE_BUILD_TYPE} ) elseif(WIN32) message ("I'd install into ${CMAKE_INSTALL_PREFIX} lib bin") - install(TARGETS cilreg + install(TARGETS cilreg RUNTIME DESTINATION bin ARCHIVE DESTINATION lib - CONFIGURATIONS ${CMAKE_BUILD_TYPE} + CONFIGURATIONS ${CMAKE_BUILD_TYPE} ) endif() @@ -126,14 +127,14 @@ if (BUILD_CUDA) message ("I'd install into ${CMAKE_INSTALL_PREFIX}/lib") install(TARGETS cilregcuda LIBRARY DESTINATION lib - CONFIGURATIONS ${CMAKE_BUILD_TYPE} + CONFIGURATIONS ${CMAKE_BUILD_TYPE} ) elseif(WIN32) message ("I'd install into ${CMAKE_INSTALL_PREFIX} lib bin") install(TARGETS cilregcuda RUNTIME DESTINATION bin ARCHIVE DESTINATION lib - CONFIGURATIONS ${CMAKE_BUILD_TYPE} + CONFIGURATIONS ${CMAKE_BUILD_TYPE} ) endif() else() diff --git a/src/Core/regularisers_CPU/FGP_TV_core.c b/src/Core/regularisers_CPU/FGP_TV_core.c index a16a2e5..ff67af2 100644 --- a/src/Core/regularisers_CPU/FGP_TV_core.c +++ b/src/Core/regularisers_CPU/FGP_TV_core.c @@ -46,12 +46,12 @@ float TV_FGP_CPU_main(float *Input, float *Output, float *infovector, float lamb float tk = 1.0f; float tkp1 =1.0f; int count = 0; - + if (dimZ <= 1) { /*2D case */ float *Output_prev=NULL, *P1=NULL, *P2=NULL, *P1_prev=NULL, *P2_prev=NULL, *R1=NULL, *R2=NULL; DimTotal = (long)(dimX*dimY); - + if (epsil != 0.0f) Output_prev = calloc(DimTotal, sizeof(float)); P1 = calloc(DimTotal, sizeof(float)); P2 = calloc(DimTotal, sizeof(float)); @@ -59,32 +59,32 @@ float TV_FGP_CPU_main(float *Input, float *Output, float *infovector, float lamb P2_prev = calloc(DimTotal, sizeof(float)); R1 = calloc(DimTotal, sizeof(float)); R2 = calloc(DimTotal, sizeof(float)); - + /* begin iterations */ for(ll=0; ll 3) break; } - + /*storing old values*/ copyIm(P1, P1_prev, (long)(dimX), (long)(dimY), (long)(dimZ)); copyIm(P2, P2_prev, (long)(dimX), (long)(dimY), (long)(dimZ)); copyIm(P3, P3_prev, (long)(dimX), (long)(dimY), (long)(dimZ)); tk = tkp1; } - + if (epsil != 0.0f) free(Output_prev); free(P1); free(P2); free(P3); free(P1_prev); free(P2_prev); free(P3_prev); free(R1); free(R2); free(R3); } - + /*adding info into info_vector */ infovector[0] = (float)(ll); /*iterations number (if stopped earlier based on tolerance)*/ infovector[1] = re; /* reached tolerance */ - + return 0; } @@ -202,36 +202,6 @@ float Grad_func2D(float *P1, float *P2, float *D, float *R1, float *R2, float la }} return 1; } -float Proj_func2D(float *P1, float *P2, int methTV, long DimTotal) -{ - float val1, val2, denom, sq_denom; - long i; - if (methTV == 0) { - /* isotropic TV*/ -#pragma omp parallel for shared(P1,P2) private(i,denom,sq_denom) - for(i=0; i 1.0f) { - 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(i,val1,val2) - for(i=0; i 3) break; + } + /*get updated solution*/ + + getX2D(U, U_old, theta, DimTotal); + } + free(P1); free(P2); free(U_old); + } + else { + /*3D case*/ + } + /*adding info into info_vector */ + infovector[0] = (float)(ll); /*iterations number (if stopped earlier based on tolerance)*/ + infovector[1] = re; /* reached tolerance */ + + return 0; +} + +/*****************************************************************/ +/************************2D-case related Functions */ +/*****************************************************************/ + +/*Calculating dual variable (using forward differences)*/ +float DualP2D(float *U, float *P1, float *P2, long dimX, long dimY, float sigma) +{ + long i,j,index; + #pragma omp parallel for shared(U,P1,P2) private(index,i,j) + for(j=0; j +#include +#include +#include +#include +#include "omp.h" +#include "utils.h" +#include "CCPiDefines.h" + +/* C-OMP implementation of Primal-Dual TV [1] by Chambolle Pock denoising/regularization model (2D/3D case) + * + * Input Parameters: + * 1. Noisy image/volume + * 2. lambdaPar - regularization parameter + * 3. Number of iterations + * 4. eplsilon: tolerance constant + * 5. lipschitz_const: convergence related parameter + * 6. TV-type: methodTV - 'iso' (0) or 'l1' (1) + * 7. nonneg: 'nonnegativity (0 is OFF by default, 1 is ON) + + * Output: + * [1] TV - Filtered/regularized image/volume + * [2] Information vector which contains [iteration no., reached tolerance] + * + * [1] Antonin Chambolle, Thomas Pock. "A First-Order Primal-Dual Algorithm for Convex Problems with Applications to Imaging", 2010 + */ + +#ifdef __cplusplus +extern "C" { +#endif +CCPI_EXPORT float PDTV_CPU_main(float *Input, float *U, float *infovector, float lambdaPar, int iterationsNumb, float epsil, float lipschitz_const, int methodTV, int nonneg, int dimX, int dimY, int dimZ); + +CCPI_EXPORT float DualP2D(float *U, float *P1, float *P2, long dimX, long dimY, float sigma); +CCPI_EXPORT float DivProj2D(float *U, float *Input, float *P1, float *P2, long dimX, long dimY, float lt, float tau); +CCPI_EXPORT float getX2D(float *U, float *U_old, float theta, long DimTotal); +#ifdef __cplusplus +} +#endif diff --git a/src/Core/regularisers_CPU/utils.c b/src/Core/regularisers_CPU/utils.c index 5bb3a5c..cf4ad72 100644 --- a/src/Core/regularisers_CPU/utils.c +++ b/src/Core/regularisers_CPU/utils.c @@ -65,7 +65,7 @@ float TV_energy2D(float *U, float *U0, float *E_val, float lambda, int type, int { int i, j, i1, j1, index; float NOMx_2, NOMy_2, E_Grad=0.0f, E_Data=0.0f; - + /* first calculate \grad U_xy*/ for(j=0; j 1.0f) { + 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(i,val1,val2) + for(i=0; i>>>>>>'); +delete PD_TV_core* +fprintf('%s \n', '<<<<<<< CPU regularisers were successfully compiled! >>>>>>>'); pathA2 = sprintf(['..' fsep '..' fsep '..' fsep '..' fsep 'demos' fsep 'Matlab_demos'], 1i); cd(pathA2); \ No newline at end of file diff --git a/src/Matlab/mex_compile/regularisers_CPU/PD_TV.c b/src/Matlab/mex_compile/regularisers_CPU/PD_TV.c new file mode 100644 index 0000000..eac2d18 --- /dev/null +++ b/src/Matlab/mex_compile/regularisers_CPU/PD_TV.c @@ -0,0 +1,98 @@ +/* + * This work is part of the Core Imaging Library developed by + * Visual Analytics and Imaging System Group of the Science Technology + * Facilities Council, STFC + * + * Copyright 2019 Daniil Kazantsev + * Copyright 2019 Srikanth Nagella, Edoardo Pasca + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * http://www.apache.org/licenses/LICENSE-2.0 + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +#include "matrix.h" +#include "mex.h" +#include "PD_TV_core.h" + +/* C-OMP implementation of Primal-Dual TV [1] by Chambolle Pock denoising/regularization model (2D/3D case) + * + * Input Parameters: + * 1. Noisy image/volume + * 2. lambdaPar - regularization parameter + * 3. Number of iterations + * 4. eplsilon: tolerance constant + * 5. TV-type: methodTV - 'iso' (0) or 'l1' (1) + * 6. nonneg: 'nonnegativity (0 is OFF by default, 1 is ON) + * 7. lipschitz_const: convergence related parameter + + * Output: + * [1] TV - Filtered/regularized image/volume + * [2] Information vector which contains [iteration no., reached tolerance] + * + * [1] Antonin Chambolle, Thomas Pock. "A First-Order Primal-Dual Algorithm for Convex Problems with Applications to Imaging", 2010 + */ +void mexFunction( + int nlhs, mxArray *plhs[], + int nrhs, const mxArray *prhs[]) + +{ + int number_of_dims, iter, methTV, nonneg; + mwSize dimX, dimY, dimZ; + const mwSize *dim_array; + float *Input, *infovec=NULL, *Output=NULL, lambda, epsil, lipschitz_const; + + number_of_dims = mxGetNumberOfDimensions(prhs[0]); + dim_array = mxGetDimensions(prhs[0]); + + /*Handling Matlab input data*/ + if ((nrhs < 2) || (nrhs > 7)) mexErrMsgTxt("At least 2 parameters is required, all parameters are: Image(2D/3D), Regularization parameter, iterations number, tolerance, penalty type ('iso' or 'l1'), nonnegativity switch, lipschitz_const"); + + Input = (float *) mxGetData(prhs[0]); /*noisy image (2D/3D) */ + lambda = (float) mxGetScalar(prhs[1]); /* regularization parameter */ + iter = 400; /* default iterations number */ + epsil = 1.0e-06; /* default tolerance constant */ + methTV = 0; /* default isotropic TV penalty */ + nonneg = 0; /* default nonnegativity switch, off - 0 */ + lipschitz_const = 12.0; /* lipschitz_const */ + + if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); } + + if ((nrhs == 3) || (nrhs == 4) || (nrhs == 5) || (nrhs == 6) || (nrhs == 7)) iter = (int) mxGetScalar(prhs[2]); /* iterations number */ + if ((nrhs == 4) || (nrhs == 5) || (nrhs == 6) || (nrhs == 7)) epsil = (float) mxGetScalar(prhs[3]); /* tolerance constant */ + if ((nrhs == 5) || (nrhs == 6) || (nrhs == 7)) { + 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 ((nrhs == 6) || (nrhs == 7)) { + nonneg = (int) mxGetScalar(prhs[5]); + if ((nonneg != 0) && (nonneg != 1)) mexErrMsgTxt("Nonnegativity constraint can be enabled by choosing 1 or off - 0"); + } + if (nrhs == 7) lipschitz_const = (float) mxGetScalar(prhs[6]); + + + /*Handling Matlab output data*/ + dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2]; + + if (number_of_dims == 2) { + dimZ = 1; /*2D case*/ + Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); + } + if (number_of_dims == 3) { + Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); + } + mwSize vecdim[1]; + vecdim[0] = 2; + infovec = (float*)mxGetPr(plhs[1] = mxCreateNumericArray(1, vecdim, mxSINGLE_CLASS, mxREAL)); + + /* running the function */ + PDTV_CPU_main(Input, Output, infovec, lambda, iter, epsil, lipschitz_const, methTV, nonneg, dimX, dimY, dimZ); +} diff --git a/src/Python/ccpi/filters/regularisers.py b/src/Python/ccpi/filters/regularisers.py index 0b5b2ee..d65c0b9 100644 --- a/src/Python/ccpi/filters/regularisers.py +++ b/src/Python/ccpi/filters/regularisers.py @@ -2,7 +2,7 @@ script which assigns a proper device core function based on a flag ('cpu' or 'gpu') """ -from ccpi.filters.cpu_regularisers import TV_ROF_CPU, TV_FGP_CPU, TV_SB_CPU, dTV_FGP_CPU, TNV_CPU, NDF_CPU, Diff4th_CPU, TGV_CPU, LLT_ROF_CPU, PATCHSEL_CPU, NLTV_CPU +from ccpi.filters.cpu_regularisers import TV_ROF_CPU, TV_FGP_CPU, TV_PD_CPU, TV_SB_CPU, dTV_FGP_CPU, TNV_CPU, NDF_CPU, Diff4th_CPU, TGV_CPU, LLT_ROF_CPU, PATCHSEL_CPU, NLTV_CPU try: from ccpi.filters.gpu_regularisers import TV_ROF_GPU, TV_FGP_GPU, TV_SB_GPU, dTV_FGP_GPU, NDF_GPU, Diff4th_GPU, TGV_GPU, LLT_ROF_GPU, PATCHSEL_GPU gpu_enabled = True @@ -51,6 +51,31 @@ def FGP_TV(inputData, regularisation_parameter,iterations, raise ValueError ('GPU is not available') raise ValueError('Unknown device {0}. Expecting gpu or cpu'\ .format(device)) + +def PD_TV(inputData, regularisation_parameter, iterations, + tolerance_param, methodTV, nonneg, lipschitz_const, device='cpu'): + if device == 'cpu': + return TV_PD_CPU(inputData, + regularisation_parameter, + iterations, + tolerance_param, + methodTV, + nonneg, + lipschitz_const) + elif device == 'gpu' and gpu_enabled: + return TV_PD_CPU(inputData, + regularisation_parameter, + iterations, + tolerance_param, + methodTV, + nonneg, + lipschitz_const) + else: + if not gpu_enabled and device == 'gpu': + raise ValueError ('GPU is not available') + raise ValueError('Unknown device {0}. Expecting gpu or cpu'\ + .format(device)) + def SB_TV(inputData, regularisation_parameter, iterations, tolerance_param, methodTV, device='cpu'): if device == 'cpu': @@ -212,4 +237,3 @@ def NDF_INP(inputData, maskData, regularisation_parameter, edge_parameter, itera def NVM_INP(inputData, maskData, SW_increment, iterations): return NVM_INPAINT_CPU(inputData, maskData, SW_increment, iterations) - diff --git a/src/Python/setup-regularisers.py.in b/src/Python/setup-regularisers.py.in index 4c578e3..9bcd46d 100644 --- a/src/Python/setup-regularisers.py.in +++ b/src/Python/setup-regularisers.py.in @@ -8,13 +8,13 @@ from Cython.Distutils import build_ext import os import sys import numpy -import platform +import platform cil_version=os.environ['CIL_VERSION'] if cil_version == '': print("Please set the environmental variable CIL_VERSION") sys.exit(1) - + library_include_path = "" library_lib_path = "" try: @@ -23,7 +23,7 @@ try: except: library_include_path = os.environ['PREFIX']+'/include' pass - + extra_include_dirs = [numpy.get_include(), library_include_path] #extra_library_dirs = [os.path.join(library_include_path, "..", "lib")] extra_compile_args = [] @@ -38,6 +38,7 @@ extra_include_dirs += [os.path.join(".." , "Core"), os.path.join(".." , "Core", "regularisers_CPU"), os.path.join(".." , "Core", "inpainters_CPU"), os.path.join(".." , "Core", "regularisers_GPU" , "TV_FGP" ) , + os.path.join(".." , "Core", "regularisers_GPU" , "TV_PD" ) , os.path.join(".." , "Core", "regularisers_GPU" , "TV_ROF" ) , os.path.join(".." , "Core", "regularisers_GPU" , "TV_SB" ) , os.path.join(".." , "Core", "regularisers_GPU" , "TGV" ) , @@ -48,12 +49,12 @@ extra_include_dirs += [os.path.join(".." , "Core"), os.path.join(".." , "Core", "regularisers_GPU" , "PatchSelect" ) , "."] -if platform.system() == 'Windows': - extra_compile_args[0:] = ['/DWIN32','/EHsc','/DBOOST_ALL_NO_LIB' , '/openmp' ] +if platform.system() == 'Windows': + extra_compile_args[0:] = ['/DWIN32','/EHsc','/DBOOST_ALL_NO_LIB' , '/openmp' ] else: extra_compile_args = ['-fopenmp','-O2', '-funsigned-char', '-Wall', '-std=c++0x'] extra_libraries += [@EXTRA_OMP_LIB@] - + setup( name='ccpi', description='CCPi Core Imaging Library - Image regularisers', @@ -61,13 +62,13 @@ setup( cmdclass = {'build_ext': build_ext}, ext_modules = [Extension("ccpi.filters.cpu_regularisers", sources=[os.path.join("." , "src", "cpu_regularisers.pyx" ) ], - include_dirs=extra_include_dirs, - library_dirs=extra_library_dirs, - extra_compile_args=extra_compile_args, - libraries=extra_libraries ), - + include_dirs=extra_include_dirs, + library_dirs=extra_library_dirs, + extra_compile_args=extra_compile_args, + libraries=extra_libraries ), + ], - zip_safe = False, + zip_safe = False, packages = {'ccpi', 'ccpi.filters', 'ccpi.supp'}, ) diff --git a/src/Python/src/cpu_regularisers.pyx b/src/Python/src/cpu_regularisers.pyx index 4917d06..724634b 100644 --- a/src/Python/src/cpu_regularisers.pyx +++ b/src/Python/src/cpu_regularisers.pyx @@ -20,6 +20,7 @@ cimport numpy as np cdef extern float TV_ROF_CPU_main(float *Input, float *Output, float *infovector, float *lambdaPar, int lambda_is_arr, int iterationsNumb, float tau, float epsil, int dimX, int dimY, int dimZ); cdef extern float TV_FGP_CPU_main(float *Input, float *Output, float *infovector, float lambdaPar, int iterationsNumb, float epsil, int methodTV, int nonneg, int dimX, int dimY, int dimZ); +cdef extern float PDTV_CPU_main(float *Input, float *U, float *infovector, float lambdaPar, int iterationsNumb, float epsil, float lipschitz_const, int methodTV, int nonneg, int dimX, int dimY, int dimZ); cdef extern float SB_TV_CPU_main(float *Input, float *Output, float *infovector, float mu, int iter, float epsil, int methodTV, int dimX, int dimY, int dimZ); cdef extern float LLT_ROF_CPU_main(float *Input, float *Output, float *infovector, float lambdaROF, float lambdaLLT, int iterationsNumb, float tau, float epsil, int dimX, int dimY, int dimZ); cdef extern float TGV_main(float *Input, float *Output, float *infovector, float lambdaPar, float alpha1, float alpha0, int iterationsNumb, float L2, float epsil, int dimX, int dimY, int dimZ); @@ -58,9 +59,6 @@ def TV_ROF_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, cdef np.ndarray[np.float32_t, ndim=1, mode="c"] infovec = \ np.ones([2], dtype='float32') - # Run ROF iterations for 2D data - # TV_ROF_CPU_main(&inputData[0,0], &outputData[0,0], &infovec[0], regularisation_parameter, iterationsNumb, marching_step_parameter, tolerance_param, dims[1], dims[0], 1) - # Run ROF iterations for 2D data if isinstance (regularisation_parameter, np.ndarray): reg = regularisation_parameter.copy() TV_ROF_CPU_main(&inputData[0,0], &outputData[0,0], &infovec[0], ®[0,0], 1, iterationsNumb, marching_step_parameter, tolerance_param, dims[1], dims[0], 1) @@ -158,6 +156,44 @@ def TV_FGP_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, dims[2], dims[1], dims[0]) return (outputData,infovec) +#****************************************************************# +#****************** Total-variation Primal-dual *****************# +#****************************************************************# +def TV_PD_CPU(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, lipschitz_const): + if inputData.ndim == 2: + return TV_PD_2D(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, lipschitz_const) + elif inputData.ndim == 3: + return 0 + +def TV_PD_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, + float regularisation_parameter, + int iterationsNumb, + float tolerance_param, + int methodTV, + int nonneg, + float lipschitz_const): + + 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') + + cdef np.ndarray[np.float32_t, ndim=1, mode="c"] infovec = \ + np.ones([2], dtype='float32') + + #/* Run FGP-TV iterations for 2D data */ + PDTV_CPU_main(&inputData[0,0], &outputData[0,0], &infovec[0], regularisation_parameter, + iterationsNumb, + tolerance_param, + lipschitz_const, + methodTV, + nonneg, + dims[1],dims[0], 1) + + return (outputData,infovec) + #***************************************************************# #********************** Total-variation SB *********************# #***************************************************************# -- cgit v1.2.3