summaryrefslogtreecommitdiffstats
path: root/Wrappers/Python
diff options
context:
space:
mode:
Diffstat (limited to 'Wrappers/Python')
-rw-r--r--Wrappers/Python/src/cpu_regularizers.cpp546
-rw-r--r--Wrappers/Python/src/gpu_regularizers.pyx69
-rw-r--r--Wrappers/Python/test/test_cpu_vs_gpu.py10
-rw-r--r--Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py8
4 files changed, 353 insertions, 280 deletions
diff --git a/Wrappers/Python/src/cpu_regularizers.cpp b/Wrappers/Python/src/cpu_regularizers.cpp
index e311570..43d5d11 100644
--- a/Wrappers/Python/src/cpu_regularizers.cpp
+++ b/Wrappers/Python/src/cpu_regularizers.cpp
@@ -27,7 +27,7 @@ limitations under the License.
#include "boost/tuple/tuple.hpp"
#include "SplitBregman_TV_core.h"
-#include "FGP_TV_core.h"
+//#include "FGP_TV_core.h"
#include "LLT_model_core.h"
#include "PatchBased_Regul_core.h"
#include "TGV_PD_core.h"
@@ -305,289 +305,289 @@ bp::list SplitBregman_TV(np::ndarray input, double d_mu, int iter, double d_epsi
-bp::list FGP_TV(np::ndarray input, double d_mu, int iter, double d_epsil, int methTV) {
+//bp::list FGP_TV(np::ndarray input, double d_mu, int iter, double d_epsil, int methTV) {
- // the result is in the following list
- bp::list result;
+ //// the result is in the following list
+ //bp::list result;
- int number_of_dims, dimX, dimY, dimZ, ll, j, count;
- float *A, *D = NULL, *D_old = NULL, *P1 = NULL, *P2 = NULL, *P3 = NULL, *P1_old = NULL, *P2_old = NULL, *P3_old = NULL, *R1 = NULL, *R2 = NULL, *R3 = NULL;
- float lambda, tk, tkp1, re, re1, re_old, epsil, funcval;
+ //int number_of_dims, dimX, dimY, dimZ, ll, j, count;
+ //float *A, *D = NULL, *D_old = NULL, *P1 = NULL, *P2 = NULL, *P3 = NULL, *P1_old = NULL, *P2_old = NULL, *P3_old = NULL, *R1 = NULL, *R2 = NULL, *R3 = NULL;
+ //float lambda, tk, tkp1, re, re1, re_old, epsil, funcval;
- //number_of_dims = mxGetNumberOfDimensions(prhs[0]);
- //dim_array = mxGetDimensions(prhs[0]);
-
- number_of_dims = input.get_nd();
- int dim_array[3];
+ ////number_of_dims = mxGetNumberOfDimensions(prhs[0]);
+ ////dim_array = mxGetDimensions(prhs[0]);
- dim_array[0] = input.shape(0);
- dim_array[1] = input.shape(1);
- if (number_of_dims == 2) {
- dim_array[2] = -1;
- }
- else {
- dim_array[2] = input.shape(2);
- }
- // Parameter handling is be done in Python
- ///*Handling Matlab input data*/
- //if ((nrhs < 2) || (nrhs > 5)) mexErrMsgTxt("At least 2 parameters is required: Image(2D/3D), Regularization parameter. The full list of parameters: Image(2D/3D), Regularization parameter, iterations number, tolerance, penalty type ('iso' or 'l1')");
-
- ///*Handling Matlab input data*/
- //A = (float *)mxGetData(prhs[0]); /*noisy image (2D/3D) */
- A = reinterpret_cast<float *>(input.get_data());
-
- //mu = (float)mxGetScalar(prhs[1]); /* regularization parameter */
- lambda = (float)d_mu;
-
- //iter = 35; /* default iterations number */
+ //number_of_dims = input.get_nd();
+ //int dim_array[3];
- //epsil = 0.0001; /* default tolerance constant */
- epsil = (float)d_epsil;
- //methTV = 0; /* default isotropic TV penalty */
- //if ((nrhs == 3) || (nrhs == 4) || (nrhs == 5)) iter = (int)mxGetScalar(prhs[2]); /* iterations number */
- //if ((nrhs == 4) || (nrhs == 5)) epsil = (float)mxGetScalar(prhs[3]); /* tolerance constant */
- //if (nrhs == 5) {
- // 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);
+ //dim_array[0] = input.shape(0);
+ //dim_array[1] = input.shape(1);
+ //if (number_of_dims == 2) {
+ //dim_array[2] = -1;
//}
- //if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) { mexErrMsgTxt("The input image must be in a single precision"); }
-
- //plhs[1] = mxCreateNumericMatrix(1, 1, mxSINGLE_CLASS, mxREAL);
- bp::tuple shape1 = bp::make_tuple(dim_array[0], dim_array[1]);
- np::dtype dtype = np::dtype::get_builtin<float>();
- np::ndarray out1 = np::zeros(shape1, dtype);
+ //else {
+ //dim_array[2] = input.shape(2);
+ //}
+ //// Parameter handling is be done in Python
+ /////*Handling Matlab input data*/
+ ////if ((nrhs < 2) || (nrhs > 5)) mexErrMsgTxt("At least 2 parameters is required: Image(2D/3D), Regularization parameter. The full list of parameters: Image(2D/3D), Regularization parameter, iterations number, tolerance, penalty type ('iso' or 'l1')");
+
+ /////*Handling Matlab input data*/
+ ////A = (float *)mxGetData(prhs[0]); /*noisy image (2D/3D) */
+ //A = reinterpret_cast<float *>(input.get_data());
+
+ ////mu = (float)mxGetScalar(prhs[1]); /* regularization parameter */
+ //lambda = (float)d_mu;
+
+ ////iter = 35; /* default iterations number */
+
+ ////epsil = 0.0001; /* default tolerance constant */
+ //epsil = (float)d_epsil;
+ ////methTV = 0; /* default isotropic TV penalty */
+ ////if ((nrhs == 3) || (nrhs == 4) || (nrhs == 5)) iter = (int)mxGetScalar(prhs[2]); /* iterations number */
+ ////if ((nrhs == 4) || (nrhs == 5)) epsil = (float)mxGetScalar(prhs[3]); /* tolerance constant */
+ ////if (nrhs == 5) {
+ //// 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 (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) { mexErrMsgTxt("The input image must be in a single precision"); }
+
+ ////plhs[1] = mxCreateNumericMatrix(1, 1, mxSINGLE_CLASS, mxREAL);
+ //bp::tuple shape1 = bp::make_tuple(dim_array[0], dim_array[1]);
+ //np::dtype dtype = np::dtype::get_builtin<float>();
+ //np::ndarray out1 = np::zeros(shape1, dtype);
- //float *funcvalA = (float *)mxGetData(plhs[1]);
- float * funcvalA = reinterpret_cast<float *>(out1.get_data());
- //if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) { mexErrMsgTxt("The input image must be in a single precision"); }
-
- /*Handling Matlab output data*/
- dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2];
-
- tk = 1.0f;
- tkp1 = 1.0f;
- count = 1;
- re_old = 0.0f;
-
- if (number_of_dims == 2) {
- dimZ = 1; /*2D case*/
- /*D = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
- D_old = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
- P1 = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
- P2 = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
- P1_old = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
- P2_old = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
- R1 = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
- R2 = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));*/
-
- bp::tuple shape = bp::make_tuple(dim_array[0], dim_array[1]);
- np::dtype dtype = np::dtype::get_builtin<float>();
-
-
- np::ndarray npD = np::zeros(shape, dtype);
- np::ndarray npD_old = np::zeros(shape, dtype);
- np::ndarray npP1 = np::zeros(shape, dtype);
- np::ndarray npP2 = np::zeros(shape, dtype);
- np::ndarray npP1_old = np::zeros(shape, dtype);
- np::ndarray npP2_old = np::zeros(shape, dtype);
- np::ndarray npR1 = np::zeros(shape, dtype);
- np::ndarray npR2 = np::zeros(shape, dtype);
-
- D = reinterpret_cast<float *>(npD.get_data());
- D_old = reinterpret_cast<float *>(npD_old.get_data());
- P1 = reinterpret_cast<float *>(npP1.get_data());
- P2 = reinterpret_cast<float *>(npP2.get_data());
- P1_old = reinterpret_cast<float *>(npP1_old.get_data());
- P2_old = reinterpret_cast<float *>(npP2_old.get_data());
- R1 = reinterpret_cast<float *>(npR1.get_data());
- R2 = reinterpret_cast<float *>(npR2.get_data());
-
- /* begin iterations */
- for (ll = 0; ll<iter; ll++) {
- /* computing the gradient of the objective function */
- Obj_func2D(A, D, R1, R2, lambda, dimX, dimY);
-
- /*Taking a step towards minus of the gradient*/
- Grad_func2D(P1, P2, D, R1, R2, lambda, dimX, dimY);
+ ////float *funcvalA = (float *)mxGetData(plhs[1]);
+ //float * funcvalA = reinterpret_cast<float *>(out1.get_data());
+ ////if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) { mexErrMsgTxt("The input image must be in a single precision"); }
+
+ ///*Handling Matlab output data*/
+ //dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2];
+
+ //tk = 1.0f;
+ //tkp1 = 1.0f;
+ //count = 1;
+ //re_old = 0.0f;
+
+ //if (number_of_dims == 2) {
+ //dimZ = 1; /*2D case*/
+ ///*D = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
+ //D_old = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
+ //P1 = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
+ //P2 = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
+ //P1_old = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
+ //P2_old = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
+ //R1 = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
+ //R2 = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));*/
+
+ //bp::tuple shape = bp::make_tuple(dim_array[0], dim_array[1]);
+ //np::dtype dtype = np::dtype::get_builtin<float>();
+
+
+ //np::ndarray npD = np::zeros(shape, dtype);
+ //np::ndarray npD_old = np::zeros(shape, dtype);
+ //np::ndarray npP1 = np::zeros(shape, dtype);
+ //np::ndarray npP2 = np::zeros(shape, dtype);
+ //np::ndarray npP1_old = np::zeros(shape, dtype);
+ //np::ndarray npP2_old = np::zeros(shape, dtype);
+ //np::ndarray npR1 = np::zeros(shape, dtype);
+ //np::ndarray npR2 = np::zeros(shape, dtype);
+
+ //D = reinterpret_cast<float *>(npD.get_data());
+ //D_old = reinterpret_cast<float *>(npD_old.get_data());
+ //P1 = reinterpret_cast<float *>(npP1.get_data());
+ //P2 = reinterpret_cast<float *>(npP2.get_data());
+ //P1_old = reinterpret_cast<float *>(npP1_old.get_data());
+ //P2_old = reinterpret_cast<float *>(npP2_old.get_data());
+ //R1 = reinterpret_cast<float *>(npR1.get_data());
+ //R2 = reinterpret_cast<float *>(npR2.get_data());
+
+ ///* begin iterations */
+ //for (ll = 0; ll<iter; ll++) {
+ ///* computing the gradient of the objective function */
+ //Obj_func2D(A, D, R1, R2, lambda, dimX, dimY);
+
+ ///*Taking a step towards minus of the gradient*/
+ //Grad_func2D(P1, P2, D, R1, R2, lambda, dimX, dimY);
- /* projection step */
- Proj_func2D(P1, P2, methTV, dimX, dimY);
+ ///* projection step */
+ //Proj_func2D(P1, P2, methTV, dimX, dimY);
- /*updating R and t*/
- tkp1 = (1.0f + sqrt(1.0f + 4.0f*tk*tk))*0.5f;
- Rupd_func2D(P1, P1_old, P2, P2_old, R1, R2, tkp1, tk, dimX, dimY);
-
- /* calculate norm */
- re = 0.0f; re1 = 0.0f;
- for (j = 0; j<dimX*dimY*dimZ; j++)
- {
- re += pow(D[j] - D_old[j], 2);
- re1 += pow(D[j], 2);
- }
- re = sqrt(re) / sqrt(re1);
- if (re < epsil) count++;
- if (count > 4) {
- Obj_func2D(A, D, P1, P2, lambda, dimX, dimY);
- funcval = 0.0f;
- for (j = 0; j<dimX*dimY*dimZ; j++) funcval += pow(D[j], 2);
- //funcvalA[0] = sqrt(funcval);
- float fv = sqrt(funcval);
- std::memcpy(funcvalA, &fv, sizeof(float));
- break;
- }
-
- /* check that the residual norm is decreasing */
- if (ll > 2) {
- if (re > re_old) {
- Obj_func2D(A, D, P1, P2, lambda, dimX, dimY);
- funcval = 0.0f;
- for (j = 0; j<dimX*dimY*dimZ; j++) funcval += pow(D[j], 2);
- //funcvalA[0] = sqrt(funcval);
- float fv = sqrt(funcval);
- std::memcpy(funcvalA, &fv, sizeof(float));
- break;
- }
- }
- re_old = re;
- /*printf("%f %i %i \n", re, ll, count); */
-
- /*storing old values*/
- copyIm(D, D_old, dimX, dimY, dimZ);
- copyIm(P1, P1_old, dimX, dimY, dimZ);
- copyIm(P2, P2_old, dimX, dimY, dimZ);
- tk = tkp1;
-
- /* calculating the objective function value */
- if (ll == (iter - 1)) {
- Obj_func2D(A, D, P1, P2, lambda, dimX, dimY);
- funcval = 0.0f;
- for (j = 0; j<dimX*dimY*dimZ; j++) funcval += pow(D[j], 2);
- //funcvalA[0] = sqrt(funcval);
- float fv = sqrt(funcval);
- std::memcpy(funcvalA, &fv, sizeof(float));
- }
- }
- //printf("FGP-TV iterations stopped at iteration %i with the function value %f \n", ll, funcvalA[0]);
- result.append<np::ndarray>(npD);
- result.append<np::ndarray>(out1);
- result.append<int>(ll);
- }
- if (number_of_dims == 3) {
- /*D = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
- D_old = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
- P1 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
- P2 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
- P3 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
- P1_old = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
- P2_old = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
- P3_old = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
- R1 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
- R2 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
- R3 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));*/
- bp::tuple shape = bp::make_tuple(dim_array[0], dim_array[1], dim_array[2]);
- np::dtype dtype = np::dtype::get_builtin<float>();
+ ///*updating R and t*/
+ //tkp1 = (1.0f + sqrt(1.0f + 4.0f*tk*tk))*0.5f;
+ //Rupd_func2D(P1, P1_old, P2, P2_old, R1, R2, tkp1, tk, dimX, dimY);
+
+ ///* calculate norm */
+ //re = 0.0f; re1 = 0.0f;
+ //for (j = 0; j<dimX*dimY*dimZ; j++)
+ //{
+ //re += pow(D[j] - D_old[j], 2);
+ //re1 += pow(D[j], 2);
+ //}
+ //re = sqrt(re) / sqrt(re1);
+ //if (re < epsil) count++;
+ //if (count > 4) {
+ //Obj_func2D(A, D, P1, P2, lambda, dimX, dimY);
+ //funcval = 0.0f;
+ //for (j = 0; j<dimX*dimY*dimZ; j++) funcval += pow(D[j], 2);
+ ////funcvalA[0] = sqrt(funcval);
+ //float fv = sqrt(funcval);
+ //std::memcpy(funcvalA, &fv, sizeof(float));
+ //break;
+ //}
+
+ ///* check that the residual norm is decreasing */
+ //if (ll > 2) {
+ //if (re > re_old) {
+ //Obj_func2D(A, D, P1, P2, lambda, dimX, dimY);
+ //funcval = 0.0f;
+ //for (j = 0; j<dimX*dimY*dimZ; j++) funcval += pow(D[j], 2);
+ ////funcvalA[0] = sqrt(funcval);
+ //float fv = sqrt(funcval);
+ //std::memcpy(funcvalA, &fv, sizeof(float));
+ //break;
+ //}
+ //}
+ //re_old = re;
+ ///*printf("%f %i %i \n", re, ll, count); */
+
+ ///*storing old values*/
+ //copyIm(D, D_old, dimX, dimY, dimZ);
+ //copyIm(P1, P1_old, dimX, dimY, dimZ);
+ //copyIm(P2, P2_old, dimX, dimY, dimZ);
+ //tk = tkp1;
+
+ ///* calculating the objective function value */
+ //if (ll == (iter - 1)) {
+ //Obj_func2D(A, D, P1, P2, lambda, dimX, dimY);
+ //funcval = 0.0f;
+ //for (j = 0; j<dimX*dimY*dimZ; j++) funcval += pow(D[j], 2);
+ ////funcvalA[0] = sqrt(funcval);
+ //float fv = sqrt(funcval);
+ //std::memcpy(funcvalA, &fv, sizeof(float));
+ //}
+ //}
+ ////printf("FGP-TV iterations stopped at iteration %i with the function value %f \n", ll, funcvalA[0]);
+ //result.append<np::ndarray>(npD);
+ //result.append<np::ndarray>(out1);
+ //result.append<int>(ll);
+ //}
+ //if (number_of_dims == 3) {
+ ///*D = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
+ //D_old = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
+ //P1 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
+ //P2 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
+ //P3 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
+ //P1_old = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
+ //P2_old = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
+ //P3_old = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
+ //R1 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
+ //R2 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
+ //R3 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));*/
+ //bp::tuple shape = bp::make_tuple(dim_array[0], dim_array[1], dim_array[2]);
+ //np::dtype dtype = np::dtype::get_builtin<float>();
- np::ndarray npD = np::zeros(shape, dtype);
- np::ndarray npD_old = np::zeros(shape, dtype);
- np::ndarray npP1 = np::zeros(shape, dtype);
- np::ndarray npP2 = np::zeros(shape, dtype);
- np::ndarray npP3 = np::zeros(shape, dtype);
- np::ndarray npP1_old = np::zeros(shape, dtype);
- np::ndarray npP2_old = np::zeros(shape, dtype);
- np::ndarray npP3_old = np::zeros(shape, dtype);
- np::ndarray npR1 = np::zeros(shape, dtype);
- np::ndarray npR2 = np::zeros(shape, dtype);
- np::ndarray npR3 = np::zeros(shape, dtype);
-
- D = reinterpret_cast<float *>(npD.get_data());
- D_old = reinterpret_cast<float *>(npD_old.get_data());
- P1 = reinterpret_cast<float *>(npP1.get_data());
- P2 = reinterpret_cast<float *>(npP2.get_data());
- P3 = reinterpret_cast<float *>(npP3.get_data());
- P1_old = reinterpret_cast<float *>(npP1_old.get_data());
- P2_old = reinterpret_cast<float *>(npP2_old.get_data());
- P3_old = reinterpret_cast<float *>(npP3_old.get_data());
- R1 = reinterpret_cast<float *>(npR1.get_data());
- R2 = reinterpret_cast<float *>(npR2.get_data());
- R3 = reinterpret_cast<float *>(npR3.get_data());
- /* begin iterations */
- for (ll = 0; ll<iter; ll++) {
- /* computing the gradient of the objective function */
- Obj_func3D(A, D, R1, R2, R3, lambda, dimX, dimY, dimZ);
- /*Taking a step towards minus of the gradient*/
- Grad_func3D(P1, P2, P3, D, R1, R2, R3, lambda, dimX, dimY, dimZ);
-
- /* projection step */
- Proj_func3D(P1, P2, P3, dimX, dimY, dimZ);
-
- /*updating R and t*/
- tkp1 = (1.0f + sqrt(1.0f + 4.0f*tk*tk))*0.5f;
- Rupd_func3D(P1, P1_old, P2, P2_old, P3, P3_old, R1, R2, R3, tkp1, tk, dimX, dimY, dimZ);
-
- /* calculate norm - stopping rules*/
- re = 0.0f; re1 = 0.0f;
- for (j = 0; j<dimX*dimY*dimZ; j++)
- {
- re += pow(D[j] - D_old[j], 2);
- re1 += pow(D[j], 2);
- }
- re = sqrt(re) / sqrt(re1);
- /* stop if the norm residual is less than the tolerance EPS */
- if (re < epsil) count++;
- if (count > 3) {
- Obj_func3D(A, D, P1, P2, P3, lambda, dimX, dimY, dimZ);
- funcval = 0.0f;
- for (j = 0; j<dimX*dimY*dimZ; j++) funcval += pow(D[j], 2);
- //funcvalA[0] = sqrt(funcval);
- float fv = sqrt(funcval);
- std::memcpy(funcvalA, &fv, sizeof(float));
- break;
- }
-
- /* check that the residual norm is decreasing */
- if (ll > 2) {
- if (re > re_old) {
- Obj_func3D(A, D, P1, P2, P3, lambda, dimX, dimY, dimZ);
- funcval = 0.0f;
- for (j = 0; j<dimX*dimY*dimZ; j++) funcval += pow(D[j], 2);
- //funcvalA[0] = sqrt(funcval);
- float fv = sqrt(funcval);
- std::memcpy(funcvalA, &fv, sizeof(float));
- break;
- }
- }
-
- re_old = re;
- /*printf("%f %i %i \n", re, ll, count); */
-
- /*storing old values*/
- copyIm(D, D_old, dimX, dimY, dimZ);
- copyIm(P1, P1_old, dimX, dimY, dimZ);
- copyIm(P2, P2_old, dimX, dimY, dimZ);
- copyIm(P3, P3_old, dimX, dimY, dimZ);
- tk = tkp1;
-
- if (ll == (iter - 1)) {
- Obj_func3D(A, D, P1, P2, P3, lambda, dimX, dimY, dimZ);
- funcval = 0.0f;
- for (j = 0; j<dimX*dimY*dimZ; j++) funcval += pow(D[j], 2);
- //funcvalA[0] = sqrt(funcval);
- float fv = sqrt(funcval);
- std::memcpy(funcvalA, &fv, sizeof(float));
- }
-
- }
- //printf("FGP-TV iterations stopped at iteration %i with the function value %f \n", ll, funcvalA[0]);
- result.append<np::ndarray>(npD);
- result.append<np::ndarray>(out1);
- result.append<int>(ll);
- }
+ //np::ndarray npD = np::zeros(shape, dtype);
+ //np::ndarray npD_old = np::zeros(shape, dtype);
+ //np::ndarray npP1 = np::zeros(shape, dtype);
+ //np::ndarray npP2 = np::zeros(shape, dtype);
+ //np::ndarray npP3 = np::zeros(shape, dtype);
+ //np::ndarray npP1_old = np::zeros(shape, dtype);
+ //np::ndarray npP2_old = np::zeros(shape, dtype);
+ //np::ndarray npP3_old = np::zeros(shape, dtype);
+ //np::ndarray npR1 = np::zeros(shape, dtype);
+ //np::ndarray npR2 = np::zeros(shape, dtype);
+ //np::ndarray npR3 = np::zeros(shape, dtype);
+
+ //D = reinterpret_cast<float *>(npD.get_data());
+ //D_old = reinterpret_cast<float *>(npD_old.get_data());
+ //P1 = reinterpret_cast<float *>(npP1.get_data());
+ //P2 = reinterpret_cast<float *>(npP2.get_data());
+ //P3 = reinterpret_cast<float *>(npP3.get_data());
+ //P1_old = reinterpret_cast<float *>(npP1_old.get_data());
+ //P2_old = reinterpret_cast<float *>(npP2_old.get_data());
+ //P3_old = reinterpret_cast<float *>(npP3_old.get_data());
+ //R1 = reinterpret_cast<float *>(npR1.get_data());
+ //R2 = reinterpret_cast<float *>(npR2.get_data());
+ //R3 = reinterpret_cast<float *>(npR3.get_data());
+ ///* begin iterations */
+ //for (ll = 0; ll<iter; ll++) {
+ ///* computing the gradient of the objective function */
+ //Obj_func3D(A, D, R1, R2, R3, lambda, dimX, dimY, dimZ);
+ ///*Taking a step towards minus of the gradient*/
+ //Grad_func3D(P1, P2, P3, D, R1, R2, R3, lambda, dimX, dimY, dimZ);
+
+ ///* projection step */
+ //Proj_func3D(P1, P2, P3, dimX, dimY, dimZ);
+
+ ///*updating R and t*/
+ //tkp1 = (1.0f + sqrt(1.0f + 4.0f*tk*tk))*0.5f;
+ //Rupd_func3D(P1, P1_old, P2, P2_old, P3, P3_old, R1, R2, R3, tkp1, tk, dimX, dimY, dimZ);
+
+ ///* calculate norm - stopping rules*/
+ //re = 0.0f; re1 = 0.0f;
+ //for (j = 0; j<dimX*dimY*dimZ; j++)
+ //{
+ //re += pow(D[j] - D_old[j], 2);
+ //re1 += pow(D[j], 2);
+ //}
+ //re = sqrt(re) / sqrt(re1);
+ ///* stop if the norm residual is less than the tolerance EPS */
+ //if (re < epsil) count++;
+ //if (count > 3) {
+ //Obj_func3D(A, D, P1, P2, P3, lambda, dimX, dimY, dimZ);
+ //funcval = 0.0f;
+ //for (j = 0; j<dimX*dimY*dimZ; j++) funcval += pow(D[j], 2);
+ ////funcvalA[0] = sqrt(funcval);
+ //float fv = sqrt(funcval);
+ //std::memcpy(funcvalA, &fv, sizeof(float));
+ //break;
+ //}
+
+ ///* check that the residual norm is decreasing */
+ //if (ll > 2) {
+ //if (re > re_old) {
+ //Obj_func3D(A, D, P1, P2, P3, lambda, dimX, dimY, dimZ);
+ //funcval = 0.0f;
+ //for (j = 0; j<dimX*dimY*dimZ; j++) funcval += pow(D[j], 2);
+ ////funcvalA[0] = sqrt(funcval);
+ //float fv = sqrt(funcval);
+ //std::memcpy(funcvalA, &fv, sizeof(float));
+ //break;
+ //}
+ //}
+
+ //re_old = re;
+ ///*printf("%f %i %i \n", re, ll, count); */
+
+ ///*storing old values*/
+ //copyIm(D, D_old, dimX, dimY, dimZ);
+ //copyIm(P1, P1_old, dimX, dimY, dimZ);
+ //copyIm(P2, P2_old, dimX, dimY, dimZ);
+ //copyIm(P3, P3_old, dimX, dimY, dimZ);
+ //tk = tkp1;
+
+ //if (ll == (iter - 1)) {
+ //Obj_func3D(A, D, P1, P2, P3, lambda, dimX, dimY, dimZ);
+ //funcval = 0.0f;
+ //for (j = 0; j<dimX*dimY*dimZ; j++) funcval += pow(D[j], 2);
+ ////funcvalA[0] = sqrt(funcval);
+ //float fv = sqrt(funcval);
+ //std::memcpy(funcvalA, &fv, sizeof(float));
+ //}
+
+ //}
+ ////printf("FGP-TV iterations stopped at iteration %i with the function value %f \n", ll, funcvalA[0]);
+ //result.append<np::ndarray>(npD);
+ //result.append<np::ndarray>(out1);
+ //result.append<int>(ll);
+ //}
- return result;
-}
+ //return result;
+//}
bp::list LLT_model(np::ndarray input, double d_lambda, double d_tau, int iter, double d_epsil, int switcher) {
// the result is in the following list
diff --git a/Wrappers/Python/src/gpu_regularizers.pyx b/Wrappers/Python/src/gpu_regularizers.pyx
index c724471..263fa4a 100644
--- a/Wrappers/Python/src/gpu_regularizers.pyx
+++ b/Wrappers/Python/src/gpu_regularizers.pyx
@@ -25,7 +25,9 @@ cdef extern void NLM_GPU_kernel(float *A, float* B, float *Eucl_Vec,
int N, int M, int Z, int dimension,
int SearchW, int SimilW,
int SearchW_real, float denh2, float lambdaf);
-cdef extern void TV_ROF_GPU_kernel(float* A, float* B, int N, int M, int Z, int iter, float tau, float lambdaf);
+cdef extern void TV_ROF_GPU(float* Input, float* Output, int N, int M, int Z, int iter, float tau, float lambdaf);
+cdef extern void TV_FGP_GPU(float *Input, float *Output, float lambda, int iter, float epsil, int methodTV, int nonneg, int printM, int N, int M, int Z);
+
cdef extern float pad_crop(float *A, float *Ap,
int OldSizeX, int OldSizeY, int OldSizeZ,
int NewSizeX, int NewSizeY, int NewSizeZ,
@@ -343,7 +345,7 @@ def ROFTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
np.zeros([dims[0],dims[1]], dtype='float32')
# Running CUDA code here
- TV_ROF_GPU_kernel(
+ TV_ROF_GPU(
&inputData[0,0], &B[0,0],
dims[0], dims[1], 1,
iterations ,
@@ -366,7 +368,7 @@ def ROFTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
np.zeros([dims[0],dims[1],dims[2]], dtype='float32')
# Running CUDA code here
- TV_ROF_GPU_kernel(
+ TV_ROF_GPU(
&inputData[0,0,0], &B[0,0,0],
dims[0], dims[1], dims[2],
iterations ,
@@ -374,3 +376,64 @@ def ROFTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
regularization_parameter);
return B
+
+
+def TVFGP2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
+ float regularization_parameter,
+ int iterations,
+ float tolerance_param,
+ int methodTV,
+ int nonneg,
+ int printM):
+
+ cdef long dims[2]
+ dims[0] = inputData.shape[0]
+ dims[1] = inputData.shape[1]
+
+ cdef np.ndarray[np.float32_t, ndim=2, mode="c"] B = \
+ np.zeros([dims[0],dims[1]], dtype='float32')
+
+ # Running CUDA code here
+ TV_FGP_GPU(
+ &inputData[0,0], &B[0,0],
+ regularization_parameter ,
+ iterations,
+ tolerance_param,
+ methodTV,
+ nonneg,
+ printM,
+ dims[0], dims[1], 1);
+
+ return B
+
+def TVFGP3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
+ float regularization_parameter,
+ int iterations,
+ float tolerance_param,
+ int methodTV,
+ int nonneg,
+ int printM):
+
+ cdef long dims[3]
+ dims[0] = inputData.shape[0]
+ dims[1] = inputData.shape[1]
+ dims[2] = inputData.shape[2]
+
+ cdef np.ndarray[np.float32_t, ndim=3, mode="c"] B = \
+ np.zeros([dims[0],dims[1],dims[2]], dtype='float32')
+
+ # Running CUDA code here
+ TV_FGP_GPU(
+ &inputData[0,0,0], &B[0,0,0],
+ regularization_parameter ,
+ iterations,
+ tolerance_param,
+ methodTV,
+ nonneg,
+ printM,
+ dims[0], dims[1], dims[2]);
+
+ return B
+
+
+
diff --git a/Wrappers/Python/test/test_cpu_vs_gpu.py b/Wrappers/Python/test/test_cpu_vs_gpu.py
new file mode 100644
index 0000000..74d65dd
--- /dev/null
+++ b/Wrappers/Python/test/test_cpu_vs_gpu.py
@@ -0,0 +1,10 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Feb 21 12:12:22 2018
+
+# CPU vs GPU comparison tests
+
+@author: algol
+"""
+
diff --git a/Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py b/Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py
index d742a7f..6344021 100644
--- a/Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py
+++ b/Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py
@@ -12,8 +12,8 @@ import matplotlib.pyplot as plt
import numpy as np
import os
import timeit
-from ccpi.filters.gpu_regularizers import Diff4thHajiaboli, NML, GPU_ROF_TV
-from ccpi.filters.cpu_regularizers_cython import ROF_TV
+from ccpi.filters.gpu_regularizers import Diff4thHajiaboli, NML, TV_ROF_GPU
+from ccpi.filters.cpu_regularizers_cython import TV_ROF_CPU
###############################################################################
def printParametersToString(pars):
txt = r''
@@ -64,7 +64,7 @@ pars = {'algorithm': ROF_TV , \
}
print ("#################ROF TV CPU#####################")
start_time = timeit.default_timer()
-rof_cpu = ROF_TV(pars['input'],
+rof_cpu = TV_ROF_CPU(pars['input'],
pars['number_of_iterations'],
pars['regularization_parameter'],
pars['time_marching_parameter']
@@ -89,7 +89,7 @@ plt.title('{}'.format('CPU results'))
print ("#################ROF TV GPU#####################")
start_time = timeit.default_timer()
-rof_gpu = GPU_ROF_TV(pars['input'],
+rof_gpu = TV_ROF_GPU(pars['input'],
pars['number_of_iterations'],
pars['time_marching_parameter'],
pars['regularization_parameter'])
generated by cgit v1.2.3 (git 2.25.1) at 2024-12-29 09:40:58 +0000