summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Core/regularizers_CPU/FGP_TV_core.c139
-rw-r--r--Core/regularizers_CPU/FGP_TV_core.h8
-rw-r--r--Core/regularizers_CPU/utils.c13
-rw-r--r--Wrappers/Matlab/mex_compile/regularizers_CPU/FGP_TV.c142
-rw-r--r--Wrappers/Python/demo/test_cpu_regularizers.py10
-rw-r--r--Wrappers/Python/src/cpu_regularizers.pyx5
6 files changed, 89 insertions, 228 deletions
diff --git a/Core/regularizers_CPU/FGP_TV_core.c b/Core/regularizers_CPU/FGP_TV_core.c
index 9c0fcfc..7a8ce48 100644
--- a/Core/regularizers_CPU/FGP_TV_core.c
+++ b/Core/regularizers_CPU/FGP_TV_core.c
@@ -46,7 +46,7 @@ float TV_FGP_CPU_main(float *Input, float *Output, float lambdaPar, int iteratio
int count = 0;
if (dimZ <= 1) {
- /*2D case */
+ /*2D case */
float *Output_prev=NULL, *P1=NULL, *P2=NULL, *P1_prev=NULL, *P2_prev=NULL, *R1=NULL, *R2=NULL;
DimTotal = dimX*dimY;
@@ -65,28 +65,28 @@ float TV_FGP_CPU_main(float *Input, float *Output, float lambdaPar, int iteratio
Obj_func2D(Input, Output, R1, R2, lambdaPar, dimX, dimY);
/* apply nonnegativity */
- if (nonneg == 1) for(j=0; j<dimX*dimY; j++) {if (Output[j] < 0.0f) Output[j] = 0.0f;}
+ if (nonneg == 1) for(j=0; j<DimTotal; j++) {if (Output[j] < 0.0f) Output[j] = 0.0f;}
/*Taking a step towards minus of the gradient*/
Grad_func2D(P1, P2, Output, R1, R2, lambdaPar, dimX, dimY);
/* projection step */
- Proj_func2D(P1, P2, methodTV, dimX, dimY);
+ Proj_func2D(P1, P2, methodTV, DimTotal);
/*updating R and t*/
tkp1 = (1.0f + sqrt(1.0f + 4.0f*tk*tk))*0.5f;
- Rupd_func2D(P1, P1_prev, P2, P2_prev, R1, R2, tkp1, tk, dimX, dimY);
+ Rupd_func2D(P1, P1_prev, P2, P2_prev, R1, R2, tkp1, tk, DimTotal);
/* check early stopping criteria */
re = 0.0f; re1 = 0.0f;
- for(j=0; j<dimX*dimY; j++)
+ for(j=0; j<DimTotal; j++)
{
re += pow(Output[j] - Output_prev[j],2);
re1 += pow(Output[j],2);
}
re = sqrt(re)/sqrt(re1);
if (re < epsil) count++;
- if (count > 4) break;
+ if (count > 4) break;
/*storing old values*/
copyIm(Output, Output_prev, dimX, dimY, 1);
@@ -120,21 +120,21 @@ float TV_FGP_CPU_main(float *Input, float *Output, float lambdaPar, int iteratio
Obj_func3D(Input, Output, R1, R2, R3, lambdaPar, dimX, dimY, dimZ);
/* apply nonnegativity */
- if (nonneg == 1) for(j=0; j<dimX*dimY*dimZ; j++) {if (Output[j] < 0.0f) Output[j] = 0.0f;}
+ if (nonneg == 1) for(j=0; j<DimTotal; j++) {if (Output[j] < 0.0f) Output[j] = 0.0f;}
/*Taking a step towards minus of the gradient*/
Grad_func3D(P1, P2, P3, Output, R1, R2, R3, lambdaPar, dimX, dimY, dimZ);
/* projection step */
- Proj_func3D(P1, P2, P3, methodTV, dimX, dimY, dimZ);
+ Proj_func3D(P1, P2, P3, methodTV, DimTotal);
/*updating R and t*/
tkp1 = (1.0f + sqrt(1.0f + 4.0f*tk*tk))*0.5f;
- Rupd_func3D(P1, P1_prev, P2, P2_prev, P3, P3_prev, R1, R2, R3, tkp1, tk, dimX, dimY, dimZ);
+ Rupd_func3D(P1, P1_prev, P2, P2_prev, P3, P3_prev, R1, R2, R3, tkp1, tk, DimTotal);
/* calculate norm - stopping rules*/
re = 0.0f; re1 = 0.0f;
- for(j=0; j<dimX*dimY*dimZ; j++)
+ for(j=0; j<DimTotal; j++)
{
re += pow(Output[j] - Output_prev[j],2);
re1 += pow(Output[j],2);
@@ -189,52 +189,46 @@ float Grad_func2D(float *P1, float *P2, float *D, float *R1, float *R2, float la
}}
return 1;
}
-float Proj_func2D(float *P1, float *P2, int methTV, int dimX, int dimY)
+float Proj_func2D(float *P1, float *P2, int methTV, int DimTotal)
{
float val1, val2, denom, sq_denom;
- int i,j,index;
+ int i;
if (methTV == 0) {
/* isotropic TV*/
-#pragma omp parallel for shared(P1,P2) private(index,i,j,denom,sq_denom)
- for(i=0; i<dimX; i++) {
- for(j=0; j<dimY; j++) {
- index = j*dimX+i;
- denom = pow(P1[index],2) + pow(P2[index],2);
+#pragma omp parallel for shared(P1,P2) private(i,denom,sq_denom)
+ for(i=0; i<DimTotal; i++) {
+ denom = powf(P1[i],2) + powf(P2[i],2);
if (denom > 1.0f) {
- sq_denom = 1.0f/sqrt(denom);
- P1[index] = P1[index]*sq_denom;
- P2[index] = P2[index]*sq_denom;
+ sq_denom = 1.0f/sqrtf(denom);
+ P1[i] = P1[i]*sq_denom;
+ P2[i] = P2[i]*sq_denom;
}
- }}
+ }
}
else {
/* anisotropic TV*/
-#pragma omp parallel for shared(P1,P2) private(index,i,j,val1,val2)
- for(i=0; i<dimX; i++) {
- for(j=0; j<dimY; j++) {
- index = j*dimX+i;
- val1 = fabs(P1[index]);
- val2 = fabs(P2[index]);
+#pragma omp parallel for shared(P1,P2) private(i,val1,val2)
+ for(i=0; i<DimTotal; i++) {
+ val1 = fabs(P1[i]);
+ val2 = fabs(P2[i]);
if (val1 < 1.0f) {val1 = 1.0f;}
if (val2 < 1.0f) {val2 = 1.0f;}
- P1[index] = P1[index]/val1;
- P2[index] = P2[index]/val2;
- }}
+ P1[i] = P1[i]/val1;
+ P2[i] = P2[i]/val2;
+ }
}
return 1;
}
-float Rupd_func2D(float *P1, float *P1_old, float *P2, float *P2_old, float *R1, float *R2, float tkp1, float tk, int dimX, int dimY)
+float Rupd_func2D(float *P1, float *P1_old, float *P2, float *P2_old, float *R1, float *R2, float tkp1, float tk, int DimTotal)
{
- int i,j,index;
+ int i;
float multip;
multip = ((tk-1.0f)/tkp1);
-#pragma omp parallel for shared(P1,P2,P1_old,P2_old,R1,R2,multip) private(index,i,j)
- for(i=0; i<dimX; i++) {
- for(j=0; j<dimY; j++) {
- index = j*dimX+i;
- R1[index] = P1[index] + multip*(P1[index] - P1_old[index]);
- R2[index] = P2[index] + multip*(P2[index] - P2_old[index]);
- }}
+#pragma omp parallel for shared(P1,P2,P1_old,P2_old,R1,R2,multip) private(i)
+ for(i=0; i<DimTotal; i++) {
+ R1[i] = P1[i] + multip*(P1[i] - P1_old[i]);
+ R2[i] = P2[i] + multip*(P2[i] - P2_old[i]);
+ }
return 1;
}
@@ -248,7 +242,7 @@ float Obj_func3D(float *A, float *D, float *R1, float *R2, float *R3, float lamb
for(i=0; i<dimX; i++) {
for(j=0; j<dimY; j++) {
for(k=0; k<dimZ; k++) {
- index = (dimX*dimY)*k + j*dimX+i;
+ index = (dimX*dimY)*k + j*dimX+i;
/* boundary conditions */
if (i == 0) {val1 = 0.0f;} else {val1 = R1[(dimX*dimY)*k + j*dimX + (i-1)];}
if (j == 0) {val2 = 0.0f;} else {val2 = R2[(dimX*dimY)*k + (j-1)*dimX + i];}
@@ -277,59 +271,50 @@ float Grad_func3D(float *P1, float *P2, float *P3, float *D, float *R1, float *R
}}}
return 1;
}
-float Proj_func3D(float *P1, float *P2, float *P3, int methTV, int dimX, int dimY, int dimZ)
+float Proj_func3D(float *P1, float *P2, float *P3, int methTV, int DimTotal)
{
float val1, val2, val3, denom, sq_denom;
- int i,j,k, index;
+ int i;
if (methTV == 0) {
/* isotropic TV*/
- #pragma omp parallel for shared(P1,P2,P3) private(index,i,j,k,val1,val2,val3,sq_denom)
- for(i=0; i<dimX; i++) {
- for(j=0; j<dimY; j++) {
- for(k=0; k<dimZ; k++) {
- index = (dimX*dimY)*k + j*dimX+i;
- denom = pow(P1[index],2) + pow(P2[index],2) + pow(P3[index],2);
+ #pragma omp parallel for shared(P1,P2,P3) private(i,val1,val2,val3,sq_denom)
+ for(i=0; i<DimTotal; i++) {
+ denom = powf(P1[i],2) + powf(P2[i],2) + powf(P3[i],2);
if (denom > 1.0f) {
- sq_denom = 1.0f/sqrt(denom);
- P1[index] = P1[index]*sq_denom;
- P2[index] = P2[index]*sq_denom;
- P3[index] = P3[index]*sq_denom;
- }
- }}}
+ sq_denom = 1.0f/sqrtf(denom);
+ P1[i] = P1[i]*sq_denom;
+ P2[i] = P2[i]*sq_denom;
+ P3[i] = P3[i]*sq_denom;
+ }
+ }
}
else {
/* anisotropic TV*/
-#pragma omp parallel for shared(P1,P2,P3) private(index,i,j,k,val1,val2,val3)
- for(i=0; i<dimX; i++) {
- for(j=0; j<dimY; j++) {
- for(k=0; k<dimZ; k++) {
- index = (dimX*dimY)*k + j*dimX+i;
- val1 = fabs(P1[index]);
- val2 = fabs(P2[index]);
- val3 = fabs(P3[index]);
+#pragma omp parallel for shared(P1,P2,P3) private(i,val1,val2,val3)
+ for(i=0; i<DimTotal; i++) {
+ val1 = fabs(P1[i]);
+ val2 = fabs(P2[i]);
+ val3 = fabs(P3[i]);
if (val1 < 1.0f) {val1 = 1.0f;}
if (val2 < 1.0f) {val2 = 1.0f;}
if (val3 < 1.0f) {val3 = 1.0f;}
- P1[index] = P1[index]/val1;
- P2[index] = P2[index]/val2;
- P3[index] = P3[index]/val3;
- }}}
+ P1[i] = P1[i]/val1;
+ P2[i] = P2[i]/val2;
+ P3[i] = P3[i]/val3;
+ }
}
return 1;
}
-float Rupd_func3D(float *P1, float *P1_old, float *P2, float *P2_old, float *P3, float *P3_old, float *R1, float *R2, float *R3, float tkp1, float tk, int dimX, int dimY, int dimZ)
+float Rupd_func3D(float *P1, float *P1_old, float *P2, float *P2_old, float *P3, float *P3_old, float *R1, float *R2, float *R3, float tkp1, float tk, int DimTotal)
{
- int i,j,k,index;
+ int i;
float multip;
multip = ((tk-1.0f)/tkp1);
-#pragma omp parallel for shared(P1,P2,P3,P1_old,P2_old,P3_old,R1,R2,R3,multip) private(index,i,j,k)
- for(i=0; i<dimX; i++) {
- for(j=0; j<dimY; j++) {
- for(k=0; k<dimZ; k++) {
- index = (dimX*dimY)*k + j*dimX+i;
- R1[index] = P1[index] + multip*(P1[index] - P1_old[index]);
- R2[index] = P2[index] + multip*(P2[index] - P2_old[index]);
- R3[index] = P3[index] + multip*(P3[index] - P3_old[index]);
- }}}
+#pragma omp parallel for shared(P1,P2,P3,P1_old,P2_old,P3_old,R1,R2,R3,multip) private(i)
+ for(i=0; i<DimTotal; i++) {
+ R1[i] = P1[i] + multip*(P1[i] - P1_old[i]);
+ R2[i] = P2[i] + multip*(P2[i] - P2_old[i]);
+ R3[i] = P3[i] + multip*(P3[i] - P3_old[i]);
+ }
return 1;
}
diff --git a/Core/regularizers_CPU/FGP_TV_core.h b/Core/regularizers_CPU/FGP_TV_core.h
index 43a8519..37e32f7 100644
--- a/Core/regularizers_CPU/FGP_TV_core.h
+++ b/Core/regularizers_CPU/FGP_TV_core.h
@@ -51,13 +51,13 @@ CCPI_EXPORT float TV_FGP_CPU_main(float *Input, float *Output, float lambdaPar,
CCPI_EXPORT float Obj_func2D(float *A, float *D, float *R1, float *R2, float lambda, int dimX, int dimY);
CCPI_EXPORT float Grad_func2D(float *P1, float *P2, float *D, float *R1, float *R2, float lambda, int dimX, int dimY);
-CCPI_EXPORT float Proj_func2D(float *P1, float *P2, int methTV, int dimX, int dimY);
-CCPI_EXPORT float Rupd_func2D(float *P1, float *P1_old, float *P2, float *P2_old, float *R1, float *R2, float tkp1, float tk, int dimX, int dimY);
+CCPI_EXPORT float Proj_func2D(float *P1, float *P2, int methTV, int DimTotal);
+CCPI_EXPORT float Rupd_func2D(float *P1, float *P1_old, float *P2, float *P2_old, float *R1, float *R2, float tkp1, float tk, int DimTotal);
CCPI_EXPORT float Obj_func3D(float *A, float *D, float *R1, float *R2, float *R3, float lambda, int dimX, int dimY, int dimZ);
CCPI_EXPORT float Grad_func3D(float *P1, float *P2, float *P3, float *D, float *R1, float *R2, float *R3, float lambda, int dimX, int dimY, int dimZ);
-CCPI_EXPORT float Proj_func3D(float *P1, float *P2, float *P3, int methTV, int dimX, int dimY, int dimZ);
-CCPI_EXPORT float Rupd_func3D(float *P1, float *P1_old, float *P2, float *P2_old, float *P3, float *P3_old, float *R1, float *R2, float *R3, float tkp1, float tk, int dimX, int dimY, int dimZ);
+CCPI_EXPORT float Proj_func3D(float *P1, float *P2, float *P3, int methTV, int DimTotal);
+CCPI_EXPORT float Rupd_func3D(float *P1, float *P1_old, float *P2, float *P2_old, float *P3, float *P3_old, float *R1, float *R2, float *R3, float tkp1, float tk, int DimTotal);
#ifdef __cplusplus
}
#endif
diff --git a/Core/regularizers_CPU/utils.c b/Core/regularizers_CPU/utils.c
index cdf3d0e..0c02c45 100644
--- a/Core/regularizers_CPU/utils.c
+++ b/Core/regularizers_CPU/utils.c
@@ -18,6 +18,7 @@ limitations under the License.
*/
#include "utils.h"
+#include <math.h>
/* Copy Image */
float copyIm(float *A, float *U, int dimX, int dimY, int dimZ)
@@ -34,7 +35,7 @@ float copyIm(float *A, float *U, int dimX, int dimY, int dimZ)
float TV_energy2D(float *U, float *U0, float *E_val, float lambda, int dimX, int dimY)
{
int i, j, i1, j1, index;
- float NOMx_2, NOMy_2, E_Grad, E_Data;
+ float NOMx_2, NOMy_2, E_Grad=0.0f, E_Data=0.0f;
/* first calculate \grad U_xy*/
for(j=0; j<dimY; j++) {
@@ -44,11 +45,11 @@ float TV_energy2D(float *U, float *U0, float *E_val, float lambda, int dimX, int
i1 = i + 1; if (i == dimX-1) i1 = i;
j1 = j + 1; if (j == dimY-1) j1 = j;
- /* Forward differences */
- NOMx_2 = powf(U[j1*dimX + i] - U[index],2); /* x+ */
- NOMy_2 = powf(U[j*dimX + i1] - U[index],2); /* y+ */
- E_Grad += sqrtf(NOMx_2 + NOMy_2); /* gradient term energy */
- E_Data += 0.5f * lambda*(powf((U[index]-U0[index]),2)); /* fidelity term energy */
+ /* Forward differences */
+ NOMx_2 = powf((float)(U[j1*dimX + i] - U[index]),2); /* x+ */
+ NOMy_2 = powf((float)(U[j*dimX + i1] - U[index]),2); /* y+ */
+ E_Grad += sqrtf((float)(NOMx_2) + (float)(NOMy_2)); /* gradient term energy */
+ E_Data += 0.5f * lambda*(powf((float)(U[index]-U0[index]),2)); /* fidelity term energy */
}
}
E_val[0] = E_Grad + E_Data;
diff --git a/Wrappers/Matlab/mex_compile/regularizers_CPU/FGP_TV.c b/Wrappers/Matlab/mex_compile/regularizers_CPU/FGP_TV.c
index 30cea1a..7cc861a 100644
--- a/Wrappers/Matlab/mex_compile/regularizers_CPU/FGP_TV.c
+++ b/Wrappers/Matlab/mex_compile/regularizers_CPU/FGP_TV.c
@@ -53,9 +53,9 @@ void mexFunction(
int nrhs, const mxArray *prhs[])
{
- int number_of_dims, iter, dimX, dimY, dimZ, ll, j, count, methTV;
+ int number_of_dims, iter, dimX, dimY, dimZ, methTV;
const int *dim_array;
- 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, lambda, tk, tkp1, re, re1, re_old, epsil;
+ float *Input, *Output, lambda, epsil;
number_of_dims = mxGetNumberOfDimensions(prhs[0]);
dim_array = mxGetDimensions(prhs[0]);
@@ -63,9 +63,9 @@ void mexFunction(
/*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')");
- A = (float *) mxGetData(prhs[0]); /*noisy image (2D/3D) */
+ Input = (float *) mxGetData(prhs[0]); /*noisy image (2D/3D) */
lambda = (float) mxGetScalar(prhs[1]); /* regularization parameter */
- iter = 50; /* default iterations number */
+ iter = 300; /* default iterations number */
epsil = 0.0001; /* default tolerance constant */
methTV = 0; /* default isotropic TV penalty */
@@ -78,139 +78,17 @@ void mexFunction(
if (strcmp(penalty_type, "l1") == 0) methTV = 1; /* enable 'l1' penalty */
mxFree(penalty_type);
}
- /*output function value (last iteration) */
- plhs[1] = mxCreateNumericMatrix(1, 1, mxSINGLE_CLASS, mxREAL);
- float *funcvalA = (float *) mxGetData(plhs[1]);
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 = 0;
- re_old = 0.0f;
+ dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2];
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));
-
- /* 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_func_CALC2D(A, D, funcvalA, lambda, dimX, dimY);
- break; }
-
- /* check that the residual norm is decreasing */
- if (ll > 2) {
- if (re > re_old) {
- Obj_func_CALC2D(A, D, funcvalA, lambda, dimX, dimY);
- 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_func_CALC2D(A, D, funcvalA, lambda, dimX, dimY);
- }
- printf("FGP-TV iterations stopped at iteration %i with the function value %f \n", ll, funcvalA[0]);
- }
- 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));
-
- /* 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_func_CALC3D(A, D, funcvalA, lambda, dimX, dimY, dimZ);
- break;}
-
- /* check that the residual norm is decreasing */
- if (ll > 2) {
- if (re > re_old) {
- Obj_func_CALC3D(A, D, funcvalA, lambda, dimX, dimY, dimZ);
- }}
- 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_func_CALC3D(A, D, funcvalA, lambda, dimX, dimY, dimZ);
- }
- printf("FGP-TV iterations stopped at iteration %i with the function value %f \n", ll, funcvalA[0]);
- }
+ Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
+ TV_FGP_CPU_main(Input, Output, lambda, iter, epsil, methTV, 0, 0, dimX, dimY, dimZ);
+ }
+ if (number_of_dims == 3) {
+ }
}
diff --git a/Wrappers/Python/demo/test_cpu_regularizers.py b/Wrappers/Python/demo/test_cpu_regularizers.py
index f1eb3c3..7f08605 100644
--- a/Wrappers/Python/demo/test_cpu_regularizers.py
+++ b/Wrappers/Python/demo/test_cpu_regularizers.py
@@ -111,12 +111,10 @@ rms = rmse(Im, splitbregman)
pars['rmse'] = rms
txtstr = printParametersToString(pars)
txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-
+print (txtstr)
a=fig.add_subplot(2,4,2)
-
# these are matplotlib.patch.Patch properties
props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
# place a text box in upper left in axes coords
@@ -130,14 +128,14 @@ imgplot = plt.imshow(splitbregman,\
start_time = timeit.default_timer()
pars = {'algorithm' : TV_FGP_CPU , \
- 'input' : u0,
+ 'input' : u0,\
'regularization_parameter':0.07, \
'number_of_iterations' :300 ,\
'tolerance_constant':0.00001,\
'methodTV': 0 ,\
'nonneg': 0 ,\
- 'printingOut': 0
-}
+ 'printingOut': 0
+ }
out = TV_FGP_CPU (pars['input'],
pars['regularization_parameter'],
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
generated by cgit v1.2.3 (git 2.25.1) at 2024-12-04 07:30:34 +0000