diff options
Diffstat (limited to 'Wrappers/Matlab')
7 files changed, 211 insertions, 3 deletions
diff --git a/Wrappers/Matlab/demos/demoMatlab_3Ddenoise.m b/Wrappers/Matlab/demos/demoMatlab_3Ddenoise.m index fb55097..502b6bd 100644 --- a/Wrappers/Matlab/demos/demoMatlab_3Ddenoise.m +++ b/Wrappers/Matlab/demos/demoMatlab_3Ddenoise.m @@ -53,6 +53,15 @@ figure; imshow(u_sb(:,:,15), [0 1]); title('SB-TV denoised volume (CPU)'); % tic; u_sbG = SB_TV_GPU(single(vol3D), lambda_reg, iter_sb, epsil_tol); toc; % figure; imshow(u_sbG(:,:,15), [0 1]); title('SB-TV denoised volume (GPU)'); %% +%% +fprintf('Denoise a volume using Nonlinear-Diffusion model (CPU) \n'); +iter_diff = 300; % number of diffusion iterations +lambda_regDiff = 0.06; % regularisation for the diffusivity +sigmaPar = 0.04; % edge-preserving parameter +tau_param = 0.025; % time-marching constant +tic; u_diff = NonlDiff(single(vol3D), lambda_regDiff, sigmaPar, iter_diff, tau_param, 'Huber'); toc; +figure; imshow(u_diff(:,:,15), [0 1]); title('Diffusion denoised volume (CPU)'); +%% %>>>>>>>>>>>>>> MULTI-CHANNEL priors <<<<<<<<<<<<<<< % fprintf('Denoise a volume using the FGP-dTV model (CPU) \n'); diff --git a/Wrappers/Matlab/demos/demoMatlab_denoise.m b/Wrappers/Matlab/demos/demoMatlab_denoise.m index dab98dc..4a0a19a 100644 --- a/Wrappers/Matlab/demos/demoMatlab_denoise.m +++ b/Wrappers/Matlab/demos/demoMatlab_denoise.m @@ -46,6 +46,22 @@ figure; imshow(u_sb, [0 1]); title('SB-TV denoised image (CPU)'); % tic; u_sbG = SB_TV_GPU(single(u0), lambda_reg, iter_sb, epsil_tol); toc; % figure; imshow(u_sbG, [0 1]); title('SB-TV denoised image (GPU)'); %% +fprintf('Denoise using Nonlinear-Diffusion model (CPU) \n'); +iter_diff = 800; % number of diffusion iterations +lambda_regDiff = 0.06; % regularisation for the diffusivity +sigmaPar = 0.04; % edge-preserving parameter +tau_param = 0.025; % time-marching constant +tic; u_diff = NonlDiff(single(u0), lambda_regDiff, sigmaPar, iter_diff, tau_param, 'Huber'); toc; +figure; imshow(u_diff, [0 1]); title('Diffusion denoised image (CPU)'); +%% +% fprintf('Denoise using Nonlinear-Diffusion model (GPU) \n'); +% iter_diff = 800; % number of diffusion iterations +% lambda_regDiff = 0.06; % regularisation for the diffusivity +% sigmaPar = 0.04; % edge-preserving parameter +% tau_param = 0.025; % time-marching constant +% tic; u_diff_g = NonlDiff_GPU(single(u0), lambda_regDiff, sigmaPar, iter_diff, tau_param, 'Huber'); toc; +% figure; imshow(u_diff_g, [0 1]); title('Diffusion denoised image (GPU)'); +%% %>>>>>>>>>>>>>> MULTI-CHANNEL priors <<<<<<<<<<<<<<< % fprintf('Denoise using the FGP-dTV model (CPU) \n'); diff --git a/Wrappers/Matlab/mex_compile/compileCPU_mex.m b/Wrappers/Matlab/mex_compile/compileCPU_mex.m index 9892d73..ec799bd 100644 --- a/Wrappers/Matlab/mex_compile/compileCPU_mex.m +++ b/Wrappers/Matlab/mex_compile/compileCPU_mex.m @@ -20,7 +20,10 @@ movefile FGP_dTV.mex* ../installed/ mex TNV.c TNV_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp" movefile TNV.mex* ../installed/ -delete SB_TV_core* ROF_TV_core* FGP_TV_core* FGP_dTV_core* TNV_core* utils* CCPiDefines.h +mex NonlDiff.c Diffusion_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp" +movefile NonlDiff.mex* ../installed/ + +delete SB_TV_core* ROF_TV_core* FGP_TV_core* FGP_dTV_core* TNV_core* utils* Diffusion_core* CCPiDefines.h fprintf('%s \n', 'All successfully compiled!'); diff --git a/Wrappers/Matlab/mex_compile/compileGPU_mex.m b/Wrappers/Matlab/mex_compile/compileGPU_mex.m index 3dbeb8a..55b51eb 100644 --- a/Wrappers/Matlab/mex_compile/compileGPU_mex.m +++ b/Wrappers/Matlab/mex_compile/compileGPU_mex.m @@ -31,7 +31,11 @@ movefile SB_TV_GPU.mex* ../installed/ mex -g -I/usr/local/cuda-7.5/include -L/usr/local/cuda-7.5/lib64 -lcudart -lcufft -lmwgpu FGP_dTV_GPU.cpp dTV_FGP_GPU_core.o movefile FGP_dTV_GPU.mex* ../installed/ -delete TV_ROF_GPU_core* TV_FGP_GPU_core* TV_SB_GPU_core* dTV_FGP_GPU_core* CCPiDefines.h +!/usr/local/cuda/bin/nvcc -O0 -c NonlDiff_GPU_core.cu -Xcompiler -fPIC -I~/SOFT/MATLAB9/extern/include/ +mex -g -I/usr/local/cuda-7.5/include -L/usr/local/cuda-7.5/lib64 -lcudart -lcufft -lmwgpu NonlDiff_GPU.cpp NonlDiff_GPU_core.o +movefile NonlDiff_GPU.mex* ../installed/ + +delete TV_ROF_GPU_core* TV_FGP_GPU_core* TV_SB_GPU_core* dTV_FGP_GPU_core* NonlDiff_GPU_core* CCPiDefines.h fprintf('%s \n', 'All successfully compiled!'); cd ../../ diff --git a/Wrappers/Matlab/mex_compile/regularisers_CPU/NonlDiff.c b/Wrappers/Matlab/mex_compile/regularisers_CPU/NonlDiff.c new file mode 100644 index 0000000..e05f5d4 --- /dev/null +++ b/Wrappers/Matlab/mex_compile/regularisers_CPU/NonlDiff.c @@ -0,0 +1,87 @@ +/* + * 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 2017 Daniil Kazantsev + * Copyright 2017 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 "Diffusion_core.h" + +/* C-OMP implementation of linear and nonlinear diffusion with the regularisation model [1] (2D/3D case) + * The minimisation is performed using explicit scheme. + * + * Input Parameters: + * 1. Noisy image/volume + * 2. lambda - regularization parameter + * 3. Edge-preserving parameter (sigma), when sigma equals to zero nonlinear diffusion -> linear diffusion + * 4. Number of iterations, for explicit scheme >= 150 is recommended [OPTIONAL parameter] + * 5. tau - time-marching step for explicit scheme [OPTIONAL parameter] + * 6. Penalty type: 1 - Huber, 2 - Perona-Malik, 3 - Tukey Biweight [OPTIONAL parameter] + * + * Output: + * [1] Regularized image/volume + * + * This function is based on the paper by + * [1] Perona, P. and Malik, J., 1990. Scale-space and edge detection using anisotropic diffusion. IEEE Transactions on pattern analysis and machine intelligence, 12(7), pp.629-639. + */ + +void mexFunction( + int nlhs, mxArray *plhs[], + int nrhs, const mxArray *prhs[]) + +{ + int number_of_dims, iter_numb, dimX, dimY, dimZ, penaltytype; + const int *dim_array; + float *Input, *Output=NULL, lambda, tau, sigma; + + dim_array = mxGetDimensions(prhs[0]); + number_of_dims = mxGetNumberOfDimensions(prhs[0]); + + /*Handling Matlab input data*/ + Input = (float *) mxGetData(prhs[0]); + lambda = (float) mxGetScalar(prhs[1]); /* regularization parameter */ + sigma = (float) mxGetScalar(prhs[2]); /* Edge-preserving parameter */ + iter_numb = 300; /* iterations number */ + tau = 0.025; /* marching step parameter */ + penaltytype = 1; /* Huber penalty by default */ + + if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); } + if ((nrhs < 3) || (nrhs > 6)) mexErrMsgTxt("At least 3 parameters is required, all parameters are: Image(2D/3D), Regularisation parameter, Edge-preserving parameter, iterations number, time-marching constant, penalty type - Huber, PM or Tukey"); + if ((nrhs == 4) || (nrhs == 5) || (nrhs == 6)) iter_numb = (int) mxGetScalar(prhs[3]); /* iterations number */ + if ((nrhs == 5) || (nrhs == 6)) tau = (float) mxGetScalar(prhs[4]); /* marching step parameter */ + if (nrhs == 6) { + char *penalty_type; + penalty_type = mxArrayToString(prhs[5]); /* Huber, PM or Tukey 'Huber' is the default */ + if ((strcmp(penalty_type, "Huber") != 0) && (strcmp(penalty_type, "PM") != 0) && (strcmp(penalty_type, "Tukey") != 0)) mexErrMsgTxt("Choose penalty: 'Huber', 'PM' or 'Tukey',"); + if (strcmp(penalty_type, "Huber") == 0) penaltytype = 1; /* enable 'Huber' penalty */ + if (strcmp(penalty_type, "PM") == 0) penaltytype = 2; /* enable Perona-Malik penalty */ + if (strcmp(penalty_type, "Tukey") == 0) penaltytype = 3; /* enable Tikey Biweight penalty */ + mxFree(penalty_type); + } + + /*Handling Matlab output data*/ + dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2]; + + /* output arrays*/ + if (number_of_dims == 2) { + dimZ = 1; /*2D case*/ + /* output image/volume */ + 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)); + + Diffusion_CPU_main(Input, Output, lambda, sigma, iter_numb, tau, penaltytype, dimX, dimY, dimZ); +}
\ No newline at end of file diff --git a/Wrappers/Matlab/mex_compile/regularisers_GPU/NonlDiff_GPU.cpp b/Wrappers/Matlab/mex_compile/regularisers_GPU/NonlDiff_GPU.cpp new file mode 100644 index 0000000..bfba9ea --- /dev/null +++ b/Wrappers/Matlab/mex_compile/regularisers_GPU/NonlDiff_GPU.cpp @@ -0,0 +1,90 @@ +/* + * 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 2017 Daniil Kazantsev + * Copyright 2017 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 <stdio.h> +#include <string.h> +#include "NonlDiff_GPU_core.h" + +/* CUDA implementation of linear and nonlinear diffusion with the regularisation model [1,2] (2D/3D case) + * The minimisation is performed using explicit scheme. + * + * Input Parameters: + * 1. Noisy image/volume + * 2. lambda - regularization parameter + * 3. Edge-preserving parameter (sigma), when sigma equals to zero nonlinear diffusion -> linear diffusion + * 4. Number of iterations, for explicit scheme >= 150 is recommended + * 5. tau - time-marching step for explicit scheme + * 6. Penalty type: 1 - Huber, 2 - Perona-Malik, 3 - Tukey Biweight + * + * Output: + * [1] Regularized image/volume + * + * This function is based on the paper by + * [1] Perona, P. and Malik, J., 1990. Scale-space and edge detection using anisotropic diffusion. IEEE Transactions on pattern analysis and machine intelligence, 12(7), pp.629-639. + * [2] Black, M.J., Sapiro, G., Marimont, D.H. and Heeger, D., 1998. Robust anisotropic diffusion. IEEE Transactions on image processing, 7(3), pp.421-432. + */ + +void mexFunction( + int nlhs, mxArray *plhs[], + int nrhs, const mxArray *prhs[]) + +{ + int number_of_dims, iter_numb, dimX, dimY, dimZ, penaltytype; + const int *dim_array; + float *Input, *Output=NULL, lambda, tau, sigma; + + dim_array = mxGetDimensions(prhs[0]); + number_of_dims = mxGetNumberOfDimensions(prhs[0]); + + /*Handling Matlab input data*/ + Input = (float *) mxGetData(prhs[0]); + lambda = (float) mxGetScalar(prhs[1]); /* regularization parameter */ + sigma = (float) mxGetScalar(prhs[2]); /* Edge-preserving parameter */ + iter_numb = 300; /* iterations number */ + tau = 0.025; /* marching step parameter */ + penaltytype = 1; /* Huber penalty by default */ + + if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); } + if ((nrhs < 3) || (nrhs > 6)) mexErrMsgTxt("At least 3 parameters is required, all parameters are: Image(2D/3D), Regularisation parameter, Edge-preserving parameter, iterations number, time-marching constant, penalty type - Huber, PM or Tukey"); + if ((nrhs == 4) || (nrhs == 5) || (nrhs == 6)) iter_numb = (int) mxGetScalar(prhs[3]); /* iterations number */ + if ((nrhs == 5) || (nrhs == 6)) tau = (float) mxGetScalar(prhs[4]); /* marching step parameter */ + if (nrhs == 6) { + char *penalty_type; + penalty_type = mxArrayToString(prhs[5]); /* Huber, PM or Tukey 'Huber' is the default */ + if ((strcmp(penalty_type, "Huber") != 0) && (strcmp(penalty_type, "PM") != 0) && (strcmp(penalty_type, "Tukey") != 0)) mexErrMsgTxt("Choose penalty: 'Huber', 'PM' or 'Tukey',"); + if (strcmp(penalty_type, "Huber") == 0) penaltytype = 1; /* enable 'Huber' penalty */ + if (strcmp(penalty_type, "PM") == 0) penaltytype = 2; /* enable Perona-Malik penalty */ + if (strcmp(penalty_type, "Tukey") == 0) penaltytype = 3; /* enable Tikey Biweight penalty */ + mxFree(penalty_type); + } + + /*Handling Matlab output data*/ + dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2]; + + /* output arrays*/ + if (number_of_dims == 2) { + dimZ = 1; /*2D case*/ + /* output image/volume */ + 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)); + + NonlDiff_GPU_main(Input, Output, lambda, sigma, iter_numb, tau, penaltytype, dimX, dimY, dimZ); +}
\ No newline at end of file diff --git a/Wrappers/Matlab/mex_compile/regularisers_GPU/ROF_TV_GPU.cpp b/Wrappers/Matlab/mex_compile/regularisers_GPU/ROF_TV_GPU.cpp index 7bbe3af..f60ba7b 100644 --- a/Wrappers/Matlab/mex_compile/regularisers_GPU/ROF_TV_GPU.cpp +++ b/Wrappers/Matlab/mex_compile/regularisers_GPU/ROF_TV_GPU.cpp @@ -37,7 +37,6 @@ * * D. Kazantsev, 2016-18 */ - void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) |