summaryrefslogtreecommitdiffstats
path: root/src/Matlab
diff options
context:
space:
mode:
Diffstat (limited to 'src/Matlab')
-rw-r--r--src/Matlab/mex_compile/compileGPU_mex.m6
-rw-r--r--src/Matlab/mex_compile/regularisers_CPU/PD_TV.c22
-rw-r--r--src/Matlab/mex_compile/regularisers_GPU/PD_TV_GPU.cpp101
3 files changed, 119 insertions, 10 deletions
diff --git a/src/Matlab/mex_compile/compileGPU_mex.m b/src/Matlab/mex_compile/compileGPU_mex.m
index 56fcd38..577630f 100644
--- a/src/Matlab/mex_compile/compileGPU_mex.m
+++ b/src/Matlab/mex_compile/compileGPU_mex.m
@@ -41,6 +41,11 @@ fprintf('%s \n', 'Compiling SB-TV...');
mex -g -I/usr/local/cuda-10.0/include -L/usr/local/cuda-10.0/lib64 -lcudart -lcufft -lmwgpu SB_TV_GPU.cpp TV_SB_GPU_core.o
movefile('SB_TV_GPU.mex*',Pathmove);
+fprintf('%s \n', 'Compiling PD-TV...');
+!/usr/local/cuda/bin/nvcc -O0 -c TV_PD_GPU_core.cu -Xcompiler -fPIC -I~/SOFT/MATLAB9/extern/include/
+mex -g -I/usr/local/cuda-10.0/include -L/usr/local/cuda-10.0/lib64 -lcudart -lcufft -lmwgpu PD_TV_GPU.cpp TV_PD_GPU_core.o
+movefile('PD_TV_GPU.mex*',Pathmove);
+
fprintf('%s \n', 'Compiling TGV...');
!/usr/local/cuda/bin/nvcc -O0 -c TGV_GPU_core.cu -Xcompiler -fPIC -I~/SOFT/MATLAB9/extern/include/
mex -g -I/usr/local/cuda-10.0/include -L/usr/local/cuda-10.0/lib64 -lcudart -lcufft -lmwgpu TGV_GPU.cpp TGV_GPU_core.o
@@ -72,6 +77,7 @@ mex -g -I/usr/local/cuda-10.0/include -L/usr/local/cuda-10.0/lib64 -lcudart -lcu
movefile('PatchSelect_GPU.mex*',Pathmove);
delete TV_ROF_GPU_core* TV_FGP_GPU_core* TV_SB_GPU_core* dTV_FGP_GPU_core* NonlDiff_GPU_core* Diffus_4thO_GPU_core* TGV_GPU_core* LLT_ROF_GPU_core* CCPiDefines.h
+delete TV_PD_GPU_core*
delete PatchSelect_GPU_core* Nonlocal_TV_core* shared.h
fprintf('%s \n', 'All successfully compiled!');
diff --git a/src/Matlab/mex_compile/regularisers_CPU/PD_TV.c b/src/Matlab/mex_compile/regularisers_CPU/PD_TV.c
index eac2d18..e5ab1e4 100644
--- a/src/Matlab/mex_compile/regularisers_CPU/PD_TV.c
+++ b/src/Matlab/mex_compile/regularisers_CPU/PD_TV.c
@@ -30,6 +30,7 @@
* 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
+ * 8. tau: convergence related parameter
* Output:
* [1] TV - Filtered/regularized image/volume
@@ -45,7 +46,7 @@ void mexFunction(
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;
+ float *Input, *infovec=NULL, *Output=NULL, lambda, epsil, lipschitz_const, tau;
number_of_dims = mxGetNumberOfDimensions(prhs[0]);
dim_array = mxGetDimensions(prhs[0]);
@@ -55,29 +56,30 @@ void mexFunction(
Input = (float *) mxGetData(prhs[0]); /*noisy image (2D/3D) */
lambda = (float) mxGetScalar(prhs[1]); /* regularization parameter */
- iter = 400; /* default iterations number */
+ iter = 500; /* 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 */
+ lipschitz_const = 8.0; /* lipschitz_const */
+ tau = 0.0025; /* tau convergence 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)) {
+ if ((nrhs == 3) || (nrhs == 4) || (nrhs == 5) || (nrhs == 6) || (nrhs == 7) || (nrhs == 8)) iter = (int) mxGetScalar(prhs[2]); /* iterations number */
+ if ((nrhs == 4) || (nrhs == 5) || (nrhs == 6) || (nrhs == 7) || (nrhs == 8)) epsil = (float) mxGetScalar(prhs[3]); /* tolerance constant */
+ if ((nrhs == 5) || (nrhs == 6) || (nrhs == 7) || (nrhs == 8)) {
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) || (nrhs == 7) || (nrhs == 8)) {
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]);
-
+ if ((nrhs == 7) || (nrhs == 8)) lipschitz_const = (float) mxGetScalar(prhs[6]);
+ if (nrhs == 8) tau = (float) mxGetScalar(prhs[7]);
/*Handling Matlab output data*/
dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2];
@@ -94,5 +96,5 @@ void mexFunction(
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);
+ PDTV_CPU_main(Input, Output, infovec, lambda, iter, epsil, lipschitz_const, methTV, nonneg, tau, dimX, dimY, dimZ);
}
diff --git a/src/Matlab/mex_compile/regularisers_GPU/PD_TV_GPU.cpp b/src/Matlab/mex_compile/regularisers_GPU/PD_TV_GPU.cpp
new file mode 100644
index 0000000..e853dd3
--- /dev/null
+++ b/src/Matlab/mex_compile/regularisers_GPU/PD_TV_GPU.cpp
@@ -0,0 +1,101 @@
+/*
+ * 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 "TV_PD_GPU_core.h"
+
+/* GPU 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
+ * 8. tau: 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, tau;
+
+ 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 = 500; /* 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 = 8.0; /* lipschitz_const */
+ tau = 0.0025; /* tau convergence 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) || (nrhs == 8)) iter = (int) mxGetScalar(prhs[2]); /* iterations number */
+ if ((nrhs == 4) || (nrhs == 5) || (nrhs == 6) || (nrhs == 7) || (nrhs == 8)) epsil = (float) mxGetScalar(prhs[3]); /* tolerance constant */
+ if ((nrhs == 5) || (nrhs == 6) || (nrhs == 7) || (nrhs == 8)) {
+ 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) || (nrhs == 8)) {
+ 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) || (nrhs == 8)) lipschitz_const = (float) mxGetScalar(prhs[6]);
+ if (nrhs == 8) tau = (float) mxGetScalar(prhs[7]);
+
+ /*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 */
+ TV_PD_GPU_main(Input, Output, infovec, lambda, iter, epsil, lipschitz_const, methTV, nonneg, tau, dimX, dimY, dimZ);
+}
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-09 23:33:27 +0000