summaryrefslogtreecommitdiffstats
path: root/src/Core/regularisers_CPU/FGP_TV_core.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/Core/regularisers_CPU/FGP_TV_core.c')
-rw-r--r--src/Core/regularisers_CPU/FGP_TV_core.c258
1 files changed, 129 insertions, 129 deletions
diff --git a/src/Core/regularisers_CPU/FGP_TV_core.c b/src/Core/regularisers_CPU/FGP_TV_core.c
index 69c92bc..a16a2e5 100644
--- a/src/Core/regularisers_CPU/FGP_TV_core.c
+++ b/src/Core/regularisers_CPU/FGP_TV_core.c
@@ -1,21 +1,21 @@
/*
-This work is part of the Core Imaging Library developed by
-Visual Analytics and Imaging System Group of the Science Technology
-Facilities Council, STFC
-
-Copyright 2019 Daniil Kazantsev
-Copyright 2019 Srikanth Nagella, Edoardo Pasca
-
-Licensed under the Apache License, Version 2.0 (the "License");
-you may not use this file except in compliance with the License.
-You may obtain a copy of the License at
-http://www.apache.org/licenses/LICENSE-2.0
-Unless required by applicable law or agreed to in writing, software
-distributed under the License is distributed on an "AS IS" BASIS,
-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.
-*/
+ * This work is part of the Core Imaging Library developed by
+ * Visual Analytics and Imaging System Group of the Science Technology
+ * Facilities Council, STFC
+ *
+ * Copyright 2019 Daniil Kazantsev
+ * Copyright 2019 Srikanth Nagella, Edoardo Pasca
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ * http://www.apache.org/licenses/LICENSE-2.0
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * 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.
+ */
#include "FGP_TV_core.h"
@@ -46,66 +46,66 @@ float TV_FGP_CPU_main(float *Input, float *Output, float *infovector, float lamb
float tk = 1.0f;
float tkp1 =1.0f;
int count = 0;
-
- if (dimZ <= 1) {
- /*2D case */
- float *Output_prev=NULL, *P1=NULL, *P2=NULL, *P1_prev=NULL, *P2_prev=NULL, *R1=NULL, *R2=NULL;
- DimTotal = (long)(dimX*dimY);
-
- if (epsil != 0.0f) Output_prev = calloc(DimTotal, sizeof(float));
+
+ if (dimZ <= 1) {
+ /*2D case */
+ float *Output_prev=NULL, *P1=NULL, *P2=NULL, *P1_prev=NULL, *P2_prev=NULL, *R1=NULL, *R2=NULL;
+ DimTotal = (long)(dimX*dimY);
+
+ if (epsil != 0.0f) Output_prev = calloc(DimTotal, sizeof(float));
P1 = calloc(DimTotal, sizeof(float));
P2 = calloc(DimTotal, sizeof(float));
P1_prev = calloc(DimTotal, sizeof(float));
P2_prev = calloc(DimTotal, sizeof(float));
R1 = calloc(DimTotal, sizeof(float));
R2 = calloc(DimTotal, sizeof(float));
-
- /* begin iterations */
+
+ /* 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 */
Obj_func2D(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_func2D(P1, P2, Output, R1, R2, lambdaPar, (long)(dimX), (long)(dimY));
-
+
/* projection step */
Proj_func2D(P1, P2, methodTV, DimTotal);
-
+
/*updating R and t*/
tkp1 = (1.0f + sqrtf(1.0f + 4.0f*tk*tk))*0.5f;
Rupd_func2D(P1, P1_prev, P2, P2_prev, R1, R2, tkp1, tk, DimTotal);
-
+
/*storing old values*/
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;
- for(j=0; j<DimTotal; j++)
- {
- re += powf(Output[j] - Output_prev[j],2);
- re1 += powf(Output[j],2);
- }
- re = sqrtf(re)/sqrtf(re1);
- if (re < epsil) count++;
- if (count > 3) break;
+ re = 0.0f; re1 = 0.0f;
+ for(j=0; j<DimTotal; j++)
+ {
+ re += powf(Output[j] - Output_prev[j],2);
+ re1 += powf(Output[j],2);
+ }
+ re = sqrtf(re)/sqrtf(re1);
+ if (re < epsil) count++;
+ if (count > 3) break;
}
- }
- if (epsil != 0.0f) free(Output_prev);
- free(P1); free(P2); free(P1_prev); free(P2_prev); free(R1); free(R2);
- }
- else {
- /*3D case*/
- float *Output_prev=NULL, *P1=NULL, *P2=NULL, *P3=NULL, *P1_prev=NULL, *P2_prev=NULL, *P3_prev=NULL, *R1=NULL, *R2=NULL, *R3=NULL;
+ }
+ if (epsil != 0.0f) free(Output_prev);
+ free(P1); free(P2); free(P1_prev); free(P2_prev); free(R1); free(R2);
+ }
+ else {
+ /*3D case*/
+ float *Output_prev=NULL, *P1=NULL, *P2=NULL, *P3=NULL, *P1_prev=NULL, *P2_prev=NULL, *P3_prev=NULL, *R1=NULL, *R2=NULL, *R3=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));
@@ -116,58 +116,58 @@ float TV_FGP_CPU_main(float *Input, float *Output, float *infovector, float lamb
R1 = calloc(DimTotal, sizeof(float));
R2 = calloc(DimTotal, sizeof(float));
R3 = calloc(DimTotal, sizeof(float));
-
- /* begin iterations */
+
+ /* 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));
-
+
/* computing the gradient of the objective function */
Obj_func3D(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_func3D(P1, P2, P3, Output, R1, R2, R3, lambdaPar, (long)(dimX), (long)(dimY), (long)(dimZ));
-
+
/* projection step */
Proj_func3D(P1, P2, P3, methodTV, DimTotal);
-
+
/*updating R and t*/
tkp1 = (1.0f + sqrtf(1.0f + 4.0f*tk*tk))*0.5f;
Rupd_func3D(P1, P1_prev, P2, P2_prev, P3, P3_prev, R1, R2, R3, tkp1, tk, DimTotal);
-
+
/* calculate norm - stopping rules*/
if ((epsil != 0.0f) && (ll % 5 == 0)) {
- re = 0.0f; re1 = 0.0f;
- for(j=0; j<DimTotal; j++)
- {
- re += powf(Output[j] - Output_prev[j],2);
- re1 += powf(Output[j],2);
- }
- re = sqrtf(re)/sqrtf(re1);
- /* stop if the norm residual is less than the tolerance EPS */
- if (re < epsil) count++;
- if (count > 3) break;
+ re = 0.0f; re1 = 0.0f;
+ for(j=0; j<DimTotal; j++)
+ {
+ re += powf(Output[j] - Output_prev[j],2);
+ re1 += powf(Output[j],2);
+ }
+ re = sqrtf(re)/sqrtf(re1);
+ /* stop if the norm residual is less than the tolerance EPS */
+ if (re < epsil) count++;
+ if (count > 3) break;
}
-
+
/*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;
}
-
- if (epsil != 0.0f) free(Output_prev);
- free(P1); free(P2); free(P3); free(P1_prev); free(P2_prev); free(P3_prev); free(R1); free(R2); free(R3);
- }
-
- /*adding info into info_vector */
- infovector[0] = (float)(ll); /*iterations number (if stopped earlier based on tolerance)*/
- infovector[1] = re; /* reached tolerance */
-
- return 0;
+
+ if (epsil != 0.0f) free(Output_prev);
+ free(P1); free(P2); free(P3); free(P1_prev); free(P2_prev); free(P3_prev); free(R1); free(R2); free(R3);
+ }
+
+ /*adding info into info_vector */
+ infovector[0] = (float)(ll); /*iterations number (if stopped earlier based on tolerance)*/
+ infovector[1] = re; /* reached tolerance */
+
+ return 0;
}
float Obj_func2D(float *A, float *D, float *R1, float *R2, float lambda, long dimX, long dimY)
@@ -177,7 +177,7 @@ float Obj_func2D(float *A, float *D, float *R1, float *R2, float lambda, long di
#pragma omp parallel for shared(A,D,R1,R2) private(index,i,j,val1,val2)
for(j=0; j<dimY; j++) {
for(i=0; i<dimX; i++) {
- index = j*dimX+i;
+ index = j*dimX+i;
/* boundary conditions */
if (i == 0) {val1 = 0.0f;} else {val1 = R1[j*dimX + (i-1)];}
if (j == 0) {val2 = 0.0f;} else {val2 = R2[(j-1)*dimX + i];}
@@ -193,7 +193,7 @@ float Grad_func2D(float *P1, float *P2, float *D, float *R1, float *R2, float la
#pragma omp parallel for shared(P1,P2,D,R1,R2,multip) private(index,i,j,val1,val2)
for(j=0; j<dimY; j++) {
for(i=0; i<dimX; i++) {
- index = j*dimX+i;
+ index = j*dimX+i;
/* 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];
@@ -210,25 +210,25 @@ float Proj_func2D(float *P1, float *P2, int methTV, long DimTotal)
/* 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;
- }
+ 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;
- }
+ 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;
}
@@ -239,9 +239,9 @@ float Rupd_func2D(float *P1, float *P1_old, float *P2, float *P2_old, float *R1,
multip = ((tk-1.0f)/tkp1);
#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]);
- }
+ R1[i] = P1[i] + multip*(P1[i] - P1_old[i]);
+ R2[i] = P2[i] + multip*(P2[i] - P2_old[i]);
+ }
return 1;
}
@@ -255,7 +255,7 @@ float Obj_func3D(float *A, float *D, float *R1, float *R2, float *R3, float lamb
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;
+ 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];}
@@ -273,7 +273,7 @@ float Grad_func3D(float *P1, float *P2, float *P3, float *D, float *R1, float *R
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;
+ index = (dimX*dimY)*k + j*dimX+i;
/* boundary conditions */
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];
@@ -289,33 +289,33 @@ 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;
- }
- }
- }
+ /* 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*/
+ /* 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;
- }
- }
+ 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)
@@ -325,9 +325,9 @@ float Rupd_func3D(float *P1, float *P1_old, float *P2, float *P2_old, float *P3,
multip = ((tk-1.0f)/tkp1);
#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]);
- }
+ 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;
}
generated by cgit v1.2.1 (git 2.18.0) at 2024-05-23 21:01:38 +0000