diff options
| -rw-r--r-- | demos/demo_cpu_regularisers.py | 8 | ||||
| -rw-r--r-- | demos/demo_cpu_regularisers3D.py | 6 | ||||
| -rw-r--r-- | demos/demo_cpu_vs_gpu_regularisers.py | 12 | ||||
| -rw-r--r-- | demos/demo_gpu_regularisers.py | 6 | ||||
| -rw-r--r-- | demos/demo_gpu_regularisers3D.py | 6 | ||||
| -rw-r--r-- | src/Core/regularisers_CPU/PD_TV_core.c | 6 | ||||
| -rw-r--r-- | src/Core/regularisers_CPU/PD_TV_core.h | 2 | ||||
| -rw-r--r-- | src/Core/regularisers_GPU/TV_PD_GPU_core.cu | 13 | ||||
| -rw-r--r-- | src/Core/regularisers_GPU/TV_PD_GPU_core.h | 2 | ||||
| -rw-r--r-- | src/Matlab/mex_compile/regularisers_CPU/PD_TV.c | 50 | ||||
| -rw-r--r-- | src/Matlab/mex_compile/regularisers_GPU/PD_TV_GPU.cpp | 49 | ||||
| -rw-r--r-- | src/Python/ccpi/filters/regularisers.py | 8 | ||||
| -rw-r--r-- | src/Python/src/cpu_regularisers.pyx | 18 | ||||
| -rw-r--r-- | src/Python/src/gpu_regularisers.pyx | 16 | ||||
| -rw-r--r-- | test/test_CPU_regularisers.py | 2 | ||||
| -rwxr-xr-x | test/test_run_test.py | 9 | 
16 files changed, 90 insertions, 123 deletions
diff --git a/demos/demo_cpu_regularisers.py b/demos/demo_cpu_regularisers.py index 7d66b7f..50a5065 100644 --- a/demos/demo_cpu_regularisers.py +++ b/demos/demo_cpu_regularisers.py @@ -176,11 +176,10 @@ pars = {'algorithm' : PD_TV, \          'input' : u0,\          'regularisation_parameter':0.02, \          'number_of_iterations' :1500 ,\ -        'tolerance_constant':1e-08,\ +        'tolerance_constant':1e-06,\          'methodTV': 0 ,\          'nonneg': 1, -        'lipschitz_const' : 8, -        'tau' : 0.0025} +        'lipschitz_const' : 8}  print ("#############PD TV CPU####################")  start_time = timeit.default_timer() @@ -190,8 +189,7 @@ start_time = timeit.default_timer()                pars['tolerance_constant'],                 pars['methodTV'],                pars['nonneg'], -              pars['lipschitz_const'], -              pars['tau'],'cpu') +              pars['lipschitz_const'],'cpu')  Qtools = QualityTools(Im, pd_cpu)  pars['rmse'] = Qtools.rmse() diff --git a/demos/demo_cpu_regularisers3D.py b/demos/demo_cpu_regularisers3D.py index cfdd2d4..0e7e9be 100644 --- a/demos/demo_cpu_regularisers3D.py +++ b/demos/demo_cpu_regularisers3D.py @@ -188,8 +188,7 @@ pars = {'algorithm' : PD_TV, \          'tolerance_constant':1e-06,\          'methodTV': 0 ,\          'nonneg': 0, -        'lipschitz_const' : 8, -        'tau' : 0.0025} +        'lipschitz_const' : 8}  print ("#############FGP TV GPU####################")  start_time = timeit.default_timer() @@ -199,8 +198,7 @@ start_time = timeit.default_timer()                pars['tolerance_constant'],                 pars['methodTV'],                pars['nonneg'], -              pars['lipschitz_const'], -              pars['tau'],'cpu') +              pars['lipschitz_const'],'cpu')  Qtools = QualityTools(idealVol, pd_cpu3D)  pars['rmse'] = Qtools.rmse() diff --git a/demos/demo_cpu_vs_gpu_regularisers.py b/demos/demo_cpu_vs_gpu_regularisers.py index 015dfc6..a34bc19 100644 --- a/demos/demo_cpu_vs_gpu_regularisers.py +++ b/demos/demo_cpu_vs_gpu_regularisers.py @@ -241,8 +241,7 @@ pars = {'algorithm' : PD_TV, \          'tolerance_constant':0.0,\          'methodTV': 0 ,\          'nonneg': 0, -        'lipschitz_const' : 8, -        'tau' : 0.0025} +        'lipschitz_const' : 8}  print ("#############PD TV CPU####################")  start_time = timeit.default_timer() @@ -252,8 +251,7 @@ start_time = timeit.default_timer()                pars['tolerance_constant'],                 pars['methodTV'],                pars['nonneg'], -              pars['lipschitz_const'], -              pars['tau'],'cpu') +              pars['lipschitz_const'],'cpu')  Qtools = QualityTools(Im, pd_cpu)  pars['rmse'] = Qtools.rmse() @@ -279,8 +277,7 @@ pars = {'algorithm' : PD_TV, \          'tolerance_constant':0.0,\          'methodTV': 0 ,\          'nonneg': 0, -        'lipschitz_const' : 8, -        'tau' : 0.0025} +        'lipschitz_const' : 8}  print ("#############PD TV CPU####################")  start_time = timeit.default_timer() @@ -290,8 +287,7 @@ start_time = timeit.default_timer()                pars['tolerance_constant'],                 pars['methodTV'],                pars['nonneg'], -              pars['lipschitz_const'], -              pars['tau'],'gpu') +              pars['lipschitz_const'],'gpu')  Qtools = QualityTools(Im, pd_gpu)  pars['rmse'] = Qtools.rmse() diff --git a/demos/demo_gpu_regularisers.py b/demos/demo_gpu_regularisers.py index 5131847..c6114db 100644 --- a/demos/demo_gpu_regularisers.py +++ b/demos/demo_gpu_regularisers.py @@ -176,8 +176,7 @@ pars = {'algorithm' : PD_TV, \          'tolerance_constant':1e-06,\          'methodTV': 0 ,\          'nonneg': 1, -        'lipschitz_const' : 8, -        'tau' : 0.0025} +        'lipschitz_const' : 8}  print ("#############PD TV CPU####################")  start_time = timeit.default_timer() @@ -187,8 +186,7 @@ start_time = timeit.default_timer()                pars['tolerance_constant'],                 pars['methodTV'],                pars['nonneg'], -              pars['lipschitz_const'], -              pars['tau'],'gpu') +              pars['lipschitz_const'],'gpu')  Qtools = QualityTools(Im, pd_gpu)  pars['rmse'] = Qtools.rmse() diff --git a/demos/demo_gpu_regularisers3D.py b/demos/demo_gpu_regularisers3D.py index 2c25d01..18d97e5 100644 --- a/demos/demo_gpu_regularisers3D.py +++ b/demos/demo_gpu_regularisers3D.py @@ -192,8 +192,7 @@ pars = {'algorithm' : PD_TV, \          'tolerance_constant':1e-06,\          'methodTV': 0 ,\          'nonneg': 0, -        'lipschitz_const' : 8, -        'tau' : 0.0025} +        'lipschitz_const' : 8}  print ("#############PD TV GPU####################")  start_time = timeit.default_timer() @@ -203,8 +202,7 @@ start_time = timeit.default_timer()                pars['tolerance_constant'],                 pars['methodTV'],                pars['nonneg'], -              pars['lipschitz_const'], -              pars['tau'],'gpu') +              pars['lipschitz_const'],'gpu')  Qtools = QualityTools(idealVol, pd_gpu3D)  pars['rmse'] = Qtools.rmse() diff --git a/src/Core/regularisers_CPU/PD_TV_core.c b/src/Core/regularisers_CPU/PD_TV_core.c index 534091b..c1b21e7 100644 --- a/src/Core/regularisers_CPU/PD_TV_core.c +++ b/src/Core/regularisers_CPU/PD_TV_core.c @@ -29,7 +29,6 @@   * 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) - * 8. tau: time marching parameter   * Output:   * [1] TV - Filtered/regularized image/volume @@ -38,16 +37,17 @@   * [1] Antonin Chambolle, Thomas Pock. "A First-Order Primal-Dual Algorithm for Convex Problems with Applications to Imaging", 2010   */ -float PDTV_CPU_main(float *Input, float *U, float *infovector, float lambdaPar, int iterationsNumb, float epsil, float lipschitz_const, int methodTV, int nonneg, float tau, int dimX, int dimY, int dimZ) +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)  {      int ll;      long j, DimTotal; -    float re, re1, sigma, theta, lt; +    float re, re1, sigma, theta, lt, tau;      re = 0.0f; re1 = 0.0f;      int count = 0;      //tau = 1.0/powf(lipschitz_const,0.5);      //sigma = 1.0/powf(lipschitz_const,0.5); +    tau = lambdaPar*0.1f;      sigma = 1.0/(lipschitz_const*tau);      theta = 1.0f;      lt = tau/lambdaPar; diff --git a/src/Core/regularisers_CPU/PD_TV_core.h b/src/Core/regularisers_CPU/PD_TV_core.h index 97edc05..294e75c 100644 --- a/src/Core/regularisers_CPU/PD_TV_core.h +++ b/src/Core/regularisers_CPU/PD_TV_core.h @@ -47,7 +47,7 @@ limitations under the License.  #ifdef __cplusplus  extern "C" {  #endif -float PDTV_CPU_main(float *Input, float *U, float *infovector, float lambdaPar, int iterationsNumb, float epsil, float lipschitz_const, int methodTV, int nonneg, float tau, int dimX, int dimY, int dimZ); +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); diff --git a/src/Core/regularisers_GPU/TV_PD_GPU_core.cu b/src/Core/regularisers_GPU/TV_PD_GPU_core.cu index e57020a..01e0ab5 100644 --- a/src/Core/regularisers_GPU/TV_PD_GPU_core.cu +++ b/src/Core/regularisers_GPU/TV_PD_GPU_core.cu @@ -33,7 +33,6 @@ limitations under the License.   * 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) - * 8. tau: time marching parameter   * Output:   * [1] TV - Filtered/regularized image/volume @@ -42,8 +41,8 @@ limitations under the License.   * [1] Antonin Chambolle, Thomas Pock. "A First-Order Primal-Dual Algorithm for Convex Problems with Applications to Imaging", 2010   */ -#define BLKXSIZE2D 8 -#define BLKYSIZE2D 8 +#define BLKXSIZE2D 16 +#define BLKYSIZE2D 16  #define BLKXSIZE 8  #define BLKYSIZE 8 @@ -322,12 +321,9 @@ __global__ void PDResidCalc3D_kernel(float *Input1, float *Input2, float* Output          Output[index] = Input1[index] - Input2[index];      }  } - -  /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/ -  ////////////MAIN HOST FUNCTION /////////////// -extern "C" int TV_PD_GPU_main(float *Input, float *Output, float *infovector, float lambdaPar, int iter, float epsil, float lipschitz_const, int methodTV, int nonneg, float tau, int dimX, int dimY, int dimZ) +extern "C" int TV_PD_GPU_main(float *Input, float *Output, float *infovector, float lambdaPar, int iter, float epsil, float lipschitz_const, int methodTV, int nonneg, int dimX, int dimY, int dimZ)  {     int deviceCount = -1; // number of devices     cudaGetDeviceCount(&deviceCount); @@ -336,9 +332,10 @@ extern "C" int TV_PD_GPU_main(float *Input, float *Output, float *infovector, fl         return -1;     }     int count = 0, i; -   float re, sigma, theta, lt; +   float re, sigma, theta, lt, tau;     re = 0.0f; +   tau = lambdaPar*0.1f;     sigma = 1.0/(lipschitz_const*tau);     theta = 1.0f;     lt = tau/lambdaPar; diff --git a/src/Core/regularisers_GPU/TV_PD_GPU_core.h b/src/Core/regularisers_GPU/TV_PD_GPU_core.h index 2b123d9..48e353e 100644 --- a/src/Core/regularisers_GPU/TV_PD_GPU_core.h +++ b/src/Core/regularisers_GPU/TV_PD_GPU_core.h @@ -4,6 +4,6 @@  #include "CCPiDefines.h"  #include <memory.h> -extern "C" CCPI_EXPORT int TV_PD_GPU_main(float *Input, float *Output, float *infovector, float lambdaPar, int iter, float epsil, float lipschitz_const, int methodTV, int nonneg, float tau, int dimX, int dimY, int dimZ); +extern "C" CCPI_EXPORT int TV_PD_GPU_main(float *Input, float *Output, float *infovector, float lambdaPar, int iter, float epsil, float lipschitz_const, int methodTV, int nonneg, int dimX, int dimY, int dimZ);  #endif diff --git a/src/Matlab/mex_compile/regularisers_CPU/PD_TV.c b/src/Matlab/mex_compile/regularisers_CPU/PD_TV.c index e5ab1e4..f8f5272 100644 --- a/src/Matlab/mex_compile/regularisers_CPU/PD_TV.c +++ b/src/Matlab/mex_compile/regularisers_CPU/PD_TV.c @@ -30,8 +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   * [2] Information vector which contains [iteration no., reached tolerance] @@ -41,19 +40,19 @@  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; -     +    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"); -     +    if ((nrhs < 2) || (nrhs > 6)) 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 */ @@ -61,40 +60,39 @@ void mexFunction(      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))  { + +    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) || (nrhs == 8))  { +    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) || (nrhs == 8)) lipschitz_const = (float) mxGetScalar(prhs[6]); -    if (nrhs == 8) tau = (float) mxGetScalar(prhs[7]);    -     +    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];     -        +    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));         +        Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));      } -    if (number_of_dims == 3) {     +    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, tau, dimX, dimY, dimZ); + +    /* running the function */ +    PDTV_CPU_main(Input, Output, infovec, lambda, iter,  epsil, lipschitz_const, methTV, nonneg, 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 index e853dd3..2c037a5 100644 --- a/src/Matlab/mex_compile/regularisers_GPU/PD_TV_GPU.cpp +++ b/src/Matlab/mex_compile/regularisers_GPU/PD_TV_GPU.cpp @@ -30,8 +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   * [2] Information vector which contains [iteration no., reached tolerance] @@ -42,19 +41,19 @@  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; -     +    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"); -     +    if ((nrhs < 2) || (nrhs > 6)) 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 */ @@ -62,40 +61,38 @@ void mexFunction(      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))  { + +    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) || (nrhs == 8))  { +    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) || (nrhs == 8)) lipschitz_const = (float) mxGetScalar(prhs[6]); -    if (nrhs == 8) tau = (float) mxGetScalar(prhs[7]);    -     +    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];     -        +    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));         +        Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));      } -    if (number_of_dims == 3) {     +    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); + +    /* running the function */ +    TV_PD_GPU_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 5f4001a..e3a984e 100644 --- a/src/Python/ccpi/filters/regularisers.py +++ b/src/Python/ccpi/filters/regularisers.py @@ -53,7 +53,7 @@ def FGP_TV(inputData, regularisation_parameter,iterations,                           .format(device))  def PD_TV(inputData, regularisation_parameter, iterations, -                     tolerance_param, methodTV, nonneg, lipschitz_const, tau, device='cpu'): +                     tolerance_param, methodTV, nonneg, lipschitz_const, device='cpu'):      if device == 'cpu':          return TV_PD_CPU(inputData,                       regularisation_parameter, @@ -61,8 +61,7 @@ def PD_TV(inputData, regularisation_parameter, iterations,                       tolerance_param,                       methodTV,                       nonneg, -                     lipschitz_const, -                     tau) +                     lipschitz_const)      elif device == 'gpu' and gpu_enabled:          return TV_PD_GPU(inputData,                       regularisation_parameter, @@ -70,8 +69,7 @@ def PD_TV(inputData, regularisation_parameter, iterations,                       tolerance_param,                       methodTV,                       nonneg, -                     lipschitz_const, -                     tau) +                     lipschitz_const)      else:          if not gpu_enabled and device == 'gpu':              raise ValueError ('GPU is not available') diff --git a/src/Python/src/cpu_regularisers.pyx b/src/Python/src/cpu_regularisers.pyx index 8de6aea..08e247c 100644 --- a/src/Python/src/cpu_regularisers.pyx +++ b/src/Python/src/cpu_regularisers.pyx @@ -20,7 +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, float tau, 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); @@ -159,11 +159,11 @@ def TV_FGP_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,  #****************************************************************#  #****************** Total-variation Primal-dual *****************#  #****************************************************************# -def TV_PD_CPU(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, lipschitz_const, tau): +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, tau) +        return TV_PD_2D(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, lipschitz_const)      elif inputData.ndim == 3: -        return TV_PD_3D(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, lipschitz_const, tau) +        return TV_PD_3D(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, lipschitz_const)  def TV_PD_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,                       float regularisation_parameter, @@ -171,8 +171,7 @@ def TV_PD_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,                       float tolerance_param,                       int methodTV,                       int nonneg, -                     float lipschitz_const, -                     float tau): +                     float lipschitz_const):      cdef long dims[2]      dims[0] = inputData.shape[0] @@ -191,7 +190,6 @@ def TV_PD_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,                         lipschitz_const,                         methodTV,                         nonneg, -                       tau,                         dims[1],dims[0], 1)      return (outputData,infovec) @@ -200,9 +198,8 @@ def TV_PD_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,                       int iterationsNumb,                       float tolerance_param,                       int methodTV, -                     int nonneg,                      -                     float lipschitz_const, -                     float tau): +                     int nonneg, +                     float lipschitz_const):      cdef long dims[3]      dims[0] = inputData.shape[0] @@ -221,7 +218,6 @@ def TV_PD_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,                         lipschitz_const,                         methodTV,                         nonneg, -                       tau,                         dims[2], dims[1], dims[0])      return (outputData,infovec) diff --git a/src/Python/src/gpu_regularisers.pyx b/src/Python/src/gpu_regularisers.pyx index b22d15e..8a4568e 100644 --- a/src/Python/src/gpu_regularisers.pyx +++ b/src/Python/src/gpu_regularisers.pyx @@ -22,7 +22,7 @@ CUDAErrorMessage = 'CUDA error'  cdef extern int TV_ROF_GPU_main(float* Input, float* Output, float *infovector, float *lambdaPar, int lambda_is_arr, int iter, float tau, float epsil, int N, int M, int Z);  cdef extern int TV_FGP_GPU_main(float *Input, float *Output, float *infovector, float lambdaPar, int iter, float epsil, int methodTV, int nonneg, int N, int M, int Z); -cdef extern int TV_PD_GPU_main(float *Input, float *Output, float *infovector, float lambdaPar, int iter, float epsil, float lipschitz_const, int methodTV, int nonneg, float tau, int dimX, int dimY, int dimZ); +cdef extern int TV_PD_GPU_main(float *Input, float *Output, float *infovector, float lambdaPar, int iter, float epsil, float lipschitz_const, int methodTV, int nonneg, int dimX, int dimY, int dimZ);  cdef extern int TV_SB_GPU_main(float *Input, float *Output, float *infovector, float lambdaPar, int iter, float epsil, int methodTV, int N, int M, int Z);  cdef extern int LLT_ROF_GPU_main(float *Input, float *Output, float *infovector, float lambdaROF, float lambdaLLT, int iterationsNumb, float tau,  float epsil, int N, int M, int Z);  cdef extern int TGV_GPU_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); @@ -72,11 +72,11 @@ def TV_FGP_GPU(inputData,                       methodTV,                       nonneg)  # Total-variation Primal-Dual (PD) -def TV_PD_GPU(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, lipschitz_const, tau): +def TV_PD_GPU(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, lipschitz_const):      if inputData.ndim == 2: -        return TVPD2D(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, lipschitz_const, tau) +        return TVPD2D(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, lipschitz_const)      elif inputData.ndim == 3: -        return TVPD3D(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, lipschitz_const, tau) +        return TVPD3D(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, lipschitz_const)  def TVPD2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,                       float regularisation_parameter, @@ -84,8 +84,7 @@ def TVPD2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,                       float tolerance_param,                       int methodTV,                       int nonneg, -                     float lipschitz_const, -                     float tau): +                     float lipschitz_const):      cdef long dims[2]      dims[0] = inputData.shape[0] @@ -103,7 +102,6 @@ def TVPD2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,                         lipschitz_const,                         methodTV,                         nonneg, -                       tau,                         dims[1],dims[0], 1) ==0):          return (outputData,infovec)      else: @@ -115,8 +113,7 @@ def TVPD3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,                       float tolerance_param,                       int methodTV,                       int nonneg, -                     float lipschitz_const, -                     float tau): +                     float lipschitz_const):      cdef long dims[3]      dims[0] = inputData.shape[0] @@ -134,7 +131,6 @@ def TVPD3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,                         lipschitz_const,                         methodTV,                         nonneg, -                       tau,                         dims[2], dims[1], dims[0]) ==0):          return (outputData,infovec)      else: diff --git a/test/test_CPU_regularisers.py b/test/test_CPU_regularisers.py index 266ca8a..1eba479 100644 --- a/test/test_CPU_regularisers.py +++ b/test/test_CPU_regularisers.py @@ -42,7 +42,7 @@ class TestRegularisers(unittest.TestCase):      def test_PD_TV_CPU(self):          Im,input,ref = self.getPars() -        pd_cpu,info = PD_TV(input, 0.02, 300, 0.0, 0, 0, 8, 0.0025, 'cpu'); +        pd_cpu,info = PD_TV(input, 0.02, 300, 0.0, 0, 0, 8, 'cpu');          rms = rmse(Im, pd_cpu) diff --git a/test/test_run_test.py b/test/test_run_test.py index 1707aec..f562593 100755 --- a/test/test_run_test.py +++ b/test/test_run_test.py @@ -200,8 +200,7 @@ class TestRegularisers(unittest.TestCase):                  'tolerance_constant':0.0,\
                  'methodTV': 0 ,\
                  'nonneg': 0,
 -                'lipschitz_const' : 8,
 -                'tau' : 0.0025}
 +                'lipschitz_const' : 8}
          print ("#############PD TV CPU####################")
          start_time = timeit.default_timer()
 @@ -211,8 +210,7 @@ class TestRegularisers(unittest.TestCase):                        pars['tolerance_constant'], 
                        pars['methodTV'],
                        pars['nonneg'],
 -                      pars['lipschitz_const'],
 -                      pars['tau'],'cpu')
 +                      pars['lipschitz_const'],'cpu')
          rms = rmse(Im, pd_cpu)
          pars['rmse'] = rms
 @@ -230,8 +228,7 @@ class TestRegularisers(unittest.TestCase):                pars['tolerance_constant'], 
                pars['methodTV'],
                pars['nonneg'],
 -              pars['lipschitz_const'],
 -              pars['tau'],'gpu')
 +              pars['lipschitz_const'],'gpu')
          except ValueError as ve:
              self.skipTest("Results not comparable. GPU computing error.")
  | 
