From 75917255b4f0b8aae2e6ce9492a75e0f749bfb3e Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <edo.paskino@gmail.com>
Date: Wed, 21 Feb 2018 12:49:59 +0000
Subject: added TV_ROF
---
Wrappers/Python/src/gpu_regularizers.pyx | 8 ++++----
1 file changed, 4 insertions(+), 4 deletions(-)
(limited to 'Wrappers/Python/src')
diff --git a/Wrappers/Python/src/gpu_regularizers.pyx b/Wrappers/Python/src/gpu_regularizers.pyx
index 5a5d274..fcb91cc 100644
--- a/Wrappers/Python/src/gpu_regularizers.pyx
+++ b/Wrappers/Python/src/gpu_regularizers.pyx
@@ -25,7 +25,7 @@ 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(float* A, float* B, int N, int M, int Z, int iter, float tau, 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 float pad_crop(float *A, float *Ap,
int OldSizeX, int OldSizeY, int OldSizeZ,
int NewSizeX, int NewSizeY, int NewSizeZ,
@@ -67,7 +67,7 @@ def NML(inputData,
h,
lambdaf)
-def ROF_TV_GPU(inputData,
+def GPU_ROF_TV(inputData,
iterations,
time_marching_parameter,
regularization_parameter):
@@ -343,7 +343,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(
+ TV_ROF_GPU_kernel(
&inputData[0,0], &B[0,0],
dims[0], dims[1], 0,
iterations ,
@@ -366,7 +366,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(
+ TV_ROF_GPU_kernel(
&inputData[0,0,0], &B[0,0,0],
dims[0], dims[1], dims[2],
iterations ,
--
cgit v1.2.3
From d4d12945f5396b10259fdeeb2b09f99b0e2c6afd Mon Sep 17 00:00:00 2001
From: algol <dkazanc@hotmail.com>
Date: Thu, 22 Feb 2018 12:34:39 +0000
Subject: GPU test fixed, cpu vs gpu test added for ROF-TV
---
Wrappers/Python/src/gpu_regularizers.pyx | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
(limited to 'Wrappers/Python/src')
diff --git a/Wrappers/Python/src/gpu_regularizers.pyx b/Wrappers/Python/src/gpu_regularizers.pyx
index fcb91cc..c724471 100644
--- a/Wrappers/Python/src/gpu_regularizers.pyx
+++ b/Wrappers/Python/src/gpu_regularizers.pyx
@@ -345,7 +345,7 @@ def ROFTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
# Running CUDA code here
TV_ROF_GPU_kernel(
&inputData[0,0], &B[0,0],
- dims[0], dims[1], 0,
+ dims[0], dims[1], 1,
iterations ,
time_marching_parameter,
regularization_parameter);
--
cgit v1.2.3
From bcdf186ccdca54a3df383512ad5a117004500a60 Mon Sep 17 00:00:00 2001
From: Daniil Kazantsev <dkazanc@hotmail.com>
Date: Mon, 26 Feb 2018 12:50:58 +0000
Subject: CPU-GPU naming consistency
---
Wrappers/Python/src/cpu_regularizers.cpp | 546 +++++++++++++++----------------
Wrappers/Python/src/gpu_regularizers.pyx | 69 +++-
2 files changed, 339 insertions(+), 276 deletions(-)
(limited to 'Wrappers/Python/src')
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
+
+
+
--
cgit v1.2.3
From 8082a76d4dfd9588590bab3fd26eae976b744a94 Mon Sep 17 00:00:00 2001
From: Daniil Kazantsev <dkazanc@hotmail.com>
Date: Mon, 26 Feb 2018 13:20:33 +0000
Subject: cpu-related cython file updated #33
---
Wrappers/Python/src/cpu_regularizers.pyx | 77 +++++++++++++++++++++++++++-----
Wrappers/Python/src/gpu_regularizers.pyx | 2 +-
2 files changed, 68 insertions(+), 11 deletions(-)
(limited to 'Wrappers/Python/src')
diff --git a/Wrappers/Python/src/cpu_regularizers.pyx b/Wrappers/Python/src/cpu_regularizers.pyx
index e7ff78c..b8089a8 100644
--- a/Wrappers/Python/src/cpu_regularizers.pyx
+++ b/Wrappers/Python/src/cpu_regularizers.pyx
@@ -18,20 +18,20 @@ import cython
import numpy as np
cimport numpy as np
-cdef extern float TV_ROF(float *A, float *B, int dimX, int dimY, int dimZ,
- int iterationsNumb, float tau, float flambda);
+cdef extern float TV_ROF_CPU(float *Input, float *Output, int dimX, int dimY, int dimZ, int iterationsNumb, float tau, float flambda);
+cdef extern float TV_FGP_CPU(float *Input, float *Output, float lambda, int iter, float epsil, int methodTV, int nonneg, int printM, int dimX, int dimY, int dimZ);
-
-def ROF_TV(inputData, iterations, regularization_parameter,
+# Can we use the same name here in "def" as the C function?
+def TV_ROF_CPU(inputData, iterations, regularization_parameter,
marching_step_parameter):
if inputData.ndim == 2:
- return ROF_TV_2D(inputData, iterations, regularization_parameter,
+ return TV_ROF_2D(inputData, iterations, regularization_parameter,
marching_step_parameter)
elif inputData.ndim == 3:
- return ROF_TV_3D(inputData, iterations, regularization_parameter,
+ return TV_ROF_3D(inputData, iterations, regularization_parameter,
marching_step_parameter)
-def ROF_TV_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
+def TV_ROF_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
int iterations,
float regularization_parameter,
float marching_step_parameter
@@ -45,13 +45,13 @@ def ROF_TV_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
np.zeros([dims[0],dims[1]], dtype='float32')
#/* Run ROF iterations for 2D data */
- TV_ROF(&inputData[0,0], &B[0,0], dims[0], dims[1], 1, iterations,
+ TV_ROF_CPU(&inputData[0,0], &B[0,0], dims[0], dims[1], 1, iterations,
marching_step_parameter, regularization_parameter)
return B
-def ROF_TV_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
+def TV_ROF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
int iterations,
float regularization_parameter,
float marching_step_parameter
@@ -65,8 +65,65 @@ def ROF_TV_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
np.zeros([dims[0],dims[1],dims[2]], dtype='float32')
#/* Run ROF iterations for 3D data */
- TV_ROF(&inputData[0,0, 0], &B[0,0, 0], dims[0], dims[1], dims[2], iterations,
+ TV_FGP_CPU(&inputData[0,0,0], &B[0,0,0], dims[0], dims[1], dims[2], iterations,
marching_step_parameter, regularization_parameter)
return B
+def TV_FGP_CPU(inputData, regularization_parameter, iterations, tolerance_param, methodTV, nonneg, printM):
+ if inputData.ndim == 2:
+ return TV_FGP_2D(inputData, regularization_parameter, iterations, tolerance_param, methodTV, nonneg, printM)
+ elif inputData.ndim == 3:
+ return TV_FGP_3D(inputData, regularization_parameter, iterations, tolerance_param, methodTV, nonneg, printM)
+
+def TV_FGP_2D(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')
+
+ #/* Run ROF iterations for 2D data */
+ TV_FGP_CPU(&inputData[0,0], &B[0,0], regularization_parameter,
+ iterations,
+ tolerance_param,
+ methodTV,
+ nonneg,
+ printM,
+ dims[0], dims[1], 1)
+
+ return B
+
+
+def TV_FGP_3D(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[2]
+ 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')
+
+ #/* Run ROF iterations for 3D data */
+ TV_FGP_CPU(&inputData[0,0, 0], &B[0,0, 0], iterations,
+ tolerance_param,
+ methodTV,
+ nonneg,
+ printM,
+ dims[0], dims[1], dims[2])
+
+ return B
diff --git a/Wrappers/Python/src/gpu_regularizers.pyx b/Wrappers/Python/src/gpu_regularizers.pyx
index 263fa4a..a14a20d 100644
--- a/Wrappers/Python/src/gpu_regularizers.pyx
+++ b/Wrappers/Python/src/gpu_regularizers.pyx
@@ -396,7 +396,7 @@ def TVFGP2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
# Running CUDA code here
TV_FGP_GPU(
&inputData[0,0], &B[0,0],
- regularization_parameter ,
+ regularization_parameter,
iterations,
tolerance_param,
methodTV,
--
cgit v1.2.3
From 74ff5b5f319b2f7ea3e078c62041ec8a1bb28335 Mon Sep 17 00:00:00 2001
From: algol <dkazanc@hotmail.com>
Date: Mon, 5 Mar 2018 18:12:01 +0000
Subject: Cmake/Cython fixes to compile ROF-FGP
---
Wrappers/Python/src/cpu_regularizers.cpp | 1 -
Wrappers/Python/src/cpu_regularizers.pyx | 25 ++++++++++++-------------
Wrappers/Python/src/gpu_regularizers.pyx | 2 +-
3 files changed, 13 insertions(+), 15 deletions(-)
(limited to 'Wrappers/Python/src')
diff --git a/Wrappers/Python/src/cpu_regularizers.cpp b/Wrappers/Python/src/cpu_regularizers.cpp
index 43d5d11..b8156d5 100644
--- a/Wrappers/Python/src/cpu_regularizers.cpp
+++ b/Wrappers/Python/src/cpu_regularizers.cpp
@@ -1040,7 +1040,6 @@ BOOST_PYTHON_MODULE(cpu_regularizers_boost)
np::dtype dt2 = np::dtype::get_builtin<uint16_t>();
def("SplitBregman_TV", SplitBregman_TV);
- def("FGP_TV", FGP_TV);
def("LLT_model", LLT_model);
def("PatchBased_Regul", PatchBased_Regul);
def("TGV_PD", TGV_PD);
diff --git a/Wrappers/Python/src/cpu_regularizers.pyx b/Wrappers/Python/src/cpu_regularizers.pyx
index b8089a8..448da31 100644
--- a/Wrappers/Python/src/cpu_regularizers.pyx
+++ b/Wrappers/Python/src/cpu_regularizers.pyx
@@ -18,8 +18,8 @@ import cython
import numpy as np
cimport numpy as np
-cdef extern float TV_ROF_CPU(float *Input, float *Output, int dimX, int dimY, int dimZ, int iterationsNumb, float tau, float flambda);
-cdef extern float TV_FGP_CPU(float *Input, float *Output, float lambda, int iter, float epsil, int methodTV, int nonneg, int printM, int dimX, int dimY, int dimZ);
+cdef extern float TV_ROF_CPU_main(float *Input, float *Output, int dimX, int dimY, int dimZ, int iterationsNumb, float tau, float flambda);
+cdef extern float TV_FGP_CPU_main(float *Input, float *Output, float lambdaPar, int iter, float epsil, int methodTV, int nonneg, int printM, int dimX, int dimY, int dimZ);
# Can we use the same name here in "def" as the C function?
def TV_ROF_CPU(inputData, iterations, regularization_parameter,
@@ -45,11 +45,10 @@ def TV_ROF_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
np.zeros([dims[0],dims[1]], dtype='float32')
#/* Run ROF iterations for 2D data */
- TV_ROF_CPU(&inputData[0,0], &B[0,0], dims[0], dims[1], 1, iterations,
+ TV_ROF_CPU_main(&inputData[0,0], &B[0,0], dims[0], dims[1], 1, iterations,
marching_step_parameter, regularization_parameter)
- return B
-
+ return B
def TV_ROF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
int iterations,
@@ -65,7 +64,7 @@ def TV_ROF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
np.zeros([dims[0],dims[1],dims[2]], dtype='float32')
#/* Run ROF iterations for 3D data */
- TV_FGP_CPU(&inputData[0,0,0], &B[0,0,0], dims[0], dims[1], dims[2], iterations,
+ TV_ROF_CPU_main(&inputData[0,0,0], &B[0,0,0], dims[0], dims[1], dims[2], iterations,
marching_step_parameter, regularization_parameter)
return B
@@ -88,11 +87,11 @@ def TV_FGP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
dims[0] = inputData.shape[0]
dims[1] = inputData.shape[1]
- cdef np.ndarray[np.float32_t, ndim=2, mode="c"] B = \
+ cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
np.zeros([dims[0],dims[1]], dtype='float32')
#/* Run ROF iterations for 2D data */
- TV_FGP_CPU(&inputData[0,0], &B[0,0], regularization_parameter,
+ TV_FGP_CPU_main(&inputData[0,0], &outputData[0,0], regularization_parameter,
iterations,
tolerance_param,
methodTV,
@@ -100,8 +99,7 @@ def TV_FGP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
printM,
dims[0], dims[1], 1)
- return B
-
+ return outputData
def TV_FGP_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
float regularization_parameter,
@@ -115,15 +113,16 @@ def TV_FGP_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
dims[1] = inputData.shape[1]
dims[2] = inputData.shape[2]
- cdef np.ndarray[np.float32_t, ndim=3, mode="c"] B = \
+ cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
np.zeros([dims[0],dims[1],dims[2]], dtype='float32')
#/* Run ROF iterations for 3D data */
- TV_FGP_CPU(&inputData[0,0, 0], &B[0,0, 0], iterations,
+ TV_FGP_CPU_main(&inputData[0,0,0], &outputData[0,0, 0], regularization_parameter,
+ iterations,
tolerance_param,
methodTV,
nonneg,
printM,
dims[0], dims[1], dims[2])
- return B
+ return outputData
diff --git a/Wrappers/Python/src/gpu_regularizers.pyx b/Wrappers/Python/src/gpu_regularizers.pyx
index a14a20d..e99bfa7 100644
--- a/Wrappers/Python/src/gpu_regularizers.pyx
+++ b/Wrappers/Python/src/gpu_regularizers.pyx
@@ -26,7 +26,7 @@ cdef extern void NLM_GPU_kernel(float *A, float* B, float *Eucl_Vec,
int SearchW, int SimilW,
int SearchW_real, float denh2, 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 void TV_FGP_GPU(float *Input, float *Output, float lambdaPar, 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,
--
cgit v1.2.3
From ccf9b61bba1004af783c6333d58ea9611c0f81f2 Mon Sep 17 00:00:00 2001
From: Daniil Kazantsev <dkazanc@hotmail.com>
Date: Tue, 6 Mar 2018 11:16:41 +0000
Subject: ROF_FGP_unified syntax
---
Wrappers/Python/src/cpu_regularizers.pyx | 59 +++++++++++++++-----------------
1 file changed, 27 insertions(+), 32 deletions(-)
(limited to 'Wrappers/Python/src')
diff --git a/Wrappers/Python/src/cpu_regularizers.pyx b/Wrappers/Python/src/cpu_regularizers.pyx
index 448da31..2654831 100644
--- a/Wrappers/Python/src/cpu_regularizers.pyx
+++ b/Wrappers/Python/src/cpu_regularizers.pyx
@@ -11,73 +11,68 @@ 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.
-Author: Edoardo Pasca
+Author: Edoardo Pasca, Daniil Kazantsev
"""
import cython
import numpy as np
cimport numpy as np
-cdef extern float TV_ROF_CPU_main(float *Input, float *Output, int dimX, int dimY, int dimZ, int iterationsNumb, float tau, float flambda);
-cdef extern float TV_FGP_CPU_main(float *Input, float *Output, float lambdaPar, int iter, float epsil, int methodTV, int nonneg, int printM, int dimX, int dimY, int dimZ);
+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);
-# Can we use the same name here in "def" as the C function?
-def TV_ROF_CPU(inputData, iterations, regularization_parameter,
+def TV_ROF_CPU(inputData, regularization_parameter, iterationsNumb
marching_step_parameter):
if inputData.ndim == 2:
- return TV_ROF_2D(inputData, iterations, regularization_parameter,
+ return TV_ROF_2D(inputData, regularization_parameter, iterationsNumb
marching_step_parameter)
elif inputData.ndim == 3:
- return TV_ROF_3D(inputData, iterations, regularization_parameter,
+ return TV_ROF_3D(inputData, regularization_parameter, iterationsNumb
marching_step_parameter)
def TV_ROF_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
- int iterations,
float regularization_parameter,
- float marching_step_parameter
- ):
-
+ int iterationsNumb,
+ float marching_step_parameter):
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 = \
+ cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
np.zeros([dims[0],dims[1]], dtype='float32')
- #/* Run ROF iterations for 2D data */
- TV_ROF_CPU_main(&inputData[0,0], &B[0,0], dims[0], dims[1], 1, iterations,
- marching_step_parameter, regularization_parameter)
+ # Run ROF iterations for 2D data
+ TV_ROF_CPU_main(&inputData[0,0], &outputData[0,0], regularization_parameter, iterationsNumb, marching_step_parameter, dims[0], dims[1], 1)
- return B
+ return outputData
def TV_ROF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
int iterations,
float regularization_parameter,
float marching_step_parameter
):
- cdef long dims[2]
+ 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 = \
+ cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
np.zeros([dims[0],dims[1],dims[2]], dtype='float32')
- #/* Run ROF iterations for 3D data */
- TV_ROF_CPU_main(&inputData[0,0,0], &B[0,0,0], dims[0], dims[1], dims[2], iterations,
- marching_step_parameter, regularization_parameter)
+ # Run ROF iterations for 3D data
+ TV_ROF_CPU_main(&inputData[0,0,0], &outputData[0,0,0], regularization_parameter, iterationsNumb, marching_step_parameter, dims[0], dims[1], dims[2])
- return B
+ return outputData
-def TV_FGP_CPU(inputData, regularization_parameter, iterations, tolerance_param, methodTV, nonneg, printM):
+def TV_FGP_CPU(inputData, regularization_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, printM):
if inputData.ndim == 2:
- return TV_FGP_2D(inputData, regularization_parameter, iterations, tolerance_param, methodTV, nonneg, printM)
+ return TV_FGP_2D(inputData, regularization_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, printM)
elif inputData.ndim == 3:
- return TV_FGP_3D(inputData, regularization_parameter, iterations, tolerance_param, methodTV, nonneg, printM)
+ return TV_FGP_3D(inputData, regularization_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, printM)
def TV_FGP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
float regularization_parameter,
- int iterations,
+ int iterationsNumb,
float tolerance_param,
int methodTV,
int nonneg,
@@ -92,7 +87,7 @@ def TV_FGP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
#/* Run ROF iterations for 2D data */
TV_FGP_CPU_main(&inputData[0,0], &outputData[0,0], regularization_parameter,
- iterations,
+ iterationsNumb,
tolerance_param,
methodTV,
nonneg,
@@ -103,22 +98,22 @@ def TV_FGP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
def TV_FGP_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
float regularization_parameter,
- int iterations,
+ int iterationsNumb,
float tolerance_param,
int methodTV,
int nonneg,
int printM):
- cdef long dims[2]
+ 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"] outputData = \
- np.zeros([dims[0],dims[1],dims[2]], dtype='float32')
+ np.zeros([dims[0], dims[1], dims[2]], dtype='float32')
#/* Run ROF iterations for 3D data */
- TV_FGP_CPU_main(&inputData[0,0,0], &outputData[0,0, 0], regularization_parameter,
- iterations,
+ TV_FGP_CPU_main(&inputData[0,0,0], &outputData[0,0,0], regularization_parameter,
+ iterationsNumb,
tolerance_param,
methodTV,
nonneg,
--
cgit v1.2.3
From 8d310478254f3cda63f3663729b416f425ad70b6 Mon Sep 17 00:00:00 2001
From: algol <dkazanc@hotmail.com>
Date: Tue, 6 Mar 2018 11:45:53 +0000
Subject: work on FGP intergration
---
Wrappers/Python/src/cpu_regularizers.pyx | 14 +++++---------
1 file changed, 5 insertions(+), 9 deletions(-)
(limited to 'Wrappers/Python/src')
diff --git a/Wrappers/Python/src/cpu_regularizers.pyx b/Wrappers/Python/src/cpu_regularizers.pyx
index 2654831..d62ca59 100644
--- a/Wrappers/Python/src/cpu_regularizers.pyx
+++ b/Wrappers/Python/src/cpu_regularizers.pyx
@@ -21,14 +21,11 @@ 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);
-def TV_ROF_CPU(inputData, regularization_parameter, iterationsNumb
- marching_step_parameter):
+def TV_ROF_CPU(inputData, regularization_parameter, iterationsNumb, marching_step_parameter):
if inputData.ndim == 2:
- return TV_ROF_2D(inputData, regularization_parameter, iterationsNumb
- marching_step_parameter)
+ return TV_ROF_2D(inputData, regularization_parameter, iterationsNumb, marching_step_parameter)
elif inputData.ndim == 3:
- return TV_ROF_3D(inputData, regularization_parameter, iterationsNumb
- marching_step_parameter)
+ return TV_ROF_3D(inputData, regularization_parameter, iterationsNumb, marching_step_parameter)
def TV_ROF_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
float regularization_parameter,
@@ -47,10 +44,9 @@ def TV_ROF_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
return outputData
def TV_ROF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
- int iterations,
+ int iterationsNumb,
float regularization_parameter,
- float marching_step_parameter
- ):
+ float marching_step_parameter):
cdef long dims[3]
dims[0] = inputData.shape[0]
dims[1] = inputData.shape[1]
--
cgit v1.2.3
From 309d84445b5ec2980db7c79832958a6481343f17 Mon Sep 17 00:00:00 2001
From: Daniil Kazantsev <dkazanc@hotmail.com>
Date: Tue, 6 Mar 2018 12:46:09 +0000
Subject: FGP/ROF updates
---
Wrappers/Python/src/cpu_regularizers.pyx | 5 ++---
1 file changed, 2 insertions(+), 3 deletions(-)
(limited to 'Wrappers/Python/src')
diff --git a/Wrappers/Python/src/cpu_regularizers.pyx b/Wrappers/Python/src/cpu_regularizers.pyx
index d62ca59..60e8627 100644
--- a/Wrappers/Python/src/cpu_regularizers.pyx
+++ b/Wrappers/Python/src/cpu_regularizers.pyx
@@ -93,7 +93,7 @@ def TV_FGP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
return outputData
def TV_FGP_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
- float regularization_parameter,
+ float regularization_parameter,
int iterationsNumb,
float tolerance_param,
int methodTV,
@@ -109,11 +109,10 @@ def TV_FGP_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
#/* Run ROF iterations for 3D data */
TV_FGP_CPU_main(&inputData[0,0,0], &outputData[0,0,0], regularization_parameter,
- iterationsNumb,
+ iterationsNumb,
tolerance_param,
methodTV,
nonneg,
printM,
dims[0], dims[1], dims[2])
-
return outputData
--
cgit v1.2.3
From 5411ebbd4165c81b398b010d6ad9d11d2e973aad Mon Sep 17 00:00:00 2001
From: Daniil Kazantsev <dkazanc@hotmail.com>
Date: Tue, 6 Mar 2018 14:40:11 +0000
Subject: work on cythonization2
---
Wrappers/Python/src/gpu_regularizers.pyx | 97 ++++++++++++++++++--------------
1 file changed, 56 insertions(+), 41 deletions(-)
(limited to 'Wrappers/Python/src')
diff --git a/Wrappers/Python/src/gpu_regularizers.pyx b/Wrappers/Python/src/gpu_regularizers.pyx
index e99bfa7..cb94e86 100644
--- a/Wrappers/Python/src/gpu_regularizers.pyx
+++ b/Wrappers/Python/src/gpu_regularizers.pyx
@@ -11,7 +11,7 @@ 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.
-Author: Edoardo Pasca
+Author: Edoardo Pasca, Daniil Kazantsev
"""
import cython
@@ -25,14 +25,16 @@ 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(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 lambdaPar, int iter, float epsil, int methodTV, int nonneg, int printM, int N, int M, int Z);
+cdef extern void TV_ROF_GPU_main(float* Input, float* Output, float lambdaPar, int iter, float tau, int N, int M, int Z);
+#cdef extern void TV_FGP_GPU(float *Input, float *Output, float lambdaPar, int iter, float epsil, int methodTV, int nonneg, int printM, int N, int M, int Z);
+# correct the function
+cdef extern void TV_FGP_GPU_main(float *Input, float *Output, float lambdaPar, 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,
int padXY, int switchpad_crop);
-
+#Diffusion 4th order regularizer
def Diff4thHajiaboli(inputData,
edge_preserv_parameter,
iterations,
@@ -50,7 +52,7 @@ def Diff4thHajiaboli(inputData,
iterations,
time_marching_parameter,
regularization_parameter)
-
+# patch-based nonlocal regularization
def NML(inputData,
SearchW_real,
SimilW,
@@ -68,23 +70,37 @@ def NML(inputData,
SimilW,
h,
lambdaf)
-
-def GPU_ROF_TV(inputData,
+# Total-variation Rudin-Osher-Fatemi (ROF)
+def TV_ROF_GPU(inputData,
+ regularization_parameter,
iterations,
- time_marching_parameter,
- regularization_parameter):
+ time_marching_parameter):
if inputData.ndim == 2:
return ROFTV2D(inputData,
- iterations,
- time_marching_parameter,
- regularization_parameter)
+ regularization_parameter,
+ iterations,
+ time_marching_parameter)
elif inputData.ndim == 3:
return ROFTV3D(inputData,
+ regularization_parameter
iterations,
- time_marching_parameter,
- regularization_parameter)
-
-
+ time_marching_parameter)
+# Total-variation Fast-Gradient-Projection (FGP)
+def TV_FGP_GPU(inputData,
+ regularization_parameter,
+ iterations,
+ time_marching_parameter):
+ if inputData.ndim == 2:
+ return FGPTV2D(inputData,
+ regularization_parameter,
+ iterations,
+ time_marching_parameter)
+ elif inputData.ndim == 3:
+ return FGPTV3D(inputData,
+ regularization_parameter
+ iterations,
+ time_marching_parameter)
+#****************************************************************#
def Diff4thHajiaboli2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
float edge_preserv_parameter,
int iterations,
@@ -333,52 +349,51 @@ def NML3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
return B
def ROFTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
+ float regularization_parameter,
int iterations,
- float time_marching_parameter,
- float regularization_parameter):
+ float time_marching_parameter):
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 = \
+ cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
np.zeros([dims[0],dims[1]], dtype='float32')
# Running CUDA code here
- TV_ROF_GPU(
- &inputData[0,0], &B[0,0],
- dims[0], dims[1], 1,
+ TV_ROF_GPU_main(
+ &inputData[0,0], &outputData[0,0],
+ regularization_parameter,
iterations ,
time_marching_parameter,
- regularization_parameter);
+ dims[0], dims[1], 1);
- return B
+ return outputData
def ROFTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
+ float regularization_parameter,
int iterations,
- float time_marching_parameter,
- float regularization_parameter):
+ float time_marching_parameter):
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 = \
+ cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
np.zeros([dims[0],dims[1],dims[2]], dtype='float32')
# Running CUDA code here
- TV_ROF_GPU(
- &inputData[0,0,0], &B[0,0,0],
- dims[0], dims[1], dims[2],
+ TV_ROF_GPU_main(
+ &inputData[0,0,0], &outputData[0,0,0],
+ regularization_parameter,
iterations ,
time_marching_parameter,
- regularization_parameter);
+ dims[0], dims[1], dims[2]);
- return B
-
+ return outputData
-def TVFGP2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
+def FGPTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
float regularization_parameter,
int iterations,
float tolerance_param,
@@ -390,12 +405,12 @@ def TVFGP2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
dims[0] = inputData.shape[0]
dims[1] = inputData.shape[1]
- cdef np.ndarray[np.float32_t, ndim=2, mode="c"] B = \
+ cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
np.zeros([dims[0],dims[1]], dtype='float32')
# Running CUDA code here
TV_FGP_GPU(
- &inputData[0,0], &B[0,0],
+ &inputData[0,0], &outputData[0,0],
regularization_parameter,
iterations,
tolerance_param,
@@ -404,9 +419,9 @@ def TVFGP2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
printM,
dims[0], dims[1], 1);
- return B
+ return outputData
-def TVFGP3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
+def FGPTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
float regularization_parameter,
int iterations,
float tolerance_param,
@@ -419,12 +434,12 @@ def TVFGP3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
dims[1] = inputData.shape[1]
dims[2] = inputData.shape[2]
- cdef np.ndarray[np.float32_t, ndim=3, mode="c"] B = \
+ cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
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],
+ &inputData[0,0,0], &outputData[0,0,0],
regularization_parameter ,
iterations,
tolerance_param,
@@ -433,7 +448,7 @@ def TVFGP3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
printM,
dims[0], dims[1], dims[2]);
- return B
+ return outputData
--
cgit v1.2.3
From 63ad8306c261a90c572d95084bf1dd8db2b3dce7 Mon Sep 17 00:00:00 2001
From: algol <dkazanc@hotmail.com>
Date: Tue, 6 Mar 2018 16:18:20 +0000
Subject: FGP/ROF-CPU/GPU fully working with new interfaces and regularizers.py
script. This closes #39 and #38
---
Wrappers/Python/src/gpu_regularizers.pyx | 43 +++++++++++++++++++-------------
1 file changed, 26 insertions(+), 17 deletions(-)
(limited to 'Wrappers/Python/src')
diff --git a/Wrappers/Python/src/gpu_regularizers.pyx b/Wrappers/Python/src/gpu_regularizers.pyx
index cb94e86..f96156a 100644
--- a/Wrappers/Python/src/gpu_regularizers.pyx
+++ b/Wrappers/Python/src/gpu_regularizers.pyx
@@ -26,14 +26,14 @@ cdef extern void NLM_GPU_kernel(float *A, float* B, float *Eucl_Vec,
int SearchW, int SimilW,
int SearchW_real, float denh2, float lambdaf);
cdef extern void TV_ROF_GPU_main(float* Input, float* Output, float lambdaPar, int iter, float tau, int N, int M, int Z);
-#cdef extern void TV_FGP_GPU(float *Input, float *Output, float lambdaPar, int iter, float epsil, int methodTV, int nonneg, int printM, int N, int M, int Z);
-# correct the function
cdef extern void TV_FGP_GPU_main(float *Input, float *Output, float lambdaPar, int iter, float epsil, int methodTV, int nonneg, int printM, int N, int M, int Z);
+# padding function
cdef extern float pad_crop(float *A, float *Ap,
int OldSizeX, int OldSizeY, int OldSizeZ,
int NewSizeX, int NewSizeY, int NewSizeZ,
int padXY, int switchpad_crop);
+
#Diffusion 4th order regularizer
def Diff4thHajiaboli(inputData,
edge_preserv_parameter,
@@ -70,6 +70,7 @@ def NML(inputData,
SimilW,
h,
lambdaf)
+
# Total-variation Rudin-Osher-Fatemi (ROF)
def TV_ROF_GPU(inputData,
regularization_parameter,
@@ -82,24 +83,34 @@ def TV_ROF_GPU(inputData,
time_marching_parameter)
elif inputData.ndim == 3:
return ROFTV3D(inputData,
- regularization_parameter
+ regularization_parameter,
iterations,
time_marching_parameter)
+
# Total-variation Fast-Gradient-Projection (FGP)
def TV_FGP_GPU(inputData,
regularization_parameter,
iterations,
- time_marching_parameter):
+ tolerance_param,
+ methodTV,
+ nonneg,
+ printM):
if inputData.ndim == 2:
- return FGPTV2D(inputData,
+ return FGPTV2D(inputData,
regularization_parameter,
- iterations,
- time_marching_parameter)
+ iterations,
+ tolerance_param,
+ methodTV,
+ nonneg,
+ printM)
elif inputData.ndim == 3:
- return FGPTV3D(inputData,
- regularization_parameter
+ return FGPTV3D(inputData,
+ regularization_parameter,
iterations,
- time_marching_parameter)
+ tolerance_param,
+ methodTV,
+ nonneg,
+ printM)
#****************************************************************#
def Diff4thHajiaboli2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
float edge_preserv_parameter,
@@ -347,7 +358,8 @@ def NML3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
switchpad_crop)
return B
-
+
+# Total-variation ROF
def ROFTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
float regularization_parameter,
int iterations,
@@ -393,6 +405,7 @@ def ROFTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
return outputData
+# Total-variation Fast-Gradient-Projection (FGP)
def FGPTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
float regularization_parameter,
int iterations,
@@ -409,8 +422,7 @@ def FGPTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
np.zeros([dims[0],dims[1]], dtype='float32')
# Running CUDA code here
- TV_FGP_GPU(
- &inputData[0,0], &outputData[0,0],
+ TV_FGP_GPU_main(&inputData[0,0], &outputData[0,0],
regularization_parameter,
iterations,
tolerance_param,
@@ -438,7 +450,7 @@ def FGPTV3D(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_FGP_GPU(
+ TV_FGP_GPU_main(
&inputData[0,0,0], &outputData[0,0,0],
regularization_parameter ,
iterations,
@@ -449,6 +461,3 @@ def FGPTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
dims[0], dims[1], dims[2]);
return outputData
-
-
-
--
cgit v1.2.3
From 4593040d228fec5ef1f1a301e511708d47f0d55e Mon Sep 17 00:00:00 2001
From: algol <dkazanc@hotmail.com>
Date: Tue, 6 Mar 2018 16:51:31 +0000
Subject: demos fixed #40
---
Wrappers/Python/src/cpu_regularizers.cpp | 289 -------------------------------
1 file changed, 289 deletions(-)
(limited to 'Wrappers/Python/src')
diff --git a/Wrappers/Python/src/cpu_regularizers.cpp b/Wrappers/Python/src/cpu_regularizers.cpp
index b8156d5..88b7c3c 100644
--- a/Wrappers/Python/src/cpu_regularizers.cpp
+++ b/Wrappers/Python/src/cpu_regularizers.cpp
@@ -303,292 +303,6 @@ 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) {
-
- //// 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;
-
- ////number_of_dims = mxGetNumberOfDimensions(prhs[0]);
- ////dim_array = mxGetDimensions(prhs[0]);
-
- //number_of_dims = input.get_nd();
- //int dim_array[3];
-
- //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 */
-
- ////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);
-
- ///* 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>();
-
- //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;
-//}
-
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
bp::list result;
@@ -1022,9 +736,6 @@ bp::list TGV_PD(np::ndarray input, double d_lambda, double d_alpha1, double d_al
result.append<np::ndarray>(npU);
}
-
-
-
return result;
}
--
cgit v1.2.3
From fbb7a12f7714978c251f02bf84ab9a66c762f428 Mon Sep 17 00:00:00 2001
From: algol <dkazanc@hotmail.com>
Date: Wed, 7 Mar 2018 10:15:13 +0000
Subject: small correct #40
---
Wrappers/Python/src/cpu_regularizers.cpp | 1 -
1 file changed, 1 deletion(-)
(limited to 'Wrappers/Python/src')
diff --git a/Wrappers/Python/src/cpu_regularizers.cpp b/Wrappers/Python/src/cpu_regularizers.cpp
index 88b7c3c..3529ebd 100644
--- a/Wrappers/Python/src/cpu_regularizers.cpp
+++ b/Wrappers/Python/src/cpu_regularizers.cpp
@@ -27,7 +27,6 @@ limitations under the License.
#include "boost/tuple/tuple.hpp"
#include "SplitBregman_TV_core.h"
-//#include "FGP_TV_core.h"
#include "LLT_model_core.h"
#include "PatchBased_Regul_core.h"
#include "TGV_PD_core.h"
--
cgit v1.2.3