From 4c6769401570415a5a543b81130e314af6b95d9a Mon Sep 17 00:00:00 2001 From: Daniil Kazantsev Date: Thu, 26 Sep 2019 23:07:17 +0100 Subject: loop CPU routines upgrades --- src/Core/regularisers_CPU/Diffus4th_order_core.c | 283 +++++++------- src/Core/regularisers_CPU/Diffusion_core.c | 361 +++++++++--------- src/Core/regularisers_CPU/FGP_TV_core.c | 258 ++++++------- src/Core/regularisers_CPU/FGP_dTV_core.c | 363 +++++++++--------- src/Core/regularisers_CPU/LLT_ROF_core.c | 448 +++++++++++------------ src/Core/regularisers_CPU/Nonlocal_TV_core.c | 167 +++++---- src/Core/regularisers_CPU/PatchSelect_core.c | 46 +-- src/Core/regularisers_CPU/ROF_TV_core.c | 124 +++---- src/Core/regularisers_CPU/SB_TV_core.c | 301 +++++++-------- src/Core/regularisers_CPU/TGV_core.c | 186 +++++----- 10 files changed, 1266 insertions(+), 1271 deletions(-) (limited to 'src') diff --git a/src/Core/regularisers_CPU/Diffus4th_order_core.c b/src/Core/regularisers_CPU/Diffus4th_order_core.c index 28ac8a9..840b736 100644 --- a/src/Core/regularisers_CPU/Diffus4th_order_core.c +++ b/src/Core/regularisers_CPU/Diffus4th_order_core.c @@ -50,35 +50,35 @@ float Diffus4th_CPU_main(float *Input, float *Output, float *infovector, float l float *W_Lapl=NULL, *Output_prev=NULL; sigmaPar2 = sigmaPar*sigmaPar; DimTotal = dimX*dimY*dimZ; - + W_Lapl = calloc(DimTotal, sizeof(float)); - + if (epsil != 0.0f) Output_prev = calloc(DimTotal, sizeof(float)); - + /* copy into output */ copyIm(Input, Output, (long)(dimX), (long)(dimY), (long)(dimZ)); - + for(i=0; i < iterationsNumb; i++) { - if ((epsil != 0.0f) && (i % 5 == 0)) copyIm(Output, Output_prev, (long)(dimX), (long)(dimY), (long)(dimZ)); - - if (dimZ == 1) { + if ((epsil != 0.0f) && (i % 5 == 0)) copyIm(Output, Output_prev, (long)(dimX), (long)(dimY), (long)(dimZ)); + + if (dimZ == 1) { /* running 2D diffusion iterations */ /* Calculating weighted Laplacian */ Weighted_Laplc2D(W_Lapl, Output, sigmaPar2, dimX, dimY); /* Perform iteration step */ Diffusion_update_step2D(Output, Input, W_Lapl, lambdaPar, sigmaPar2, tau, (long)(dimX), (long)(dimY)); - } - else { + } + else { /* running 3D diffusion iterations */ /* Calculating weighted Laplacian */ Weighted_Laplc3D(W_Lapl, Output, sigmaPar2, dimX, dimY, dimZ); /* Perform iteration step */ Diffusion_update_step3D(Output, Input, W_Lapl, lambdaPar, sigmaPar2, tau, (long)(dimX), (long)(dimY), (long)(dimZ)); - } - - /* check early stopping criteria */ - if ((epsil != 0.0f) && (i % 5 == 0)) { - re = 0.0f; re1 = 0.0f; + } + + /* check early stopping criteria */ + if ((epsil != 0.0f) && (i % 5 == 0)) { + re = 0.0f; re1 = 0.0f; for(j=0; j 3) break; - } - } - free(W_Lapl); - + } + } + free(W_Lapl); + if (epsil != 0.0f) free(Output_prev); /*adding info into info_vector */ infovector[0] = (float)(i); /*iterations number (if stopped earlier based on tolerance)*/ @@ -104,74 +104,71 @@ float Weighted_Laplc2D(float *W_Lapl, float *U0, float sigma, long dimX, long di { long i,j,i1,i2,j1,j2,index; float gradX, gradX_sq, gradY, gradY_sq, gradXX, gradYY, gradXY, xy_2, denom, V_norm, V_orth, c, c_sq; - - #pragma omp parallel for shared(W_Lapl) private(i,j,i1,i2,j1,j2,index,gradX, gradX_sq, gradY, gradY_sq, gradXX, gradYY, gradXY, xy_2, denom, V_norm, V_orth, c, c_sq) + +#pragma omp parallel for shared(W_Lapl) private(i,j,i1,i2,j1,j2,index,gradX, gradX_sq, gradY, gradY_sq, gradXX, gradYY, gradXY, xy_2, denom, V_norm, V_orth, c, c_sq) + for(j=0; j 3) break; - } - } - + re = sqrtf(re)/sqrtf(re1); + /* stop if the norm residual is less than the tolerance EPS */ + if (re < epsil) count++; + if (count > 3) break; + } + } + free(Output_prev); - /*adding info into info_vector */ + /*adding info into info_vector */ infovector[0] = (float)(i); /*iterations number (if stopped earlier based on tolerance)*/ infovector[1] = re; /* reached tolerance */ return 0; } - /********************************************************************/ /***************************2D Functions*****************************/ /********************************************************************/ /* linear diffusion (heat equation) */ float LinearDiff2D(float *Input, float *Output, float lambdaPar, float tau, long dimX, long dimY) { - long i,j,i1,i2,j1,j2,index; - float e,w,n,s,e1,w1,n1,s1; - + long i,j,i1,i2,j1,j2,index; + float e,w,n,s,e1,w1,n1,s1; + #pragma omp parallel for shared(Input) private(index,i,j,i1,i2,j1,j2,e,w,n,s,e1,w1,n1,s1) - for(i=0; i sigmaPar) e1 = signNDFc(e1); - else e1 = e1/sigmaPar; - - if (fabs(w1) > sigmaPar) w1 = signNDFc(w1); - else w1 = w1/sigmaPar; - - if (fabs(n1) > sigmaPar) n1 = signNDFc(n1); - else n1 = n1/sigmaPar; - - if (fabs(s1) > sigmaPar) s1 = signNDFc(s1); - else s1 = s1/sigmaPar; + /* Huber penalty */ + if (fabs(e1) > sigmaPar) e1 = signNDFc(e1); + else e1 = e1/sigmaPar; + + if (fabs(w1) > sigmaPar) w1 = signNDFc(w1); + else w1 = w1/sigmaPar; + + if (fabs(n1) > sigmaPar) n1 = signNDFc(n1); + else n1 = n1/sigmaPar; + + if (fabs(s1) > sigmaPar) s1 = signNDFc(s1); + else s1 = s1/sigmaPar; } else if (penaltytype == 2) { - /* Perona-Malik */ - e1 = (e1)/(1.0f + powf((e1/sigmaPar),2)); - w1 = (w1)/(1.0f + powf((w1/sigmaPar),2)); - n1 = (n1)/(1.0f + powf((n1/sigmaPar),2)); - s1 = (s1)/(1.0f + powf((s1/sigmaPar),2)); + /* Perona-Malik */ + e1 = (e1)/(1.0f + powf((e1/sigmaPar),2)); + w1 = (w1)/(1.0f + powf((w1/sigmaPar),2)); + n1 = (n1)/(1.0f + powf((n1/sigmaPar),2)); + s1 = (s1)/(1.0f + powf((s1/sigmaPar),2)); } else if (penaltytype == 3) { - /* Tukey Biweight */ - if (fabs(e1) <= sigmaPar) e1 = e1*powf((1.0f - powf((e1/sigmaPar),2)), 2); - else e1 = 0.0f; - if (fabs(w1) <= sigmaPar) w1 = w1*powf((1.0f - powf((w1/sigmaPar),2)), 2); - else w1 = 0.0f; - if (fabs(n1) <= sigmaPar) n1 = n1*powf((1.0f - powf((n1/sigmaPar),2)), 2); - else n1 = 0.0f; - if (fabs(s1) <= sigmaPar) s1 = s1*powf((1.0f - powf((s1/sigmaPar),2)), 2); - else s1 = 0.0f; + /* Tukey Biweight */ + if (fabs(e1) <= sigmaPar) e1 = e1*powf((1.0f - powf((e1/sigmaPar),2)), 2); + else e1 = 0.0f; + if (fabs(w1) <= sigmaPar) w1 = w1*powf((1.0f - powf((w1/sigmaPar),2)), 2); + else w1 = 0.0f; + if (fabs(n1) <= sigmaPar) n1 = n1*powf((1.0f - powf((n1/sigmaPar),2)), 2); + else n1 = 0.0f; + if (fabs(s1) <= sigmaPar) s1 = s1*powf((1.0f - powf((s1/sigmaPar),2)), 2); + else s1 = 0.0f; } else { - printf("%s \n", "No penalty function selected! Use 1,2 or 3."); - break; - } - Output[index] += tau*(lambdaPar*(e1 + w1 + n1 + s1) - (Output[index] - Input[index])); - }} - return *Output; + printf("%s \n", "No penalty function selected! Use 1,2 or 3."); + break; + } + Output[index] += tau*(lambdaPar*(e1 + w1 + n1 + s1) - (Output[index] - Input[index])); + }} + return *Output; } /********************************************************************/ /***************************3D Functions*****************************/ @@ -210,125 +209,125 @@ float NonLinearDiff2D(float *Input, float *Output, float lambdaPar, float sigmaP /* linear diffusion (heat equation) */ float LinearDiff3D(float *Input, float *Output, float lambdaPar, float tau, long dimX, long dimY, long dimZ) { - long i,j,k,i1,i2,j1,j2,k1,k2,index; - float e,w,n,s,u,d,e1,w1,n1,s1,u1,d1; - + long i,j,k,i1,i2,j1,j2,k1,k2,index; + float e,w,n,s,u,d,e1,w1,n1,s1,u1,d1; + #pragma omp parallel for shared(Input) private(index,i,j,i1,i2,j1,j2,e,w,n,s,e1,w1,n1,s1,k,k1,k2,u1,d1,u,d) -for(k=0; k sigmaPar) e1 = signNDFc(e1); - else e1 = e1/sigmaPar; - - if (fabs(w1) > sigmaPar) w1 = signNDFc(w1); - else w1 = w1/sigmaPar; - - if (fabs(n1) > sigmaPar) n1 = signNDFc(n1); - else n1 = n1/sigmaPar; - - if (fabs(s1) > sigmaPar) s1 = signNDFc(s1); - else s1 = s1/sigmaPar; - - if (fabs(u1) > sigmaPar) u1 = signNDFc(u1); - else u1 = u1/sigmaPar; - - if (fabs(d1) > sigmaPar) d1 = signNDFc(d1); - else d1 = d1/sigmaPar; - } - else if (penaltytype == 2) { - /* Perona-Malik */ - e1 = (e1)/(1.0f + powf((e1/sigmaPar),2)); - w1 = (w1)/(1.0f + powf((w1/sigmaPar),2)); - n1 = (n1)/(1.0f + powf((n1/sigmaPar),2)); - s1 = (s1)/(1.0f + powf((s1/sigmaPar),2)); - u1 = (u1)/(1.0f + powf((u1/sigmaPar),2)); - d1 = (d1)/(1.0f + powf((d1/sigmaPar),2)); - } - else if (penaltytype == 3) { - /* Tukey Biweight */ - if (fabs(e1) <= sigmaPar) e1 = e1*powf((1.0f - powf((e1/sigmaPar),2)), 2); - else e1 = 0.0f; - if (fabs(w1) <= sigmaPar) w1 = w1*powf((1.0f - powf((w1/sigmaPar),2)), 2); - else w1 = 0.0f; - if (fabs(n1) <= sigmaPar) n1 = n1*powf((1.0f - powf((n1/sigmaPar),2)), 2); - else n1 = 0.0f; - if (fabs(s1) <= sigmaPar) s1 = s1*powf((1.0f - powf((s1/sigmaPar),2)), 2); - else s1 = 0.0f; - if (fabs(u1) <= sigmaPar) u1 = u1*powf((1.0f - powf((u1/sigmaPar),2)), 2); - else u1 = 0.0f; - if (fabs(d1) <= sigmaPar) d1 = d1*powf((1.0f - powf((d1/sigmaPar),2)), 2); - else d1 = 0.0f; - } - else { - printf("%s \n", "No penalty function selected! Use 1,2 or 3."); - break; - } - + + if (penaltytype == 1){ + /* Huber penalty */ + if (fabs(e1) > sigmaPar) e1 = signNDFc(e1); + else e1 = e1/sigmaPar; + + if (fabs(w1) > sigmaPar) w1 = signNDFc(w1); + else w1 = w1/sigmaPar; + + if (fabs(n1) > sigmaPar) n1 = signNDFc(n1); + else n1 = n1/sigmaPar; + + if (fabs(s1) > sigmaPar) s1 = signNDFc(s1); + else s1 = s1/sigmaPar; + + if (fabs(u1) > sigmaPar) u1 = signNDFc(u1); + else u1 = u1/sigmaPar; + + if (fabs(d1) > sigmaPar) d1 = signNDFc(d1); + else d1 = d1/sigmaPar; + } + else if (penaltytype == 2) { + /* Perona-Malik */ + e1 = (e1)/(1.0f + powf((e1/sigmaPar),2)); + w1 = (w1)/(1.0f + powf((w1/sigmaPar),2)); + n1 = (n1)/(1.0f + powf((n1/sigmaPar),2)); + s1 = (s1)/(1.0f + powf((s1/sigmaPar),2)); + u1 = (u1)/(1.0f + powf((u1/sigmaPar),2)); + d1 = (d1)/(1.0f + powf((d1/sigmaPar),2)); + } + else if (penaltytype == 3) { + /* Tukey Biweight */ + if (fabs(e1) <= sigmaPar) e1 = e1*powf((1.0f - powf((e1/sigmaPar),2)), 2); + else e1 = 0.0f; + if (fabs(w1) <= sigmaPar) w1 = w1*powf((1.0f - powf((w1/sigmaPar),2)), 2); + else w1 = 0.0f; + if (fabs(n1) <= sigmaPar) n1 = n1*powf((1.0f - powf((n1/sigmaPar),2)), 2); + else n1 = 0.0f; + if (fabs(s1) <= sigmaPar) s1 = s1*powf((1.0f - powf((s1/sigmaPar),2)), 2); + else s1 = 0.0f; + if (fabs(u1) <= sigmaPar) u1 = u1*powf((1.0f - powf((u1/sigmaPar),2)), 2); + else u1 = 0.0f; + if (fabs(d1) <= sigmaPar) d1 = d1*powf((1.0f - powf((d1/sigmaPar),2)), 2); + else d1 = 0.0f; + } + else { + printf("%s \n", "No penalty function selected! Use 1,2 or 3."); + break; + } + Output[index] += tau*(lambdaPar*(e1 + w1 + n1 + s1 + u1 + d1) - (Output[index] - Input[index])); - }}} - return *Output; + }}} + return *Output; } 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 3) break; + re = 0.0f; re1 = 0.0f; + for(j=0; j 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 3) break; + re = 0.0f; re1 = 0.0f; + for(j=0; j 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 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 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 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 3) break; + re = 0.0f; re1 = 0.0f; + for(j=0; j 3) break; } } - } - else { - /*3D case*/ - float *P3=NULL, *P3_prev=NULL, *R3=NULL, *InputRef_z=NULL; - + } + 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 */ + + /* 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 3) break; + re = 0.0f; re1 = 0.0f; + for(j=0; j 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; + + 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; } @@ -191,9 +191,9 @@ float GradNorm_func2D(float *B, float *B_x, float *B_y, float eta, long dimX, lo long i,j,index; float val1, val2, gradX, gradY, magn; #pragma omp parallel for shared(B, B_x, B_y) private(i,j,index,val1,val2,gradX,gradY,magn) - for(i=0; i 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 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 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 3) break; - } - + } + + /* check early stopping criteria */ + if ((epsil != 0.0f) && (ll % 5 == 0)) { + re = 0.0f; re1 = 0.0f; + for(j=0; j 3) break; + } + } /*end of iterations*/ free(D1_LLT);free(D2_LLT);free(D3_LLT); free(D1_ROF);free(D2_ROF);free(D3_ROF); if (epsil != 0.0f) free(Output_prev); - + /*adding info into info_vector */ infovector[0] = (float)(ll); /*iterations number (if stopped earlier based on tolerance)*/ infovector[1] = re; /* reached tolerance */ - return 0; + return 0; } /*************************************************************************/ @@ -130,65 +130,65 @@ float LLT_ROF_CPU_main(float *Input, float *Output, float *infovector, float lam /*************************************************************************/ float der2D_LLT(float *U, float *D1, float *D2, long dimX, long dimY, long dimZ) { - long i, j, index, i_p, i_m, j_m, j_p; - float dxx, dyy, denom_xx, denom_yy; + long i, j, index, i_p, i_m, j_m, j_p; + float dxx, dyy, denom_xx, denom_yy; #pragma omp parallel for shared(U,D1,D2) private(i, j, index, i_p, i_m, j_m, j_p, denom_xx, denom_yy, dxx, dyy) - for (i = 0; i 1) { #pragma omp parallel for shared (A, D1, dimX, dimY, dimZ) private(index, i, j, k, i1, j1, k1, i2, j2, k2, NOMx_1,NOMy_1,NOMy_0,NOMz_1,NOMz_0,denom1,denom2,denom3,T1) - for(j=0; j= dimX) i1 = i-1; i2 = i - 1; if (i2 < 0) i2 = i+1; @@ -213,17 +213,17 @@ float D1_func_ROF(float *A, float *D1, long dimX, long dimY, long dimZ) j2 = j - 1; if (j2 < 0) j2 = j+1; k1 = k + 1; if (k1 >= dimZ) k1 = k-1; k2 = k - 1; if (k2 < 0) k2 = k+1; - + /* Forward-backward differences */ NOMx_1 = A[(dimX*dimY)*k + j1*dimX + i] - A[index]; /* x+ */ NOMy_1 = A[(dimX*dimY)*k + j*dimX + i1] - A[index]; /* y+ */ /*NOMx_0 = (A[(i)*dimY + j] - A[(i2)*dimY + j]); */ /* x- */ NOMy_0 = A[index] - A[(dimX*dimY)*k + j*dimX + i2]; /* y- */ - + NOMz_1 = A[(dimX*dimY)*k1 + j*dimX + i] - A[index]; /* z+ */ NOMz_0 = A[index] - A[(dimX*dimY)*k2 + j*dimX + i]; /* z- */ - - + + denom1 = NOMx_1*NOMx_1; denom2 = 0.5f*(signLLT(NOMy_1) + signLLT(NOMy_0))*(MIN(fabs(NOMy_1),fabs(NOMy_0))); denom2 = denom2*denom2; @@ -237,19 +237,19 @@ float D1_func_ROF(float *A, float *D1, long dimX, long dimY, long dimZ) #pragma omp parallel for shared (A, D1, dimX, dimY) private(i, j, i1, j1, i2, j2,NOMx_1,NOMy_1,NOMy_0,denom1,denom2,T1,index) for(j=0; j= dimX) i1 = i-1; i2 = i - 1; if (i2 < 0) i2 = i+1; j1 = j + 1; if (j1 >= dimY) j1 = j-1; j2 = j - 1; if (j2 < 0) j2 = j+1; - + /* Forward-backward differences */ NOMx_1 = A[j1*dimX + i] - A[index]; /* x+ */ NOMy_1 = A[j*dimX + i1] - A[index]; /* y+ */ /*NOMx_0 = (A[(i)*dimY + j] - A[(i2)*dimY + j]); */ /* x- */ NOMy_0 = A[index] - A[(j)*dimX + i2]; /* y- */ - + denom1 = NOMx_1*NOMx_1; denom2 = 0.5f*(signLLT(NOMy_1) + signLLT(NOMy_0))*(MIN(fabs(NOMy_1),fabs(NOMy_0))); denom2 = denom2*denom2; @@ -264,13 +264,13 @@ float D2_func_ROF(float *A, float *D2, long dimX, long dimY, long dimZ) { float NOMx_1, NOMy_1, NOMx_0, NOMz_1, NOMz_0, denom1, denom2, denom3, T2; long i,j,k,i1,i2,k1,j1,j2,k2,index; - + if (dimZ > 1) { #pragma omp parallel for shared (A, D2, dimX, dimY, dimZ) private(index, i, j, k, i1, j1, k1, i2, j2, k2, NOMx_1, NOMy_1, NOMx_0, NOMz_1, NOMz_0, denom1, denom2, denom3, T2) - for(j=0; j= dimX) i1 = i-1; i2 = i - 1; if (i2 < 0) i2 = i+1; @@ -278,16 +278,16 @@ float D2_func_ROF(float *A, float *D2, long dimX, long dimY, long dimZ) j2 = j - 1; if (j2 < 0) j2 = j+1; k1 = k + 1; if (k1 >= dimZ) k1 = k-1; k2 = k - 1; if (k2 < 0) k2 = k+1; - - + + /* Forward-backward differences */ NOMx_1 = A[(dimX*dimY)*k + (j1)*dimX + i] - A[index]; /* x+ */ NOMy_1 = A[(dimX*dimY)*k + (j)*dimX + i1] - A[index]; /* y+ */ NOMx_0 = A[index] - A[(dimX*dimY)*k + (j2)*dimX + i]; /* x- */ NOMz_1 = A[(dimX*dimY)*k1 + j*dimX + i] - A[index]; /* z+ */ NOMz_0 = A[index] - A[(dimX*dimY)*k2 + (j)*dimX + i]; /* z- */ - - + + denom1 = NOMy_1*NOMy_1; denom2 = 0.5f*(signLLT(NOMx_1) + signLLT(NOMx_0))*(MIN(fabs(NOMx_1),fabs(NOMx_0))); denom2 = denom2*denom2; @@ -301,19 +301,19 @@ float D2_func_ROF(float *A, float *D2, long dimX, long dimY, long dimZ) #pragma omp parallel for shared (A, D2, dimX, dimY) private(i, j, i1, j1, i2, j2, NOMx_1,NOMy_1,NOMx_0,denom1,denom2,T2,index) for(j=0; j= dimX) i1 = i-1; i2 = i - 1; if (i2 < 0) i2 = i+1; j1 = j + 1; if (j1 >= dimY) j1 = j-1; j2 = j - 1; if (j2 < 0) j2 = j+1; - + /* Forward-backward differences */ NOMx_1 = A[j1*dimX + i] - A[index]; /* x+ */ NOMy_1 = A[j*dimX + i1] - A[index]; /* y+ */ NOMx_0 = A[index] - A[j2*dimX + i]; /* x- */ /*NOMy_0 = A[(i)*dimY + j] - A[(i)*dimY + j2]; */ /* y- */ - + denom1 = NOMy_1*NOMy_1; denom2 = 0.5f*(signLLT(NOMx_1) + signLLT(NOMx_0))*(MIN(fabs(NOMx_1),fabs(NOMx_0))); denom2 = denom2*denom2; @@ -329,12 +329,12 @@ float D3_func_ROF(float *A, float *D3, long dimX, long dimY, long dimZ) { float NOMx_1, NOMy_1, NOMx_0, NOMy_0, NOMz_1, denom1, denom2, denom3, T3; long index,i,j,k,i1,i2,k1,j1,j2,k2; - + #pragma omp parallel for shared (A, D3, dimX, dimY, dimZ) private(index, i, j, k, i1, j1, k1, i2, j2, k2, NOMx_1, NOMy_1, NOMy_0, NOMx_0, NOMz_1, denom1, denom2, denom3, T3) - for(j=0; j= dimX) i1 = i-1; i2 = i - 1; if (i2 < 0) i2 = i+1; @@ -342,7 +342,7 @@ float D3_func_ROF(float *A, float *D3, long dimX, long dimY, long dimZ) j2 = j - 1; if (j2 < 0) j2 = j+1; k1 = k + 1; if (k1 >= dimZ) k1 = k-1; k2 = k - 1; if (k2 < 0) k2 = k+1; - + /* Forward-backward differences */ NOMx_1 = A[(dimX*dimY)*k + (j1)*dimX + i] - A[index]; /* x+ */ NOMy_1 = A[(dimX*dimY)*k + (j)*dimX + i1] - A[index]; /* y+ */ @@ -350,7 +350,7 @@ float D3_func_ROF(float *A, float *D3, long dimX, long dimY, long dimZ) NOMx_0 = A[index] - A[(dimX*dimY)*k + (j2)*dimX + i]; /* x- */ NOMz_1 = A[(dimX*dimY)*k1 + j*dimX + i] - A[index]; /* z+ */ /*NOMz_0 = A[(dimX*dimY)*k + (i)*dimY + j] - A[(dimX*dimY)*k2 + (i)*dimY + j]; */ /* z- */ - + denom1 = NOMz_1*NOMz_1; denom2 = 0.5f*(signLLT(NOMx_1) + signLLT(NOMx_0))*(MIN(fabs(NOMx_1),fabs(NOMx_0))); denom2 = denom2*denom2; @@ -368,69 +368,69 @@ float D3_func_ROF(float *A, float *D3, long dimX, long dimY, long dimZ) float Update2D_LLT_ROF(float *U0, float *U, float *D1_LLT, float *D2_LLT, float *D1_ROF, float *D2_ROF, float lambdaROF, float lambdaLLT, float tau, long dimX, long dimY, long dimZ) { - long i, j, index, i_p, i_m, j_m, j_p; - float div, laplc, dxx, dyy, dv1, dv2; + long i, j, index, i_p, i_m, j_m, j_p; + float div, laplc, dxx, dyy, dv1, dv2; #pragma omp parallel for shared(U,U0) private(i, j, index, i_p, i_m, j_m, j_p, laplc, div, dxx, dyy, dv1, dv2) - for (i = 0; i>>>**********/ float NLM_H1_2D(float *A, float *A_orig, unsigned short *H_i, unsigned short *H_j, float *Weights, long i, long j, long dimX, long dimY, int NumNeighb, float lambdaReg) { - long x, i1, j1, index, index_m; - float value = 0.0f, normweight = 0.0f; - - index_m = j*dimX+i; - for(x=0; x < NumNeighb; x++) { - index = (dimX*dimY*x) + j*dimX+i; - i1 = H_i[index]; - j1 = H_j[index]; - value += A[j1*dimX+i1]*Weights[index]; - normweight += Weights[index]; - } - A[index_m] = (lambdaReg*A_orig[index_m] + value)/(lambdaReg + normweight); + long x, i1, j1, index, index_m; + float value = 0.0f, normweight = 0.0f; + + index_m = j*dimX+i; + for(x=0; x < NumNeighb; x++) { + index = (dimX*dimY*x) + j*dimX+i; + i1 = H_i[index]; + j1 = H_j[index]; + value += A[j1*dimX+i1]*Weights[index]; + normweight += Weights[index]; + } + A[index_m] = (lambdaReg*A_orig[index_m] + value)/(lambdaReg + normweight); return *A; } /*3D version*/ float NLM_H1_3D(float *A, float *A_orig, unsigned short *H_i, unsigned short *H_j, unsigned short *H_k, float *Weights, long i, long j, long k, long dimX, long dimY, long dimZ, int NumNeighb, float lambdaReg) { - long x, i1, j1, k1, index; - float value = 0.0f, normweight = 0.0f; - - for(x=0; x < NumNeighb; x++) { - index = dimX*dimY*dimZ*x + (dimX*dimY*k) + j*dimX+i; - i1 = H_i[index]; - j1 = H_j[index]; - k1 = H_k[index]; - value += A[(dimX*dimY*k1) + j1*dimX+i1]*Weights[index]; - normweight += Weights[index]; - } + long x, i1, j1, k1, index; + float value = 0.0f, normweight = 0.0f; + + for(x=0; x < NumNeighb; x++) { + index = dimX*dimY*dimZ*x + (dimX*dimY*k) + j*dimX+i; + i1 = H_i[index]; + j1 = H_j[index]; + k1 = H_k[index]; + value += A[(dimX*dimY*k1) + j1*dimX+i1]*Weights[index]; + normweight += Weights[index]; + } A[(dimX*dimY*k) + j*dimX+i] = (lambdaReg*A_orig[(dimX*dimY*k) + j*dimX+i] + value)/(lambdaReg + normweight); return *A; } - /***********<<<
>>>**********/ float NLM_TV_2D(float *A, float *A_orig, unsigned short *H_i, unsigned short *H_j, float *Weights, long i, long j, long dimX, long dimY, int NumNeighb, float lambdaReg) { - long x, i1, j1, index, index_m; - float value = 0.0f, normweight = 0.0f, NLgrad_magn = 0.0f, NLCoeff; - - index_m = j*dimX+i; - - for(x=0; x < NumNeighb; x++) { - index = (dimX*dimY*x) + j*dimX+i; /*c*/ - i1 = H_i[index]; - j1 = H_j[index]; - NLgrad_magn += powf((A[j1*dimX+i1] - A[index_m]),2)*Weights[index]; - } - + long x, i1, j1, index, index_m; + float value = 0.0f, normweight = 0.0f, NLgrad_magn = 0.0f, NLCoeff; + + index_m = j*dimX+i; + + for(x=0; x < NumNeighb; x++) { + index = (dimX*dimY*x) + j*dimX+i; /*c*/ + i1 = H_i[index]; + j1 = H_j[index]; + NLgrad_magn += powf((A[j1*dimX+i1] - A[index_m]),2)*Weights[index]; + } + NLgrad_magn = sqrtf(NLgrad_magn); /*Non Local Gradients Magnitude */ NLCoeff = 2.0f*(1.0f/(NLgrad_magn + EPS)); - + for(x=0; x < NumNeighb; x++) { - index = (dimX*dimY*x) + j*dimX+i; /*c*/ - i1 = H_i[index]; - j1 = H_j[index]; + index = (dimX*dimY*x) + j*dimX+i; /*c*/ + i1 = H_i[index]; + j1 = H_j[index]; value += A[j1*dimX+i1]*NLCoeff*Weights[index]; normweight += Weights[index]*NLCoeff; } @@ -151,25 +150,25 @@ float NLM_TV_2D(float *A, float *A_orig, unsigned short *H_i, unsigned short *H_ /*3D version*/ float NLM_TV_3D(float *A, float *A_orig, unsigned short *H_i, unsigned short *H_j, unsigned short *H_k, float *Weights, long i, long j, long k, long dimX, long dimY, long dimZ, int NumNeighb, float lambdaReg) { - long x, i1, j1, k1, index; - float value = 0.0f, normweight = 0.0f, NLgrad_magn = 0.0f, NLCoeff; - - for(x=0; x < NumNeighb; x++) { - index = dimX*dimY*dimZ*x + (dimX*dimY*k) + j*dimX+i; - i1 = H_i[index]; - j1 = H_j[index]; - k1 = H_k[index]; - NLgrad_magn += powf((A[(dimX*dimY*k1) + j1*dimX+i1] - A[(dimX*dimY*k1) + j*dimX+i]),2)*Weights[index]; - } - + long x, i1, j1, k1, index; + float value = 0.0f, normweight = 0.0f, NLgrad_magn = 0.0f, NLCoeff; + + for(x=0; x < NumNeighb; x++) { + index = dimX*dimY*dimZ*x + (dimX*dimY*k) + j*dimX+i; + i1 = H_i[index]; + j1 = H_j[index]; + k1 = H_k[index]; + NLgrad_magn += powf((A[(dimX*dimY*k1) + j1*dimX+i1] - A[(dimX*dimY*k1) + j*dimX+i]),2)*Weights[index]; + } + NLgrad_magn = sqrtf(NLgrad_magn); /*Non Local Gradients Magnitude */ NLCoeff = 2.0f*(1.0f/(NLgrad_magn + EPS)); - + for(x=0; x < NumNeighb; x++) { - index = dimX*dimY*dimZ*x + (dimX*dimY*k) + j*dimX+i; - i1 = H_i[index]; - j1 = H_j[index]; - k1 = H_k[index]; + index = dimX*dimY*dimZ*x + (dimX*dimY*k) + j*dimX+i; + i1 = H_i[index]; + j1 = H_j[index]; + k1 = H_k[index]; value += A[(dimX*dimY*k1) + j1*dimX+i1]*NLCoeff*Weights[index]; normweight += Weights[index]*NLCoeff; } diff --git a/src/Core/regularisers_CPU/PatchSelect_core.c b/src/Core/regularisers_CPU/PatchSelect_core.c index 8581797..d80353b 100644 --- a/src/Core/regularisers_CPU/PatchSelect_core.c +++ b/src/Core/regularisers_CPU/PatchSelect_core.c @@ -77,13 +77,13 @@ float PatchSelect_CPU_main(float *A, unsigned short *H_i, unsigned short *H_j, u }} /*main neighb loop */ /* for each pixel store indeces of the most similar neighbours (patches) */ #pragma omp parallel for shared (A, Weights, H_i, H_j) private(i,j) - for(i=0; i<(long)(dimX); i++) { - for(j=0; j<(long)(dimY); j++) { + for(j=0; j<(long)(dimY); j++) { + for(i=0; i<(long)(dimX); i++) { Indeces2D(A, H_i, H_j, Weights, i, j, (long)(dimX), (long)(dimY), Eucl_Vec, NumNeighb, SearchWindow, SimilarWin, h2); }} } else { - /****************3D INPUT ***************/ + /****************3D INPUT ***************/ /* generate a 3D Gaussian kernel for NLM procedure */ Eucl_Vec = (float*) calloc ((2*SimilarWin+1)*(2*SimilarWin+1)*(2*SimilarWin+1),sizeof(float)); counterG = 0; @@ -93,12 +93,12 @@ float PatchSelect_CPU_main(float *A, unsigned short *H_i, unsigned short *H_j, u Eucl_Vec[counterG] = (float)exp(-(pow(((float) i), 2) + pow(((float) j), 2) + pow(((float) k), 2))/(2*SimilarWin*SimilarWin*SimilarWin)); counterG++; }}} /*main neighb loop */ - + /* for each voxel store indeces of the most similar neighbours (patches) */ #pragma omp parallel for shared (A, Weights, H_i, H_j, H_k) private(i,j,k) - for(i=0; i 1) D3_func(Output, D3, (long)(dimX), (long)(dimY), (long)(dimZ)); - TV_kernel(D1, D2, D3, Output, Input, lambdaPar, lambda_is_arr, tau, (long)(dimX), (long)(dimY), (long)(dimZ)); - - /* check early stopping criteria */ - if ((epsil != 0.0f) && (i % 5 == 0)) { + if ((epsil != 0.0f) && (i % 5 == 0)) copyIm(Output, Output_prev, (long)(dimX), (long)(dimY), (long)(dimZ)); + + /* calculate differences */ + D1_func(Output, D1, (long)(dimX), (long)(dimY), (long)(dimZ)); + D2_func(Output, D2, (long)(dimX), (long)(dimY), (long)(dimZ)); + if (dimZ > 1) D3_func(Output, D3, (long)(dimX), (long)(dimY), (long)(dimZ)); + TV_kernel(D1, D2, D3, Output, Input, lambdaPar, lambda_is_arr, tau, (long)(dimX), (long)(dimY), (long)(dimZ)); + + /* check early stopping criteria */ + if ((epsil != 0.0f) && (i % 5 == 0)) { re = 0.0f; re1 = 0.0f; - for(j=0; j 3) break; + for(j=0; j 3) break; + } + } free(D1);free(D2); free(D3); if (epsil != 0.0f) free(Output_prev); - + /*adding info into info_vector */ infovector[0] = (float)(i); /*iterations number (if stopped earlier based on tolerance)*/ infovector[1] = re; /* reached tolerance */ - - return 0; + + return 0; } /* calculate differences 1 */ @@ -103,13 +103,13 @@ float D1_func(float *A, float *D1, long dimX, long dimY, long dimZ) { float NOMx_1, NOMy_1, NOMy_0, NOMz_1, NOMz_0, denom1, denom2,denom3, T1; long i,j,k,i1,i2,k1,j1,j2,k2,index; - + if (dimZ > 1) { #pragma omp parallel for shared (A, D1, dimX, dimY, dimZ) private(index, i, j, k, i1, j1, k1, i2, j2, k2, NOMx_1,NOMy_1,NOMy_0,NOMz_1,NOMz_0,denom1,denom2,denom3,T1) - for(j=0; j= dimX) i1 = i-1; i2 = i - 1; if (i2 < 0) i2 = i+1; @@ -117,17 +117,17 @@ float D1_func(float *A, float *D1, long dimX, long dimY, long dimZ) j2 = j - 1; if (j2 < 0) j2 = j+1; k1 = k + 1; if (k1 >= dimZ) k1 = k-1; k2 = k - 1; if (k2 < 0) k2 = k+1; - + /* Forward-backward differences */ NOMx_1 = A[(dimX*dimY)*k + j1*dimX + i] - A[index]; /* x+ */ NOMy_1 = A[(dimX*dimY)*k + j*dimX + i1] - A[index]; /* y+ */ /*NOMx_0 = (A[(i)*dimY + j] - A[(i2)*dimY + j]); */ /* x- */ NOMy_0 = A[index] - A[(dimX*dimY)*k + j*dimX + i2]; /* y- */ - + NOMz_1 = A[(dimX*dimY)*k1 + j*dimX + i] - A[index]; /* z+ */ NOMz_0 = A[index] - A[(dimX*dimY)*k2 + j*dimX + i]; /* z- */ - - + + denom1 = NOMx_1*NOMx_1; denom2 = 0.5f*(sign(NOMy_1) + sign(NOMy_0))*(MIN(fabs(NOMy_1),fabs(NOMy_0))); denom2 = denom2*denom2; @@ -141,19 +141,19 @@ float D1_func(float *A, float *D1, long dimX, long dimY, long dimZ) #pragma omp parallel for shared (A, D1, dimX, dimY) private(i, j, i1, j1, i2, j2,NOMx_1,NOMy_1,NOMy_0,denom1,denom2,T1,index) for(j=0; j= dimX) i1 = i-1; i2 = i - 1; if (i2 < 0) i2 = i+1; j1 = j + 1; if (j1 >= dimY) j1 = j-1; j2 = j - 1; if (j2 < 0) j2 = j+1; - + /* Forward-backward differences */ NOMx_1 = A[j1*dimX + i] - A[index]; /* x+ */ NOMy_1 = A[j*dimX + i1] - A[index]; /* y+ */ /*NOMx_0 = (A[(i)*dimY + j] - A[(i2)*dimY + j]); */ /* x- */ NOMy_0 = A[index] - A[(j)*dimX + i2]; /* y- */ - + denom1 = NOMx_1*NOMx_1; denom2 = 0.5f*(sign(NOMy_1) + sign(NOMy_0))*(MIN(fabs(NOMy_1),fabs(NOMy_0))); denom2 = denom2*denom2; @@ -168,12 +168,12 @@ float D2_func(float *A, float *D2, long dimX, long dimY, long dimZ) { float NOMx_1, NOMy_1, NOMx_0, NOMz_1, NOMz_0, denom1, denom2, denom3, T2; long i,j,k,i1,i2,k1,j1,j2,k2,index; - + if (dimZ > 1) { #pragma omp parallel for shared (A, D2, dimX, dimY, dimZ) private(index, i, j, k, i1, j1, k1, i2, j2, k2, NOMx_1, NOMy_1, NOMx_0, NOMz_1, NOMz_0, denom1, denom2, denom3, T2) - for(j=0; j= dimX) i1 = i-1; @@ -182,15 +182,15 @@ float D2_func(float *A, float *D2, long dimX, long dimY, long dimZ) j2 = j - 1; if (j2 < 0) j2 = j+1; k1 = k + 1; if (k1 >= dimZ) k1 = k-1; k2 = k - 1; if (k2 < 0) k2 = k+1; - + /* Forward-backward differences */ NOMx_1 = A[(dimX*dimY)*k + (j1)*dimX + i] - A[index]; /* x+ */ NOMy_1 = A[(dimX*dimY)*k + (j)*dimX + i1] - A[index]; /* y+ */ NOMx_0 = A[index] - A[(dimX*dimY)*k + (j2)*dimX + i]; /* x- */ NOMz_1 = A[(dimX*dimY)*k1 + j*dimX + i] - A[index]; /* z+ */ NOMz_0 = A[index] - A[(dimX*dimY)*k2 + (j)*dimX + i]; /* z- */ - - + + denom1 = NOMy_1*NOMy_1; denom2 = 0.5f*(sign(NOMx_1) + sign(NOMx_0))*(MIN(fabs(NOMx_1),fabs(NOMx_0))); denom2 = denom2*denom2; @@ -204,19 +204,19 @@ float D2_func(float *A, float *D2, long dimX, long dimY, long dimZ) #pragma omp parallel for shared (A, D2, dimX, dimY) private(i, j, i1, j1, i2, j2, NOMx_1,NOMy_1,NOMx_0,denom1,denom2,T2,index) for(j=0; j= dimX) i1 = i-1; i2 = i - 1; if (i2 < 0) i2 = i+1; j1 = j + 1; if (j1 >= dimY) j1 = j-1; j2 = j - 1; if (j2 < 0) j2 = j+1; - + /* Forward-backward differences */ NOMx_1 = A[j1*dimX + i] - A[index]; /* x+ */ NOMy_1 = A[j*dimX + i1] - A[index]; /* y+ */ NOMx_0 = A[index] - A[j2*dimX + i]; /* x- */ /*NOMy_0 = A[(i)*dimY + j] - A[(i)*dimY + j2]; */ /* y- */ - + denom1 = NOMy_1*NOMy_1; denom2 = 0.5f*(sign(NOMx_1) + sign(NOMx_0))*(MIN(fabs(NOMx_1),fabs(NOMx_0))); denom2 = denom2*denom2; @@ -232,12 +232,12 @@ float D3_func(float *A, float *D3, long dimX, long dimY, long dimZ) { float NOMx_1, NOMy_1, NOMx_0, NOMy_0, NOMz_1, denom1, denom2, denom3, T3; long index,i,j,k,i1,i2,k1,j1,j2,k2; - + #pragma omp parallel for shared (A, D3, dimX, dimY, dimZ) private(index, i, j, k, i1, j1, k1, i2, j2, k2, NOMx_1, NOMy_1, NOMy_0, NOMx_0, NOMz_1, denom1, denom2, denom3, T3) - for(j=0; j= dimX) i1 = i-1; i2 = i - 1; if (i2 < 0) i2 = i+1; @@ -245,7 +245,7 @@ float D3_func(float *A, float *D3, long dimX, long dimY, long dimZ) j2 = j - 1; if (j2 < 0) j2 = j+1; k1 = k + 1; if (k1 >= dimZ) k1 = k-1; k2 = k - 1; if (k2 < 0) k2 = k+1; - + /* Forward-backward differences */ NOMx_1 = A[(dimX*dimY)*k + (j1)*dimX + i] - A[index]; /* x+ */ NOMy_1 = A[(dimX*dimY)*k + (j)*dimX + i1] - A[index]; /* y+ */ @@ -253,7 +253,7 @@ float D3_func(float *A, float *D3, long dimX, long dimY, long dimZ) NOMx_0 = A[index] - A[(dimX*dimY)*k + (j2)*dimX + i]; /* x- */ NOMz_1 = A[(dimX*dimY)*k1 + j*dimX + i] - A[index]; /* z+ */ /*NOMz_0 = A[(dimX*dimY)*k + (i)*dimY + j] - A[(dimX*dimY)*k2 + (i)*dimY + j]; */ /* z- */ - + denom1 = NOMz_1*NOMz_1; denom2 = 0.5f*(sign(NOMx_1) + sign(NOMx_0))*(MIN(fabs(NOMx_1),fabs(NOMx_0))); denom2 = denom2*denom2; @@ -270,12 +270,12 @@ float TV_kernel(float *D1, float *D2, float *D3, float *B, float *A, float *lamb { float dv1, dv2, dv3, lambda_val; long index,i,j,k,i1,i2,k1,j1,j2,k2; - + if (dimZ > 1) { #pragma omp parallel for shared (D1, D2, D3, B, dimX, dimY, dimZ) private(index, i, j, k, i1, j1, k1, i2, j2, k2, dv1,dv2,dv3,lambda_val) - for(j=0; j= dimZ) k1 = k-1; k2 = k - 1; if (k2 < 0) k2 = k+1; - + /*divergence components */ dv1 = D1[index] - D1[(dimX*dimY)*k + j2*dimX+i]; dv2 = D2[index] - D2[(dimX*dimY)*k + j*dimX+i2]; dv3 = D3[index] - D3[(dimX*dimY)*k2 + j*dimX+i]; - + B[index] += tau*(lambda_val*(dv1 + dv2 + dv3) - (B[index] - A[index])); }}} } @@ -305,11 +305,11 @@ float TV_kernel(float *D1, float *D2, float *D3, float *B, float *A, float *lamb i2 = i - 1; if (i2 < 0) i2 = i+1; j1 = j + 1; if (j1 >= dimY) j1 = j-1; j2 = j - 1; if (j2 < 0) j2 = j+1; - + /* divergence components */ dv1 = D1[index] - D1[j2*dimX + i]; dv2 = D2[index] - D2[j*dimX + i2]; - + B[index] += tau*(lambda_val*(dv1 + dv2) - (B[index] - A[index])); }} } diff --git a/src/Core/regularisers_CPU/SB_TV_core.c b/src/Core/regularisers_CPU/SB_TV_core.c index 8d80787..1215baf 100755 --- a/src/Core/regularisers_CPU/SB_TV_core.c +++ b/src/Core/regularisers_CPU/SB_TV_core.c @@ -1,146 +1,146 @@ /* -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 2017 Daniil Kazantsev -Copyright 2017 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 2017 Daniil Kazantsev + * Copyright 2017 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 "SB_TV_core.h" /* C-OMP implementation of Split Bregman - TV denoising-regularisation model (2D/3D) [1] -* -* Input Parameters: -* 1. Noisy image/volume -* 2. lambda - regularisation parameter -* 3. Number of iterations [OPTIONAL parameter] -* 4. eplsilon - tolerance constant [OPTIONAL parameter] -* 5. TV-type: 'iso' or 'l1' [OPTIONAL parameter] -* -* Output: -* [1] Filtered/regularized image/volume -* [2] Information vector which contains [iteration no., reached tolerance] -* -* [1]. Goldstein, T. and Osher, S., 2009. The split Bregman method for L1-regularized problems. SIAM journal on imaging sciences, 2(2), pp.323-343. -*/ + * + * Input Parameters: + * 1. Noisy image/volume + * 2. lambda - regularisation parameter + * 3. Number of iterations [OPTIONAL parameter] + * 4. eplsilon - tolerance constant [OPTIONAL parameter] + * 5. TV-type: 'iso' or 'l1' [OPTIONAL parameter] + * + * Output: + * [1] Filtered/regularized image/volume + * [2] Information vector which contains [iteration no., reached tolerance] + * + * [1]. Goldstein, T. and Osher, S., 2009. The split Bregman method for L1-regularized problems. SIAM journal on imaging sciences, 2(2), pp.323-343. + */ float SB_TV_CPU_main(float *Input, float *Output, float *infovector, float mu, int iter, float epsil, int methodTV, int dimX, int dimY, int dimZ) { - int ll; + int ll; long j, DimTotal; - float re, re1, lambda; + float re, re1, lambda; re = 0.0f; re1 = 0.0f; int count = 0; mu = 1.0f/mu; lambda = 2.0f*mu; - + float *Output_prev=NULL, *Dx=NULL, *Dy=NULL, *Bx=NULL, *By=NULL; DimTotal = (long)(dimX*dimY*dimZ); Output_prev = calloc(DimTotal, sizeof(float)); - Dx = calloc(DimTotal, sizeof(float)); - Dy = calloc(DimTotal, sizeof(float)); - Bx = calloc(DimTotal, sizeof(float)); - By = calloc(DimTotal, sizeof(float)); - - if (dimZ == 1) { - /* 2D case */ + Dx = calloc(DimTotal, sizeof(float)); + Dy = calloc(DimTotal, sizeof(float)); + Bx = calloc(DimTotal, sizeof(float)); + By = calloc(DimTotal, sizeof(float)); + + if (dimZ == 1) { + /* 2D case */ copyIm(Input, Output, (long)(dimX), (long)(dimY), 1l); /*initialize */ - + /* begin outer SB iterations */ for(ll=0; ll 3) break; + re = 0.0f; re1 = 0.0f; + for(j=0; j 3) break; } } - } - else { - /* 3D case */ - float *Dz=NULL, *Bz=NULL; - - Dz = calloc(DimTotal, sizeof(float)); - Bz = calloc(DimTotal, sizeof(float)); - + } + else { + /* 3D case */ + float *Dz=NULL, *Bz=NULL; + + Dz = calloc(DimTotal, sizeof(float)); + Bz = calloc(DimTotal, sizeof(float)); + copyIm(Input, Output, (long)(dimX), (long)(dimY), (long)(dimZ)); /*initialize */ - + /* begin outer SB iterations */ for(ll=0; ll 3) break; + re = 0.0f; re1 = 0.0f; + for(j=0; j 3) break; } } - free(Dz); free(Bz); - } - - free(Output_prev); free(Dx); free(Dy); free(Bx); free(By); - /*adding info into info_vector */ - infovector[0] = (float)(ll); /*iterations number (if stopped earlier based on tolerance)*/ - infovector[1] = re; /* reached tolerance */ - return 0; + free(Dz); free(Bz); + } + + free(Output_prev); free(Dx); free(Dy); free(Bx); free(By); + /*adding info into info_vector */ + infovector[0] = (float)(ll); /*iterations number (if stopped earlier based on tolerance)*/ + infovector[1] = re; /* reached tolerance */ + return 0; } /********************************************************************/ @@ -151,18 +151,18 @@ float gauss_seidel2D(float *U, float *A, float *U_prev, float *Dx, float *Dy, fl float sum, normConst; long i,j,i1,i2,j1,j2,index; normConst = 1.0f/(mu + 4.0f*lambda); - + #pragma omp parallel for shared(U) private(index,i,j,i1,i2,j1,j2,sum) - for(i=0; i divergence and projection of P*/ DivProjP_2D(U, U0, P1, P2, (long)(dimX), (long)(dimY), lambda, tau); - + /*get updated solution U*/ newU(U, U_old, (long)(dimX), (long)(dimY)); - + /*saving V into V_old*/ copyIm(V1, V1_old, (long)(dimX), (long)(dimY), 1l); copyIm(V2, V2_old, (long)(dimX), (long)(dimY), 1l); - + /* upd V*/ UpdV_2D(V1, V2, P1, P2, Q1, Q2, Q3, (long)(dimX), (long)(dimY), tau); - + /*get new V*/ newU(V1, V1_old, (long)(dimX), (long)(dimY)); newU(V2, V2_old, (long)(dimX), (long)(dimY)); - + /* check early stopping criteria */ if ((epsil != 0.0f) && (ll % 5 == 0)) { - re = 0.0f; re1 = 0.0f; - for(j=0; j 3) break; - } + re = 0.0f; re1 = 0.0f; + for(j=0; j 3) break; + } } /*end of iterations*/ } else { /*3D case*/ float *P3, *Q4, *Q5, *Q6, *V3, *V3_old; - + P3 = calloc(DimTotal, sizeof(float)); Q4 = calloc(DimTotal, sizeof(float)); Q5 = calloc(DimTotal, sizeof(float)); Q6 = calloc(DimTotal, sizeof(float)); V3 = calloc(DimTotal, sizeof(float)); V3_old = calloc(DimTotal, sizeof(float)); - + /* Primal-dual iterations begin here */ for(ll = 0; ll < iter; ll++) { - + /* Calculate Dual Variable P */ DualP_3D(U, V1, V2, V3, P1, P2, P3, (long)(dimX), (long)(dimY), (long)(dimZ), sigma); - + /*Projection onto convex set for P*/ ProjP_3D(P1, P2, P3, (long)(dimX), (long)(dimY), (long)(dimZ), alpha1); - + /* Calculate Dual Variable Q */ DualQ_3D(V1, V2, V3, Q1, Q2, Q3, Q4, Q5, Q6, (long)(dimX), (long)(dimY), (long)(dimZ), sigma); - + /*Projection onto convex set for Q*/ ProjQ_3D(Q1, Q2, Q3, Q4, Q5, Q6, (long)(dimX), (long)(dimY), (long)(dimZ), alpha0); - + /*saving U into U_old*/ copyIm(U, U_old, (long)(dimX), (long)(dimY), (long)(dimZ)); - + /*adjoint operation -> divergence and projection of P*/ DivProjP_3D(U, U0, P1, P2, P3, (long)(dimX), (long)(dimY), (long)(dimZ), lambda, tau); - + /*get updated solution U*/ newU3D(U, U_old, (long)(dimX), (long)(dimY), (long)(dimZ)); - + /*saving V into V_old*/ copyIm_3Ar(V1, V2, V3, V1_old, V2_old, V3_old, (long)(dimX), (long)(dimY), (long)(dimZ)); - + /* upd V*/ UpdV_3D(V1, V2, V3, P1, P2, P3, Q1, Q2, Q3, Q4, Q5, Q6, (long)(dimX), (long)(dimY), (long)(dimZ), tau); - + /*get new V*/ newU3D_3Ar(V1, V2, V3, V1_old, V2_old, V3_old, (long)(dimX), (long)(dimY), (long)(dimZ)); - + /* check early stopping criteria */ if ((epsil != 0.0f) && (ll % 5 == 0)) { - re = 0.0f; re1 = 0.0f; - for(j=0; j 3) break; - } - + re = 0.0f; re1 = 0.0f; + for(j=0; j 3) break; + } + } /*end of iterations*/ free(P3);free(Q4);free(Q5);free(Q6);free(V3);free(V3_old); } - + /*freeing*/ free(P1);free(P2);free(Q1);free(Q2);free(Q3);free(U_old); free(V1);free(V2);free(V1_old);free(V2_old); - + /*adding info into info_vector */ infovector[0] = (float)(ll); /*iterations number (if stopped earlier based on tolerance)*/ infovector[1] = re; /* reached tolerance */ - + return 0; } @@ -201,15 +201,15 @@ float DualP_2D(float *U, float *V1, float *V2, float *P1, float *P2, long dimX, { long i,j, index; #pragma omp parallel for shared(U,V1,V2,P1,P2) private(i,j,index) - for(i=0; i 1.0f) { @@ -236,8 +236,8 @@ float DualQ_2D(float *V1, float *V2, float *Q1, float *Q2, float *Q3, long dimX, long i,j,index; float q1, q2, q11, q22; #pragma omp parallel for shared(Q1,Q2,Q3,V1,V2) private(i,j,index,q1,q2,q11,q22) - for(i=0; i 1.0f) { @@ -391,9 +391,9 @@ float DualQ_3D(float *V1, float *V2, float *V3, float *Q1, float *Q2, float *Q3, long i,j,k,index; float q1, q2, q3, q11, q22, q33, q44, q55, q66; #pragma omp parallel for shared(Q1,Q2,Q3,Q4,Q5,Q6,V1,V2,V3) private(i,j,k,index,q1,q2,q3,q11,q22,q33,q44,q55,q66) - for(i=0; i