diff options
| author | Kazantsev <kjy41806@ws097.diamond.ac.uk> | 2018-09-14 12:29:26 +0100 |
|---|---|---|
| committer | Kazantsev <kjy41806@ws097.diamond.ac.uk> | 2018-09-14 12:29:26 +0100 |
| commit | c16582003fb3a52e528f0a415e932c2e869e9952 (patch) | |
| tree | b4c14324001bb21e7d595619db6cffd477cb5af4 | |
| parent | a7ce37a74e0f7917d78a3ecdbfcd0adda95d42a0 (diff) | |
| download | regularization-c16582003fb3a52e528f0a415e932c2e869e9952.tar.gz regularization-c16582003fb3a52e528f0a415e932c2e869e9952.tar.bz2 regularization-c16582003fb3a52e528f0a415e932c2e869e9952.tar.xz regularization-c16582003fb3a52e528f0a415e932c2e869e9952.zip | |
CPU modules fixed for big data
24 files changed, 339 insertions, 333 deletions
diff --git a/Core/inpainters_CPU/Diffusion_Inpaint_core.c b/Core/inpainters_CPU/Diffusion_Inpaint_core.c index 56e0421..08b168a 100644 --- a/Core/inpainters_CPU/Diffusion_Inpaint_core.c +++ b/Core/inpainters_CPU/Diffusion_Inpaint_core.c @@ -47,12 +47,12 @@ int signNDF_inc(float x) { float Diffusion_Inpaint_CPU_main(float *Input, unsigned char *Mask, float *Output, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int penaltytype, int dimX, int dimY, int dimZ) { - int i,pointsone; + long i, pointsone; float sigmaPar2; sigmaPar2 = sigmaPar/sqrt(2.0f); /* copy into output */ - copyIm(Input, Output, dimX, dimY, dimZ); + copyIm(Input, Output, (long)(dimX), (long)(dimY), (long)(dimZ)); pointsone = 0; for (i=0; i<dimY*dimX*dimZ; i++) if (Mask[i] == 1) pointsone++; @@ -63,17 +63,17 @@ float Diffusion_Inpaint_CPU_main(float *Input, unsigned char *Mask, float *Outpu if (dimZ == 1) { /* running 2D diffusion iterations */ for(i=0; i < iterationsNumb; i++) { - if (sigmaPar == 0.0f) LinearDiff_Inp_2D(Input, Mask, Output, lambdaPar, tau, dimX, dimY); /* linear diffusion (heat equation) */ - else NonLinearDiff_Inp_2D(Input, Mask, Output, lambdaPar, sigmaPar2, tau, penaltytype, dimX, dimY); /* nonlinear diffusion */ + if (sigmaPar == 0.0f) LinearDiff_Inp_2D(Input, Mask, Output, lambdaPar, tau, (long)(dimX), (long)(dimY)); /* linear diffusion (heat equation) */ + else NonLinearDiff_Inp_2D(Input, Mask, Output, lambdaPar, sigmaPar2, tau, penaltytype, (long)(dimX), (long)(dimY)); /* nonlinear diffusion */ } } else { /* running 3D diffusion iterations */ for(i=0; i < iterationsNumb; i++) { - if (sigmaPar == 0.0f) LinearDiff_Inp_3D(Input, Mask, Output, lambdaPar, tau, dimX, dimY, dimZ); - else NonLinearDiff_Inp_3D(Input, Mask, Output, lambdaPar, sigmaPar2, tau, penaltytype, dimX, dimY, dimZ); - } - } + if (sigmaPar == 0.0f) LinearDiff_Inp_3D(Input, Mask, Output, lambdaPar, tau, (long)(dimX), (long)(dimY), (long)(dimZ)); + else NonLinearDiff_Inp_3D(Input, Mask, Output, lambdaPar, sigmaPar2, tau, penaltytype, (long)(dimX), (long)(dimY), (long)(dimZ)); + } + } } return *Output; } @@ -81,9 +81,9 @@ float Diffusion_Inpaint_CPU_main(float *Input, unsigned char *Mask, float *Outpu /***************************2D Functions*****************************/ /********************************************************************/ /* linear diffusion (heat equation) */ -float LinearDiff_Inp_2D(float *Input, unsigned char *Mask, float *Output, float lambdaPar, float tau, int dimX, int dimY) +float LinearDiff_Inp_2D(float *Input, unsigned char *Mask, float *Output, float lambdaPar, float tau, long dimX, long dimY) { - int i,j,i1,i2,j1,j2,index; + long i,j,i1,i2,j1,j2,index; float e,w,n,s,e1,w1,n1,s1; #pragma omp parallel for shared(Input,Mask) private(index,i,j,i1,i2,j1,j2,e,w,n,s,e1,w1,n1,s1) @@ -116,9 +116,9 @@ float LinearDiff_Inp_2D(float *Input, unsigned char *Mask, float *Output, float } /* nonlinear diffusion */ -float NonLinearDiff_Inp_2D(float *Input, unsigned char *Mask, float *Output, float lambdaPar, float sigmaPar, float tau, int penaltytype, int dimX, int dimY) +float NonLinearDiff_Inp_2D(float *Input, unsigned char *Mask, float *Output, float lambdaPar, float sigmaPar, float tau, int penaltytype, long dimX, long dimY) { - int i,j,i1,i2,j1,j2,index; + long i,j,i1,i2,j1,j2,index; float e,w,n,s,e1,w1,n1,s1; #pragma omp parallel for shared(Input,Mask) private(index,i,j,i1,i2,j1,j2,e,w,n,s,e1,w1,n1,s1) @@ -189,9 +189,9 @@ float NonLinearDiff_Inp_2D(float *Input, unsigned char *Mask, float *Output, flo /***************************3D Functions*****************************/ /********************************************************************/ /* linear diffusion (heat equation) */ -float LinearDiff_Inp_3D(float *Input, unsigned char *Mask, float *Output, float lambdaPar, float tau, int dimX, int dimY, int dimZ) +float LinearDiff_Inp_3D(float *Input, unsigned char *Mask, float *Output, float lambdaPar, float tau, long dimX, long dimY, long dimZ) { - int i,j,k,i1,i2,j1,j2,k1,k2,index; + 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,Mask) private(index,i,j,i1,i2,j1,j2,e,w,n,s,e1,w1,n1,s1,k,k1,k2,u1,d1,u,d) @@ -231,9 +231,9 @@ for(k=0; k<dimZ; k++) { return *Output; } -float NonLinearDiff_Inp_3D(float *Input, unsigned char *Mask, float *Output, float lambdaPar, float sigmaPar, float tau, int penaltytype, int dimX, int dimY, int dimZ) +float NonLinearDiff_Inp_3D(float *Input, unsigned char *Mask, float *Output, float lambdaPar, float sigmaPar, float tau, int penaltytype, long dimX, long dimY, long dimZ) { - int i,j,k,i1,i2,j1,j2,k1,k2,index; + 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,Mask) private(index,i,j,i1,i2,j1,j2,e,w,n,s,e1,w1,n1,s1,k,k1,k2,u1,d1,u,d) diff --git a/Core/inpainters_CPU/Diffusion_Inpaint_core.h b/Core/inpainters_CPU/Diffusion_Inpaint_core.h index 0fc2608..a96fe79 100644 --- a/Core/inpainters_CPU/Diffusion_Inpaint_core.h +++ b/Core/inpainters_CPU/Diffusion_Inpaint_core.h @@ -51,10 +51,11 @@ limitations under the License. extern "C" { #endif CCPI_EXPORT float Diffusion_Inpaint_CPU_main(float *Input, unsigned char *Mask, float *Output, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int penaltytype, int dimX, int dimY, int dimZ); -CCPI_EXPORT float LinearDiff_Inp_2D(float *Input, unsigned char *Mask, float *Output, float lambdaPar, float tau, int dimX, int dimY); -CCPI_EXPORT float NonLinearDiff_Inp_2D(float *Input, unsigned char *Mask, float *Output, float lambdaPar, float sigmaPar, float tau, int penaltytype, int dimX, int dimY); -CCPI_EXPORT float LinearDiff_Inp_3D(float *Input, unsigned char *Mask, float *Output, float lambdaPar, float tau, int dimX, int dimY, int dimZ); -CCPI_EXPORT float NonLinearDiff_Inp_3D(float *Input, unsigned char *Mask, float *Output, float lambdaPar, float sigmaPar, float tau, int penaltytype, int dimX, int dimY, int dimZ); + +CCPI_EXPORT float LinearDiff_Inp_2D(float *Input, unsigned char *Mask, float *Output, float lambdaPar, float tau, long dimX, long dimY); +CCPI_EXPORT float NonLinearDiff_Inp_2D(float *Input, unsigned char *Mask, float *Output, float lambdaPar, float sigmaPar, float tau, int penaltytype, long dimX, long dimY); +CCPI_EXPORT float LinearDiff_Inp_3D(float *Input, unsigned char *Mask, float *Output, float lambdaPar, float tau, long dimX, long dimY, long dimZ); +CCPI_EXPORT float NonLinearDiff_Inp_3D(float *Input, unsigned char *Mask, float *Output, float lambdaPar, float sigmaPar, float tau, int penaltytype, long dimX, long dimY, long dimZ); #ifdef __cplusplus } #endif diff --git a/Core/inpainters_CPU/NonlocalMarching_Inpaint_core.c b/Core/inpainters_CPU/NonlocalMarching_Inpaint_core.c index 66be15c..b488ca4 100644 --- a/Core/inpainters_CPU/NonlocalMarching_Inpaint_core.c +++ b/Core/inpainters_CPU/NonlocalMarching_Inpaint_core.c @@ -35,7 +35,7 @@ * 1. Inpainted image or a sinogram * 2. updated mask * - * Reference: TBA + * Reference: D. Kazantsev (paper in preparation) */ float NonlocalMarching_Inpaint_main(float *Input, unsigned char *M, float *Output, unsigned char *M_upd, int SW_increment, int iterationsNumb, int trigger, int dimX, int dimY, int dimZ) diff --git a/Core/regularisers_CPU/Diffus4th_order_core.c b/Core/regularisers_CPU/Diffus4th_order_core.c index 069a42e..01f4f64 100644 --- a/Core/regularisers_CPU/Diffus4th_order_core.c +++ b/Core/regularisers_CPU/Diffus4th_order_core.c @@ -50,7 +50,7 @@ float Diffus4th_CPU_main(float *Input, float *Output, float lambdaPar, float sig W_Lapl = calloc(DimTotal, sizeof(float)); /* copy into output */ - copyIm(Input, Output, dimX, dimY, dimZ); + copyIm(Input, Output, (long)(dimX), (long)(dimY), (long)(dimZ)); if (dimZ == 1) { /* running 2D diffusion iterations */ @@ -58,7 +58,7 @@ float Diffus4th_CPU_main(float *Input, float *Output, float lambdaPar, float sig /* Calculating weighted Laplacian */ Weighted_Laplc2D(W_Lapl, Output, sigmaPar2, dimX, dimY); /* Perform iteration step */ - Diffusion_update_step2D(Output, Input, W_Lapl, lambdaPar, sigmaPar2, tau, dimX, dimY); + Diffusion_update_step2D(Output, Input, W_Lapl, lambdaPar, sigmaPar2, tau, (long)(dimX), (long)(dimY)); } } else { @@ -67,7 +67,7 @@ float Diffus4th_CPU_main(float *Input, float *Output, float lambdaPar, float sig /* 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, dimX, dimY, dimZ); + Diffusion_update_step3D(Output, Input, W_Lapl, lambdaPar, sigmaPar2, tau, (long)(dimX), (long)(dimY), (long)(dimZ)); } } free(W_Lapl); @@ -76,9 +76,9 @@ float Diffus4th_CPU_main(float *Input, float *Output, float lambdaPar, float sig /********************************************************************/ /***************************2D Functions*****************************/ /********************************************************************/ -float Weighted_Laplc2D(float *W_Lapl, float *U0, float sigma, int dimX, int dimY) +float Weighted_Laplc2D(float *W_Lapl, float *U0, float sigma, long dimX, long dimY) { - int i,j,i1,i2,j1,j2,index; + 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) @@ -125,9 +125,9 @@ float Weighted_Laplc2D(float *W_Lapl, float *U0, float sigma, int dimX, int dimY return *W_Lapl; } -float Diffusion_update_step2D(float *Output, float *Input, float *W_Lapl, float lambdaPar, float sigmaPar2, float tau, int dimX, int dimY) +float Diffusion_update_step2D(float *Output, float *Input, float *W_Lapl, float lambdaPar, float sigmaPar2, float tau, long dimX, long dimY) { - int i,j,i1,i2,j1,j2,index; + long i,j,i1,i2,j1,j2,index; float gradXXc, gradYYc; #pragma omp parallel for shared(Output, Input, W_Lapl) private(i,j,i1,i2,j1,j2,index,gradXXc,gradYYc) @@ -152,9 +152,9 @@ float Diffusion_update_step2D(float *Output, float *Input, float *W_Lapl, float /********************************************************************/ /***************************3D Functions*****************************/ /********************************************************************/ -float Weighted_Laplc3D(float *W_Lapl, float *U0, float sigma, int dimX, int dimY, int dimZ) +float Weighted_Laplc3D(float *W_Lapl, float *U0, float sigma, long dimX, long dimY, long dimZ) { - int i,j,k,i1,i2,j1,j2,k1,k2,index; + long i,j,k,i1,i2,j1,j2,k1,k2,index; float gradX, gradX_sq, gradY, gradY_sq, gradXX, gradYY, gradXY, xy_2, denom, V_norm, V_orth, c, c_sq, gradZ, gradZ_sq, gradZZ, gradXZ, gradYZ, xyz_1, xyz_2; #pragma omp parallel for shared(W_Lapl) private(i,j,k,i1,i2,j1,j2,k1,k2,index,gradX, gradX_sq, gradY, gradY_sq, gradXX, gradYY, gradXY, xy_2, denom, V_norm, V_orth, c, c_sq, gradZ, gradZ_sq, gradZZ, gradXZ, gradYZ, xyz_1, xyz_2) @@ -216,9 +216,9 @@ float Weighted_Laplc3D(float *W_Lapl, float *U0, float sigma, int dimX, int dimY return *W_Lapl; } -float Diffusion_update_step3D(float *Output, float *Input, float *W_Lapl, float lambdaPar, float sigmaPar2, float tau, int dimX, int dimY, int dimZ) +float Diffusion_update_step3D(float *Output, float *Input, float *W_Lapl, float lambdaPar, float sigmaPar2, float tau, long dimX, long dimY, long dimZ) { - int i,j,i1,i2,j1,j2,index, k,k1,k2; + long i,j,i1,i2,j1,j2,index,k,k1,k2; float gradXXc, gradYYc, gradZZc; #pragma omp parallel for shared(Output, Input, W_Lapl) private(i,j,i1,i2,j1,j2,k,k1,k2,index,gradXXc,gradYYc,gradZZc) diff --git a/Core/regularisers_CPU/Diffus4th_order_core.h b/Core/regularisers_CPU/Diffus4th_order_core.h index 7b8fc08..d81afcb 100644 --- a/Core/regularisers_CPU/Diffus4th_order_core.h +++ b/Core/regularisers_CPU/Diffus4th_order_core.h @@ -46,10 +46,10 @@ limitations under the License. extern "C" { #endif CCPI_EXPORT float Diffus4th_CPU_main(float *Input, float *Output, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int dimX, int dimY, int dimZ); -CCPI_EXPORT float Weighted_Laplc2D(float *W_Lapl, float *U0, float sigma, int dimX, int dimY); -CCPI_EXPORT float Diffusion_update_step2D(float *Output, float *Input, float *W_Lapl, float lambdaPar, float sigmaPar2, float tau, int dimX, int dimY); -CCPI_EXPORT float Weighted_Laplc3D(float *W_Lapl, float *U0, float sigma, int dimX, int dimY, int dimZ); -CCPI_EXPORT float Diffusion_update_step3D(float *Output, float *Input, float *W_Lapl, float lambdaPar, float sigmaPar2, float tau, int dimX, int dimY, int dimZ); +CCPI_EXPORT float Weighted_Laplc2D(float *W_Lapl, float *U0, float sigma, long dimX, long dimY); +CCPI_EXPORT float Diffusion_update_step2D(float *Output, float *Input, float *W_Lapl, float lambdaPar, float sigmaPar2, float tau, long dimX, long dimY); +CCPI_EXPORT float Weighted_Laplc3D(float *W_Lapl, float *U0, float sigma, long dimX, long dimY, long dimZ); +CCPI_EXPORT float Diffusion_update_step3D(float *Output, float *Input, float *W_Lapl, float lambdaPar, float sigmaPar2, float tau, long dimX, long dimY, long dimZ); #ifdef __cplusplus } #endif diff --git a/Core/regularisers_CPU/Diffusion_core.c b/Core/regularisers_CPU/Diffusion_core.c index 51d0a57..b765796 100644 --- a/Core/regularisers_CPU/Diffusion_core.c +++ b/Core/regularisers_CPU/Diffusion_core.c @@ -55,20 +55,20 @@ float Diffusion_CPU_main(float *Input, float *Output, float lambdaPar, float sig sigmaPar2 = sigmaPar/sqrt(2.0f); /* copy into output */ - copyIm(Input, Output, dimX, dimY, dimZ); + copyIm(Input, Output, (long)(dimX), (long)(dimY), (long)(dimZ)); if (dimZ == 1) { /* running 2D diffusion iterations */ for(i=0; i < iterationsNumb; i++) { - if (sigmaPar == 0.0f) LinearDiff2D(Input, Output, lambdaPar, tau, dimX, dimY); /* linear diffusion (heat equation) */ - else NonLinearDiff2D(Input, Output, lambdaPar, sigmaPar2, tau, penaltytype, dimX, dimY); /* nonlinear diffusion */ + if (sigmaPar == 0.0f) LinearDiff2D(Input, Output, lambdaPar, tau, (long)(dimX), (long)(dimY)); /* linear diffusion (heat equation) */ + else NonLinearDiff2D(Input, Output, lambdaPar, sigmaPar2, tau, penaltytype, (long)(dimX), (long)(dimY)); /* nonlinear diffusion */ } } else { /* running 3D diffusion iterations */ for(i=0; i < iterationsNumb; i++) { - if (sigmaPar == 0.0f) LinearDiff3D(Input, Output, lambdaPar, tau, dimX, dimY, dimZ); - else NonLinearDiff3D(Input, Output, lambdaPar, sigmaPar2, tau, penaltytype, dimX, dimY, dimZ); + if (sigmaPar == 0.0f) LinearDiff3D(Input, Output, lambdaPar, tau, (long)(dimX), (long)(dimY), (long)(dimZ)); + else NonLinearDiff3D(Input, Output, lambdaPar, sigmaPar2, tau, penaltytype, (long)(dimX), (long)(dimY), (long)(dimZ)); } } return *Output; @@ -79,9 +79,9 @@ float Diffusion_CPU_main(float *Input, float *Output, float lambdaPar, float sig /***************************2D Functions*****************************/ /********************************************************************/ /* linear diffusion (heat equation) */ -float LinearDiff2D(float *Input, float *Output, float lambdaPar, float tau, int dimX, int dimY) +float LinearDiff2D(float *Input, float *Output, float lambdaPar, float tau, long dimX, long dimY) { - int i,j,i1,i2,j1,j2,index; + 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) @@ -111,9 +111,9 @@ float LinearDiff2D(float *Input, float *Output, float lambdaPar, float tau, int } /* nonlinear diffusion */ -float NonLinearDiff2D(float *Input, float *Output, float lambdaPar, float sigmaPar, float tau, int penaltytype, int dimX, int dimY) +float NonLinearDiff2D(float *Input, float *Output, float lambdaPar, float sigmaPar, float tau, int penaltytype, long dimX, long dimY) { - int i,j,i1,i2,j1,j2,index; + 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) @@ -181,9 +181,9 @@ float NonLinearDiff2D(float *Input, float *Output, float lambdaPar, float sigmaP /***************************3D Functions*****************************/ /********************************************************************/ /* linear diffusion (heat equation) */ -float LinearDiff3D(float *Input, float *Output, float lambdaPar, float tau, int dimX, int dimY, int dimZ) +float LinearDiff3D(float *Input, float *Output, float lambdaPar, float tau, long dimX, long dimY, long dimZ) { - int i,j,k,i1,i2,j1,j2,k1,k2,index; + 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) @@ -219,9 +219,9 @@ for(k=0; k<dimZ; k++) { return *Output; } -float NonLinearDiff3D(float *Input, float *Output, float lambdaPar, float sigmaPar, float tau, int penaltytype, int dimX, int dimY, int dimZ) +float NonLinearDiff3D(float *Input, float *Output, float lambdaPar, float sigmaPar, float tau, int penaltytype, long dimX, long dimY, long dimZ) { - int i,j,k,i1,i2,j1,j2,k1,k2,index; + 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) diff --git a/Core/regularisers_CPU/Diffusion_core.h b/Core/regularisers_CPU/Diffusion_core.h index 0b4149a..cc36dad 100644 --- a/Core/regularisers_CPU/Diffusion_core.h +++ b/Core/regularisers_CPU/Diffusion_core.h @@ -50,10 +50,10 @@ limitations under the License. extern "C" { #endif CCPI_EXPORT float Diffusion_CPU_main(float *Input, float *Output, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int penaltytype, int dimX, int dimY, int dimZ); -CCPI_EXPORT float LinearDiff2D(float *Input, float *Output, float lambdaPar, float tau, int dimX, int dimY); -CCPI_EXPORT float NonLinearDiff2D(float *Input, float *Output, float lambdaPar, float sigmaPar, float tau, int penaltytype, int dimX, int dimY); -CCPI_EXPORT float LinearDiff3D(float *Input, float *Output, float lambdaPar, float tau, int dimX, int dimY, int dimZ); -CCPI_EXPORT float NonLinearDiff3D(float *Input, float *Output, float lambdaPar, float sigmaPar, float tau, int penaltytype, int dimX, int dimY, int dimZ); +CCPI_EXPORT float LinearDiff2D(float *Input, float *Output, float lambdaPar, float tau, long dimX, long dimY); +CCPI_EXPORT float NonLinearDiff2D(float *Input, float *Output, float lambdaPar, float sigmaPar, float tau, int penaltytype, long dimX, long dimY); +CCPI_EXPORT float LinearDiff3D(float *Input, float *Output, float lambdaPar, float tau, long dimX, long dimY, long dimZ); +CCPI_EXPORT float NonLinearDiff3D(float *Input, float *Output, float lambdaPar, float sigmaPar, float tau, int penaltytype, long dimX, long dimY, long dimZ); #ifdef __cplusplus } #endif diff --git a/Core/regularisers_CPU/FGP_TV_core.c b/Core/regularisers_CPU/FGP_TV_core.c index 15bb45c..d828d48 100644 --- a/Core/regularisers_CPU/FGP_TV_core.c +++ b/Core/regularisers_CPU/FGP_TV_core.c @@ -39,7 +39,8 @@ limitations under the License. float TV_FGP_CPU_main(float *Input, float *Output, float lambdaPar, int iterationsNumb, float epsil, int methodTV, int nonneg, int printM, int dimX, int dimY, int dimZ) { - int ll, j, DimTotal; + int ll; + long j, DimTotal; float re, re1; float tk = 1.0f; float tkp1=1.0f; @@ -48,7 +49,7 @@ float TV_FGP_CPU_main(float *Input, float *Output, float lambdaPar, int iteratio if (dimZ <= 1) { /*2D case */ float *Output_prev=NULL, *P1=NULL, *P2=NULL, *P1_prev=NULL, *P2_prev=NULL, *R1=NULL, *R2=NULL; - DimTotal = dimX*dimY; + DimTotal = (long)(dimX*dimY); Output_prev = calloc(DimTotal, sizeof(float)); P1 = calloc(DimTotal, sizeof(float)); @@ -62,13 +63,13 @@ float TV_FGP_CPU_main(float *Input, float *Output, float lambdaPar, int iteratio for(ll=0; ll<iterationsNumb; ll++) { /* computing the gradient of the objective function */ - Obj_func2D(Input, Output, R1, R2, lambdaPar, dimX, dimY); + Obj_func2D(Input, Output, R1, R2, lambdaPar, (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_func2D(P1, P2, Output, R1, R2, lambdaPar, dimX, dimY); + Grad_func2D(P1, P2, Output, R1, R2, lambdaPar, (long)(dimY), (long)(dimZ)); /* projection step */ Proj_func2D(P1, P2, methodTV, DimTotal); @@ -89,9 +90,9 @@ float TV_FGP_CPU_main(float *Input, float *Output, float lambdaPar, int iteratio if (count > 4) break; /*storing old values*/ - copyIm(Output, Output_prev, dimX, dimY, 1); - copyIm(P1, P1_prev, dimX, dimY, 1); - copyIm(P2, P2_prev, dimX, dimY, 1); + copyIm(Output, Output_prev, (long)(dimY), (long)(dimZ), 1l); + copyIm(P1, P1_prev, (long)(dimY), (long)(dimZ), 1l); + copyIm(P2, P2_prev, (long)(dimY), (long)(dimZ), 1l); tk = tkp1; } if (printM == 1) printf("FGP-TV iterations stopped at iteration %i \n", ll); @@ -100,7 +101,7 @@ float TV_FGP_CPU_main(float *Input, float *Output, float lambdaPar, int iteratio 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 = dimX*dimY*dimZ; + DimTotal = (long)(dimX*dimY*dimZ); Output_prev = calloc(DimTotal, sizeof(float)); P1 = calloc(DimTotal, sizeof(float)); @@ -117,13 +118,13 @@ float TV_FGP_CPU_main(float *Input, float *Output, float lambdaPar, int iteratio for(ll=0; ll<iterationsNumb; ll++) { /* computing the gradient of the objective function */ - Obj_func3D(Input, Output, R1, R2, R3, lambdaPar, dimX, dimY, dimZ); + 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, dimX, dimY, dimZ); + 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); @@ -145,10 +146,10 @@ float TV_FGP_CPU_main(float *Input, float *Output, float lambdaPar, int iteratio if (count > 4) break; /*storing old values*/ - copyIm(Output, Output_prev, dimX, dimY, dimZ); - copyIm(P1, P1_prev, dimX, dimY, dimZ); - copyIm(P2, P2_prev, dimX, dimY, dimZ); - copyIm(P3, P3_prev, dimX, dimY, dimZ); + copyIm(Output, Output_prev, (long)(dimX), (long)(dimY), (long)(dimZ)); + 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 (printM == 1) printf("FGP-TV iterations stopped at iteration %i \n", ll); @@ -157,10 +158,10 @@ float TV_FGP_CPU_main(float *Input, float *Output, float lambdaPar, int iteratio return *Output; } -float Obj_func2D(float *A, float *D, float *R1, float *R2, float lambda, int dimX, int dimY) +float Obj_func2D(float *A, float *D, float *R1, float *R2, float lambda, long dimX, long dimY) { float val1, val2; - int i,j,index; + long i,j,index; #pragma omp parallel for shared(A,D,R1,R2) private(index,i,j,val1,val2) for(i=0; i<dimX; i++) { for(j=0; j<dimY; j++) { @@ -172,10 +173,10 @@ float Obj_func2D(float *A, float *D, float *R1, float *R2, float lambda, int dim }} return *D; } -float Grad_func2D(float *P1, float *P2, float *D, float *R1, float *R2, float lambda, int dimX, int dimY) +float Grad_func2D(float *P1, float *P2, float *D, float *R1, float *R2, float lambda, long dimX, long dimY) { float val1, val2, multip; - int i,j,index; + long i,j,index; multip = (1.0f/(8.0f*lambda)); #pragma omp parallel for shared(P1,P2,D,R1,R2,multip) private(index,i,j,val1,val2) for(i=0; i<dimX; i++) { @@ -189,10 +190,10 @@ 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 DimTotal) +float Proj_func2D(float *P1, float *P2, int methTV, long DimTotal) { float val1, val2, denom, sq_denom; - int i; + long i; if (methTV == 0) { /* isotropic TV*/ #pragma omp parallel for shared(P1,P2) private(i,denom,sq_denom) @@ -219,9 +220,9 @@ float Proj_func2D(float *P1, float *P2, int methTV, int DimTotal) } return 1; } -float Rupd_func2D(float *P1, float *P1_old, float *P2, float *P2_old, float *R1, float *R2, float tkp1, float tk, int DimTotal) +float Rupd_func2D(float *P1, float *P1_old, float *P2, float *P2_old, float *R1, float *R2, float tkp1, float tk, long DimTotal) { - int i; + long i; float multip; multip = ((tk-1.0f)/tkp1); #pragma omp parallel for shared(P1,P2,P1_old,P2_old,R1,R2,multip) private(i) @@ -234,10 +235,10 @@ float Rupd_func2D(float *P1, float *P1_old, float *P2, float *P2_old, float *R1, /* 3D-case related Functions */ /*****************************************************************/ -float Obj_func3D(float *A, float *D, float *R1, float *R2, float *R3, float lambda, int dimX, int dimY, int dimZ) +float Obj_func3D(float *A, float *D, float *R1, float *R2, float *R3, float lambda, long dimX, long dimY, long dimZ) { float val1, val2, val3; - int i,j,k,index; + long i,j,k,index; #pragma omp parallel for shared(A,D,R1,R2,R3) private(index,i,j,k,val1,val2,val3) for(i=0; i<dimX; i++) { for(j=0; j<dimY; j++) { @@ -251,10 +252,10 @@ float Obj_func3D(float *A, float *D, float *R1, float *R2, float *R3, float lamb }}} return *D; } -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) +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) { float val1, val2, val3, multip; - int i,j,k, index; + long i,j,k, index; multip = (1.0f/(26.0f*lambda)); #pragma omp parallel for shared(P1,P2,P3,D,R1,R2,R3,multip) private(index,i,j,k,val1,val2,val3) for(i=0; i<dimX; i++) { @@ -271,10 +272,10 @@ 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 DimTotal) +float Proj_func3D(float *P1, float *P2, float *P3, int methTV, long DimTotal) { float val1, val2, val3, denom, sq_denom; - int i; + long i; if (methTV == 0) { /* isotropic TV*/ #pragma omp parallel for shared(P1,P2,P3) private(i,val1,val2,val3,sq_denom) @@ -305,9 +306,9 @@ float Proj_func3D(float *P1, float *P2, float *P3, int methTV, int DimTotal) } 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 DimTotal) +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) { - int i; + long 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(i) diff --git a/Core/regularisers_CPU/FGP_TV_core.h b/Core/regularisers_CPU/FGP_TV_core.h index 37e32f7..3418604 100644 --- a/Core/regularisers_CPU/FGP_TV_core.h +++ b/Core/regularisers_CPU/FGP_TV_core.h @@ -49,15 +49,15 @@ extern "C" { #endif CCPI_EXPORT float TV_FGP_CPU_main(float *Input, float *Output, float lambdaPar, int iterationsNumb, float epsil, int methodTV, int nonneg, int printM, int dimX, int dimY, int dimZ); -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 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_func2D(float *A, float *D, float *R1, float *R2, float lambda, long dimX, long dimY); +CCPI_EXPORT float Grad_func2D(float *P1, float *P2, float *D, float *R1, float *R2, float lambda, long dimX, long dimY); +CCPI_EXPORT float Proj_func2D(float *P1, float *P2, int methTV, long DimTotal); +CCPI_EXPORT float Rupd_func2D(float *P1, float *P1_old, float *P2, float *P2_old, float *R1, float *R2, float tkp1, float tk, long 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 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); +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 } #endif diff --git a/Core/regularisers_CPU/FGP_dTV_core.c b/Core/regularisers_CPU/FGP_dTV_core.c index f6b4f79..17b75ff 100644 --- a/Core/regularisers_CPU/FGP_dTV_core.c +++ b/Core/regularisers_CPU/FGP_dTV_core.c @@ -44,7 +44,8 @@ limitations under the License. float dTV_FGP_CPU_main(float *Input, float *InputRef, float *Output, float lambdaPar, int iterationsNumb, float epsil, float eta, int methodTV, int nonneg, int printM, int dimX, int dimY, int dimZ) { - int ll, j, DimTotal; + int ll; + long j, DimTotal; float re, re1; float tk = 1.0f; float tkp1=1.0f; @@ -53,7 +54,7 @@ float dTV_FGP_CPU_main(float *Input, float *InputRef, float *Output, float lambd if (dimZ <= 1) { /*2D case */ 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 = dimX*dimY; + DimTotal = (long)(dimX*dimY); Output_prev = calloc(DimTotal, sizeof(float)); P1 = calloc(DimTotal, sizeof(float)); @@ -66,22 +67,22 @@ float dTV_FGP_CPU_main(float *Input, float *InputRef, float *Output, float lambd InputRef_y = calloc(DimTotal, sizeof(float)); /* calculate gradient field (smoothed) for the reference image */ - GradNorm_func2D(InputRef, InputRef_x, InputRef_y, eta, dimX, dimY); + GradNorm_func2D(InputRef, InputRef_x, InputRef_y, eta, (long)(dimX), (long)(dimY)); /* begin iterations */ for(ll=0; ll<iterationsNumb; ll++) { /*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, dimX, dimY); + 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, dimX, dimY); + 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, dimX, dimY); + Grad_dfunc2D(P1, P2, Output, R1, R2, InputRef_x, InputRef_y, lambdaPar, (long)(dimX), (long)(dimY)); /* projection step */ Proj_dfunc2D(P1, P2, methodTV, DimTotal); @@ -102,9 +103,9 @@ float dTV_FGP_CPU_main(float *Input, float *InputRef, float *Output, float lambd if (count > 4) break; /*storing old values*/ - copyIm(Output, Output_prev, dimX, dimY, 1); - copyIm(P1, P1_prev, dimX, dimY, 1); - copyIm(P2, P2_prev, dimX, dimY, 1); + copyIm(Output, Output_prev, (long)(dimX), (long)(dimY), 1l); + copyIm(P1, P1_prev, (long)(dimX), (long)(dimY), 1l); + copyIm(P2, P2_prev, (long)(dimX), (long)(dimY), 1l); tk = tkp1; } if (printM == 1) printf("FGP-dTV iterations stopped at iteration %i \n", ll); @@ -113,7 +114,7 @@ float dTV_FGP_CPU_main(float *Input, float *InputRef, float *Output, float lambd 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, *InputRef_x=NULL, *InputRef_y=NULL, *InputRef_z=NULL; - DimTotal = dimX*dimY*dimZ; + DimTotal = (long)(dimX*dimY*dimZ); Output_prev = calloc(DimTotal, sizeof(float)); P1 = calloc(DimTotal, sizeof(float)); @@ -130,22 +131,22 @@ float dTV_FGP_CPU_main(float *Input, float *InputRef, float *Output, float lambd InputRef_z = calloc(DimTotal, sizeof(float)); /* calculate gradient field (smoothed) for the reference volume */ - GradNorm_func3D(InputRef, InputRef_x, InputRef_y, InputRef_z, eta, dimX, dimY, dimZ); + GradNorm_func3D(InputRef, InputRef_x, InputRef_y, InputRef_z, eta, (long)(dimX), (long)(dimY), (long)(dimZ)); /* begin iterations */ for(ll=0; ll<iterationsNumb; ll++) { /*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, dimX, dimY, dimZ); + 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, dimX, dimY, dimZ); + 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, dimX, dimY, dimZ); + 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); @@ -167,10 +168,10 @@ float dTV_FGP_CPU_main(float *Input, float *InputRef, float *Output, float lambd if (count > 4) break; /*storing old values*/ - copyIm(Output, Output_prev, dimX, dimY, dimZ); - copyIm(P1, P1_prev, dimX, dimY, dimZ); - copyIm(P2, P2_prev, dimX, dimY, dimZ); - copyIm(P3, P3_prev, dimX, dimY, dimZ); + copyIm(Output, Output_prev, (long)(dimX), (long)(dimY), (long)(dimZ)); + 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 (printM == 1) printf("FGP-dTV iterations stopped at iteration %i \n", ll); @@ -184,9 +185,9 @@ float dTV_FGP_CPU_main(float *Input, float *InputRef, float *Output, float lambd /***************************2D Functions*****************************/ /********************************************************************/ -float GradNorm_func2D(float *B, float *B_x, float *B_y, float eta, int dimX, int dimY) +float GradNorm_func2D(float *B, float *B_x, float *B_y, float eta, long dimX, long dimY) { - int i,j,index; + 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<dimX; i++) { @@ -205,9 +206,9 @@ float GradNorm_func2D(float *B, float *B_x, float *B_y, float eta, int dimX, int return 1; } -float ProjectVect_func2D(float *R1, float *R2, float *B_x, float *B_y, int dimX, int dimY) +float ProjectVect_func2D(float *R1, float *R2, float *B_x, float *B_y, long dimX, long dimY) { - int i,j,index; + long i,j,index; float in_prod; #pragma omp parallel for shared(R1, R2, B_x, B_y) private(index,i,j,in_prod) for(i=0; i<dimX; i++) { @@ -220,10 +221,10 @@ float ProjectVect_func2D(float *R1, float *R2, float *B_x, float *B_y, int dimX, return 1; } -float Obj_dfunc2D(float *A, float *D, float *R1, float *R2, float lambda, int dimX, int dimY) +float Obj_dfunc2D(float *A, float *D, float *R1, float *R2, float lambda, long dimX, long dimY) { float val1, val2; - int i,j,index; + long i,j,index; #pragma omp parallel for shared(A,D,R1,R2) private(index,i,j,val1,val2) for(i=0; i<dimX; i++) { for(j=0; j<dimY; j++) { @@ -235,10 +236,10 @@ float Obj_dfunc2D(float *A, float *D, float *R1, float *R2, float lambda, int di }} return *D; } -float Grad_dfunc2D(float *P1, float *P2, float *D, float *R1, float *R2, float *B_x, float *B_y, float lambda, int dimX, int dimY) +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) { float val1, val2, multip, in_prod; - int i,j,index; + long i,j,index; multip = (1.0f/(8.0f*lambda)); #pragma omp parallel for shared(P1,P2,D,R1,R2,B_x,B_y,multip) private(i,j,index,val1,val2,in_prod) for(i=0; i<dimX; i++) { @@ -258,10 +259,10 @@ float Grad_dfunc2D(float *P1, float *P2, float *D, float *R1, float *R2, float * }} return 1; } -float Proj_dfunc2D(float *P1, float *P2, int methTV, int DimTotal) +float Proj_dfunc2D(float *P1, float *P2, int methTV, long DimTotal) { float val1, val2, denom, sq_denom; - int i; + long i; if (methTV == 0) { /* isotropic TV*/ #pragma omp parallel for shared(P1,P2) private(i,denom,sq_denom) @@ -288,9 +289,9 @@ float Proj_dfunc2D(float *P1, float *P2, int methTV, int DimTotal) } return 1; } -float Rupd_dfunc2D(float *P1, float *P1_old, float *P2, float *P2_old, float *R1, float *R2, float tkp1, float tk, int DimTotal) +float Rupd_dfunc2D(float *P1, float *P1_old, float *P2, float *P2_old, float *R1, float *R2, float tkp1, float tk, long DimTotal) { - int i; + long i; float multip; multip = ((tk-1.0f)/tkp1); #pragma omp parallel for shared(P1,P2,P1_old,P2_old,R1,R2,multip) private(i) @@ -304,9 +305,9 @@ float Rupd_dfunc2D(float *P1, float *P1_old, float *P2, float *P2_old, float *R1 /********************************************************************/ /***************************3D Functions*****************************/ /********************************************************************/ -float GradNorm_func3D(float *B, float *B_x, float *B_y, float *B_z, float eta, int dimX, int dimY, int dimZ) +float GradNorm_func3D(float *B, float *B_x, float *B_y, float *B_z, float eta, long dimX, long dimY, long dimZ) { - int i, j, k, index; + long i, j, k, index; float val1, val2, val3, gradX, gradY, gradZ, magn; #pragma omp parallel for shared(B, B_x, B_y, B_z) private(i,j,k,index,val1,val2,val3,gradX,gradY,gradZ,magn) for(i=0; i<dimX; i++) { @@ -331,9 +332,9 @@ float GradNorm_func3D(float *B, float *B_x, float *B_y, float *B_z, float eta, i return 1; } -float ProjectVect_func3D(float *R1, float *R2, float *R3, float *B_x, float *B_y, float *B_z, int dimX, int dimY, int dimZ) +float ProjectVect_func3D(float *R1, float *R2, float *R3, float *B_x, float *B_y, float *B_z, long dimX, long dimY, long dimZ) { - int i,j,k,index; + long i,j,k,index; float in_prod; #pragma omp parallel for shared(R1, R2, R3, B_x, B_y, B_z) private(index,i,j,k,in_prod) for(i=0; i<dimX; i++) { @@ -348,10 +349,10 @@ float ProjectVect_func3D(float *R1, float *R2, float *R3, float *B_x, float *B_y return 1; } -float Obj_dfunc3D(float *A, float *D, float *R1, float *R2, float *R3, float lambda, int dimX, int dimY, int dimZ) +float Obj_dfunc3D(float *A, float *D, float *R1, float *R2, float *R3, float lambda, long dimX, long dimY, long dimZ) { float val1, val2, val3; - int i,j,k,index; + long i,j,k,index; #pragma omp parallel for shared(A,D,R1,R2,R3) private(index,i,j,k,val1,val2,val3) for(i=0; i<dimX; i++) { for(j=0; j<dimY; j++) { @@ -365,10 +366,10 @@ float Obj_dfunc3D(float *A, float *D, float *R1, float *R2, float *R3, float lam }}} return *D; } -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, int dimX, int dimY, int dimZ) +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) { float val1, val2, val3, multip, in_prod; - int i,j,k, index; + long i,j,k, index; multip = (1.0f/(26.0f*lambda)); #pragma omp parallel for shared(P1,P2,P3,D,R1,R2,R3,multip) private(index,i,j,k,val1,val2,val3,in_prod) for(i=0; i<dimX; i++) { @@ -391,10 +392,10 @@ float Grad_dfunc3D(float *P1, float *P2, float *P3, float *D, float *R1, float * }}} return 1; } -float Proj_dfunc3D(float *P1, float *P2, float *P3, int methTV, int DimTotal) +float Proj_dfunc3D(float *P1, float *P2, float *P3, int methTV, long DimTotal) { float val1, val2, val3, denom, sq_denom; - int i; + long i; if (methTV == 0) { /* isotropic TV*/ #pragma omp parallel for shared(P1,P2,P3) private(i,val1,val2,val3,sq_denom) @@ -425,9 +426,9 @@ float Proj_dfunc3D(float *P1, float *P2, float *P3, int methTV, int DimTotal) } 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, int DimTotal) +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) { - int i; + long 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(i) diff --git a/Core/regularisers_CPU/FGP_dTV_core.h b/Core/regularisers_CPU/FGP_dTV_core.h index 8d721d5..442dd30 100644 --- a/Core/regularisers_CPU/FGP_dTV_core.h +++ b/Core/regularisers_CPU/FGP_dTV_core.h @@ -54,19 +54,19 @@ extern "C" { #endif CCPI_EXPORT float dTV_FGP_CPU_main(float *Input, float *InputRef, float *Output, float lambdaPar, int iterationsNumb, float epsil, float eta, int methodTV, int nonneg, int printM, int dimX, int dimY, int dimZ); -CCPI_EXPORT float GradNorm_func2D(float *B, float *B_x, float *B_y, float eta, int dimX, int dimY); -CCPI_EXPORT float ProjectVect_func2D(float *R1, float *R2, float *B_x, float *B_y, int dimX, int dimY); -CCPI_EXPORT float Obj_dfunc2D(float *A, float *D, float *R1, float *R2, float lambda, int dimX, int dimY); -CCPI_EXPORT float Grad_dfunc2D(float *P1, float *P2, float *D, float *R1, float *R2, float *B_x, float *B_y, float lambda, int dimX, int dimY); -CCPI_EXPORT float Proj_dfunc2D(float *P1, float *P2, int methTV, int DimTotal); -CCPI_EXPORT float Rupd_dfunc2D(float *P1, float *P1_old, float *P2, float *P2_old, float *R1, float *R2, float tkp1, float tk, int DimTotal); +CCPI_EXPORT float GradNorm_func2D(float *B, float *B_x, float *B_y, float eta, long dimX, long dimY); +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, int dimX, int dimY, int dimZ); -CCPI_EXPORT float ProjectVect_func3D(float *R1, float *R2, float *R3, float *B_x, float *B_y, float *B_z, int dimX, int dimY, int dimZ); -CCPI_EXPORT float Obj_dfunc3D(float *A, float *D, float *R1, float *R2, float *R3, float lambda, int dimX, int dimY, int 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, int dimX, int dimY, int dimZ); -CCPI_EXPORT float Proj_dfunc3D(float *P1, float *P2, float *P3, int methTV, int 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, int 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 } #endif diff --git a/Core/regularisers_CPU/LLT_ROF_core.c b/Core/regularisers_CPU/LLT_ROF_core.c index 1584a29..8416a14 100644 --- a/Core/regularisers_CPU/LLT_ROF_core.c +++ b/Core/regularisers_CPU/LLT_ROF_core.c @@ -51,10 +51,11 @@ int signLLT(float x) { float LLT_ROF_CPU_main(float *Input, float *Output, float lambdaROF, float lambdaLLT, int iterationsNumb, float tau, int dimX, int dimY, int dimZ) { - int DimTotal, ll; + long DimTotal; + int ll; float *D1_LLT=NULL, *D2_LLT=NULL, *D3_LLT=NULL, *D1_ROF=NULL, *D2_ROF=NULL, *D3_ROF=NULL; - DimTotal = dimX*dimY*dimZ; + DimTotal = (long)(dimX*dimY*dimZ); D1_ROF = calloc(DimTotal, sizeof(float)); D2_ROF = calloc(DimTotal, sizeof(float)); @@ -64,32 +65,32 @@ float LLT_ROF_CPU_main(float *Input, float *Output, float lambdaROF, float lambd D2_LLT = calloc(DimTotal, sizeof(float)); D3_LLT = calloc(DimTotal, sizeof(float)); - copyIm(Input, Output, dimX, dimY, dimZ); /* initialize */ + copyIm(Input, Output, (long)(dimX), (long)(dimY), (long)(dimZ)); /* initialize */ for(ll = 0; ll < iterationsNumb; ll++) { if (dimZ == 1) { /* 2D case */ /****************ROF******************/ /* calculate first-order differences */ - D1_func_ROF(Output, D1_ROF, dimX, dimY, dimZ); - D2_func_ROF(Output, D2_ROF, dimX, dimY, dimZ); + D1_func_ROF(Output, D1_ROF, (long)(dimX), (long)(dimY), 1l); + D2_func_ROF(Output, D2_ROF, (long)(dimX), (long)(dimY), 1l); /****************LLT******************/ /* estimate second-order derrivatives */ - der2D_LLT(Output, D1_LLT, D2_LLT, dimX, dimY, dimZ); + der2D_LLT(Output, D1_LLT, D2_LLT, (long)(dimX), (long)(dimY), 1l); /* Joint update for ROF and LLT models */ - Update2D_LLT_ROF(Input, Output, D1_LLT, D2_LLT, D1_ROF, D2_ROF, lambdaROF, lambdaLLT, tau, dimX, dimY, dimZ); + Update2D_LLT_ROF(Input, Output, D1_LLT, D2_LLT, D1_ROF, D2_ROF, lambdaROF, lambdaLLT, tau, (long)(dimX), (long)(dimY), 1l); } else { /* 3D case */ /* calculate first-order differences */ - D1_func_ROF(Output, D1_ROF, dimX, dimY, dimZ); - D2_func_ROF(Output, D2_ROF, dimX, dimY, dimZ); - D3_func_ROF(Output, D3_ROF, dimX, dimY, dimZ); + D1_func_ROF(Output, D1_ROF, (long)(dimX), (long)(dimY), (long)(dimZ)); + D2_func_ROF(Output, D2_ROF, (long)(dimX), (long)(dimY), (long)(dimZ)); + D3_func_ROF(Output, D3_ROF, (long)(dimX), (long)(dimY), (long)(dimZ)); /****************LLT******************/ /* estimate second-order derrivatives */ - der3D_LLT(Output, D1_LLT, D2_LLT, D3_LLT, dimX, dimY, dimZ); + der3D_LLT(Output, D1_LLT, D2_LLT, D3_LLT,(long)(dimX), (long)(dimY), (long)(dimZ)); /* Joint update for ROF and LLT models */ - Update3D_LLT_ROF(Input, Output, D1_LLT, D2_LLT, D3_LLT, D1_ROF, D2_ROF, D3_ROF, lambdaROF, lambdaLLT, tau, dimX, dimY, dimZ); + Update3D_LLT_ROF(Input, Output, D1_LLT, D2_LLT, D3_LLT, D1_ROF, D2_ROF, D3_ROF, lambdaROF, lambdaLLT, tau, (long)(dimX), (long)(dimY), (long)(dimZ)); } } /*end of iterations*/ free(D1_LLT);free(D2_LLT);free(D3_LLT); @@ -100,9 +101,9 @@ float LLT_ROF_CPU_main(float *Input, float *Output, float lambdaROF, float lambd /*************************************************************************/ /**********************LLT-related functions *****************************/ /*************************************************************************/ -float der2D_LLT(float *U, float *D1, float *D2, int dimX, int dimY, int dimZ) +float der2D_LLT(float *U, float *D1, float *D2, long dimX, long dimY, long dimZ) { - int i, j, index, i_p, i_m, j_m, j_p; + 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<dimX; i++) { @@ -127,9 +128,9 @@ float der2D_LLT(float *U, float *D1, float *D2, int dimX, int dimY, int dimZ) return 1; } -float der3D_LLT(float *U, float *D1, float *D2, float *D3, int dimX, int dimY, int dimZ) +float der3D_LLT(float *U, float *D1, float *D2, float *D3, long dimX, long dimY, long dimZ) { - int i, j, k, i_p, i_m, j_m, j_p, k_p, k_m, index; + long i, j, k, i_p, i_m, j_m, j_p, k_p, k_m, index; float dxx, dyy, dzz, denom_xx, denom_yy, denom_zz; #pragma omp parallel for shared(U,D1,D2,D3) private(i, j, index, k, i_p, i_m, j_m, j_p, k_p, k_m, denom_xx, denom_yy, denom_zz, dxx, dyy, dzz) for (i = 0; i<dimX; i++) { @@ -167,10 +168,10 @@ float der3D_LLT(float *U, float *D1, float *D2, float *D3, int dimX, int dimY, i /*************************************************************************/ /* calculate differences 1 */ -float D1_func_ROF(float *A, float *D1, int dimX, int dimY, int dimZ) +float D1_func_ROF(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; - int i,j,k,i1,i2,k1,j1,j2,k2,index; + 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) @@ -232,10 +233,10 @@ float D1_func_ROF(float *A, float *D1, int dimX, int dimY, int dimZ) return *D1; } /* calculate differences 2 */ -float D2_func_ROF(float *A, float *D2, int dimX, int dimY, int dimZ) +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; - int i,j,k,i1,i2,k1,j1,j2,k2,index; + 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) @@ -297,10 +298,10 @@ float D2_func_ROF(float *A, float *D2, int dimX, int dimY, int dimZ) } /* calculate differences 3 */ -float D3_func_ROF(float *A, float *D3, int dimY, int dimX, int dimZ) +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; - int index,i,j,k,i1,i2,k1,j1,j2,k2; + 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<dimY; j++) { @@ -338,9 +339,9 @@ float D3_func_ROF(float *A, float *D3, int dimY, int dimX, int dimZ) /**********************ROF-LLT-related functions *************************/ /*************************************************************************/ -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, int dimX, int dimY, int 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) { - int i, j, index, i_p, i_m, j_m, j_p; + 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<dimX; i++) { @@ -369,9 +370,9 @@ float Update2D_LLT_ROF(float *U0, float *U, float *D1_LLT, float *D2_LLT, float return *U; } -float Update3D_LLT_ROF(float *U0, float *U, float *D1_LLT, float *D2_LLT, float *D3_LLT, float *D1_ROF, float *D2_ROF, float *D3_ROF, float lambdaROF, float lambdaLLT, float tau, int dimX, int dimY, int dimZ) +float Update3D_LLT_ROF(float *U0, float *U, float *D1_LLT, float *D2_LLT, float *D3_LLT, float *D1_ROF, float *D2_ROF, float *D3_ROF, float lambdaROF, float lambdaLLT, float tau, long dimX, long dimY, long dimZ) { - int i, j, k, i_p, i_m, j_m, j_p, k_p, k_m, index; + long i, j, k, i_p, i_m, j_m, j_p, k_p, k_m, index; float div, laplc, dxx, dyy, dzz, dv1, dv2, dv3; #pragma omp parallel for shared(U,U0) private(i, j, k, index, i_p, i_m, j_m, j_p, k_p, k_m, laplc, div, dxx, dyy, dzz, dv1, dv2, dv3) for (i = 0; i<dimX; i++) { diff --git a/Core/regularisers_CPU/LLT_ROF_core.h b/Core/regularisers_CPU/LLT_ROF_core.h index 65d25cd..8e6591e 100644 --- a/Core/regularisers_CPU/LLT_ROF_core.h +++ b/Core/regularisers_CPU/LLT_ROF_core.h @@ -51,15 +51,15 @@ extern "C" { #endif CCPI_EXPORT float LLT_ROF_CPU_main(float *Input, float *Output, float lambdaROF, float lambdaLLT, int iterationsNumb, float tau, int dimX, int dimY, int dimZ); -CCPI_EXPORT float der2D_LLT(float *U, float *D1, float *D2, int dimX, int dimY, int dimZ); -CCPI_EXPORT float der3D_LLT(float *U, float *D1, float *D2, float *D3, int dimX, int dimY, int dimZ); +CCPI_EXPORT float der2D_LLT(float *U, float *D1, float *D2, long dimX, long dimY, long dimZ); +CCPI_EXPORT float der3D_LLT(float *U, float *D1, float *D2, float *D3, long dimX, long dimY, long dimZ); -CCPI_EXPORT float D1_func_ROF(float *A, float *D1, int dimY, int dimX, int dimZ); -CCPI_EXPORT float D2_func_ROF(float *A, float *D2, int dimY, int dimX, int dimZ); -CCPI_EXPORT float D3_func_ROF(float *A, float *D3, int dimY, int dimX, int dimZ); +CCPI_EXPORT float D1_func_ROF(float *A, float *D1, long dimX, long dimY, long dimZ); +CCPI_EXPORT float D2_func_ROF(float *A, float *D2, long dimX, long dimY, long dimZ); +CCPI_EXPORT float D3_func_ROF(float *A, float *D3, long dimX, long dimY, long dimZ); -CCPI_EXPORT 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, int dimX, int dimY, int dimZ); -CCPI_EXPORT float Update3D_LLT_ROF(float *U0, float *U, float *D1_LLT, float *D2_LLT, float *D3_LLT, float *D1_ROF, float *D2_ROF, float *D3_ROF, float lambdaROF, float lambdaLLT, float tau, int dimX, int dimY, int dimZ); +CCPI_EXPORT 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); +CCPI_EXPORT float Update3D_LLT_ROF(float *U0, float *U, float *D1_LLT, float *D2_LLT, float *D3_LLT, float *D1_ROF, float *D2_ROF, float *D3_ROF, float lambdaROF, float lambdaLLT, float tau, long dimX, long dimY, long dimZ); #ifdef __cplusplus } #endif diff --git a/Core/regularisers_CPU/ROF_TV_core.c b/Core/regularisers_CPU/ROF_TV_core.c index fb3bc7c..1858442 100644 --- a/Core/regularisers_CPU/ROF_TV_core.c +++ b/Core/regularisers_CPU/ROF_TV_core.c @@ -49,34 +49,34 @@ int sign(float x) { float TV_ROF_CPU_main(float *Input, float *Output, float lambdaPar, int iterationsNumb, float tau, int dimX, int dimY, int dimZ) { float *D1, *D2, *D3; - int i, DimTotal; - DimTotal = dimX*dimY*dimZ; + int i; + long DimTotal; + DimTotal = (long)(dimX*dimY*dimZ); D1 = calloc(DimTotal, sizeof(float)); D2 = calloc(DimTotal, sizeof(float)); D3 = calloc(DimTotal, sizeof(float)); /* copy into output */ - copyIm(Input, Output, dimX, dimY, dimZ); + copyIm(Input, Output, (long)(dimX), (long)(dimY), (long)(dimZ)); /* start TV iterations */ - for(i=0; i < iterationsNumb; i++) { - + for(i=0; i < iterationsNumb; i++) { /* calculate differences */ - D1_func(Output, D1, dimX, dimY, dimZ); - D2_func(Output, D2, dimX, dimY, dimZ); - if (dimZ > 1) D3_func(Output, D3, dimX, dimY, dimZ); - TV_kernel(D1, D2, D3, Output, Input, lambdaPar, tau, dimX, dimY, dimZ); + 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, tau, (long)(dimX), (long)(dimY), (long)(dimZ)); } free(D1);free(D2); free(D3); return *Output; } /* calculate differences 1 */ -float D1_func(float *A, float *D1, int dimX, int dimY, int dimZ) +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; - int i,j,k,i1,i2,k1,j1,j2,k2,index; + 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) @@ -138,10 +138,10 @@ float D1_func(float *A, float *D1, int dimX, int dimY, int dimZ) return *D1; } /* calculate differences 2 */ -float D2_func(float *A, float *D2, int dimX, int dimY, int dimZ) +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; - int i,j,k,i1,i2,k1,j1,j2,k2,index; + 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) @@ -155,8 +155,7 @@ float D2_func(float *A, float *D2, int dimX, int dimY, int dimZ) j1 = j + 1; if (j1 >= dimY) j1 = j-1; 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; - + k2 = k - 1; if (k2 < 0) k2 = k+1; /* Forward-backward differences */ NOMx_1 = A[(dimX*dimY)*k + (j1)*dimX + i] - A[index]; /* x+ */ @@ -203,10 +202,10 @@ float D2_func(float *A, float *D2, int dimX, int dimY, int dimZ) } /* calculate differences 3 */ -float D3_func(float *A, float *D3, int dimY, int dimX, int dimZ) +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; - int index,i,j,k,i1,i2,k1,j1,j2,k2; + 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<dimY; j++) { @@ -241,10 +240,10 @@ float D3_func(float *A, float *D3, int dimY, int dimX, int dimZ) } /* calculate divergence */ -float TV_kernel(float *D1, float *D2, float *D3, float *B, float *A, float lambda, float tau, int dimX, int dimY, int dimZ) +float TV_kernel(float *D1, float *D2, float *D3, float *B, float *A, float lambda, float tau, long dimX, long dimY, long dimZ) { float dv1, dv2, dv3; - int index,i,j,k,i1,i2,k1,j1,j2,k2; + 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) diff --git a/Core/regularisers_CPU/ROF_TV_core.h b/Core/regularisers_CPU/ROF_TV_core.h index 14daf58..4e320e9 100644 --- a/Core/regularisers_CPU/ROF_TV_core.h +++ b/Core/regularisers_CPU/ROF_TV_core.h @@ -46,11 +46,12 @@ limitations under the License. #ifdef __cplusplus extern "C" { #endif -CCPI_EXPORT float TV_kernel(float *D1, float *D2, float *D3, float *B, float *A, float lambda, float tau, int dimY, int dimX, int dimZ); CCPI_EXPORT float TV_ROF_CPU_main(float *Input, float *Output, float lambdaPar, int iterationsNumb, float tau, int dimX, int dimY, int dimZ); -CCPI_EXPORT float D1_func(float *A, float *D1, int dimY, int dimX, int dimZ); -CCPI_EXPORT float D2_func(float *A, float *D2, int dimY, int dimX, int dimZ); -CCPI_EXPORT float D3_func(float *A, float *D3, int dimY, int dimX, int dimZ); + +CCPI_EXPORT float TV_kernel(float *D1, float *D2, float *D3, float *B, float *A, float lambda, float tau, long dimX, long dimY, long dimZ); +CCPI_EXPORT float D1_func(float *A, float *D1, long dimX, long dimY, long dimZ); +CCPI_EXPORT float D2_func(float *A, float *D2, long dimX, long dimY, long dimZ); +CCPI_EXPORT float D3_func(float *A, float *D3, long dimX, long dimY, long dimZ); #ifdef __cplusplus } -#endif +#endif
\ No newline at end of file diff --git a/Core/regularisers_CPU/SB_TV_core.c b/Core/regularisers_CPU/SB_TV_core.c index 06d3de3..769ea67 100755 --- a/Core/regularisers_CPU/SB_TV_core.c +++ b/Core/regularisers_CPU/SB_TV_core.c @@ -37,16 +37,17 @@ limitations under the License. float SB_TV_CPU_main(float *Input, float *Output, float mu, int iter, float epsil, int methodTV, int printM, int dimX, int dimY, int dimZ) { - int ll, j, DimTotal; + int ll; + long j, DimTotal; float re, re1, lambda; int count = 0; mu = 1.0f/mu; - lambda = 2.0f*mu; + lambda = 2.0f*mu; if (dimZ <= 1) { /* 2D case */ float *Output_prev=NULL, *Dx=NULL, *Dy=NULL, *Bx=NULL, *By=NULL; - DimTotal = dimX*dimY; + DimTotal = (long)(dimX*dimY); Output_prev = calloc(DimTotal, sizeof(float)); Dx = calloc(DimTotal, sizeof(float)); @@ -54,26 +55,26 @@ float SB_TV_CPU_main(float *Input, float *Output, float mu, int iter, float epsi Bx = calloc(DimTotal, sizeof(float)); By = calloc(DimTotal, sizeof(float)); - copyIm(Input, Output, dimX, dimY, 1); /*initialize */ + copyIm(Input, Output, (long)(dimX), (long)(dimY), 1l); /*initialize */ /* begin outer SB iterations */ for(ll=0; ll<iter; ll++) { /* storing old estimate */ - copyIm(Output, Output_prev, dimX, dimY, 1); + copyIm(Output, Output_prev, (long)(dimX), (long)(dimY), 1l); /* perform two GS iterations (normally 2 is enough for the convergence) */ - gauss_seidel2D(Output, Input, Output_prev, Dx, Dy, Bx, By, dimX, dimY, lambda, mu); - copyIm(Output, Output_prev, dimX, dimY, 1); + gauss_seidel2D(Output, Input, Output_prev, Dx, Dy, Bx, By, (long)(dimX), (long)(dimY), lambda, mu); + copyIm(Output, Output_prev, (long)(dimX), (long)(dimY), 1l); /*GS iteration */ - gauss_seidel2D(Output, Input, Output_prev, Dx, Dy, Bx, By, dimX, dimY, lambda, mu); + gauss_seidel2D(Output, Input, Output_prev, Dx, Dy, Bx, By, (long)(dimX), (long)(dimY), lambda, mu); /* TV-related step */ - if (methodTV == 1) updDxDy_shrinkAniso2D(Output, Dx, Dy, Bx, By, dimX, dimY, lambda); - else updDxDy_shrinkIso2D(Output, Dx, Dy, Bx, By, dimX, dimY, lambda); + if (methodTV == 1) updDxDy_shrinkAniso2D(Output, Dx, Dy, Bx, By, (long)(dimX), (long)(dimY), lambda); + else updDxDy_shrinkIso2D(Output, Dx, Dy, Bx, By, (long)(dimX), (long)(dimY), lambda); /* update for Bregman variables */ - updBxBy2D(Output, Dx, Dy, Bx, By, dimX, dimY); + updBxBy2D(Output, Dx, Dy, Bx, By, (long)(dimX), (long)(dimY)); /* check early stopping criteria if epsilon not equal zero */ if (epsil != 0) { @@ -94,7 +95,7 @@ float SB_TV_CPU_main(float *Input, float *Output, float mu, int iter, float epsi else { /* 3D case */ float *Output_prev=NULL, *Dx=NULL, *Dy=NULL, *Dz=NULL, *Bx=NULL, *By=NULL, *Bz=NULL; - DimTotal = dimX*dimY*dimZ; + DimTotal = (long)(dimX*dimY*dimZ); Output_prev = calloc(DimTotal, sizeof(float)); Dx = calloc(DimTotal, sizeof(float)); @@ -104,26 +105,26 @@ float SB_TV_CPU_main(float *Input, float *Output, float mu, int iter, float epsi By = calloc(DimTotal, sizeof(float)); Bz = calloc(DimTotal, sizeof(float)); - copyIm(Input, Output, dimX, dimY, dimZ); /*initialize */ + copyIm(Input, Output, (long)(dimX), (long)(dimY), (long)(dimZ)); /*initialize */ /* begin outer SB iterations */ for(ll=0; ll<iter; ll++) { /* storing old estimate */ - copyIm(Output, Output_prev, dimX, dimY, dimZ); + copyIm(Output, Output_prev, (long)(dimX), (long)(dimY), (long)(dimZ)); /* perform two GS iterations (normally 2 is enough for the convergence) */ - gauss_seidel3D(Output, Input, Output_prev, Dx, Dy, Dz, Bx, By, Bz, dimX, dimY, dimZ, lambda, mu); - copyIm(Output, Output_prev, dimX, dimY, dimZ); + gauss_seidel3D(Output, Input, Output_prev, Dx, Dy, Dz, Bx, By, Bz, (long)(dimX), (long)(dimY), (long)(dimZ), lambda, mu); + copyIm(Output, Output_prev, (long)(dimX), (long)(dimY), (long)(dimZ)); /*GS iteration */ - gauss_seidel3D(Output, Input, Output_prev, Dx, Dy, Dz, Bx, By, Bz, dimX, dimY, dimZ, lambda, mu); + gauss_seidel3D(Output, Input, Output_prev, Dx, Dy, Dz, Bx, By, Bz, (long)(dimX), (long)(dimY), (long)(dimZ), lambda, mu); /* TV-related step */ - if (methodTV == 1) updDxDyDz_shrinkAniso3D(Output, Dx, Dy, Dz, Bx, By, Bz, dimX, dimY, dimZ, lambda); - else updDxDyDz_shrinkIso3D(Output, Dx, Dy, Dz, Bx, By, Bz, dimX, dimY, dimZ, lambda); + if (methodTV == 1) updDxDyDz_shrinkAniso3D(Output, Dx, Dy, Dz, Bx, By, Bz, (long)(dimX), (long)(dimY), (long)(dimZ), lambda); + else updDxDyDz_shrinkIso3D(Output, Dx, Dy, Dz, Bx, By, Bz, (long)(dimX), (long)(dimY), (long)(dimZ), lambda); /* update for Bregman variables */ - updBxByBz3D(Output, Dx, Dy, Dz, Bx, By, Bz, dimX, dimY, dimZ); + updBxByBz3D(Output, Dx, Dy, Dz, Bx, By, Bz, (long)(dimX), (long)(dimY), (long)(dimZ)); /* check early stopping criteria if epsilon not equal zero */ if (epsil != 0) { @@ -147,10 +148,10 @@ float SB_TV_CPU_main(float *Input, float *Output, float mu, int iter, float epsi /********************************************************************/ /***************************2D Functions*****************************/ /********************************************************************/ -float gauss_seidel2D(float *U, float *A, float *U_prev, float *Dx, float *Dy, float *Bx, float *By, int dimX, int dimY, float lambda, float mu) +float gauss_seidel2D(float *U, float *A, float *U_prev, float *Dx, float *Dy, float *Bx, float *By, long dimX, long dimY, float lambda, float mu) { float sum, normConst; - int i,j,i1,i2,j1,j2,index; + 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) @@ -173,9 +174,9 @@ float gauss_seidel2D(float *U, float *A, float *U_prev, float *Dx, float *Dy, fl return *U; } -float updDxDy_shrinkAniso2D(float *U, float *Dx, float *Dy, float *Bx, float *By, int dimX, int dimY, float lambda) +float updDxDy_shrinkAniso2D(float *U, float *Dx, float *Dy, float *Bx, float *By, long dimX, long dimY, float lambda) { - int i,j,i1,j1,index; + long i,j,i1,j1,index; float val1, val11, val2, val22, denom_lam; denom_lam = 1.0f/lambda; #pragma omp parallel for shared(U,denom_lam) private(index,i,j,i1,j1,val1,val11,val2,val22) @@ -198,9 +199,9 @@ float updDxDy_shrinkAniso2D(float *U, float *Dx, float *Dy, float *Bx, float *By }} return 1; } -float updDxDy_shrinkIso2D(float *U, float *Dx, float *Dy, float *Bx, float *By, int dimX, int dimY, float lambda) +float updDxDy_shrinkIso2D(float *U, float *Dx, float *Dy, float *Bx, float *By, long dimX, long dimY, float lambda) { - int i,j,i1,j1,index; + long i,j,i1,j1,index; float val1, val11, val2, denom, denom_lam; denom_lam = 1.0f/lambda; @@ -230,9 +231,9 @@ float updDxDy_shrinkIso2D(float *U, float *Dx, float *Dy, float *Bx, float *By, }} return 1; } -float updBxBy2D(float *U, float *Dx, float *Dy, float *Bx, float *By, int dimX, int dimY) +float updBxBy2D(float *U, float *Dx, float *Dy, float *Bx, float *By, long dimX, long dimY) { - int i,j,i1,j1,index; + long i,j,i1,j1,index; #pragma omp parallel for shared(U) private(index,i,j,i1,j1) for(i=0; i<dimX; i++) { for(j=0; j<dimY; j++) { @@ -251,10 +252,10 @@ float updBxBy2D(float *U, float *Dx, float *Dy, float *Bx, float *By, int dimX, /***************************3D Functions*****************************/ /********************************************************************/ /*****************************************************************/ -float gauss_seidel3D(float *U, float *A, float *U_prev, float *Dx, float *Dy, float *Dz, float *Bx, float *By, float *Bz, int dimX, int dimY, int dimZ, float lambda, float mu) +float gauss_seidel3D(float *U, float *A, float *U_prev, float *Dx, float *Dy, float *Dz, float *Bx, float *By, float *Bz, long dimX, long dimY, long dimZ, float lambda, float mu) { float normConst, d_val, b_val, sum; - int i,j,i1,i2,j1,j2,k,k1,k2,index; + long i,j,i1,i2,j1,j2,k,k1,k2,index; normConst = 1.0f/(mu + 6.0f*lambda); #pragma omp parallel for shared(U) private(index,i,j,i1,i2,j1,j2,k,k1,k2,d_val,b_val,sum) for(i=0; i<dimX; i++) { @@ -280,9 +281,9 @@ float gauss_seidel3D(float *U, float *A, float *U_prev, float *Dx, float *Dy, fl return *U; } -float updDxDyDz_shrinkAniso3D(float *U, float *Dx, float *Dy, float *Dz, float *Bx, float *By, float *Bz, int dimX, int dimY, int dimZ, float lambda) +float updDxDyDz_shrinkAniso3D(float *U, float *Dx, float *Dy, float *Dz, float *Bx, float *By, float *Bz, long dimX, long dimY, long dimZ, float lambda) { - int i,j,i1,j1,k,k1,index; + long i,j,i1,j1,k,k1,index; float val1, val11, val2, val22, val3, val33, denom_lam; denom_lam = 1.0f/lambda; #pragma omp parallel for shared(U,denom_lam) private(index,i,j,i1,j1,k,k1,val1,val11,val2,val22,val3,val33) @@ -310,9 +311,9 @@ float updDxDyDz_shrinkAniso3D(float *U, float *Dx, float *Dy, float *Dz, float * }}} return 1; } -float updDxDyDz_shrinkIso3D(float *U, float *Dx, float *Dy, float *Dz, float *Bx, float *By, float *Bz, int dimX, int dimY, int dimZ, float lambda) +float updDxDyDz_shrinkIso3D(float *U, float *Dx, float *Dy, float *Dz, float *Bx, float *By, float *Bz, long dimX, long dimY, long dimZ, float lambda) { - int i,j,i1,j1,k,k1,index; + long i,j,i1,j1,k,k1,index; float val1, val11, val2, val3, denom, denom_lam; denom_lam = 1.0f/lambda; #pragma omp parallel for shared(U,denom_lam) private(index,denom,i,j,i1,j1,k,k1,val1,val11,val2,val3) @@ -346,9 +347,9 @@ float updDxDyDz_shrinkIso3D(float *U, float *Dx, float *Dy, float *Dz, float *Bx }}} return 1; } -float updBxByBz3D(float *U, float *Dx, float *Dy, float *Dz, float *Bx, float *By, float *Bz, int dimX, int dimY, int dimZ) +float updBxByBz3D(float *U, float *Dx, float *Dy, float *Dz, float *Bx, float *By, float *Bz, long dimX, long dimY, long dimZ) { - int i,j,k,i1,j1,k1,index; + long i,j,k,i1,j1,k1,index; #pragma omp parallel for shared(U) private(index,i,j,k,i1,j1,k1) for(i=0; i<dimX; i++) { for(j=0; j<dimY; j++) { diff --git a/Core/regularisers_CPU/SB_TV_core.h b/Core/regularisers_CPU/SB_TV_core.h index ad4f529..7485e3b 100644 --- a/Core/regularisers_CPU/SB_TV_core.h +++ b/Core/regularisers_CPU/SB_TV_core.h @@ -47,15 +47,15 @@ extern "C" { #endif CCPI_EXPORT float SB_TV_CPU_main(float *Input, float *Output, float mu, int iter, float epsil, int methodTV, int printM, int dimX, int dimY, int dimZ); -CCPI_EXPORT float gauss_seidel2D(float *U, float *A, float *U_prev, float *Dx, float *Dy, float *Bx, float *By, int dimX, int dimY, float lambda, float mu); -CCPI_EXPORT float updDxDy_shrinkAniso2D(float *U, float *Dx, float *Dy, float *Bx, float *By, int dimX, int dimY, float lambda); -CCPI_EXPORT float updDxDy_shrinkIso2D(float *U, float *Dx, float *Dy, float *Bx, float *By, int dimX, int dimY, float lambda); -CCPI_EXPORT float updBxBy2D(float *U, float *Dx, float *Dy, float *Bx, float *By, int dimX, int dimY); - -CCPI_EXPORT float gauss_seidel3D(float *U, float *A, float *U_prev, float *Dx, float *Dy, float *Dz, float *Bx, float *By, float *Bz, int dimX, int dimY, int dimZ, float lambda, float mu); -CCPI_EXPORT float updDxDyDz_shrinkAniso3D(float *U, float *Dx, float *Dy, float *Dz, float *Bx, float *By, float *Bz, int dimX, int dimY, int dimZ, float lambda); -CCPI_EXPORT float updDxDyDz_shrinkIso3D(float *U, float *Dx, float *Dy, float *Dz, float *Bx, float *By, float *Bz, int dimX, int dimY, int dimZ, float lambda); -CCPI_EXPORT float updBxByBz3D(float *U, float *Dx, float *Dy, float *Dz, float *Bx, float *By, float *Bz, int dimX, int dimY, int dimZ); +CCPI_EXPORT float gauss_seidel2D(float *U, float *A, float *U_prev, float *Dx, float *Dy, float *Bx, float *By, long dimX, long dimY, float lambda, float mu); +CCPI_EXPORT float updDxDy_shrinkAniso2D(float *U, float *Dx, float *Dy, float *Bx, float *By, long dimX, long dimY, float lambda); +CCPI_EXPORT float updDxDy_shrinkIso2D(float *U, float *Dx, float *Dy, float *Bx, float *By, long dimX, long dimY, float lambda); +CCPI_EXPORT float updBxBy2D(float *U, float *Dx, float *Dy, float *Bx, float *By, long dimX, long dimY); + +CCPI_EXPORT float gauss_seidel3D(float *U, float *A, float *U_prev, float *Dx, float *Dy, float *Dz, float *Bx, float *By, float *Bz, long dimX, long dimY, long dimZ, float lambda, float mu); +CCPI_EXPORT float updDxDyDz_shrinkAniso3D(float *U, float *Dx, float *Dy, float *Dz, float *Bx, float *By, float *Bz, long dimX, long dimY, long dimZ, float lambda); +CCPI_EXPORT float updDxDyDz_shrinkIso3D(float *U, float *Dx, float *Dy, float *Dz, float *Bx, float *By, float *Bz, long dimX, long dimY, long dimZ, float lambda); +CCPI_EXPORT float updBxByBz3D(float *U, float *Dx, float *Dy, float *Dz, float *Bx, float *By, float *Bz, long dimX, long dimY, long dimZ); #ifdef __cplusplus } #endif diff --git a/Core/regularisers_CPU/TGV_core.c b/Core/regularisers_CPU/TGV_core.c index db88614..edd77b2 100644 --- a/Core/regularisers_CPU/TGV_core.c +++ b/Core/regularisers_CPU/TGV_core.c @@ -39,10 +39,11 @@ limitations under the License. float TGV_main(float *U0, float *U, float lambda, float alpha1, float alpha0, int iter, float L2, int dimX, int dimY) { - int DimTotal, ll; + long DimTotal; + int ll; float *U_old, *P1, *P2, *Q1, *Q2, *Q3, *V1, *V1_old, *V2, *V2_old, tau, sigma; - DimTotal = dimX*dimY; + DimTotal = (long)(dimX*dimY); /* dual variables */ P1 = calloc(DimTotal, sizeof(float)); @@ -59,7 +60,7 @@ float TGV_main(float *U0, float *U, float lambda, float alpha1, float alpha0, in V2 = calloc(DimTotal, sizeof(float)); V2_old = calloc(DimTotal, sizeof(float)); - copyIm(U0, U, dimX, dimY, 1); /* initialize */ + copyIm(U0, U, (long)(dimX), (long)(dimY), 1l); /* initialize */ tau = pow(L2,-0.5); sigma = pow(L2,-0.5); @@ -68,36 +69,36 @@ float TGV_main(float *U0, float *U, float lambda, float alpha1, float alpha0, in for(ll = 0; ll < iter; ll++) { /* Calculate Dual Variable P */ - DualP_2D(U, V1, V2, P1, P2, dimX, dimY, sigma); + DualP_2D(U, V1, V2, P1, P2, (long)(dimX), (long)(dimY), sigma); /*Projection onto convex set for P*/ - ProjP_2D(P1, P2, dimX, dimY, alpha1); + ProjP_2D(P1, P2, (long)(dimX), (long)(dimY), alpha1); /* Calculate Dual Variable Q */ - DualQ_2D(V1, V2, Q1, Q2, Q3, dimX, dimY, sigma); + DualQ_2D(V1, V2, Q1, Q2, Q3, (long)(dimX), (long)(dimY), sigma); /*Projection onto convex set for Q*/ - ProjQ_2D(Q1, Q2, Q3, dimX, dimY, alpha0); + ProjQ_2D(Q1, Q2, Q3, (long)(dimX), (long)(dimY), alpha0); /*saving U into U_old*/ - copyIm(U, U_old, dimX, dimY, 1); + copyIm(U, U_old, (long)(dimX), (long)(dimY), 1l); /*adjoint operation -> divergence and projection of P*/ - DivProjP_2D(U, U0, P1, P2, dimX, dimY, lambda, tau); + DivProjP_2D(U, U0, P1, P2, (long)(dimX), (long)(dimY), lambda, tau); /*get updated solution U*/ - newU(U, U_old, dimX, dimY); + newU(U, U_old, (long)(dimX), (long)(dimY)); /*saving V into V_old*/ - copyIm(V1, V1_old, dimX, dimY, 1); - copyIm(V2, V2_old, dimX, dimY, 1); + 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, dimX, dimY, tau); + UpdV_2D(V1, V2, P1, P2, Q1, Q2, Q3, (long)(dimX), (long)(dimY), tau); /*get new V*/ - newU(V1, V1_old, dimX, dimY); - newU(V2, V2_old, dimX, dimY); + newU(V1, V1_old, (long)(dimX), (long)(dimY)); + newU(V2, V2_old, (long)(dimX), (long)(dimY)); } /*end of iterations*/ return *U; } @@ -107,9 +108,9 @@ float TGV_main(float *U0, float *U, float lambda, float alpha1, float alpha0, in /********************************************************************/ /*Calculating dual variable P (using forward differences)*/ -float DualP_2D(float *U, float *V1, float *V2, float *P1, float *P2, int dimX, int dimY, float sigma) +float DualP_2D(float *U, float *V1, float *V2, float *P1, float *P2, long dimX, long dimY, float sigma) { - int i,j, index; + long i,j, index; #pragma omp parallel for shared(U,V1,V2,P1,P2) private(i,j,index) for(i=0; i<dimX; i++) { for(j=0; j<dimY; j++) { @@ -123,10 +124,10 @@ float DualP_2D(float *U, float *V1, float *V2, float *P1, float *P2, int dimX, i return 1; } /*Projection onto convex set for P*/ -float ProjP_2D(float *P1, float *P2, int dimX, int dimY, float alpha1) +float ProjP_2D(float *P1, float *P2, long dimX, long dimY, float alpha1) { float grad_magn; - int i,j,index; + long i,j,index; #pragma omp parallel for shared(P1,P2) private(i,j,index,grad_magn) for(i=0; i<dimX; i++) { for(j=0; j<dimY; j++) { @@ -141,9 +142,9 @@ float ProjP_2D(float *P1, float *P2, int dimX, int dimY, float alpha1) return 1; } /*Calculating dual variable Q (using forward differences)*/ -float DualQ_2D(float *V1, float *V2, float *Q1, float *Q2, float *Q3, int dimX, int dimY, float sigma) +float DualQ_2D(float *V1, float *V2, float *Q1, float *Q2, float *Q3, long dimX, long dimY, float sigma) { - int i,j,index; + 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<dimX; i++) { @@ -172,10 +173,10 @@ float DualQ_2D(float *V1, float *V2, float *Q1, float *Q2, float *Q3, int dimX, }} return 1; } -float ProjQ_2D(float *Q1, float *Q2, float *Q3, int dimX, int dimY, float alpha0) +float ProjQ_2D(float *Q1, float *Q2, float *Q3, long dimX, long dimY, float alpha0) { float grad_magn; - int i,j,index; + long i,j,index; #pragma omp parallel for shared(Q1,Q2,Q3) private(i,j,index,grad_magn) for(i=0; i<dimX; i++) { for(j=0; j<dimY; j++) { @@ -191,9 +192,9 @@ float ProjQ_2D(float *Q1, float *Q2, float *Q3, int dimX, int dimY, float alpha0 return 1; } /* Divergence and projection for P*/ -float DivProjP_2D(float *U, float *U0, float *P1, float *P2, int dimX, int dimY, float lambda, float tau) +float DivProjP_2D(float *U, float *U0, float *P1, float *P2, long dimX, long dimY, float lambda, float tau) { - int i,j,index; + long i,j,index; float P_v1, P_v2, div; #pragma omp parallel for shared(U,U0,P1,P2) private(i,j,index,P_v1,P_v2,div) for(i=0; i<dimX; i++) { @@ -209,17 +210,17 @@ float DivProjP_2D(float *U, float *U0, float *P1, float *P2, int dimX, int dimY, return *U; } /*get updated solution U*/ -float newU(float *U, float *U_old, int dimX, int dimY) +float newU(float *U, float *U_old, long dimX, long dimY) { - int i; + long i; #pragma omp parallel for shared(U,U_old) private(i) for(i=0; i<dimX*dimY; i++) U[i] = 2*U[i] - U_old[i]; return *U; } /*get update for V*/ -float UpdV_2D(float *V1, float *V2, float *P1, float *P2, float *Q1, float *Q2, float *Q3, int dimX, int dimY, float tau) +float UpdV_2D(float *V1, float *V2, float *P1, float *P2, float *Q1, float *Q2, float *Q3, long dimX, long dimY, float tau) { - int i, j, index; + long i, j, index; float q1, q11, q2, q22, div1, div2; #pragma omp parallel for shared(V1,V2,P1,P2,Q1,Q2,Q3) private(i, j, index, q1, q11, q2, q22, div1, div2) for(i=0; i<dimX; i++) { diff --git a/Core/regularisers_CPU/TGV_core.h b/Core/regularisers_CPU/TGV_core.h index d5632b9..29482ad 100644 --- a/Core/regularisers_CPU/TGV_core.h +++ b/Core/regularisers_CPU/TGV_core.h @@ -49,13 +49,14 @@ extern "C" { #endif /* 2D functions */ CCPI_EXPORT float TGV_main(float *U0, float *U, float lambda, float alpha1, float alpha0, int iter, float L2, int dimX, int dimY); -CCPI_EXPORT float DualP_2D(float *U, float *V1, float *V2, float *P1, float *P2, int dimX, int dimY, float sigma); -CCPI_EXPORT float ProjP_2D(float *P1, float *P2, int dimX, int dimY, float alpha1); -CCPI_EXPORT float DualQ_2D(float *V1, float *V2, float *Q1, float *Q2, float *Q3, int dimX, int dimY, float sigma); -CCPI_EXPORT float ProjQ_2D(float *Q1, float *Q2, float *Q3, int dimX, int dimY, float alpha0); -CCPI_EXPORT float DivProjP_2D(float *U, float *U0, float *P1, float *P2, int dimX, int dimY, float lambda, float tau); -CCPI_EXPORT float UpdV_2D(float *V1, float *V2, float *P1, float *P2, float *Q1, float *Q2, float *Q3, int dimX, int dimY, float tau); -CCPI_EXPORT float newU(float *U, float *U_old, int dimX, int dimY); + +CCPI_EXPORT float DualP_2D(float *U, float *V1, float *V2, float *P1, float *P2, long dimX, long dimY, float sigma); +CCPI_EXPORT float ProjP_2D(float *P1, float *P2, long dimX, long dimY, float alpha1); +CCPI_EXPORT float DualQ_2D(float *V1, float *V2, float *Q1, float *Q2, float *Q3, long dimX, long dimY, float sigma); +CCPI_EXPORT float ProjQ_2D(float *Q1, float *Q2, float *Q3, long dimX, long dimY, float alpha0); +CCPI_EXPORT float DivProjP_2D(float *U, float *U0, float *P1, float *P2, long dimX, long dimY, float lambda, float tau); +CCPI_EXPORT float UpdV_2D(float *V1, float *V2, float *P1, float *P2, float *Q1, float *Q2, float *Q3, long dimX, long dimY, float tau); +CCPI_EXPORT float newU(float *U, float *U_old, long dimX, long dimY); #ifdef __cplusplus } #endif diff --git a/Core/regularisers_CPU/TNV_core.c b/Core/regularisers_CPU/TNV_core.c index c51d6cd..753cc5f 100755 --- a/Core/regularisers_CPU/TNV_core.c +++ b/Core/regularisers_CPU/TNV_core.c @@ -39,16 +39,16 @@ float TNV_CPU_main(float *Input, float *u, float lambda, int maxIter, float tol, int dimX, int dimY, int dimZ) { - int k, p, q, r, DimTotal; + long k, p, q, r, DimTotal; float taulambda; float *u_upd, *gx, *gy, *gx_upd, *gy_upd, *qx, *qy, *qx_upd, *qy_upd, *v, *vx, *vy, *gradx, *grady, *gradx_upd, *grady_upd, *gradx_ubar, *grady_ubar, *div, *div_upd; - p = 1; - q = 1; - r = 0; + p = 1l; + q = 1l; + r = 0l; lambda = 1.0f/(2.0f*lambda); - DimTotal = dimX*dimY*dimZ; + DimTotal = (long)(dimX*dimY*dimZ); /* PDHG algorithm parameters*/ float tau = 0.5f; float sigma = 0.5f; @@ -106,10 +106,10 @@ float TNV_CPU_main(float *Input, float *u, float lambda, int maxIter, float tol, for(k=0; k<dimX*dimY*dimZ; k++) {v[k] = u[k] + tau*div[k];} // Proximal solution of fidelity term -proxG(u_upd, v, Input, taulambda, dimX, dimY, dimZ); +proxG(u_upd, v, Input, taulambda, (long)(dimX), (long)(dimY), (long)(dimZ)); // Gradient of updated primal variable -gradient(u_upd, gradx_upd, grady_upd, dimX, dimY, dimZ); +gradient(u_upd, gradx_upd, grady_upd, (long)(dimX), (long)(dimY), (long)(dimZ)); // Argument of proximal mapping of regularization term #pragma omp parallel for shared(gradx_upd, grady_upd, gradx, grady) private(k, ubarx, ubary) @@ -122,7 +122,7 @@ for(k=0; k<dimX*dimY*dimZ; k++) { grady_ubar[k] = ubary; } -proxF(gx_upd, gy_upd, vx, vy, sigma, p, q, r, dimX, dimY, dimZ); +proxF(gx_upd, gy_upd, vx, vy, sigma, p, q, r, (long)(dimX), (long)(dimY), (long)(dimZ)); // Update dual variable #pragma omp parallel for shared(qx_upd, qy_upd) private(k) @@ -174,14 +174,14 @@ if(b > 1) sigma = (beta * sigma) / b; alpha = alpha0; - copyIm(u, u_upd, dimX, dimY, dimZ); - copyIm(gx, gx_upd, dimX, dimY, dimZ); - copyIm(gy, gy_upd, dimX, dimY, dimZ); - copyIm(qx, qx_upd, dimX, dimY, dimZ); - copyIm(qy, qy_upd, dimX, dimY, dimZ); - copyIm(gradx, gradx_upd, dimX, dimY, dimZ); - copyIm(grady, grady_upd, dimX, dimY, dimZ); - copyIm(div, div_upd, dimX, dimY, dimZ); + copyIm(u, u_upd, (long)(dimX), (long)(dimY), (long)(dimZ)); + copyIm(gx, gx_upd, (long)(dimX), (long)(dimY), (long)(dimZ)); + copyIm(gy, gy_upd, (long)(dimX), (long)(dimY), (long)(dimZ)); + copyIm(qx, qx_upd, (long)(dimX), (long)(dimY), (long)(dimZ)); + copyIm(qy, qy_upd, (long)(dimX), (long)(dimY), (long)(dimZ)); + copyIm(gradx, gradx_upd, (long)(dimX), (long)(dimY), (long)(dimZ)); + copyIm(grady, grady_upd, (long)(dimX), (long)(dimY), (long)(dimZ)); + copyIm(div, div_upd, (long)(dimX), (long)(dimY), (long)(dimZ)); } else if(resprimal > dual_dot_delta) { @@ -205,14 +205,14 @@ taulambda = tau * lambda; divsigma = 1.0f / sigma; divtau = 1.0f / tau; -copyIm(u_upd, u, dimX, dimY, dimZ); -copyIm(gx_upd, gx, dimX, dimY, dimZ); -copyIm(gy_upd, gy, dimX, dimY, dimZ); -copyIm(qx_upd, qx, dimX, dimY, dimZ); -copyIm(qy_upd, qy, dimX, dimY, dimZ); -copyIm(gradx_upd, gradx, dimX, dimY, dimZ); -copyIm(grady_upd, grady, dimX, dimY, dimZ); -copyIm(div_upd, div, dimX, dimY, dimZ); +copyIm(u_upd, u, (long)(dimX), (long)(dimY), (long)(dimZ)); +copyIm(gx_upd, gx, (long)(dimX), (long)(dimY), (long)(dimZ)); +copyIm(gy_upd, gy, (long)(dimX), (long)(dimY), (long)(dimZ)); +copyIm(qx_upd, qx, (long)(dimX), (long)(dimY), (long)(dimZ)); +copyIm(qy_upd, qy, (long)(dimX), (long)(dimY), (long)(dimZ)); +copyIm(gradx_upd, gradx, (long)(dimX), (long)(dimY), (long)(dimZ)); +copyIm(grady_upd, grady, (long)(dimX), (long)(dimY), (long)(dimZ)); +copyIm(div_upd, div, (long)(dimX), (long)(dimY), (long)(dimZ)); // Compute residual at current iteration residual = (resprimal + resdual) / ((float) (dimX*dimY*dimZ)); @@ -227,15 +227,14 @@ if (residual < tol) { free (u_upd); free(gx); free(gy); free(gx_upd); free(gy_upd); free(qx); free(qy); free(qx_upd); free(qy_upd); free(v); free(vx); free(vy); free(gradx); free(grady); free(gradx_upd); free(grady_upd); free(gradx_ubar); - free(grady_ubar); free(div); free(div_upd); - + free(grady_ubar); free(div); free(div_upd); return *u; } -float proxG(float *u_upd, float *v, float *f, float taulambda, int dimX, int dimY, int dimZ) +float proxG(float *u_upd, float *v, float *f, float taulambda, long dimX, long dimY, long dimZ) { float constant; - int k; + long k; constant = 1.0f + taulambda; #pragma omp parallel for shared(v, f, u_upd) private(k) for(k=0; k<dimZ*dimX*dimY; k++) { @@ -245,9 +244,9 @@ float proxG(float *u_upd, float *v, float *f, float taulambda, int dimX, int dim return *u_upd; } -float gradient(float *u_upd, float *gradx_upd, float *grady_upd, int dimX, int dimY, int dimZ) +float gradient(float *u_upd, float *gradx_upd, float *grady_upd, long dimX, long dimY, long dimZ) { - int i, j, k, l; + long i, j, k, l; // Compute discrete gradient using forward differences #pragma omp parallel for shared(gradx_upd,grady_upd,u_upd) private(i, j, k, l) for(k = 0; k < dimZ; k++) { @@ -270,7 +269,7 @@ float gradient(float *u_upd, float *gradx_upd, float *grady_upd, int dimX, int d return 1; } -float proxF(float *gx, float *gy, float *vx, float *vy, float sigma, int p, int q, int r, int dimX, int dimY, int dimZ) +float proxF(float *gx, float *gy, float *vx, float *vy, float sigma, int p, int q, int r, long dimX, long dimY, long dimZ) { // (S^p, \ell^1) norm decouples at each pixel // Spl1(gx, gy, vx, vy, sigma, p, num_channels, dim); @@ -289,7 +288,7 @@ float proxF(float *gx, float *gy, float *vx, float *vy, float sigma, int p, int // Auxiliar vector float *proj, sum, shrinkfactor ; float M1,M2,M3,valuex,valuey,T,D,det,eig1,eig2,sig1,sig2,V1, V2, V3, V4, v0,v1,v2, mu1,mu2,sig1_upd,sig2_upd,t1,t2,t3; - int i,j,k, ii, num; + long i,j,k, ii, num; #pragma omp parallel for shared (gx,gy,vx,vy,p) private(i,ii,j,k,proj,num, sum, shrinkfactor, M1,M2,M3,valuex,valuey,T,D,det,eig1,eig2,sig1,sig2,V1, V2, V3, V4,v0,v1,v2,mu1,mu2,sig1_upd,sig2_upd,t1,t2,t3) for(i=0; i<dimX; i++) { for(j=0; j<dimY; j++) { @@ -372,7 +371,7 @@ float proxF(float *gx, float *gy, float *vx, float *vy, float sigma, int p, int /*l1 projection part */ sum = fLarge; - num = 0; + num = 0l; shrinkfactor = 0.0f; while(sum > 1.0f) { @@ -424,9 +423,9 @@ float proxF(float *gx, float *gy, float *vx, float *vy, float sigma, int p, int return 1; } -float divergence(float *qx_upd, float *qy_upd, float *div_upd, int dimX, int dimY, int dimZ) +float divergence(float *qx_upd, float *qy_upd, float *div_upd, long dimX, long dimY, long dimZ) { - int i, j, k, l; + long i, j, k, l; #pragma omp parallel for shared(qx_upd,qy_upd,div_upd) private(i, j, k, l) for(k = 0; k < dimZ; k++) { for(j = 0; j < dimY; j++) { diff --git a/Core/regularisers_CPU/TNV_core.h b/Core/regularisers_CPU/TNV_core.h index c082694..aa050a4 100644 --- a/Core/regularisers_CPU/TNV_core.h +++ b/Core/regularisers_CPU/TNV_core.h @@ -38,10 +38,10 @@ extern "C" { CCPI_EXPORT float TNV_CPU_main(float *Input, float *u, float lambda, int maxIter, float tol, int dimX, int dimY, int dimZ); /*float PDHG(float *A, float *B, float tau, float sigma, float theta, float lambda, int p, int q, int r, float tol, int maxIter, int d_c, int d_w, int d_h);*/ -CCPI_EXPORT float proxG(float *u_upd, float *v, float *f, float taulambda, int dimX, int dimY, int dimZ); -CCPI_EXPORT float gradient(float *u_upd, float *gradx_upd, float *grady_upd, int dimX, int dimY, int dimZ); -CCPI_EXPORT float proxF(float *gx, float *gy, float *vx, float *vy, float sigma, int p, int q, int r, int dimX, int dimY, int dimZ); -CCPI_EXPORT float divergence(float *qx_upd, float *qy_upd, float *div_upd, int dimX, int dimY, int dimZ); +CCPI_EXPORT float proxG(float *u_upd, float *v, float *f, float taulambda, long dimX, long dimY, long dimZ); +CCPI_EXPORT float gradient(float *u_upd, float *gradx_upd, float *grady_upd, long dimX, long dimY, long dimZ); +CCPI_EXPORT float proxF(float *gx, float *gy, float *vx, float *vy, float sigma, int p, int q, int r, long dimX, long dimY, long dimZ); +CCPI_EXPORT float divergence(float *qx_upd, float *qy_upd, float *div_upd, long dimX, long dimY, long dimZ); #ifdef __cplusplus } #endif
\ No newline at end of file diff --git a/Core/regularisers_CPU/utils.c b/Core/regularisers_CPU/utils.c index 3eebdaa..7a4e80b 100644 --- a/Core/regularisers_CPU/utils.c +++ b/Core/regularisers_CPU/utils.c @@ -21,9 +21,9 @@ limitations under the License. #include <math.h> /* Copy Image (float) */ -float copyIm(float *A, float *U, int dimX, int dimY, int dimZ) +float copyIm(float *A, float *U, long dimX, long dimY, long dimZ) { - int j; + long j; #pragma omp parallel for shared(A, U) private(j) for (j = 0; j<dimX*dimY*dimZ; j++) U[j] = A[j]; return *U; @@ -88,18 +88,18 @@ float TV_energy2D(float *U, float *U0, float *E_val, float lambda, int type, int float TV_energy3D(float *U, float *U0, float *E_val, float lambda, int type, int dimX, int dimY, int dimZ) { - int i, j, k, i1, j1, k1, index; + long i, j, k, i1, j1, k1, index; float NOMx_2, NOMy_2, NOMz_2, E_Grad=0.0f, E_Data=0.0f; /* first calculate \grad U_xy*/ - for(j=0; j<dimY; j++) { - for(i=0; i<dimX; i++) { - for(k=0; k<dimZ; k++) { + for(j=0; j<(long)(dimY); j++) { + for(i=0; i<(long)(dimX); i++) { + for(k=0; k<(long)(dimZ); k++) { index = (dimX*dimY)*k + j*dimX+i; /* boundary conditions */ - i1 = i + 1; if (i == dimX-1) i1 = i; - j1 = j + 1; if (j == dimY-1) j1 = j; - k1 = k + 1; if (k == dimZ-1) k1 = k; + i1 = i + 1; if (i == (long)(dimX-1)) i1 = i; + j1 = j + 1; if (j == (long)(dimY-1)) j1 = j; + k1 = k + 1; if (k == (long)(dimZ-1)) k1 = k; /* Forward differences */ NOMx_2 = powf((float)(U[(dimX*dimY)*k + j1*dimX+i] - U[index]),2); /* x+ */ diff --git a/Core/regularisers_CPU/utils.h b/Core/regularisers_CPU/utils.h index 7ad83ab..cfaf6d7 100644 --- a/Core/regularisers_CPU/utils.h +++ b/Core/regularisers_CPU/utils.h @@ -24,7 +24,7 @@ limitations under the License. #ifdef __cplusplus extern "C" { #endif -CCPI_EXPORT float copyIm(float *A, float *U, int dimX, int dimY, int dimZ); +CCPI_EXPORT float copyIm(float *A, float *U, long dimX, long dimY, long dimZ); CCPI_EXPORT unsigned char copyIm_unchar(unsigned char *A, unsigned char *U, int dimX, int dimY, int dimZ); CCPI_EXPORT float copyIm_roll(float *A, float *U, int dimX, int dimY, int roll_value, int switcher); CCPI_EXPORT float TV_energy2D(float *U, float *U0, float *E_val, float lambda, int type, int dimX, int dimY); diff --git a/Wrappers/Matlab/mex_compile/compileCPU_mex_Linux.m b/Wrappers/Matlab/mex_compile/compileCPU_mex_Linux.m index 064b416..767d811 100644 --- a/Wrappers/Matlab/mex_compile/compileCPU_mex_Linux.m +++ b/Wrappers/Matlab/mex_compile/compileCPU_mex_Linux.m @@ -61,7 +61,7 @@ fprintf('%s \n', 'Compiling Nonlinear/Linear diffusion inpainting...'); mex NonlDiff_Inp.c Diffusion_Inpaint_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp" movefile('NonlDiff_Inp.mex*',Pathmove); -fprintf('%s \n', 'Compiling Nonlocal marching method for inpaiting...'); +fprintf('%s \n', 'Compiling Nonlocal marching method for inpainting...'); mex NonlocalMarching_Inpaint.c NonlocalMarching_Inpaint_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp" movefile('NonlocalMarching_Inpaint.mex*',Pathmove); |