From 39baef90c4b209090f006e5308653cb0a3348c4e Mon Sep 17 00:00:00 2001 From: dkazanc Date: Wed, 6 Mar 2019 15:13:58 +0000 Subject: tol work --- Readme.md | 2 +- build/run.sh | 24 ++++--- demos/SoftwareX_supp/Demo_VolumeDenoise.py | 62 ++++++++++++++++- demos/demoMatlab_denoise.m | 12 ++-- demos/demo_cpu_regularisers.py | 16 ++--- run.sh | 26 +++++++ src/Core/regularisers_CPU/FGP_TV_core.c | 88 ++++++++++++++---------- src/Core/regularisers_CPU/FGP_TV_core.h | 7 +- src/Matlab/mex_compile/compileCPU_mex_Linux.m | 13 ++-- src/Matlab/mex_compile/regularisers_CPU/FGP_TV.c | 28 ++++---- src/Python/ccpi/filters/regularisers.py | 7 +- src/Python/setup-regularisers.py.in | 2 +- src/Python/src/cpu_regularisers.pyx | 38 +++++----- 13 files changed, 217 insertions(+), 108 deletions(-) create mode 100755 run.sh diff --git a/Readme.md b/Readme.md index 92b4273..cc96b2c 100644 --- a/Readme.md +++ b/Readme.md @@ -109,7 +109,7 @@ conda install ccpi-regulariser -c ccpi -c conda-forge #### Python (conda-build) ``` - export CIL_VERSION=19.02 + export CIL_VERSION=19.03 conda build recipe/ --numpy 1.12 --python 3.5 conda install ccpi-regulariser=${CIL_VERSION} --use-local --force cd demos/ diff --git a/build/run.sh b/build/run.sh index 332d660..d450299 100644 --- a/build/run.sh +++ b/build/run.sh @@ -1,22 +1,24 @@ #!/bin/bash echo "Building CCPi-regularisation Toolkit using CMake" -rm -r build +rm -r build_proj # Requires Cython, install it first: # pip install cython -mkdir build -cd build/ +mkdir build_proj +cd build_proj/ make clean +export CIL_VERSION=19.03 # 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 # install Matlab modules with CUDA -cmake ../ -DBUILD_PYTHON_WRAPPER=OFF -DMatlab_ROOT_DIR=/dls_sw/apps/matlab/r2014a/ -DBUILD_MATLAB_WRAPPER=ON -DBUILD_CUDA=ON -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=./install +# cmake ../ -DBUILD_PYTHON_WRAPPER=OFF -DMatlab_ROOT_DIR=/dls_sw/apps/matlab/r2014a/ -DBUILD_MATLAB_WRAPPER=ON -DBUILD_CUDA=ON -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=./install make install -# cp install/lib/libcilreg.so install/python/ccpi/filters -#cd install/python -#export LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:../lib -# spyder -# one can also run Matlab in Linux as: +#### Python +#cp install/lib/libcilreg.so install/python/ccpi/filters +cd install/python +export LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:../lib +spyder +##### one can also run Matlab in Linux as: #PATH="/path/to/mex/:$PATH" LD_LIBRARY_PATH="/path/to/library:$LD_LIBRARY_PATH" matlab -PATH="/home/kjy41806/Documents/SOFT/CCPi-Regularisation-Toolkit/build/install/matlab/:$PATH" LD_LIBRARY_PATH="/home/kjy41806/Documents/SOFT/CCPi-Regularisation-Toolkit/build/install/lib:$LD_LIBRARY_PATH" matlab +#PATH="/home/kjy41806/Documents/SOFT/CCPi-Regularisation-Toolkit/build/install/matlab/:$PATH" LD_LIBRARY_PATH="/home/kjy41806/Documents/SOFT/CCPi-Regularisation-Toolkit/build/install/lib:$LD_LIBRARY_PATH" matlab diff --git a/demos/SoftwareX_supp/Demo_VolumeDenoise.py b/demos/SoftwareX_supp/Demo_VolumeDenoise.py index 2387e94..17cdf4d 100644 --- a/demos/SoftwareX_supp/Demo_VolumeDenoise.py +++ b/demos/SoftwareX_supp/Demo_VolumeDenoise.py @@ -60,13 +60,14 @@ plt.title('3D Phantom, sagittal view') plt.show() #%% print ("____________________Applying regularisers_______________________") +print ("________________________________________________________________") print ("#############ROF TV CPU####################") # set parameters pars = {'algorithm': ROF_TV, \ 'input' : phantom_noise,\ 'regularisation_parameter':0.04,\ - 'number_of_iterations': 100,\ + 'number_of_iterations': 600,\ 'time_marching_parameter': 0.0025 } @@ -116,6 +117,65 @@ win2d = win * (win.T) ssim_rof = Qtools.ssim(win2d) print("ROF-TV (gpu) ____ RMSE: {}, MMSIM: {}, run time: {} sec".format(RMSE_rof,ssim_rof[0],Run_time_rof)) +#%% +print ("#############FGP TV CPU####################") +# set parameters +pars = {'algorithm':FGP_TV, \ + 'input' : phantom_noise,\ + 'regularisation_parameter':0.04,\ + 'number_of_iterations': 300,\ + 'time_marching_parameter': 0.0025,\ + 'tolerance_constant':1e-05,\ + } + +tic=timeit.default_timer() +fgp_cpu3D = FGP_TV(pars['input'], + pars['regularisation_parameter'], + pars['number_of_iterations'], + pars['time_marching_parameter'],'cpu') +toc=timeit.default_timer() + +Run_time_fgp = toc - tic +Qtools = QualityTools(phantom_tm, fgp_cpu3D) +RMSE_rof = Qtools.rmse() + +# SSIM measure +Qtools = QualityTools(phantom_tm[128,:,:]*255, fgp_cpu3D[128,:,:]*235) +win = np.array([gaussian(11, 1.5)]) +win2d = win * (win.T) +ssim_fgp = Qtools.ssim(win2d) + +print("FGP-TV (cpu) ____ RMSE: {}, MMSIM: {}, run time: {} sec".format(RMSE_rof,ssim_fgp[0],Run_time_fgp)) +#%% +print ("#############FGP TV GPU####################") +# set parameters +pars = {'algorithm':FGP_TV, \ + 'input' : phantom_noise,\ + 'regularisation_parameter':0.04,\ + 'number_of_iterations': 300,\ + 'time_marching_parameter': 0.0025,\ + 'tolerance_constant':1e-05,\ + } +tic=timeit.default_timer() +fgp_gpu3D = FGP_TV(pars['input'], + pars['regularisation_parameter'], + pars['number_of_iterations'], + pars['time_marching_parameter'],'gpu') +toc=timeit.default_timer() + +Run_time_fgp = toc - tic +Qtools = QualityTools(phantom_tm, fgp_gpu3D) +RMSE_rof = Qtools.rmse() + +# SSIM measure +Qtools = QualityTools(phantom_tm[128,:,:]*255, fgp_gpu3D[128,:,:]*235) +win = np.array([gaussian(11, 1.5)]) +win2d = win * (win.T) +ssim_fgp = Qtools.ssim(win2d) + +print("FGP-TV (gpu) ____ RMSE: {}, MMSIM: {}, run time: {} sec".format(RMSE_rof,ssim_fgp[0],Run_time_fgp)) #%% + + diff --git a/demos/demoMatlab_denoise.m b/demos/demoMatlab_denoise.m index 5e92ee1..fa81f6d 100644 --- a/demos/demoMatlab_denoise.m +++ b/demos/demoMatlab_denoise.m @@ -2,9 +2,11 @@ clear; close all fsep = '/'; -Path1 = sprintf(['..' fsep 'src' fsep 'Matlab' fsep 'mex_compile' fsep 'installed'], 1i); -Path2 = sprintf([ data' fsep], 1i); -Path3 = sprintf(['..' filesep 'src' filesep 'Matlab' filesep 'supp'], 1i); +%Path1 = sprintf(['..' fsep 'src' fsep 'Matlab' fsep 'mex_compile' fsep 'installed'], 1i); +Path1 = ('/home/kjy41806/Documents/SOFT/CCPi-Regularisation-Toolkit/src/Matlab/mex_compile/installed'); +Path2 = sprintf(['data' fsep], 1i); +%Path3 = sprintf(['..' filesep 'src' filesep 'Matlab' filesep 'supp'], 1i); +Path3 = '/home/kjy41806/Documents/SOFT/CCPi-Regularisation-Toolkit/src/Matlab/supp'; addpath(Path1); addpath(Path2); addpath(Path3); @@ -34,8 +36,8 @@ figure; imshow(u_rof, [0 1]); title('ROF-TV denoised image (CPU)'); %% fprintf('Denoise using the FGP-TV model (CPU) \n'); lambda_reg = 0.033; -iter_fgp = 300; % number of FGP iterations -epsil_tol = 1.0e-09; % tolerance +iter_fgp = 200; % number of FGP iterations +epsil_tol = 1.0e-05; % tolerance tic; u_fgp = FGP_TV(single(u0), lambda_reg, iter_fgp, epsil_tol); toc; energyfunc_val_fgp = TV_energy(single(u_fgp),single(u0),lambda_reg, 1); % get energy function value rmseFGP = (RMSE(u_fgp(:),Im(:))); diff --git a/demos/demo_cpu_regularisers.py b/demos/demo_cpu_regularisers.py index d34607a..e58ed2f 100644 --- a/demos/demo_cpu_regularisers.py +++ b/demos/demo_cpu_regularisers.py @@ -30,8 +30,9 @@ def printParametersToString(pars): txt += '\n' return txt ############################################################################### -#%% -filename = os.path.join( "data" ,"lena_gray_512.tif") + +# filename = os.path.join( "data" ,"lena_gray_512.tif") +filename = "/home/kjy41806/Documents/SOFT/CCPi-Regularisation-Toolkit/test/lena_gray_512.tif" # read image Im = plt.imread(filename) @@ -127,23 +128,20 @@ imgplot = plt.imshow(u0,cmap="gray") pars = {'algorithm' : FGP_TV, \ 'input' : u0,\ 'regularisation_parameter':0.04, \ - 'number_of_iterations' :2000 ,\ + 'number_of_iterations' :200 ,\ 'tolerance_constant':1e-06,\ 'methodTV': 0 ,\ - 'nonneg': 0 ,\ - 'printingOut': 0 + 'nonneg': 0 } print ("#############FGP TV CPU####################") start_time = timeit.default_timer() -fgp_cpu = FGP_TV(pars['input'], +fgp_cpu,info_vec = FGP_TV(pars['input'], pars['regularisation_parameter'], pars['number_of_iterations'], pars['tolerance_constant'], pars['methodTV'], - pars['nonneg'], - pars['printingOut'],'cpu') - + pars['nonneg'],'cpu') Qtools = QualityTools(Im, fgp_cpu) pars['rmse'] = Qtools.rmse() diff --git a/run.sh b/run.sh new file mode 100755 index 0000000..c5c396a --- /dev/null +++ b/run.sh @@ -0,0 +1,26 @@ +#!/bin/bash +echo "Building CCPi-regularisation Toolkit using CMake" +rm -r build_proj +# Requires Cython, install it first: +# pip install cython +mkdir build_proj +cd build_proj/ +#make clean +export CIL_VERSION=19.03 +# 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 +# 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 +# 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 +# cmake ../ -DBUILD_PYTHON_WRAPPER=OFF -DMatlab_ROOT_DIR=/dls_sw/apps/matlab/r2014a/ -DBUILD_MATLAB_WRAPPER=ON -DBUILD_CUDA=ON -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=./install +make install +#### Python +#cp install/lib/libcilreg.so install/python/ccpi/filters +cd install/python +export LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:../lib +spyder +##### Matlab (Linux) +#PATH="/path/to/mex/:$PATH" LD_LIBRARY_PATH="/path/to/library:$LD_LIBRARY_PATH" matlab +#PATH="/home/kjy41806/Documents/SOFT/CCPi-Regularisation-Toolkit/build_proj/install/matlab/:$PATH" LD_LIBRARY_PATH="/home/kjy41806/Documents/SOFT/CCPi-Regularisation-Toolkit/build_proj/install/lib:$LD_LIBRARY_PATH" matlab diff --git a/src/Core/regularisers_CPU/FGP_TV_core.c b/src/Core/regularisers_CPU/FGP_TV_core.c index 68d58b7..f613f0e 100644 --- a/src/Core/regularisers_CPU/FGP_TV_core.c +++ b/src/Core/regularisers_CPU/FGP_TV_core.c @@ -28,28 +28,33 @@ limitations under the License. * 4. eplsilon: tolerance constant * 5. TV-type: methodTV - 'iso' (0) or 'l1' (1) * 6. nonneg: 'nonnegativity (0 is OFF by default) - * 7. print information: 0 (off) or 1 (on) * * Output: - * [1] Filtered/regularized image + * [1] Filtered/regularized image/volume + * [2] Information vector which contains [iteration no., reached tolerance] * * This function is based on the Matlab's code and paper by * [1] Amir Beck and Marc Teboulle, "Fast Gradient-Based Algorithms for Constrained Total Variation Image Denoising and Deblurring Problems" */ -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) +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) { - int ll; + int ll; long j, DimTotal; - float re, re1; - float tk = 1.0f; - float tkp1=1.0f; + float re, re1; + re = 0.0f; re1 = 0.0f; + float tk = 1.0f; + float tkp1 =1.0f; int count = 0; + + /*adding info into info_vector */ +// infovector[0] = 50.1f; /*iterations number (if stopped earlier based on tolerance)*/ +// infovector[1] = 0.55f; /* reached tolerance */ 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); + /*2D case */ + float *Output_prev=NULL, *P1=NULL, *P2=NULL, *P1_prev=NULL, *P2_prev=NULL, *R1=NULL, *R2=NULL; + DimTotal = (long)(dimX*dimY); Output_prev = calloc(DimTotal, sizeof(float)); P1 = calloc(DimTotal, sizeof(float)); @@ -59,7 +64,7 @@ float TV_FGP_CPU_main(float *Input, float *Output, float lambdaPar, int iteratio R1 = calloc(DimTotal, sizeof(float)); R2 = calloc(DimTotal, sizeof(float)); - /* begin iterations */ + /* begin iterations */ for(ll=0; ll 4) break; + + copyIm(Output, Output_prev, (long)(dimX), (long)(dimY), 1l); } - re = sqrt(re)/sqrt(re1); - if (re < epsil) count++; - if (count > 4) break; - /*storing old values*/ - copyIm(Output, Output_prev, (long)(dimX), (long)(dimY), 1l); copyIm(P1, P1_prev, (long)(dimX), (long)(dimY), 1l); copyIm(P2, P2_prev, (long)(dimX), (long)(dimY), 1l); tk = tkp1; } - if (printM == 1) printf("FGP-TV iterations stopped at iteration %i \n", ll); - free(Output_prev); free(P1); free(P2); free(P1_prev); free(P2_prev); free(R1); free(R2); + /*adding info into info_vector */ + //info_vector[0] = (float)(ll); /*iterations number (if stopped earlier based on tolerance)*/ +// info_vector[1] = 0.55f; /* reached tolerance */ + + copyIm(Input, infovector, (long)(dimX), (long)(dimY), 1l); + +// printf("%f\n", infovector[128]); + + free(Output_prev); free(P1); free(P2); free(P1_prev); free(P2_prev); free(R1); free(R2); } else { /*3D case*/ @@ -134,26 +147,31 @@ float TV_FGP_CPU_main(float *Input, float *Output, float lambdaPar, int iteratio Rupd_func3D(P1, P1_prev, P2, P2_prev, P3, P3_prev, R1, R2, R3, tkp1, tk, DimTotal); /* calculate norm - stopping rules*/ - re = 0.0f; re1 = 0.0f; - for(j=0; j 4) break; - - /*storing old values*/ - copyIm(Output, Output_prev, (long)(dimX), (long)(dimY), (long)(dimZ)); + if (count > 4) break; + + copyIm(Output, Output_prev, (long)(dimX), (long)(dimY), (long)(dimZ)); + } + + /*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 (printM == 1) printf("FGP-TV iterations stopped at iteration %i \n", ll); - 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 */ + + free(Output_prev); free(P1); free(P2); free(P3); free(P1_prev); free(P2_prev); free(P3_prev); free(R1); free(R2); free(R3); } return *Output; } diff --git a/src/Core/regularisers_CPU/FGP_TV_core.h b/src/Core/regularisers_CPU/FGP_TV_core.h index 3418604..04e6e80 100644 --- a/src/Core/regularisers_CPU/FGP_TV_core.h +++ b/src/Core/regularisers_CPU/FGP_TV_core.h @@ -35,10 +35,11 @@ limitations under the License. * 4. eplsilon: tolerance constant * 5. TV-type: methodTV - 'iso' (0) or 'l1' (1) * 6. nonneg: 'nonnegativity (0 is OFF by default) - * 7. print information: 0 (off) or 1 (on) * * Output: - * [1] Filtered/regularized image + * [1] Filtered/regularized image/volume + * [2] Information vector which contains [iteration no., reached tolerance] + * * This function is based on the Matlab's code and paper by * [1] Amir Beck and Marc Teboulle, "Fast Gradient-Based Algorithms for Constrained Total Variation Image Denoising and Deblurring Problems" @@ -47,7 +48,7 @@ limitations under the License. #ifdef __cplusplus extern "C" { #endif -CCPI_EXPORT 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); +CCPI_EXPORT 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); CCPI_EXPORT float Obj_func2D(float *A, float *D, float *R1, float *R2, float lambda, long dimX, long dimY); CCPI_EXPORT float Grad_func2D(float *P1, float *P2, float *D, float *R1, float *R2, float lambda, long dimX, long dimY); diff --git a/src/Matlab/mex_compile/compileCPU_mex_Linux.m b/src/Matlab/mex_compile/compileCPU_mex_Linux.m index 72a828e..f3d9ce1 100644 --- a/src/Matlab/mex_compile/compileCPU_mex_Linux.m +++ b/src/Matlab/mex_compile/compileCPU_mex_Linux.m @@ -2,9 +2,9 @@ fsep = '/'; -pathcopyFrom = sprintf(['..' fsep '..' fsep '..' fsep 'Core' fsep 'regularisers_CPU'], 1i); -pathcopyFrom1 = sprintf(['..' fsep '..' fsep '..' fsep 'Core' fsep 'CCPiDefines.h'], 1i); -pathcopyFrom2 = sprintf(['..' fsep '..' fsep '..' fsep 'Core' fsep 'inpainters_CPU'], 1i); +pathcopyFrom = sprintf(['..' fsep '..' fsep 'Core' fsep 'regularisers_CPU'], 1i); +pathcopyFrom1 = sprintf(['..' fsep '..' fsep 'Core' fsep 'CCPiDefines.h'], 1i); +pathcopyFrom2 = sprintf(['..' fsep '..' fsep 'Core' fsep 'inpainters_CPU'], 1i); copyfile(pathcopyFrom, 'regularisers_CPU'); copyfile(pathcopyFrom1, 'regularisers_CPU'); @@ -76,6 +76,7 @@ delete PatchSelect_core* Nonlocal_TV_core* delete Diffusion_Inpaint_core* NonlocalMarching_Inpaint_core* fprintf('%s \n', '<<<<<<< Regularisers successfully compiled! >>>>>>>'); -pathA2 = sprintf(['..' fsep '..' fsep], 1i); -cd(pathA2); -cd demos +%pathA2 = sprintf(['..' fsep '..' fsep], 1i); +% cd(pathA2); +cd('/home/kjy41806/Documents/SOFT/CCPi-Regularisation-Toolkit/demos') +%cd demos diff --git a/src/Matlab/mex_compile/regularisers_CPU/FGP_TV.c b/src/Matlab/mex_compile/regularisers_CPU/FGP_TV.c index 642362f..603e0f4 100644 --- a/src/Matlab/mex_compile/regularisers_CPU/FGP_TV.c +++ b/src/Matlab/mex_compile/regularisers_CPU/FGP_TV.c @@ -47,13 +47,13 @@ void mexFunction( int number_of_dims, iter, methTV, printswitch, nonneg; mwSize dimX, dimY, dimZ; const mwSize *dim_array; - float *Input, *Output=NULL, lambda, epsil; + float *Input, *Output=NULL, *Output2=NULL, lambda, epsil; 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, Regularization parameter, iterations number, tolerance, penalty type ('iso' or 'l1'), nonnegativity switch, print switch"); + if ((nrhs < 2) || (nrhs > 6)) mexErrMsgTxt("At least 2 parameters is required, all parameters are: Image(2D/3D), Regularization parameter, Regularization parameter, iterations number, tolerance, penalty type ('iso' or 'l1'), nonnegativity switch"); Input = (float *) mxGetData(prhs[0]); /*noisy image (2D/3D) */ lambda = (float) mxGetScalar(prhs[1]); /* regularization parameter */ @@ -61,27 +61,24 @@ void mexFunction( epsil = 0.0001; /* default tolerance constant */ methTV = 0; /* default isotropic TV penalty */ nonneg = 0; /* default nonnegativity switch, off - 0 */ - printswitch = 0; /*default print is switched, off - 0 */ +// printswitch = 0; /*default print is switched, off - 0 */ 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)) { + if ((nrhs == 3) || (nrhs == 4) || (nrhs == 5) || (nrhs == 6)) iter = (int) mxGetScalar(prhs[2]); /* iterations number */ + if ((nrhs == 4) || (nrhs == 5) || (nrhs == 6)) epsil = (float) mxGetScalar(prhs[3]); /* tolerance constant */ + if ((nrhs == 5) || (nrhs == 6)) { 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)) { + if (nrhs == 6) { 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) { - printswitch = (int) mxGetScalar(prhs[6]); - if ((printswitch != 0) && (printswitch != 1)) mexErrMsgTxt("Print can be enabled by choosing 1 or off - 0"); - } + /*Handling Matlab output data*/ dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2]; @@ -89,9 +86,12 @@ void mexFunction( if (number_of_dims == 2) { dimZ = 1; /*2D case*/ Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); + Output2 = (float*)mxGetPr(plhs[1] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); + } + if (number_of_dims == 3) { + Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); } - if (number_of_dims == 3) Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); /* running the function */ - TV_FGP_CPU_main(Input, Output, lambda, iter, epsil, methTV, nonneg, printswitch, dimX, dimY, dimZ); -} \ No newline at end of file + TV_FGP_CPU_main(Input, Output, Output2, lambda, iter, epsil, methTV, nonneg, dimX, dimY, dimZ); +} diff --git a/src/Python/ccpi/filters/regularisers.py b/src/Python/ccpi/filters/regularisers.py index 588ea32..fb2c999 100644 --- a/src/Python/ccpi/filters/regularisers.py +++ b/src/Python/ccpi/filters/regularisers.py @@ -29,15 +29,14 @@ def ROF_TV(inputData, regularisation_parameter, iterations, .format(device)) def FGP_TV(inputData, regularisation_parameter,iterations, - tolerance_param, methodTV, nonneg, printM, device='cpu'): + tolerance_param, methodTV, nonneg, device='cpu'): if device == 'cpu': return TV_FGP_CPU(inputData, regularisation_parameter, iterations, tolerance_param, methodTV, - nonneg, - printM) + nonneg) elif device == 'gpu' and gpu_enabled: return TV_FGP_GPU(inputData, regularisation_parameter, @@ -45,7 +44,7 @@ def FGP_TV(inputData, regularisation_parameter,iterations, tolerance_param, methodTV, nonneg, - printM) + 1) else: if not gpu_enabled and device == 'gpu': raise ValueError ('GPU is not available') diff --git a/src/Python/setup-regularisers.py.in b/src/Python/setup-regularisers.py.in index 82d9f9f..39b820a 100644 --- a/src/Python/setup-regularisers.py.in +++ b/src/Python/setup-regularisers.py.in @@ -44,7 +44,7 @@ extra_include_dirs += [os.path.join(".." , "Core"), os.path.join(".." , "Core", "regularisers_GPU" , "LLTROF" ) , os.path.join(".." , "Core", "regularisers_GPU" , "NDF" ) , os.path.join(".." , "Core", "regularisers_GPU" , "dTV_FGP" ) , - os.path.join(".." , "Core", "regularisers_GPU" , "DIFF4th" ) , + os.path.join(".." , "Core", "regularisers_GPU" , "Diff4th" ) , os.path.join(".." , "Core", "regularisers_GPU" , "PatchSelect" ) , "."] diff --git a/src/Python/src/cpu_regularisers.pyx b/src/Python/src/cpu_regularisers.pyx index 11a0617..b7d029d 100644 --- a/src/Python/src/cpu_regularisers.pyx +++ b/src/Python/src/cpu_regularisers.pyx @@ -19,7 +19,7 @@ import numpy as np 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); +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 SB_TV_CPU_main(float *Input, float *Output, float lambdaPar, int iterationsNumb, float epsil, int methodTV, int printM, int dimX, int dimY, int dimZ); cdef extern float LLT_ROF_CPU_main(float *Input, float *Output, float lambdaROF, float lambdaLLT, int iterationsNumb, float tau, int dimX, int dimY, int dimZ); cdef extern float TGV_main(float *Input, float *Output, float lambdaPar, float alpha1, float alpha0, int iterationsNumb, float L2, int dimX, int dimY, int dimZ); @@ -45,7 +45,7 @@ def TV_ROF_CPU(inputData, regularisation_parameter, iterationsNumb, marching_ste def TV_ROF_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, float regularisation_parameter, - int iterationsNumb, + int iterationsNumb, float marching_step_parameter): cdef long dims[2] dims[0] = inputData.shape[0] @@ -80,45 +80,46 @@ def TV_ROF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, #********************** Total-variation FGP *********************# #****************************************************************# #******** Total-variation Fast-Gradient-Projection (FGP)*********# -def TV_FGP_CPU(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, printM): +def TV_FGP_CPU(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg): if inputData.ndim == 2: - return TV_FGP_2D(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, printM) + return TV_FGP_2D(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg) elif inputData.ndim == 3: - return TV_FGP_3D(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, printM) + return TV_FGP_3D(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg) def TV_FGP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, float regularisation_parameter, int iterationsNumb, float tolerance_param, int methodTV, - int nonneg, - int printM): - + int nonneg): + 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=2, mode="c"] infovec = \ + np.ones([dims[0],dims[1]], dtype='float32') + #/* Run FGP-TV iterations for 2D data */ - TV_FGP_CPU_main(&inputData[0,0], &outputData[0,0], regularisation_parameter, + TV_FGP_CPU_main(&inputData[0,0], &outputData[0,0], &infovec[0,0], regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, - printM, dims[1],dims[0],1) - return outputData + return (outputData,infovec) def TV_FGP_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, float regularisation_parameter, int iterationsNumb, float tolerance_param, int methodTV, - int nonneg, - int printM): + int nonneg): + cdef long dims[3] dims[0] = inputData.shape[0] dims[1] = inputData.shape[1] @@ -126,16 +127,17 @@ def TV_FGP_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \ np.zeros([dims[0], dims[1], dims[2]], dtype='float32') - + cdef np.ndarray[np.float32_t, ndim=1, mode="c"] infovec = \ + np.zeros([2], dtype='float32') + #/* Run FGP-TV iterations for 3D data */ - TV_FGP_CPU_main(&inputData[0,0,0], &outputData[0,0,0], regularisation_parameter, + TV_FGP_CPU_main(&inputData[0,0,0], &outputData[0,0,0], &infovec[0], regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, - printM, dims[2], dims[1], dims[0]) - return outputData + return (outputData,infovec) #***************************************************************# #********************** Total-variation SB *********************# -- cgit v1.2.3