summaryrefslogtreecommitdiffstats
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/Core/CMakeLists.txt2
-rw-r--r--src/Core/regularisers_CPU/FGP_TV_core.c34
-rw-r--r--src/Core/regularisers_CPU/FGP_TV_core.h1
-rw-r--r--src/Core/regularisers_CPU/FGP_dTV_core.c141
-rw-r--r--src/Core/regularisers_CPU/FGP_dTV_core.h2
-rw-r--r--src/Core/regularisers_CPU/PD_TV_core.c95
-rw-r--r--src/Core/regularisers_CPU/PD_TV_core.h5
-rw-r--r--src/Core/regularisers_CPU/utils.c37
-rw-r--r--src/Core/regularisers_CPU/utils.h1
-rw-r--r--src/Python/setup-regularisers.py.in1
-rw-r--r--src/Python/src/cpu_regularisers.pyx29
11 files changed, 196 insertions, 152 deletions
diff --git a/src/Core/CMakeLists.txt b/src/Core/CMakeLists.txt
index ac7ec91..c1bd7bb 100644
--- a/src/Core/CMakeLists.txt
+++ b/src/Core/CMakeLists.txt
@@ -66,7 +66,7 @@ message("Adding regularisers as a shared library")
add_library(cilreg SHARED
${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/FGP_TV_core.c
${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/SB_TV_core.c
- ${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/PD_TV_core.c
+ ${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/PD_TV_core.c
${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/TGV_core.c
${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/Diffusion_core.c
${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/Diffus4th_order_core.c
diff --git a/src/Core/regularisers_CPU/FGP_TV_core.c b/src/Core/regularisers_CPU/FGP_TV_core.c
index ff67af2..12f2254 100644
--- a/src/Core/regularisers_CPU/FGP_TV_core.c
+++ b/src/Core/regularisers_CPU/FGP_TV_core.c
@@ -254,40 +254,6 @@ 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, long DimTotal)
-{
- float val1, val2, val3, denom, sq_denom;
- long i;
- if (methTV == 0) {
- /* isotropic TV*/
-#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/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(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[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, long DimTotal)
{
long i;
diff --git a/src/Core/regularisers_CPU/FGP_TV_core.h b/src/Core/regularisers_CPU/FGP_TV_core.h
index 4466a35..e9c3cde 100644
--- a/src/Core/regularisers_CPU/FGP_TV_core.h
+++ b/src/Core/regularisers_CPU/FGP_TV_core.h
@@ -56,7 +56,6 @@ CCPI_EXPORT float Rupd_func2D(float *P1, float *P1_old, float *P2, float *P2_old
CCPI_EXPORT float Obj_func3D(float *A, float *D, float *R1, float *R2, float *R3, float lambda, long dimX, long dimY, long dimZ);
CCPI_EXPORT float Grad_func3D(float *P1, float *P2, float *P3, float *D, float *R1, float *R2, float *R3, float lambda, long dimX, long dimY, long dimZ);
-CCPI_EXPORT float Proj_func3D(float *P1, float *P2, float *P3, int methTV, long 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, long DimTotal);
#ifdef __cplusplus
}
diff --git a/src/Core/regularisers_CPU/FGP_dTV_core.c b/src/Core/regularisers_CPU/FGP_dTV_core.c
index afd7264..e828be6 100644
--- a/src/Core/regularisers_CPU/FGP_dTV_core.c
+++ b/src/Core/regularisers_CPU/FGP_dTV_core.c
@@ -52,11 +52,11 @@ float dTV_FGP_CPU_main(float *Input, float *InputRef, float *Output, float *info
float tk = 1.0f;
float tkp1=1.0f;
int count = 0;
-
-
+
+
float *Output_prev=NULL, *P1=NULL, *P2=NULL, *P1_prev=NULL, *P2_prev=NULL, *R1=NULL, *R2=NULL, *InputRef_x=NULL, *InputRef_y=NULL;
DimTotal = (long)(dimX*dimY*dimZ);
-
+
if (epsil != 0.0f) Output_prev = calloc(DimTotal, sizeof(float));
P1 = calloc(DimTotal, sizeof(float));
P2 = calloc(DimTotal, sizeof(float));
@@ -66,39 +66,39 @@ float dTV_FGP_CPU_main(float *Input, float *InputRef, float *Output, float *info
R2 = calloc(DimTotal, sizeof(float));
InputRef_x = calloc(DimTotal, sizeof(float));
InputRef_y = calloc(DimTotal, sizeof(float));
-
+
if (dimZ <= 1) {
/*2D case */
/* calculate gradient field (smoothed) for the reference image */
GradNorm_func2D(InputRef, InputRef_x, InputRef_y, eta, (long)(dimX), (long)(dimY));
-
+
/* begin iterations */
for(ll=0; ll<iterationsNumb; ll++) {
-
+
if ((epsil != 0.0f) && (ll % 5 == 0)) copyIm(Output, Output_prev, (long)(dimX), (long)(dimY), 1l);
/*projects a 2D vector field R-1,2 onto the orthogonal complement of another 2D vector field InputRef_xy*/
ProjectVect_func2D(R1, R2, InputRef_x, InputRef_y, (long)(dimX), (long)(dimY));
-
+
/* computing the gradient of the objective function */
Obj_dfunc2D(Input, Output, R1, R2, lambdaPar, (long)(dimX), (long)(dimY));
-
+
/* apply nonnegativity */
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_dfunc2D(P1, P2, Output, R1, R2, InputRef_x, InputRef_y, lambdaPar, (long)(dimX), (long)(dimY));
-
+
/* projection step */
- Proj_dfunc2D(P1, P2, methodTV, DimTotal);
-
+ Proj_func2D(P1, P2, methodTV, DimTotal);
+
/*updating R and t*/
tkp1 = (1.0f + sqrt(1.0f + 4.0f*tk*tk))*0.5f;
Rupd_dfunc2D(P1, P1_prev, P2, P2_prev, R1, R2, tkp1, tk, DimTotal);
-
+
copyIm(P1, P1_prev, (long)(dimX), (long)(dimY), 1l);
copyIm(P2, P2_prev, (long)(dimX), (long)(dimY), 1l);
tk = tkp1;
-
+
/* check early stopping criteria */
if ((epsil != 0.0f) && (ll % 5 == 0)) {
re = 0.0f; re1 = 0.0f;
@@ -116,45 +116,45 @@ float dTV_FGP_CPU_main(float *Input, float *InputRef, float *Output, float *info
else {
/*3D case*/
float *P3=NULL, *P3_prev=NULL, *R3=NULL, *InputRef_z=NULL;
-
+
P3 = calloc(DimTotal, sizeof(float));
P3_prev = calloc(DimTotal, sizeof(float));
R3 = calloc(DimTotal, sizeof(float));
InputRef_z = calloc(DimTotal, sizeof(float));
-
+
/* calculate gradient field (smoothed) for the reference volume */
GradNorm_func3D(InputRef, InputRef_x, InputRef_y, InputRef_z, eta, (long)(dimX), (long)(dimY), (long)(dimZ));
-
+
/* begin iterations */
for(ll=0; ll<iterationsNumb; ll++) {
-
+
if ((epsil != 0.0f) && (ll % 5 == 0)) copyIm(Output, Output_prev, (long)(dimX), (long)(dimY), (long)(dimZ));
-
+
/*projects a 3D vector field R-1,2,3 onto the orthogonal complement of another 3D vector field InputRef_xyz*/
ProjectVect_func3D(R1, R2, R3, InputRef_x, InputRef_y, InputRef_z, (long)(dimX), (long)(dimY), (long)(dimZ));
-
+
/* computing the gradient of the objective function */
Obj_dfunc3D(Input, Output, R1, R2, R3, lambdaPar, (long)(dimX), (long)(dimY), (long)(dimZ));
-
+
/* apply nonnegativity */
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_dfunc3D(P1, P2, P3, Output, R1, R2, R3, InputRef_x, InputRef_y, InputRef_z, lambdaPar, (long)(dimX), (long)(dimY), (long)(dimZ));
-
+
/* projection step */
- Proj_dfunc3D(P1, P2, P3, methodTV, DimTotal);
-
+ Proj_func3D(P1, P2, P3, methodTV, DimTotal);
+
/*updating R and t*/
tkp1 = (1.0f + sqrt(1.0f + 4.0f*tk*tk))*0.5f;
Rupd_dfunc3D(P1, P1_prev, P2, P2_prev, P3, P3_prev, R1, R2, R3, tkp1, tk, DimTotal);
-
+
/*storing old values*/
copyIm(P1, P1_prev, (long)(dimX), (long)(dimY), (long)(dimZ));
copyIm(P2, P2_prev, (long)(dimX), (long)(dimY), (long)(dimZ));
copyIm(P3, P3_prev, (long)(dimX), (long)(dimY), (long)(dimZ));
tk = tkp1;
-
+
/* check early stopping criteria */
if ((epsil != 0.0f) && (ll % 5 == 0)) {
re = 0.0f; re1 = 0.0f;
@@ -168,16 +168,16 @@ float dTV_FGP_CPU_main(float *Input, float *InputRef, float *Output, float *info
if (count > 3) break;
}
}
-
+
free(P3); free(P3_prev); free(R3); free(InputRef_z);
}
if (epsil != 0.0f) free(Output_prev);
free(P1); free(P2); free(P1_prev); free(P2_prev); free(R1); free(R2); free(InputRef_x); free(InputRef_y);
-
+
/*adding info into info_vector */
infovector[0] = (float)(ll); /*iterations number (if stopped earlier based on tolerance)*/
infovector[1] = re; /* reached tolerance */
-
+
return 0;
}
@@ -185,7 +185,6 @@ float dTV_FGP_CPU_main(float *Input, float *InputRef, float *Output, float *info
/********************************************************************/
/***************************2D Functions*****************************/
/********************************************************************/
-
float GradNorm_func2D(float *B, float *B_x, float *B_y, float eta, long dimX, long dimY)
{
long i,j,index;
@@ -249,47 +248,17 @@ float Grad_dfunc2D(float *P1, float *P2, float *D, float *R1, float *R2, float *
/* boundary conditions */
if (i == dimX-1) val1 = 0.0f; else val1 = D[index] - D[j*dimX + (i+1)];
if (j == dimY-1) val2 = 0.0f; else val2 = D[index] - D[(j+1)*dimX + i];
-
+
in_prod = val1*B_x[index] + val2*B_y[index]; /* calculate inner product */
val1 = val1 - in_prod*B_x[index];
val2 = val2 - in_prod*B_y[index];
-
+
P1[index] = R1[index] + multip*val1;
P2[index] = R2[index] + multip*val2;
-
+
}}
return 1;
}
-float Proj_dfunc2D(float *P1, float *P2, int methTV, long DimTotal)
-{
- float val1, val2, denom, sq_denom;
- long i;
- if (methTV == 0) {
- /* isotropic TV*/
-#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/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(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[i] = P1[i]/val1;
- P2[i] = P2[i]/val2;
- }
- }
- return 1;
-}
float Rupd_dfunc2D(float *P1, float *P1_old, float *P2, float *P2_old, float *R1, float *R2, float tkp1, float tk, long DimTotal)
{
long i;
@@ -314,14 +283,14 @@ float GradNorm_func3D(float *B, float *B_x, float *B_y, float *B_z, float eta, l
for(k=0; k<dimZ; k++) {
for(j=0; j<dimY; j++) {
for(i=0; i<dimX; i++) {
-
+
index = (dimX*dimY)*k + j*dimX+i;
-
+
/* zero boundary conditions */
if (i == dimX-1) {val1 = 0.0f;} else {val1 = B[(dimX*dimY)*k + j*dimX+(i+1)];}
if (j == dimY-1) {val2 = 0.0f;} else {val2 = B[(dimX*dimY)*k + (j+1)*dimX+i];}
if (k == dimZ-1) {val3 = 0.0f;} else {val3 = B[(dimX*dimY)*(k+1) + (j)*dimX+i];}
-
+
gradX = val1 - B[index];
gradY = val2 - B[index];
gradZ = val3 - B[index];
@@ -382,52 +351,18 @@ float Grad_dfunc3D(float *P1, float *P2, float *P3, float *D, float *R1, float *
if (i == dimX-1) val1 = 0.0f; else val1 = D[index] - D[(dimX*dimY)*k + j*dimX + (i+1)];
if (j == dimY-1) val2 = 0.0f; else val2 = D[index] - D[(dimX*dimY)*k + (j+1)*dimX + i];
if (k == dimZ-1) val3 = 0.0f; else val3 = D[index] - D[(dimX*dimY)*(k+1) + j*dimX + i];
-
+
in_prod = val1*B_x[index] + val2*B_y[index] + val3*B_z[index]; /* calculate inner product */
val1 = val1 - in_prod*B_x[index];
val2 = val2 - in_prod*B_y[index];
val3 = val3 - in_prod*B_z[index];
-
+
P1[index] = R1[index] + multip*val1;
P2[index] = R2[index] + multip*val2;
P3[index] = R3[index] + multip*val3;
}}}
return 1;
}
-float Proj_dfunc3D(float *P1, float *P2, float *P3, int methTV, long DimTotal)
-{
- float val1, val2, val3, denom, sq_denom;
- long i;
- if (methTV == 0) {
- /* isotropic TV*/
-#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/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(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[i] = P1[i]/val1;
- P2[i] = P2[i]/val2;
- P3[i] = P3[i]/val3;
- }
- }
- return 1;
-}
float Rupd_dfunc3D(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, long DimTotal)
{
long i;
diff --git a/src/Core/regularisers_CPU/FGP_dTV_core.h b/src/Core/regularisers_CPU/FGP_dTV_core.h
index 9ace06d..8582170 100644
--- a/src/Core/regularisers_CPU/FGP_dTV_core.h
+++ b/src/Core/regularisers_CPU/FGP_dTV_core.h
@@ -59,14 +59,12 @@ CCPI_EXPORT float GradNorm_func2D(float *B, float *B_x, float *B_y, float eta, l
CCPI_EXPORT float ProjectVect_func2D(float *R1, float *R2, float *B_x, float *B_y, long dimX, long dimY);
CCPI_EXPORT float Obj_dfunc2D(float *A, float *D, float *R1, float *R2, float lambda, long dimX, long dimY);
CCPI_EXPORT float Grad_dfunc2D(float *P1, float *P2, float *D, float *R1, float *R2, float *B_x, float *B_y, float lambda, long dimX, long dimY);
-CCPI_EXPORT float Proj_dfunc2D(float *P1, float *P2, int methTV, long DimTotal);
CCPI_EXPORT float Rupd_dfunc2D(float *P1, float *P1_old, float *P2, float *P2_old, float *R1, float *R2, float tkp1, float tk, long DimTotal);
CCPI_EXPORT float GradNorm_func3D(float *B, float *B_x, float *B_y, float *B_z, float eta, long dimX, long dimY, long dimZ);
CCPI_EXPORT float ProjectVect_func3D(float *R1, float *R2, float *R3, float *B_x, float *B_y, float *B_z, long dimX, long dimY, long dimZ);
CCPI_EXPORT float Obj_dfunc3D(float *A, float *D, float *R1, float *R2, float *R3, float lambda, long dimX, long dimY, long dimZ);
CCPI_EXPORT float Grad_dfunc3D(float *P1, float *P2, float *P3, float *D, float *R1, float *R2, float *R3, float *B_x, float *B_y, float *B_z, float lambda, long dimX, long dimY, long dimZ);
-CCPI_EXPORT float Proj_dfunc3D(float *P1, float *P2, float *P3, int methTV, long DimTotal);
CCPI_EXPORT float Rupd_dfunc3D(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, long DimTotal);
#ifdef __cplusplus
}
diff --git a/src/Core/regularisers_CPU/PD_TV_core.c b/src/Core/regularisers_CPU/PD_TV_core.c
index a8c3cfb..f476a8b 100644
--- a/src/Core/regularisers_CPU/PD_TV_core.c
+++ b/src/Core/regularisers_CPU/PD_TV_core.c
@@ -50,13 +50,13 @@ float PDTV_CPU_main(float *Input, float *U, float *infovector, float lambdaPar,
theta = 1.0f;
lt = tau/lambdaPar;
ll = 0;
+ DimTotal = (long)(dimX*dimY*dimZ);
copyIm(Input, U, (long)(dimX), (long)(dimY), (long)(dimZ));
if (dimZ <= 1) {
/*2D case */
float *U_old=NULL, *P1=NULL, *P2=NULL;
- DimTotal = (long)(dimX*dimY);
U_old = calloc(DimTotal, sizeof(float));
P1 = calloc(DimTotal, sizeof(float));
@@ -65,8 +65,7 @@ float PDTV_CPU_main(float *Input, float *U, float *infovector, float lambdaPar,
/* begin iterations */
for(ll=0; ll<iterationsNumb; ll++) {
- //if ((epsil != 0.0f) && (ll % 5 == 0)) copyIm(Output, Output_prev, (long)(dimX), (long)(dimY), 1l);
- /* computing the gradient of the objective function */
+ /* computing the the dual P variable */
DualP2D(U, P1, P2, (long)(dimX), (long)(dimY), sigma);
/* apply nonnegativity */
@@ -94,12 +93,52 @@ float PDTV_CPU_main(float *Input, float *U, float *infovector, float lambdaPar,
}
/*get updated solution*/
- getX2D(U, U_old, theta, DimTotal);
+ getX(U, U_old, theta, DimTotal);
}
free(P1); free(P2); free(U_old);
}
else {
- /*3D case*/
+ /*3D case*/
+ float *U_old=NULL, *P1=NULL, *P2=NULL, *P3=NULL;
+ U_old = calloc(DimTotal, sizeof(float));
+ P1 = calloc(DimTotal, sizeof(float));
+ P2 = calloc(DimTotal, sizeof(float));
+ P3 = calloc(DimTotal, sizeof(float));
+
+ /* begin iterations */
+ for(ll=0; ll<iterationsNumb; ll++) {
+
+ /* computing the the dual P variable */
+ DualP3D(U, P1, P2, P3, (long)(dimX), (long)(dimY), (long)(dimZ), sigma);
+
+ /* apply nonnegativity */
+ if (nonneg == 1) for(j=0; j<DimTotal; j++) {if (U[j] < 0.0f) U[j] = 0.0f;}
+
+ /* projection step */
+ Proj_func3D(P1, P2, P3, methodTV, DimTotal);
+
+ /* copy U to U_old */
+ copyIm(U, U_old, (long)(dimX), (long)(dimY), (long)(dimZ));
+
+ DivProj3D(U, Input, P1, P2, P3, (long)(dimX), (long)(dimY), (long)(dimZ), lt, tau);
+
+ /* check early stopping criteria */
+ if ((epsil != 0.0f) && (ll % 5 == 0)) {
+ re = 0.0f; re1 = 0.0f;
+ for(j=0; j<DimTotal; j++)
+ {
+ re += powf(U[j] - U_old[j],2);
+ re1 += powf(U[j],2);
+ }
+ re = sqrtf(re)/sqrtf(re1);
+ if (re < epsil) count++;
+ if (count > 3) break;
+ }
+ /*get updated solution*/
+
+ getX(U, U_old, theta, DimTotal);
+ }
+ free(P1); free(P2); free(P3); free(U_old);
}
/*adding info into info_vector */
infovector[0] = (float)(ll); /*iterations number (if stopped earlier based on tolerance)*/
@@ -134,7 +173,7 @@ float DivProj2D(float *U, float *Input, float *P1, float *P2, long dimX, long di
{
long i,j,index;
float P_v1, P_v2, div_var;
- #pragma omp parallel for shared(U,Input,P1,P2) private(i, j, P_v1, P_v2, div_var)
+ #pragma omp parallel for shared(U,Input,P1,P2) private(index, i, j, P_v1, P_v2, div_var)
for(j=0; j<dimY; j++) {
for(i=0; i<dimX; i++) {
index = j*dimX+i;
@@ -150,7 +189,7 @@ float DivProj2D(float *U, float *Input, float *P1, float *P2, long dimX, long di
}
/*get the updated solution*/
-float getX2D(float *U, float *U_old, float theta, long DimTotal)
+float getX(float *U, float *U_old, float theta, long DimTotal)
{
long i;
#pragma omp parallel for shared(U,U_old) private(i)
@@ -164,3 +203,45 @@ float getX2D(float *U, float *U_old, float theta, long DimTotal)
/*****************************************************************/
/************************3D-case related Functions */
/*****************************************************************/
+/*Calculating dual variable (using forward differences)*/
+float DualP3D(float *U, float *P1, float *P2, float *P3, long dimX, long dimY, long dimZ, float sigma)
+{
+ long i,j,k,index;
+ #pragma omp parallel for shared(U,P1,P2,P3) private(index,i,j,k)
+ for(k=0; k<dimZ; k++) {
+ for(j=0; j<dimY; j++) {
+ for(i=0; i<dimX; i++) {
+ index = (dimX*dimY)*k + j*dimX+i;
+ /* symmetric boundary conditions (Neuman) */
+ if (i == dimX-1) P1[index] += sigma*(U[(dimX*dimY)*k + j*dimX+(i-1)] - U[index]);
+ else P1[index] += sigma*(U[(dimX*dimY)*k + j*dimX+(i+1)] - U[index]);
+ if (j == dimY-1) P2[index] += sigma*(U[(dimX*dimY)*k + (j-1)*dimX+i] - U[index]);
+ else P2[index] += sigma*(U[(dimX*dimY)*k + (j+1)*dimX+i] - U[index]);
+ if (k == dimZ-1) P3[index] += sigma*(U[(dimX*dimY)*(k-1) + j*dimX+i] - U[index]);
+ else P3[index] += sigma*(U[(dimX*dimY)*(k+1) + j*dimX+i] - U[index]);
+ }}}
+ return 1;
+}
+
+/* Divergence for P dual */
+float DivProj3D(float *U, float *Input, float *P1, float *P2, float *P3, long dimX, long dimY, long dimZ, float lt, float tau)
+{
+ long i,j,k,index;
+ float P_v1, P_v2, P_v3, div_var;
+ #pragma omp parallel for shared(U,Input,P1,P2) private(index, i, j, k, P_v1, P_v2, P_v3, div_var)
+ for(k=0; k<dimZ; k++) {
+ for(j=0; j<dimY; j++) {
+ for(i=0; i<dimX; i++) {
+ index = (dimX*dimY)*k + j*dimX+i;
+ /* symmetric boundary conditions (Neuman) */
+ if (i == 0) P_v1 = -P1[index];
+ else P_v1 = -(P1[index] - P1[(dimX*dimY)*k + j*dimX+(i-1)]);
+ if (j == 0) P_v2 = -P2[index];
+ else P_v2 = -(P2[index] - P2[(dimX*dimY)*k + (j-1)*dimX+i]);
+ if (k == 0) P_v3 = -P3[index];
+ else P_v3 = -(P3[index] - P3[(dimX*dimY)*(k-1) + j*dimX+i]);
+ div_var = P_v1 + P_v2 + P_v3;
+ U[index] = (U[index] - tau*div_var + lt*Input[index])/(1.0 + lt);
+ }}}
+ return *U;
+}
diff --git a/src/Core/regularisers_CPU/PD_TV_core.h b/src/Core/regularisers_CPU/PD_TV_core.h
index b52689a..b4e8a75 100644
--- a/src/Core/regularisers_CPU/PD_TV_core.h
+++ b/src/Core/regularisers_CPU/PD_TV_core.h
@@ -51,7 +51,10 @@ CCPI_EXPORT float PDTV_CPU_main(float *Input, float *U, float *infovector, float
CCPI_EXPORT float DualP2D(float *U, float *P1, float *P2, long dimX, long dimY, float sigma);
CCPI_EXPORT float DivProj2D(float *U, float *Input, float *P1, float *P2, long dimX, long dimY, float lt, float tau);
-CCPI_EXPORT float getX2D(float *U, float *U_old, float theta, long DimTotal);
+CCPI_EXPORT float getX(float *U, float *U_old, float theta, long DimTotal);
+
+CCPI_EXPORT float DualP3D(float *U, float *P1, float *P2, float *P3, long dimX, long dimY, long dimZ, float sigma);
+CCPI_EXPORT float DivProj3D(float *U, float *Input, float *P1, float *P2, float *P3, long dimX, long dimY, long dimZ, float lt, float tau);
#ifdef __cplusplus
}
#endif
diff --git a/src/Core/regularisers_CPU/utils.c b/src/Core/regularisers_CPU/utils.c
index cf4ad72..088ace9 100644
--- a/src/Core/regularisers_CPU/utils.c
+++ b/src/Core/regularisers_CPU/utils.c
@@ -145,7 +145,7 @@ float Im_scale2D(float *Input, float *Scaled, int w, int h, int w2, int h2)
return *Scaled;
}
-/*2D Projection onto convex set for P*/
+/*2D Projection onto convex set for P (called in PD_TV, FGP_dTV and FGP_TV methods)*/
float Proj_func2D(float *P1, float *P2, int methTV, long DimTotal)
{
float val1, val2, denom, sq_denom;
@@ -176,3 +176,38 @@ float Proj_func2D(float *P1, float *P2, int methTV, long DimTotal)
}
return 1;
}
+/*3D Projection onto convex set for P (called in PD_TV, FGP_TV, FGP_dTV methods)*/
+float Proj_func3D(float *P1, float *P2, float *P3, int methTV, long DimTotal)
+{
+ float val1, val2, val3, denom, sq_denom;
+ long i;
+ if (methTV == 0) {
+ /* isotropic TV*/
+#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/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(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[i] = P1[i]/val1;
+ P2[i] = P2[i]/val2;
+ P3[i] = P3[i]/val3;
+ }
+ }
+ return 1;
+}
diff --git a/src/Core/regularisers_CPU/utils.h b/src/Core/regularisers_CPU/utils.h
index e050f86..4b57f71 100644
--- a/src/Core/regularisers_CPU/utils.h
+++ b/src/Core/regularisers_CPU/utils.h
@@ -32,6 +32,7 @@ CCPI_EXPORT float TV_energy3D(float *U, float *U0, float *E_val, float lambda, i
CCPI_EXPORT float TV_energy3D(float *U, float *U0, float *E_val, float lambda, int type, int dimX, int dimY, int dimZ);
CCPI_EXPORT float Im_scale2D(float *Input, float *Scaled, int w, int h, int w2, int h2);
CCPI_EXPORT float Proj_func2D(float *P1, float *P2, int methTV, long DimTotal);
+CCPI_EXPORT float Proj_func3D(float *P1, float *P2, float *P3, int methTV, long DimTotal);
#ifdef __cplusplus
}
#endif
diff --git a/src/Python/setup-regularisers.py.in b/src/Python/setup-regularisers.py.in
index 9bcd46d..9a5b693 100644
--- a/src/Python/setup-regularisers.py.in
+++ b/src/Python/setup-regularisers.py.in
@@ -38,7 +38,6 @@ extra_include_dirs += [os.path.join(".." , "Core"),
os.path.join(".." , "Core", "regularisers_CPU"),
os.path.join(".." , "Core", "inpainters_CPU"),
os.path.join(".." , "Core", "regularisers_GPU" , "TV_FGP" ) ,
- os.path.join(".." , "Core", "regularisers_GPU" , "TV_PD" ) ,
os.path.join(".." , "Core", "regularisers_GPU" , "TV_ROF" ) ,
os.path.join(".." , "Core", "regularisers_GPU" , "TV_SB" ) ,
os.path.join(".." , "Core", "regularisers_GPU" , "TGV" ) ,
diff --git a/src/Python/src/cpu_regularisers.pyx b/src/Python/src/cpu_regularisers.pyx
index 724634b..08e247c 100644
--- a/src/Python/src/cpu_regularisers.pyx
+++ b/src/Python/src/cpu_regularisers.pyx
@@ -163,7 +163,7 @@ def TV_PD_CPU(inputData, regularisation_parameter, iterationsNumb, tolerance_par
if inputData.ndim == 2:
return TV_PD_2D(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, lipschitz_const)
elif inputData.ndim == 3:
- return 0
+ return TV_PD_3D(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, lipschitz_const)
def TV_PD_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
float regularisation_parameter,
@@ -191,7 +191,34 @@ def TV_PD_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
methodTV,
nonneg,
dims[1],dims[0], 1)
+ return (outputData,infovec)
+
+def TV_PD_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
+ float regularisation_parameter,
+ int iterationsNumb,
+ float tolerance_param,
+ int methodTV,
+ int nonneg,
+ float lipschitz_const):
+
+ 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')
+ cdef np.ndarray[np.float32_t, ndim=1, mode="c"] infovec = \
+ np.zeros([2], dtype='float32')
+
+ #/* Run FGP-TV iterations for 3D data */
+ PDTV_CPU_main(&inputData[0,0,0], &outputData[0,0,0], &infovec[0], regularisation_parameter,
+ iterationsNumb,
+ tolerance_param,
+ lipschitz_const,
+ methodTV,
+ nonneg,
+ dims[2], dims[1], dims[0])
return (outputData,infovec)
#***************************************************************#
generated by cgit v1.2.3 (git 2.25.1) at 2025-11-04 06:21:35 +0000