loop CPU routines upgrades

12 files changed, 1280 insertions, 1287 deletions

diff --git a/demos/Matlab_demos/demoMatlab_3Ddenoise.m b/demos/Matlab_demos/demoMatlab_3Ddenoise.m
index d7ff60c..f018327 100644
--- a/demos/Matlab_demos/demoMatlab_3Ddenoise.m
+++ b/demos/Matlab_demos/demoMatlab_3Ddenoise.m
@@ -182,19 +182,18 @@ eta =  0.2; % Reference image gradient smoothing constant
 tic; u_fgp_dtv = FGP_dTV(single(vol3D), single(vol3D_ref), lambda_reg, iter_fgp, epsil_tol, eta); toc; 
 figure; imshow(u_fgp_dtv(:,:,7), [0 1]); title('FGP-dTV denoised volume (CPU)');
 %%
-fprintf('Denoise a volume using the FGP-dTV model (GPU) \n');
-
-% create another volume (reference) with slightly less amount of noise
-vol3D_ref = zeros(N,N,slices, 'single');
-for i = 1:slices
-vol3D_ref(:,:,i) = Im + .01*randn(size(Im)); 
-end
-vol3D_ref(vol3D_ref < 0) = 0;
-% vol3D_ref = zeros(size(Im),'single'); % pass zero reference (dTV -> TV)
-
-iter_fgp = 300; % number of FGP iterations
-epsil_tol =  0.0; % tolerance
-eta =  0.2; % Reference image gradient smoothing constant
-tic; u_fgp_dtv_g = FGP_dTV_GPU(single(vol3D), single(vol3D_ref), lambda_reg, iter_fgp, epsil_tol, eta); toc; 
-figure; imshow(u_fgp_dtv_g(:,:,7), [0 1]); title('FGP-dTV denoised volume (GPU)');
+% fprintf('Denoise a volume using the FGP-dTV model (GPU) \n');
+% % create another volume (reference) with slightly less amount of noise
+% vol3D_ref = zeros(N,N,slices, 'single');
+% for i = 1:slices
+% vol3D_ref(:,:,i) = Im + .01*randn(size(Im)); 
+% end
+% vol3D_ref(vol3D_ref < 0) = 0;
+% % vol3D_ref = zeros(size(Im),'single'); % pass zero reference (dTV -> TV)
+% 
+% iter_fgp = 300; % number of FGP iterations
+% epsil_tol =  0.0; % tolerance
+% eta =  0.2; % Reference image gradient smoothing constant
+% tic; u_fgp_dtv_g = FGP_dTV_GPU(single(vol3D), single(vol3D_ref), lambda_reg, iter_fgp, epsil_tol, eta); toc; 
+% figure; imshow(u_fgp_dtv_g(:,:,7), [0 1]); title('FGP-dTV denoised volume (GPU)');
 %%
diff --git a/demos/Matlab_demos/demoMatlab_denoise.m b/demos/Matlab_demos/demoMatlab_denoise.m
index 12d5570..b50eaf5 100644
--- a/demos/Matlab_demos/demoMatlab_denoise.m
+++ b/demos/Matlab_demos/demoMatlab_denoise.m
@@ -149,7 +149,6 @@ fprintf('%s %f \n', 'MSSIM error for NLTV is:', ssimval);
 figure; imagesc(u_nltv, [0 1]); colormap(gray); daspect([1 1 1]); title('Non-local Total Variation denoised image (CPU)');
 %%
 %>>>>>>>>>>>>>> MULTI-CHANNEL priors <<<<<<<<<<<<<<< %
-
 fprintf('Denoise using the FGP-dTV model (CPU) \n');
 % create another image (reference) with slightly less amount of noise
 u_ref = Im + .01*randn(size(Im)); u_ref(u_ref < 0) = 0;
diff --git a/src/Core/regularisers_CPU/Diffus4th_order_core.c b/src/Core/regularisers_CPU/Diffus4th_order_core.c
index 28ac8a9..840b736 100644
--- a/src/Core/regularisers_CPU/Diffus4th_order_core.c
+++ b/src/Core/regularisers_CPU/Diffus4th_order_core.c
@@ -50,35 +50,35 @@ float Diffus4th_CPU_main(float *Input, float *Output, float *infovector, float l
     float *W_Lapl=NULL, *Output_prev=NULL;
     sigmaPar2 = sigmaPar*sigmaPar;
     DimTotal =  dimX*dimY*dimZ;
-
+    
     W_Lapl = calloc(DimTotal, sizeof(float));
-
+    
     if (epsil != 0.0f) Output_prev = calloc(DimTotal, sizeof(float));
-
+    
     /* copy into output */
     copyIm(Input, Output, (long)(dimX), (long)(dimY), (long)(dimZ));
-
+    
     for(i=0; i < iterationsNumb; i++) {
-      if ((epsil != 0.0f) && (i % 5 == 0)) copyIm(Output, Output_prev, (long)(dimX), (long)(dimY), (long)(dimZ));
-
-      if (dimZ == 1) {
+        if ((epsil != 0.0f) && (i % 5 == 0)) copyIm(Output, Output_prev, (long)(dimX), (long)(dimY), (long)(dimZ));
+        
+        if (dimZ == 1) {
             /* running 2D diffusion iterations */
             /* Calculating weighted Laplacian */
             Weighted_Laplc2D(W_Lapl, Output, sigmaPar2, dimX, dimY);
             /* Perform iteration step */
             Diffusion_update_step2D(Output, Input, W_Lapl, lambdaPar, sigmaPar2, tau, (long)(dimX), (long)(dimY));
-		        }
-      else {
+        }
+        else {
             /* running 3D diffusion iterations */
             /* Calculating weighted Laplacian */
             Weighted_Laplc3D(W_Lapl, Output, sigmaPar2, dimX, dimY, dimZ);
             /* Perform iteration step */
             Diffusion_update_step3D(Output, Input, W_Lapl, lambdaPar, sigmaPar2, tau, (long)(dimX), (long)(dimY), (long)(dimZ));
-          }
-
-          /* check early stopping criteria */
-          if ((epsil != 0.0f) && (i % 5 == 0)) {
-          re = 0.0f; re1 = 0.0f;
+        }
+        
+        /* check early stopping criteria */
+        if ((epsil != 0.0f) && (i % 5 == 0)) {
+            re = 0.0f; re1 = 0.0f;
             for(j=0; j<DimTotal; j++)
             {
                 re += powf(Output[j] - Output_prev[j],2);
@@ -87,10 +87,10 @@ float Diffus4th_CPU_main(float *Input, float *Output, float *infovector, float l
             re = sqrtf(re)/sqrtf(re1);
             if (re < epsil)  count++;
             if (count > 3) break;
-          }
-	   }
-	  free(W_Lapl);
-
+        }
+    }
+    free(W_Lapl);
+    
     if (epsil != 0.0f) free(Output_prev);
     /*adding info into info_vector */
     infovector[0] = (float)(i);  /*iterations number (if stopped earlier based on tolerance)*/
@@ -104,74 +104,71 @@ float Weighted_Laplc2D(float *W_Lapl, float *U0, float sigma, long dimX, long di
 {
     long i,j,i1,i2,j1,j2,index;
     float gradX, gradX_sq, gradY, gradY_sq, gradXX, gradYY, gradXY, xy_2, denom, V_norm, V_orth, c, c_sq;
-
-        #pragma omp parallel for shared(W_Lapl) private(i,j,i1,i2,j1,j2,index,gradX, gradX_sq, gradY, gradY_sq, gradXX, gradYY, gradXY, xy_2, denom, V_norm, V_orth, c, c_sq)
+    
+#pragma omp parallel for shared(W_Lapl) private(i,j,i1,i2,j1,j2,index,gradX, gradX_sq, gradY, gradY_sq, gradXX, gradYY, gradXY, xy_2, denom, V_norm, V_orth, c, c_sq)
+    for(j=0; j<dimY; j++) {
+        /* symmetric boundary conditions */
+        j1 = j+1; if (j1 == dimY) j1 = j-1;
+        j2 = j-1; if (j2 < 0) j2 = j+1;
         for(i=0; i<dimX; i++) {
-			 /* symmetric boundary conditions */
-			i1 = i+1; if (i1 == dimX) i1 = i-1;
-			i2 = i-1; if (i2 < 0) i2 = i+1;
-            for(j=0; j<dimY; j++) {
-				 /* symmetric boundary conditions */
-				j1 = j+1; if (j1 == dimY) j1 = j-1;
-				j2 = j-1; if (j2 < 0) j2 = j+1;
-
-				index = j*dimX+i;
-
-				gradX = 0.5f*(U0[j*dimX+i2] - U0[j*dimX+i1]);
-				gradX_sq = pow(gradX,2);
-
-				gradY = 0.5f*(U0[j2*dimX+i] - U0[j1*dimX+i]);
-                gradY_sq = pow(gradY,2);
-
-                gradXX = U0[j*dimX+i2] + U0[j*dimX+i1] - 2*U0[index];
-                gradYY = U0[j2*dimX+i] + U0[j1*dimX+i] - 2*U0[index];
-
-                gradXY = 0.25f*(U0[j2*dimX+i2] + U0[j1*dimX+i1] - U0[j1*dimX+i2] - U0[j2*dimX+i1]);
-                xy_2 = 2.0f*gradX*gradY*gradXY;
-
-                denom =  gradX_sq + gradY_sq;
-
-                if (denom <= EPS) {
-                    V_norm = (gradXX*gradX_sq + xy_2 + gradYY*gradY_sq)/EPS;
-                    V_orth = (gradXX*gradY_sq - xy_2 + gradYY*gradX_sq)/EPS;
-                    }
-                else  {
-                    V_norm = (gradXX*gradX_sq + xy_2 + gradYY*gradY_sq)/denom;
-                    V_orth = (gradXX*gradY_sq - xy_2 + gradYY*gradX_sq)/denom;
-                    }
-
-                c = 1.0f/(1.0f + denom/sigma);
-                c_sq = c*c;
-
-                W_Lapl[index] = c_sq*V_norm + c*V_orth;
+            /* symmetric boundary conditions */
+            i1 = i+1; if (i1 == dimX) i1 = i-1;
+            i2 = i-1; if (i2 < 0) i2 = i+1;
+            index = j*dimX+i;
+            
+            gradX = 0.5f*(U0[j*dimX+i2] - U0[j*dimX+i1]);
+            gradX_sq = pow(gradX,2);
+            
+            gradY = 0.5f*(U0[j2*dimX+i] - U0[j1*dimX+i]);
+            gradY_sq = pow(gradY,2);
+            
+            gradXX = U0[j*dimX+i2] + U0[j*dimX+i1] - 2*U0[index];
+            gradYY = U0[j2*dimX+i] + U0[j1*dimX+i] - 2*U0[index];
+            
+            gradXY = 0.25f*(U0[j2*dimX+i2] + U0[j1*dimX+i1] - U0[j1*dimX+i2] - U0[j2*dimX+i1]);
+            xy_2 = 2.0f*gradX*gradY*gradXY;
+            
+            denom =  gradX_sq + gradY_sq;
+            
+            if (denom <= EPS) {
+                V_norm = (gradXX*gradX_sq + xy_2 + gradYY*gradY_sq)/EPS;
+                V_orth = (gradXX*gradY_sq - xy_2 + gradYY*gradX_sq)/EPS;
             }
-        }
-        return *W_Lapl;
+            else  {
+                V_norm = (gradXX*gradX_sq + xy_2 + gradYY*gradY_sq)/denom;
+                V_orth = (gradXX*gradY_sq - xy_2 + gradYY*gradX_sq)/denom;
+            }
+            
+            c = 1.0f/(1.0f + denom/sigma);
+            c_sq = c*c;
+            
+            W_Lapl[index] = c_sq*V_norm + c*V_orth;
+        }}
+    return *W_Lapl;
 }
 
 float Diffusion_update_step2D(float *Output, float *Input, float *W_Lapl, float lambdaPar, float sigmaPar2, float tau, long dimX, long dimY)
 {
-	long 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)
+    
+#pragma omp parallel for shared(Output, Input, W_Lapl) private(i,j,i1,i2,j1,j2,index,gradXXc,gradYYc)
+    for(j=0; j<dimY; j++) {
+        /* symmetric boundary conditions */
+        j1 = j+1; if (j1 == dimY) j1 = j-1;
+        j2 = j-1; if (j2 < 0) j2 = j+1;
         for(i=0; i<dimX; i++) {
-			 /* symmetric boundary conditions */
-			i1 = i+1; if (i1 == dimX) i1 = i-1;
-			i2 = i-1; if (i2 < 0) i2 = i+1;
-            for(j=0; j<dimY; j++) {
-				 /* symmetric boundary conditions */
-				j1 = j+1; if (j1 == dimY) j1 = j-1;
-				j2 = j-1; if (j2 < 0) j2 = j+1;
-					index = j*dimX+i;
-
-                    gradXXc = W_Lapl[j*dimX+i2] + W_Lapl[j*dimX+i1] - 2*W_Lapl[index];
-                    gradYYc = W_Lapl[j2*dimX+i] + W_Lapl[j1*dimX+i] - 2*W_Lapl[index];
-
-                    Output[index] += tau*(-lambdaPar*(gradXXc + gradYYc) - (Output[index] - Input[index]));
-                }
-            }
-	return *Output;
+            /* symmetric boundary conditions */
+            i1 = i+1; if (i1 == dimX) i1 = i-1;
+            i2 = i-1; if (i2 < 0) i2 = i+1;
+            index = j*dimX+i;
+            
+            gradXXc = W_Lapl[j*dimX+i2] + W_Lapl[j*dimX+i1] - 2*W_Lapl[index];
+            gradYYc = W_Lapl[j2*dimX+i] + W_Lapl[j1*dimX+i] - 2*W_Lapl[index];
+            
+            Output[index] += tau*(-lambdaPar*(gradXXc + gradYYc) - (Output[index] - Input[index]));
+        }}
+    return *Output;
 }
 /********************************************************************/
 /***************************3D Functions*****************************/
@@ -180,95 +177,89 @@ float Weighted_Laplc3D(float *W_Lapl, float *U0, float sigma, long dimX, long di
 {
     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)
-        for(i=0; i<dimX; i++) {
-			 /* symmetric boundary conditions */
-			i1 = i+1; if (i1 == dimX) i1 = i-1;
-			i2 = i-1; if (i2 < 0) i2 = i+1;
-            for(j=0; j<dimY; j++) {
-				/* symmetric boundary conditions */
-				j1 = j+1; if (j1 == dimY) j1 = j-1;
-				j2 = j-1; if (j2 < 0) j2 = j+1;
-
-				for(k=0; k<dimZ; k++) {
-				/* symmetric boundary conditions */
-				k1 = k+1; if (k1 == dimZ) k1 = k-1;
-				k2 = k-1; if (k2 < 0) k2 = k+1;
-
-				index = (dimX*dimY)*k + j*dimX+i;
-
-				gradX = 0.5f*(U0[(dimX*dimY)*k + j*dimX+i2] - U0[(dimX*dimY)*k + j*dimX+i1]);
-				gradX_sq = pow(gradX,2);
-
-				gradY = 0.5f*(U0[(dimX*dimY)*k + j2*dimX+i] - U0[(dimX*dimY)*k + j1*dimX+i]);
+    
+#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)
+    for(k=0; k<dimZ; k++) {
+        /* symmetric boundary conditions */
+        k1 = k+1; if (k1 == dimZ) k1 = k-1;
+        k2 = k-1; if (k2 < 0) k2 = k+1;
+        for(j=0; j<dimY; j++) {
+            /* symmetric boundary conditions */
+            j1 = j+1; if (j1 == dimY) j1 = j-1;
+            j2 = j-1; if (j2 < 0) j2 = j+1;
+            for(i=0; i<dimX; i++) {
+                /* symmetric boundary conditions */
+                i1 = i+1; if (i1 == dimX) i1 = i-1;
+                i2 = i-1; if (i2 < 0) i2 = i+1;
+                
+                index = (dimX*dimY)*k + j*dimX+i;
+                
+                gradX = 0.5f*(U0[(dimX*dimY)*k + j*dimX+i2] - U0[(dimX*dimY)*k + j*dimX+i1]);
+                gradX_sq = pow(gradX,2);
+                
+                gradY = 0.5f*(U0[(dimX*dimY)*k + j2*dimX+i] - U0[(dimX*dimY)*k + j1*dimX+i]);
                 gradY_sq = pow(gradY,2);
-
+                
                 gradZ = 0.5f*(U0[(dimX*dimY)*k2 + j*dimX+i] - U0[(dimX*dimY)*k1 + j*dimX+i]);
                 gradZ_sq = pow(gradZ,2);
-
+                
                 gradXX = U0[(dimX*dimY)*k + j*dimX+i2] + U0[(dimX*dimY)*k + j*dimX+i1] - 2*U0[index];
                 gradYY = U0[(dimX*dimY)*k + j2*dimX+i] + U0[(dimX*dimY)*k + j1*dimX+i] - 2*U0[index];
                 gradZZ = U0[(dimX*dimY)*k2 + j*dimX+i] + U0[(dimX*dimY)*k1 + j*dimX+i] - 2*U0[index];
-
+                
                 gradXY = 0.25f*(U0[(dimX*dimY)*k + j2*dimX+i2] + U0[(dimX*dimY)*k + j1*dimX+i1] - U0[(dimX*dimY)*k + j1*dimX+i2] - U0[(dimX*dimY)*k + j2*dimX+i1]);
                 gradXZ = 0.25f*(U0[(dimX*dimY)*k2 + j*dimX+i2] - U0[(dimX*dimY)*k2+j*dimX+i1] - U0[(dimX*dimY)*k1+j*dimX+i2] + U0[(dimX*dimY)*k1+j*dimX+i1]);
                 gradYZ = 0.25f*(U0[(dimX*dimY)*k2 +j2*dimX+i] - U0[(dimX*dimY)*k2+j1*dimX+i] - U0[(dimX*dimY)*k1+j2*dimX+i] + U0[(dimX*dimY)*k1+j1*dimX+i]);
-
+                
                 xy_2  = 2.0f*gradX*gradY*gradXY;
                 xyz_1 = 2.0f*gradX*gradZ*gradXZ;
                 xyz_2 = 2.0f*gradY*gradZ*gradYZ;
-
+                
                 denom =  gradX_sq + gradY_sq + gradZ_sq;
-
-					if (denom <= EPS) {
-					V_norm = (gradXX*gradX_sq + gradYY*gradY_sq + gradZZ*gradZ_sq + xy_2 + xyz_1 + xyz_2)/EPS;
+                
+                if (denom <= EPS) {
+                    V_norm = (gradXX*gradX_sq + gradYY*gradY_sq + gradZZ*gradZ_sq + xy_2 + xyz_1 + xyz_2)/EPS;
                     V_orth = ((gradY_sq + gradZ_sq)*gradXX + (gradX_sq + gradZ_sq)*gradYY + (gradX_sq + gradY_sq)*gradZZ - xy_2 - xyz_1 - xyz_2)/EPS;
-					}
-					else  {
-					V_norm = (gradXX*gradX_sq + gradYY*gradY_sq + gradZZ*gradZ_sq + xy_2 + xyz_1 + xyz_2)/denom;
+                }
+                else  {
+                    V_norm = (gradXX*gradX_sq + gradYY*gradY_sq + gradZZ*gradZ_sq + xy_2 + xyz_1 + xyz_2)/denom;
                     V_orth = ((gradY_sq + gradZ_sq)*gradXX + (gradX_sq + gradZ_sq)*gradYY + (gradX_sq + gradY_sq)*gradZZ - xy_2 - xyz_1 - xyz_2)/denom;
-					}
-
+                }
+                
                 c = 1.0f/(1.0f + denom/sigma);
                 c_sq = c*c;
-
+                
                 W_Lapl[index] = c_sq*V_norm + c*V_orth;
-				}
-            }
-        }
-        return *W_Lapl;
+            }}}
+    return *W_Lapl;
 }
 
 float Diffusion_update_step3D(float *Output, float *Input, float *W_Lapl, float lambdaPar, float sigmaPar2, float tau, long dimX, long dimY, long dimZ)
 {
-	long 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)
-        for(i=0; i<dimX; i++) {
-			 /* symmetric boundary conditions */
-			i1 = i+1; if (i1 == dimX) i1 = i-1;
-			i2 = i-1; if (i2 < 0) i2 = i+1;
-            for(j=0; j<dimY; j++) {
-				/* symmetric boundary conditions */
-				j1 = j+1; if (j1 == dimY) j1 = j-1;
-				j2 = j-1; if (j2 < 0) j2 = j+1;
-
-				for(k=0; k<dimZ; k++) {
-				/* symmetric boundary conditions */
-				k1 = k+1; if (k1 == dimZ) k1 = k-1;
-				k2 = k-1; if (k2 < 0) k2 = k+1;
-
-				index = (dimX*dimY)*k + j*dimX+i;
-
-                    gradXXc = W_Lapl[(dimX*dimY)*k + j*dimX+i2] + W_Lapl[(dimX*dimY)*k + j*dimX+i1] - 2*W_Lapl[index];
-                    gradYYc = W_Lapl[(dimX*dimY)*k + j2*dimX+i] + W_Lapl[(dimX*dimY)*k + j1*dimX+i] - 2*W_Lapl[index];
-                    gradZZc = W_Lapl[(dimX*dimY)*k2 + j*dimX+i] + W_Lapl[(dimX*dimY)*k1 + j*dimX+i] - 2*W_Lapl[index];
-
-                    Output[index] += tau*(-lambdaPar*(gradXXc + gradYYc + gradZZc) - (Output[index] - Input[index]));
-                }
-            }
-		}
-	return *Output;
+    
+#pragma omp parallel for shared(Output, Input, W_Lapl) private(i,j,i1,i2,j1,j2,k,k1,k2,index,gradXXc,gradYYc,gradZZc)
+    for(k=0; k<dimZ; k++) {
+        /* symmetric boundary conditions */
+        k1 = k+1; if (k1 == dimZ) k1 = k-1;
+        k2 = k-1; if (k2 < 0) k2 = k+1;
+        for(j=0; j<dimY; j++) {
+            /* symmetric boundary conditions */
+            j1 = j+1; if (j1 == dimY) j1 = j-1;
+            j2 = j-1; if (j2 < 0) j2 = j+1;
+            for(i=0; i<dimX; i++) {
+                /* symmetric boundary conditions */
+                i1 = i+1; if (i1 == dimX) i1 = i-1;
+                i2 = i-1; if (i2 < 0) i2 = i+1;
+                
+                index = (dimX*dimY)*k + j*dimX+i;
+                
+                gradXXc = W_Lapl[(dimX*dimY)*k + j*dimX+i2] + W_Lapl[(dimX*dimY)*k + j*dimX+i1] - 2*W_Lapl[index];
+                gradYYc = W_Lapl[(dimX*dimY)*k + j2*dimX+i] + W_Lapl[(dimX*dimY)*k + j1*dimX+i] - 2*W_Lapl[index];
+                gradZZc = W_Lapl[(dimX*dimY)*k2 + j*dimX+i] + W_Lapl[(dimX*dimY)*k1 + j*dimX+i] - 2*W_Lapl[index];
+                
+                Output[index] += tau*(-lambdaPar*(gradXXc + gradYYc + gradZZc) - (Output[index] - Input[index]));
+            }}}
+    return *Output;
 }
diff --git a/src/Core/regularisers_CPU/Diffusion_core.c b/src/Core/regularisers_CPU/Diffusion_core.c
index 7f06dd8..f0039c7 100644
--- a/src/Core/regularisers_CPU/Diffusion_core.c
+++ b/src/Core/regularisers_CPU/Diffusion_core.c
@@ -40,7 +40,7 @@ int signNDFc(float x) {
  * 5. tau - time-marching step for explicit scheme
  * 6. Penalty type: 1 - Huber, 2 - Perona-Malik, 3 - Tukey Biweight
  * 7. eplsilon - tolerance constant
-
+ *
  * Output:
  * [1] Filtered/regularized image/volume
  * [2] Information vector which contains [iteration no., reached tolerance]
@@ -60,149 +60,148 @@ float Diffusion_CPU_main(float *Input, float *Output, float *infovector, float l
     re = 0.0f; re1 = 0.0f;
     int count = 0;
     DimTotal = (long)(dimX*dimY*dimZ);
-
+    
     if (epsil != 0.0f) Output_prev = calloc(DimTotal, sizeof(float));
-
+    
     /* copy into output */
     copyIm(Input, Output, (long)(dimX), (long)(dimY), (long)(dimZ));
-
+    
     for(i=0; i < iterationsNumb; i++) {
-
-      if ((epsil != 0.0f)  && (i % 5 == 0)) copyIm(Output, Output_prev, (long)(dimX), (long)(dimY), (long)(dimZ));
-      if (dimZ == 1) {
+        
+        if ((epsil != 0.0f)  && (i % 5 == 0)) copyIm(Output, Output_prev, (long)(dimX), (long)(dimY), (long)(dimZ));
+        if (dimZ == 1) {
             /* running 2D diffusion iterations */
             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 */
-        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));
-          }
-          /* check early stopping criteria if epsilon not equal zero */
-          if ((epsil != 0.0f)  && (i % 5 == 0)) {
-          re = 0.0f; re1 = 0.0f;
+        }
+        else {
+            /* running 3D diffusion iterations */
+            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));
+        }
+        /* check early stopping criteria if epsilon not equal zero */
+        if ((epsil != 0.0f)  && (i % 5 == 0)) {
+            re = 0.0f; re1 = 0.0f;
             for(j=0; j<DimTotal; j++)
             {
                 re += powf(Output[j] - Output_prev[j],2);
                 re1 += powf(Output[j],2);
             }
-          re = sqrtf(re)/sqrtf(re1);
-          /* stop if the norm residual is less than the tolerance EPS */
-          if (re < epsil)  count++;
-          if (count > 3) break;
-          }
-		}
-
+            re = sqrtf(re)/sqrtf(re1);
+            /* stop if the norm residual is less than the tolerance EPS */
+            if (re < epsil)  count++;
+            if (count > 3) break;
+        }
+    }
+    
     free(Output_prev);
-  /*adding info into info_vector */
+    /*adding info into info_vector */
     infovector[0] = (float)(i);  /*iterations number (if stopped earlier based on tolerance)*/
     infovector[1] = re;  /* reached tolerance */
     return 0;
 }
 
-
 /********************************************************************/
 /***************************2D Functions*****************************/
 /********************************************************************/
 /* linear diffusion (heat equation) */
 float LinearDiff2D(float *Input, float *Output, float lambdaPar, float tau, long dimX, long dimY)
 {
-	long i,j,i1,i2,j1,j2,index;
-	float e,w,n,s,e1,w1,n1,s1;
-
+    long i,j,i1,i2,j1,j2,index;
+    float e,w,n,s,e1,w1,n1,s1;
+    
 #pragma omp parallel for shared(Input) private(index,i,j,i1,i2,j1,j2,e,w,n,s,e1,w1,n1,s1)
-    for(i=0; i<dimX; i++) {
+    for(j=0; j<dimY; j++) {
         /* symmetric boundary conditions (Neuman) */
-        i1 = i+1; if (i1 == dimX) i1 = i-1;
-        i2 = i-1; if (i2 < 0) i2 = i+1;
-        for(j=0; j<dimY; j++) {
+        j1 = j+1; if (j1 == dimY) j1 = j-1;
+        j2 = j-1; if (j2 < 0) j2 = j+1;
+        for(i=0; i<dimX; i++) {
             /* symmetric boundary conditions (Neuman) */
-            j1 = j+1; if (j1 == dimY) j1 = j-1;
-            j2 = j-1; if (j2 < 0) j2 = j+1;
+            i1 = i+1; if (i1 == dimX) i1 = i-1;
+            i2 = i-1; if (i2 < 0) i2 = i+1;
             index = j*dimX+i;
-
-                e = Output[j*dimX+i1];
-                w = Output[j*dimX+i2];
-                n = Output[j1*dimX+i];
-                s = Output[j2*dimX+i];
-
-                e1 = e - Output[index];
-                w1 = w - Output[index];
-                n1 = n - Output[index];
-                s1 = s - Output[index];
-
-                Output[index] += tau*(lambdaPar*(e1 + w1 + n1 + s1) - (Output[index] - Input[index]));
-		}}
-	return *Output;
+            
+            e = Output[j*dimX+i1];
+            w = Output[j*dimX+i2];
+            n = Output[j1*dimX+i];
+            s = Output[j2*dimX+i];
+            
+            e1 = e - Output[index];
+            w1 = w - Output[index];
+            n1 = n - Output[index];
+            s1 = s - Output[index];
+            
+            Output[index] += tau*(lambdaPar*(e1 + w1 + n1 + s1) - (Output[index] - Input[index]));
+        }}
+    return *Output;
 }
 
 /* nonlinear diffusion */
 float NonLinearDiff2D(float *Input, float *Output, float lambdaPar, float sigmaPar, float tau, int penaltytype, long dimX, long dimY)
 {
-	long i,j,i1,i2,j1,j2,index;
-	float e,w,n,s,e1,w1,n1,s1;
-
+    long i,j,i1,i2,j1,j2,index;
+    float e,w,n,s,e1,w1,n1,s1;
+    
 #pragma omp parallel for shared(Input) private(index,i,j,i1,i2,j1,j2,e,w,n,s,e1,w1,n1,s1)
-    for(i=0; i<dimX; i++) {
+    for(j=0; j<dimY; j++) {
         /* symmetric boundary conditions (Neuman) */
-        i1 = i+1; if (i1 == dimX) i1 = i-1;
-        i2 = i-1; if (i2 < 0) i2 = i+1;
-        for(j=0; j<dimY; j++) {
+        j1 = j+1; if (j1 == dimY) j1 = j-1;
+        j2 = j-1; if (j2 < 0) j2 = j+1;
+        for(i=0; i<dimX; i++) {
             /* symmetric boundary conditions (Neuman) */
-            j1 = j+1; if (j1 == dimY) j1 = j-1;
-            j2 = j-1; if (j2 < 0) j2 = j+1;
+            i1 = i+1; if (i1 == dimX) i1 = i-1;
+            i2 = i-1; if (i2 < 0) i2 = i+1;
             index = j*dimX+i;
-
-                e = Output[j*dimX+i1];
-                w = Output[j*dimX+i2];
-                n = Output[j1*dimX+i];
-                s = Output[j2*dimX+i];
-
-                e1 = e - Output[index];
-                w1 = w - Output[index];
-                n1 = n - Output[index];
-                s1 = s - Output[index];
-
+            
+            e = Output[j*dimX+i1];
+            w = Output[j*dimX+i2];
+            n = Output[j1*dimX+i];
+            s = Output[j2*dimX+i];
+            
+            e1 = e - Output[index];
+            w1 = w - Output[index];
+            n1 = n - Output[index];
+            s1 = s - Output[index];
+            
             if (penaltytype == 1){
-            /* Huber penalty */
-            if (fabs(e1) > sigmaPar) e1 =  signNDFc(e1);
-            else e1 = e1/sigmaPar;
-
-            if (fabs(w1) > sigmaPar) w1 =  signNDFc(w1);
-            else w1 = w1/sigmaPar;
-
-            if (fabs(n1) > sigmaPar) n1 =  signNDFc(n1);
-            else n1 = n1/sigmaPar;
-
-            if (fabs(s1) > sigmaPar) s1 =  signNDFc(s1);
-            else s1 = s1/sigmaPar;
+                /* Huber penalty */
+                if (fabs(e1) > sigmaPar) e1 =  signNDFc(e1);
+                else e1 = e1/sigmaPar;
+                
+                if (fabs(w1) > sigmaPar) w1 =  signNDFc(w1);
+                else w1 = w1/sigmaPar;
+                
+                if (fabs(n1) > sigmaPar) n1 =  signNDFc(n1);
+                else n1 = n1/sigmaPar;
+                
+                if (fabs(s1) > sigmaPar) s1 =  signNDFc(s1);
+                else s1 = s1/sigmaPar;
             }
             else if (penaltytype == 2) {
-            /* Perona-Malik */
-            e1 = (e1)/(1.0f + powf((e1/sigmaPar),2));
-            w1 = (w1)/(1.0f + powf((w1/sigmaPar),2));
-            n1 = (n1)/(1.0f + powf((n1/sigmaPar),2));
-            s1 = (s1)/(1.0f + powf((s1/sigmaPar),2));
+                /* Perona-Malik */
+                e1 = (e1)/(1.0f + powf((e1/sigmaPar),2));
+                w1 = (w1)/(1.0f + powf((w1/sigmaPar),2));
+                n1 = (n1)/(1.0f + powf((n1/sigmaPar),2));
+                s1 = (s1)/(1.0f + powf((s1/sigmaPar),2));
             }
             else if (penaltytype == 3) {
-            /* Tukey Biweight */
-            if (fabs(e1) <= sigmaPar) e1 =  e1*powf((1.0f - powf((e1/sigmaPar),2)), 2);
-            else e1 = 0.0f;
-            if (fabs(w1) <= sigmaPar) w1 =  w1*powf((1.0f - powf((w1/sigmaPar),2)), 2);
-            else w1 = 0.0f;
-            if (fabs(n1) <= sigmaPar) n1 =  n1*powf((1.0f - powf((n1/sigmaPar),2)), 2);
-            else n1 = 0.0f;
-            if (fabs(s1) <= sigmaPar) s1 =  s1*powf((1.0f - powf((s1/sigmaPar),2)), 2);
-            else s1 = 0.0f;
+                /* Tukey Biweight */
+                if (fabs(e1) <= sigmaPar) e1 =  e1*powf((1.0f - powf((e1/sigmaPar),2)), 2);
+                else e1 = 0.0f;
+                if (fabs(w1) <= sigmaPar) w1 =  w1*powf((1.0f - powf((w1/sigmaPar),2)), 2);
+                else w1 = 0.0f;
+                if (fabs(n1) <= sigmaPar) n1 =  n1*powf((1.0f - powf((n1/sigmaPar),2)), 2);
+                else n1 = 0.0f;
+                if (fabs(s1) <= sigmaPar) s1 =  s1*powf((1.0f - powf((s1/sigmaPar),2)), 2);
+                else s1 = 0.0f;
             }
             else {
-				printf("%s \n", "No penalty function selected! Use 1,2 or 3.");
-				break;
-				}
-           Output[index] += tau*(lambdaPar*(e1 + w1 + n1 + s1) - (Output[index] - Input[index]));
-		}}
-	return *Output;
+                printf("%s \n", "No penalty function selected! Use 1,2 or 3.");
+                break;
+            }
+            Output[index] += tau*(lambdaPar*(e1 + w1 + n1 + s1) - (Output[index] - Input[index]));
+        }}
+    return *Output;
 }
 /********************************************************************/
 /***************************3D Functions*****************************/
@@ -210,125 +209,125 @@ float NonLinearDiff2D(float *Input, float *Output, float lambdaPar, float sigmaP
 /* linear diffusion (heat equation) */
 float LinearDiff3D(float *Input, float *Output, float lambdaPar, float tau, long dimX, long dimY, long dimZ)
 {
-	long i,j,k,i1,i2,j1,j2,k1,k2,index;
-	float e,w,n,s,u,d,e1,w1,n1,s1,u1,d1;
-
+    long i,j,k,i1,i2,j1,j2,k1,k2,index;
+    float e,w,n,s,u,d,e1,w1,n1,s1,u1,d1;
+    
 #pragma omp parallel for shared(Input) private(index,i,j,i1,i2,j1,j2,e,w,n,s,e1,w1,n1,s1,k,k1,k2,u1,d1,u,d)
-for(k=0; k<dimZ; k++) {
-	k1 = k+1; if (k1 == dimZ) k1 = k-1;
-    k2 = k-1; if (k2 < 0) k2 = k+1;
-    for(i=0; i<dimX; i++) {
-        /* symmetric boundary conditions (Neuman) */
-        i1 = i+1; if (i1 == dimX) i1 = i-1;
-        i2 = i-1; if (i2 < 0) i2 = i+1;
+    for(k=0; k<dimZ; k++) {
+        k1 = k+1; if (k1 == dimZ) k1 = k-1;
+        k2 = k-1; if (k2 < 0) k2 = k+1;
         for(j=0; j<dimY; j++) {
             /* symmetric boundary conditions (Neuman) */
             j1 = j+1; if (j1 == dimY) j1 = j-1;
             j2 = j-1; if (j2 < 0) j2 = j+1;
-            index = (dimX*dimY)*k + j*dimX+i;
-
+            for(i=0; i<dimX; i++) {
+                /* symmetric boundary conditions (Neuman) */
+                i1 = i+1; if (i1 == dimX) i1 = i-1;
+                i2 = i-1; if (i2 < 0) i2 = i+1;
+                index = (dimX*dimY)*k + j*dimX+i;
+                
                 e = Output[(dimX*dimY)*k + j*dimX+i1];
                 w = Output[(dimX*dimY)*k + j*dimX+i2];
                 n = Output[(dimX*dimY)*k + j1*dimX+i];
                 s = Output[(dimX*dimY)*k + j2*dimX+i];
                 u = Output[(dimX*dimY)*k1 + j*dimX+i];
                 d = Output[(dimX*dimY)*k2 + j*dimX+i];
-
+                
                 e1 = e - Output[index];
                 w1 = w - Output[index];
                 n1 = n - Output[index];
                 s1 = s - Output[index];
                 u1 = u - Output[index];
                 d1 = d - Output[index];
-
+                
                 Output[index] += tau*(lambdaPar*(e1 + w1 + n1 + s1 + u1 + d1) - (Output[index] - Input[index]));
-		}}}
-	return *Output;
+            }}}
+    return *Output;
 }
 
 float NonLinearDiff3D(float *Input, float *Output, float lambdaPar, float sigmaPar, float tau, int penaltytype, long dimX, long dimY, long dimZ)
 {
-	long i,j,k,i1,i2,j1,j2,k1,k2,index;
-	float e,w,n,s,u,d,e1,w1,n1,s1,u1,d1;
-
+    long i,j,k,i1,i2,j1,j2,k1,k2,index;
+    float e,w,n,s,u,d,e1,w1,n1,s1,u1,d1;
+    
 #pragma omp parallel for shared(Input) private(index,i,j,i1,i2,j1,j2,e,w,n,s,e1,w1,n1,s1,k,k1,k2,u1,d1,u,d)
-for(k=0; k<dimZ; k++) {
-	k1 = k+1; if (k1 == dimZ) k1 = k-1;
-    k2 = k-1; if (k2 < 0) k2 = k+1;
-    for(i=0; i<dimX; i++) {
-        /* symmetric boundary conditions (Neuman) */
-        i1 = i+1; if (i1 == dimX) i1 = i-1;
-        i2 = i-1; if (i2 < 0) i2 = i+1;
+    for(k=0; k<dimZ; k++) {
+        k1 = k+1; if (k1 == dimZ) k1 = k-1;
+        k2 = k-1; if (k2 < 0) k2 = k+1;
         for(j=0; j<dimY; j++) {
             /* symmetric boundary conditions (Neuman) */
             j1 = j+1; if (j1 == dimY) j1 = j-1;
             j2 = j-1; if (j2 < 0) j2 = j+1;
-            index = (dimX*dimY)*k + j*dimX+i;
-
+            for(i=0; i<dimX; i++) {
+                /* symmetric boundary conditions (Neuman) */
+                i1 = i+1; if (i1 == dimX) i1 = i-1;
+                i2 = i-1; if (i2 < 0) i2 = i+1;
+                index = (dimX*dimY)*k + j*dimX+i;
+                
                 e = Output[(dimX*dimY)*k + j*dimX+i1];
                 w = Output[(dimX*dimY)*k + j*dimX+i2];
                 n = Output[(dimX*dimY)*k + j1*dimX+i];
                 s = Output[(dimX*dimY)*k + j2*dimX+i];
                 u = Output[(dimX*dimY)*k1 + j*dimX+i];
                 d = Output[(dimX*dimY)*k2 + j*dimX+i];
-
+                
                 e1 = e - Output[index];
                 w1 = w - Output[index];
                 n1 = n - Output[index];
                 s1 = s - Output[index];
                 u1 = u - Output[index];
                 d1 = d - Output[index];
-
-             if (penaltytype == 1){
-            /* Huber penalty */
-            if (fabs(e1) > sigmaPar) e1 =  signNDFc(e1);
-            else e1 = e1/sigmaPar;
-
-            if (fabs(w1) > sigmaPar) w1 =  signNDFc(w1);
-            else w1 = w1/sigmaPar;
-
-            if (fabs(n1) > sigmaPar) n1 =  signNDFc(n1);
-            else n1 = n1/sigmaPar;
-
-            if (fabs(s1) > sigmaPar) s1 =  signNDFc(s1);
-            else s1 = s1/sigmaPar;
-
-            if (fabs(u1) > sigmaPar) u1 =  signNDFc(u1);
-            else u1 = u1/sigmaPar;
-
-            if (fabs(d1) > sigmaPar) d1 =  signNDFc(d1);
-            else d1 = d1/sigmaPar;
-            }
-            else if (penaltytype == 2) {
-            /* Perona-Malik */
-            e1 = (e1)/(1.0f + powf((e1/sigmaPar),2));
-            w1 = (w1)/(1.0f + powf((w1/sigmaPar),2));
-            n1 = (n1)/(1.0f + powf((n1/sigmaPar),2));
-            s1 = (s1)/(1.0f + powf((s1/sigmaPar),2));
-            u1 = (u1)/(1.0f + powf((u1/sigmaPar),2));
-            d1 = (d1)/(1.0f + powf((d1/sigmaPar),2));
-            }
-            else if (penaltytype == 3) {
-            /* Tukey Biweight */
-            if (fabs(e1) <= sigmaPar) e1 =  e1*powf((1.0f - powf((e1/sigmaPar),2)), 2);
-            else e1 = 0.0f;
-            if (fabs(w1) <= sigmaPar) w1 =  w1*powf((1.0f - powf((w1/sigmaPar),2)), 2);
-            else w1 = 0.0f;
-            if (fabs(n1) <= sigmaPar) n1 =  n1*powf((1.0f - powf((n1/sigmaPar),2)), 2);
-            else n1 = 0.0f;
-            if (fabs(s1) <= sigmaPar) s1 =  s1*powf((1.0f - powf((s1/sigmaPar),2)), 2);
-            else s1 = 0.0f;
-            if (fabs(u1) <= sigmaPar) u1 =  u1*powf((1.0f - powf((u1/sigmaPar),2)), 2);
-            else u1 = 0.0f;
-            if (fabs(d1) <= sigmaPar) d1 =  d1*powf((1.0f - powf((d1/sigmaPar),2)), 2);
-            else d1 = 0.0f;
-            }
-            else {
-				printf("%s \n", "No penalty function selected! Use 1,2 or 3.");
-				break;
-				}
-
+                
+                if (penaltytype == 1){
+                    /* Huber penalty */
+                    if (fabs(e1) > sigmaPar) e1 =  signNDFc(e1);
+                    else e1 = e1/sigmaPar;
+                    
+                    if (fabs(w1) > sigmaPar) w1 =  signNDFc(w1);
+                    else w1 = w1/sigmaPar;
+                    
+                    if (fabs(n1) > sigmaPar) n1 =  signNDFc(n1);
+                    else n1 = n1/sigmaPar;
+                    
+                    if (fabs(s1) > sigmaPar) s1 =  signNDFc(s1);
+                    else s1 = s1/sigmaPar;
+                    
+                    if (fabs(u1) > sigmaPar) u1 =  signNDFc(u1);
+                    else u1 = u1/sigmaPar;
+                    
+                    if (fabs(d1) > sigmaPar) d1 =  signNDFc(d1);
+                    else d1 = d1/sigmaPar;
+                }
+                else if (penaltytype == 2) {
+                    /* Perona-Malik */
+                    e1 = (e1)/(1.0f + powf((e1/sigmaPar),2));
+                    w1 = (w1)/(1.0f + powf((w1/sigmaPar),2));
+                    n1 = (n1)/(1.0f + powf((n1/sigmaPar),2));
+                    s1 = (s1)/(1.0f + powf((s1/sigmaPar),2));
+                    u1 = (u1)/(1.0f + powf((u1/sigmaPar),2));
+                    d1 = (d1)/(1.0f + powf((d1/sigmaPar),2));
+                }
+                else if (penaltytype == 3) {
+                    /* Tukey Biweight */
+                    if (fabs(e1) <= sigmaPar) e1 =  e1*powf((1.0f - powf((e1/sigmaPar),2)), 2);
+                    else e1 = 0.0f;
+                    if (fabs(w1) <= sigmaPar) w1 =  w1*powf((1.0f - powf((w1/sigmaPar),2)), 2);
+                    else w1 = 0.0f;
+                    if (fabs(n1) <= sigmaPar) n1 =  n1*powf((1.0f - powf((n1/sigmaPar),2)), 2);
+                    else n1 = 0.0f;
+                    if (fabs(s1) <= sigmaPar) s1 =  s1*powf((1.0f - powf((s1/sigmaPar),2)), 2);
+                    else s1 = 0.0f;
+                    if (fabs(u1) <= sigmaPar) u1 =  u1*powf((1.0f - powf((u1/sigmaPar),2)), 2);
+                    else u1 = 0.0f;
+                    if (fabs(d1) <= sigmaPar) d1 =  d1*powf((1.0f - powf((d1/sigmaPar),2)), 2);
+                    else d1 = 0.0f;
+                }
+                else {
+                    printf("%s \n", "No penalty function selected! Use 1,2 or 3.");
+                    break;
+                }
+                
                 Output[index] += tau*(lambdaPar*(e1 + w1 + n1 + s1 + u1 + d1) - (Output[index] - Input[index]));
-		}}}
-	return *Output;
+            }}}
+    return *Output;
 }
diff --git a/src/Core/regularisers_CPU/FGP_TV_core.c b/src/Core/regularisers_CPU/FGP_TV_core.c
index 69c92bc..a16a2e5 100644
--- a/src/Core/regularisers_CPU/FGP_TV_core.c
+++ b/src/Core/regularisers_CPU/FGP_TV_core.c
@@ -1,21 +1,21 @@
 /*
-This work is part of the Core Imaging Library developed by
-Visual Analytics and Imaging System Group of the Science Technology
-Facilities Council, STFC
-
-Copyright 2019 Daniil Kazantsev
-Copyright 2019 Srikanth Nagella, Edoardo Pasca
-
-Licensed under the Apache License, Version 2.0 (the "License");
-you may not use this file except in compliance with the License.
-You may obtain a copy of the License at
-http://www.apache.org/licenses/LICENSE-2.0
-Unless required by applicable law or agreed to in writing, software
-distributed under the License is distributed on an "AS IS" BASIS,
-WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-See the License for the specific language governing permissions and
-limitations under the License.
-*/
+ * This work is part of the Core Imaging Library developed by
+ * Visual Analytics and Imaging System Group of the Science Technology
+ * Facilities Council, STFC
+ *
+ * Copyright 2019 Daniil Kazantsev
+ * Copyright 2019 Srikanth Nagella, Edoardo Pasca
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ * http://www.apache.org/licenses/LICENSE-2.0
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
 
 #include "FGP_TV_core.h"
 
@@ -46,66 +46,66 @@ float TV_FGP_CPU_main(float *Input, float *Output, float *infovector, float lamb
     float tk = 1.0f;
     float tkp1 =1.0f;
     int count = 0;
-
-	if (dimZ <= 1) {
-	/*2D case */
-  	float *Output_prev=NULL, *P1=NULL, *P2=NULL, *P1_prev=NULL, *P2_prev=NULL, *R1=NULL, *R2=NULL;
-	DimTotal = (long)(dimX*dimY);
-
-	      if (epsil != 0.0f) Output_prev = calloc(DimTotal, sizeof(float));
+    
+    if (dimZ <= 1) {
+        /*2D case */
+        float *Output_prev=NULL, *P1=NULL, *P2=NULL, *P1_prev=NULL, *P2_prev=NULL, *R1=NULL, *R2=NULL;
+        DimTotal = (long)(dimX*dimY);
+        
+        if (epsil != 0.0f) Output_prev = calloc(DimTotal, sizeof(float));
         P1 = calloc(DimTotal, sizeof(float));
         P2 = calloc(DimTotal, sizeof(float));
         P1_prev = calloc(DimTotal, sizeof(float));
         P2_prev = calloc(DimTotal, sizeof(float));
         R1 = calloc(DimTotal, sizeof(float));
         R2 = calloc(DimTotal, sizeof(float));
-
-	/* begin iterations */
+        
+        /* begin iterations */
         for(ll=0; ll<iterationsNumb; ll++) {
-
+            
             if ((epsil != 0.0f)  && (ll % 5 == 0)) copyIm(Output, Output_prev, (long)(dimX), (long)(dimY), 1l);
             /* computing the gradient of the objective function */
             Obj_func2D(Input, Output, R1, R2, lambdaPar, (long)(dimX), (long)(dimY));
-
+            
             /* apply nonnegativity */
             if (nonneg == 1) for(j=0; j<DimTotal; j++) {if (Output[j] < 0.0f) Output[j] = 0.0f;}
-
+            
             /*Taking a step towards minus of the gradient*/
             Grad_func2D(P1, P2, Output, R1, R2, lambdaPar, (long)(dimX), (long)(dimY));
-
+            
             /* projection step */
             Proj_func2D(P1, P2, methodTV, DimTotal);
-
+            
             /*updating R and t*/
             tkp1 = (1.0f + sqrtf(1.0f + 4.0f*tk*tk))*0.5f;
             Rupd_func2D(P1, P1_prev, P2, P2_prev, R1, R2, tkp1, tk, DimTotal);
-
+            
             /*storing old values*/
             copyIm(P1, P1_prev, (long)(dimX), (long)(dimY), 1l);
             copyIm(P2, P2_prev, (long)(dimX), (long)(dimY), 1l);
             tk = tkp1;
-
+            
             /* check early stopping criteria */
             if ((epsil != 0.0f)  && (ll % 5 == 0)) {
-            re = 0.0f; re1 = 0.0f;
-	            for(j=0; j<DimTotal; j++)
-        	    {
-        	        re += powf(Output[j] - Output_prev[j],2);
-        	        re1 += powf(Output[j],2);
-        	    }
-           re = sqrtf(re)/sqrtf(re1);
-           if (re < epsil)  count++;
-           if (count > 3) break;
+                re = 0.0f; re1 = 0.0f;
+                for(j=0; j<DimTotal; j++)
+                {
+                    re += powf(Output[j] - Output_prev[j],2);
+                    re1 += powf(Output[j],2);
+                }
+                re = sqrtf(re)/sqrtf(re1);
+                if (re < epsil)  count++;
+                if (count > 3) break;
             }
-     }
-	if (epsil != 0.0f) free(Output_prev);
-	free(P1); free(P2); free(P1_prev); free(P2_prev); free(R1); free(R2);
-	}
-	else {
-		/*3D case*/
-	float *Output_prev=NULL, *P1=NULL, *P2=NULL, *P3=NULL, *P1_prev=NULL, *P2_prev=NULL, *P3_prev=NULL, *R1=NULL, *R2=NULL, *R3=NULL;
+        }
+        if (epsil != 0.0f) free(Output_prev);
+        free(P1); free(P2); free(P1_prev); free(P2_prev); free(R1); free(R2);
+    }
+    else {
+        /*3D case*/
+        float *Output_prev=NULL, *P1=NULL, *P2=NULL, *P3=NULL, *P1_prev=NULL, *P2_prev=NULL, *P3_prev=NULL, *R1=NULL, *R2=NULL, *R3=NULL;
         DimTotal = (long)(dimX*dimY*dimZ);
-
+        
         if (epsil != 0.0f) Output_prev = calloc(DimTotal, sizeof(float));
         P1 = calloc(DimTotal, sizeof(float));
         P2 = calloc(DimTotal, sizeof(float));
@@ -116,58 +116,58 @@ float TV_FGP_CPU_main(float *Input, float *Output, float *infovector, float lamb
         R1 = calloc(DimTotal, sizeof(float));
         R2 = calloc(DimTotal, sizeof(float));
         R3 = calloc(DimTotal, sizeof(float));
-
-		    /* begin iterations */
+        
+        /* begin iterations */
         for(ll=0; ll<iterationsNumb; ll++) {
-
+            
             if ((epsil != 0.0f)  && (ll % 5 == 0)) copyIm(Output, Output_prev, (long)(dimX), (long)(dimY), (long)(dimZ));
-
+            
             /* computing the gradient of the objective function */
             Obj_func3D(Input, Output, R1, R2, R3, lambdaPar, (long)(dimX), (long)(dimY), (long)(dimZ));
-
+            
             /* apply nonnegativity */
             if (nonneg == 1) for(j=0; j<DimTotal; j++) {if (Output[j] < 0.0f) Output[j] = 0.0f;}
-
+            
             /*Taking a step towards minus of the gradient*/
             Grad_func3D(P1, P2, P3, Output, R1, R2, R3, lambdaPar, (long)(dimX), (long)(dimY), (long)(dimZ));
-
+            
             /* projection step */
             Proj_func3D(P1, P2, P3, methodTV, DimTotal);
-
+            
             /*updating R and t*/
             tkp1 = (1.0f + sqrtf(1.0f + 4.0f*tk*tk))*0.5f;
             Rupd_func3D(P1, P1_prev, P2, P2_prev, P3, P3_prev, R1, R2, R3, tkp1, tk, DimTotal);
-
+            
             /* calculate norm - stopping rules*/
             if ((epsil != 0.0f)  && (ll % 5 == 0)) {
-            re = 0.0f; re1 = 0.0f;
-	            for(j=0; j<DimTotal; j++)
-        	    {
-        	        re += powf(Output[j] - Output_prev[j],2);
-        	        re1 += powf(Output[j],2);
-        	    }
-            re = sqrtf(re)/sqrtf(re1);
-            /* stop if the norm residual is less than the tolerance EPS */
-            if (re < epsil)  count++;
-            if (count > 3) break;
+                re = 0.0f; re1 = 0.0f;
+                for(j=0; j<DimTotal; j++)
+                {
+                    re += powf(Output[j] - Output_prev[j],2);
+                    re1 += powf(Output[j],2);
+                }
+                re = sqrtf(re)/sqrtf(re1);
+                /* stop if the norm residual is less than the tolerance EPS */
+                if (re < epsil)  count++;
+                if (count > 3) break;
             }
-
+            
             /*storing old values*/
             copyIm(P1, P1_prev, (long)(dimX), (long)(dimY), (long)(dimZ));
             copyIm(P2, P2_prev, (long)(dimX), (long)(dimY), (long)(dimZ));
             copyIm(P3, P3_prev, (long)(dimX), (long)(dimY), (long)(dimZ));
             tk = tkp1;
         }
-
-	if (epsil != 0.0f) free(Output_prev);
-	free(P1); free(P2); free(P3); free(P1_prev); free(P2_prev); free(P3_prev); free(R1); free(R2); free(R3);
-	}
-
-       /*adding info into info_vector */
-        infovector[0] = (float)(ll);  /*iterations number (if stopped earlier based on tolerance)*/
-        infovector[1] = re;  /* reached tolerance */
-
-	return 0;
+        
+        if (epsil != 0.0f) free(Output_prev);
+        free(P1); free(P2); free(P3); free(P1_prev); free(P2_prev); free(P3_prev); free(R1); free(R2); free(R3);
+    }
+    
+    /*adding info into info_vector */
+    infovector[0] = (float)(ll);  /*iterations number (if stopped earlier based on tolerance)*/
+    infovector[1] = re;  /* reached tolerance */
+    
+    return 0;
 }
 
 float Obj_func2D(float *A, float *D, float *R1, float *R2, float lambda, long dimX, long dimY)
@@ -177,7 +177,7 @@ float Obj_func2D(float *A, float *D, float *R1, float *R2, float lambda, long di
 #pragma omp parallel for shared(A,D,R1,R2) private(index,i,j,val1,val2)
     for(j=0; j<dimY; j++) {
         for(i=0; i<dimX; i++) {
-			index = j*dimX+i;
+            index = j*dimX+i;
             /* boundary conditions  */
             if (i == 0) {val1 = 0.0f;} else {val1 = R1[j*dimX + (i-1)];}
             if (j == 0) {val2 = 0.0f;} else {val2 = R2[(j-1)*dimX + i];}
@@ -193,7 +193,7 @@ float Grad_func2D(float *P1, float *P2, float *D, float *R1, float *R2, float la
 #pragma omp parallel for shared(P1,P2,D,R1,R2,multip) private(index,i,j,val1,val2)
     for(j=0; j<dimY; j++) {
         for(i=0; i<dimX; i++) {
-			index = j*dimX+i;
+            index = j*dimX+i;
             /* boundary conditions */
             if (i == dimX-1) val1 = 0.0f; else val1 = D[index] - D[j*dimX + (i+1)];
             if (j == dimY-1) val2 = 0.0f; else val2 = D[index] - D[(j+1)*dimX + i];
@@ -210,25 +210,25 @@ float Proj_func2D(float *P1, float *P2, int methTV, long DimTotal)
         /* isotropic TV*/
 #pragma omp parallel for shared(P1,P2) private(i,denom,sq_denom)
         for(i=0; i<DimTotal; i++) {
-                denom = powf(P1[i],2) +  powf(P2[i],2);
-                if (denom > 1.0f) {
-					sq_denom = 1.0f/sqrtf(denom);
-                    P1[i] = P1[i]*sq_denom;
-                    P2[i] = P2[i]*sq_denom;
-                }
+            denom = powf(P1[i],2) +  powf(P2[i],2);
+            if (denom > 1.0f) {
+                sq_denom = 1.0f/sqrtf(denom);
+                P1[i] = P1[i]*sq_denom;
+                P2[i] = P2[i]*sq_denom;
             }
+        }
     }
     else {
         /* anisotropic TV*/
 #pragma omp parallel for shared(P1,P2) private(i,val1,val2)
         for(i=0; i<DimTotal; i++) {
-                val1 = fabs(P1[i]);
-                val2 = fabs(P2[i]);
-                if (val1 < 1.0f) {val1 = 1.0f;}
-                if (val2 < 1.0f) {val2 = 1.0f;}
-                P1[i] = P1[i]/val1;
-                P2[i] = P2[i]/val2;
-            }
+            val1 = fabs(P1[i]);
+            val2 = fabs(P2[i]);
+            if (val1 < 1.0f) {val1 = 1.0f;}
+            if (val2 < 1.0f) {val2 = 1.0f;}
+            P1[i] = P1[i]/val1;
+            P2[i] = P2[i]/val2;
+        }
     }
     return 1;
 }
@@ -239,9 +239,9 @@ float Rupd_func2D(float *P1, float *P1_old, float *P2, float *P2_old, float *R1,
     multip = ((tk-1.0f)/tkp1);
 #pragma omp parallel for shared(P1,P2,P1_old,P2_old,R1,R2,multip) private(i)
     for(i=0; i<DimTotal; i++) {
-            R1[i] = P1[i] + multip*(P1[i] - P1_old[i]);
-            R2[i] = P2[i] + multip*(P2[i] - P2_old[i]);
-        }
+        R1[i] = P1[i] + multip*(P1[i] - P1_old[i]);
+        R2[i] = P2[i] + multip*(P2[i] - P2_old[i]);
+    }
     return 1;
 }
 
@@ -255,7 +255,7 @@ float Obj_func3D(float *A, float *D, float *R1, float *R2, float *R3, float lamb
     for(k=0; k<dimZ; k++) {
         for(j=0; j<dimY; j++) {
             for(i=0; i<dimX; i++) {
-		index = (dimX*dimY)*k + j*dimX+i;
+                index = (dimX*dimY)*k + j*dimX+i;
                 /* boundary conditions */
                 if (i == 0) {val1 = 0.0f;} else {val1 = R1[(dimX*dimY)*k + j*dimX + (i-1)];}
                 if (j == 0) {val2 = 0.0f;} else {val2 = R2[(dimX*dimY)*k + (j-1)*dimX + i];}
@@ -273,7 +273,7 @@ float Grad_func3D(float *P1, float *P2, float *P3, float *D, float *R1, float *R
     for(k=0; k<dimZ; k++) {
         for(j=0; j<dimY; j++) {
             for(i=0; i<dimX; i++) {
-		index = (dimX*dimY)*k + j*dimX+i;
+                index = (dimX*dimY)*k + j*dimX+i;
                 /* boundary conditions */
                 if (i == dimX-1) val1 = 0.0f; else val1 = D[index] - D[(dimX*dimY)*k + j*dimX + (i+1)];
                 if (j == dimY-1) val2 = 0.0f; else val2 = D[index] - D[(dimX*dimY)*k + (j+1)*dimX + i];
@@ -289,33 +289,33 @@ float Proj_func3D(float *P1, float *P2, float *P3, int methTV, long DimTotal)
     float val1, val2, val3, denom, sq_denom;
     long i;
     if (methTV == 0) {
-	/* isotropic TV*/
-	#pragma omp parallel for shared(P1,P2,P3) private(i,val1,val2,val3,sq_denom)
-    for(i=0; i<DimTotal; i++) {
-				denom = powf(P1[i],2) + powf(P2[i],2) + powf(P3[i],2);
-                if (denom > 1.0f) {
-					sq_denom = 1.0f/sqrtf(denom);
-                    P1[i] = P1[i]*sq_denom;
-                    P2[i] = P2[i]*sq_denom;
-                    P3[i] = P3[i]*sq_denom;
-                }
-			}
-	}
+        /* isotropic TV*/
+#pragma omp parallel for shared(P1,P2,P3) private(i,val1,val2,val3,sq_denom)
+        for(i=0; i<DimTotal; i++) {
+            denom = powf(P1[i],2) + powf(P2[i],2) + powf(P3[i],2);
+            if (denom > 1.0f) {
+                sq_denom = 1.0f/sqrtf(denom);
+                P1[i] = P1[i]*sq_denom;
+                P2[i] = P2[i]*sq_denom;
+                P3[i] = P3[i]*sq_denom;
+            }
+        }
+    }
     else {
-    /* anisotropic TV*/
+        /* anisotropic TV*/
 #pragma omp parallel for shared(P1,P2,P3) private(i,val1,val2,val3)
-    for(i=0; i<DimTotal; i++) {
-                val1 = fabs(P1[i]);
-                val2 = fabs(P2[i]);
-                val3 = fabs(P3[i]);
-                if (val1 < 1.0f) {val1 = 1.0f;}
-                if (val2 < 1.0f) {val2 = 1.0f;}
-                if (val3 < 1.0f) {val3 = 1.0f;}
-                P1[i] = P1[i]/val1;
-                P2[i] = P2[i]/val2;
-                P3[i] = P3[i]/val3;
-            }
-		}
+        for(i=0; i<DimTotal; i++) {
+            val1 = fabs(P1[i]);
+            val2 = fabs(P2[i]);
+            val3 = fabs(P3[i]);
+            if (val1 < 1.0f) {val1 = 1.0f;}
+            if (val2 < 1.0f) {val2 = 1.0f;}
+            if (val3 < 1.0f) {val3 = 1.0f;}
+            P1[i] = P1[i]/val1;
+            P2[i] = P2[i]/val2;
+            P3[i] = P3[i]/val3;
+        }
+    }
     return 1;
 }
 float Rupd_func3D(float *P1, float *P1_old, float *P2, float *P2_old, float *P3, float *P3_old, float *R1, float *R2, float *R3, float tkp1, float tk, long DimTotal)
@@ -325,9 +325,9 @@ float Rupd_func3D(float *P1, float *P1_old, float *P2, float *P2_old, float *P3,
     multip = ((tk-1.0f)/tkp1);
 #pragma omp parallel for shared(P1,P2,P3,P1_old,P2_old,P3_old,R1,R2,R3,multip) private(i)
     for(i=0; i<DimTotal; i++) {
-                R1[i] = P1[i] + multip*(P1[i] - P1_old[i]);
-                R2[i] = P2[i] + multip*(P2[i] - P2_old[i]);
-                R3[i] = P3[i] + multip*(P3[i] - P3_old[i]);
-            }
+        R1[i] = P1[i] + multip*(P1[i] - P1_old[i]);
+        R2[i] = P2[i] + multip*(P2[i] - P2_old[i]);
+        R3[i] = P3[i] + multip*(P3[i] - P3_old[i]);
+    }
     return 1;
 }
diff --git a/src/Core/regularisers_CPU/FGP_dTV_core.c b/src/Core/regularisers_CPU/FGP_dTV_core.c
index 4e1e38c..afd7264 100644
--- a/src/Core/regularisers_CPU/FGP_dTV_core.c
+++ b/src/Core/regularisers_CPU/FGP_dTV_core.c
@@ -1,21 +1,21 @@
 /*
-This work is part of the Core Imaging Library developed by
-Visual Analytics and Imaging System Group of the Science Technology
-Facilities Council, STFC
-
-Copyright 2019 Daniil Kazantsev
-Copyright 2019 Srikanth Nagella, Edoardo Pasca
-
-Licensed under the Apache License, Version 2.0 (the "License");
-you may not use this file except in compliance with the License.
-You may obtain a copy of the License at
-http://www.apache.org/licenses/LICENSE-2.0
-Unless required by applicable law or agreed to in writing, software
-distributed under the License is distributed on an "AS IS" BASIS,
-WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-See the License for the specific language governing permissions and
-limitations under the License.
-*/
+ * This work is part of the Core Imaging Library developed by
+ * Visual Analytics and Imaging System Group of the Science Technology
+ * Facilities Council, STFC
+ *
+ * Copyright 2019 Daniil Kazantsev
+ * Copyright 2019 Srikanth Nagella, Edoardo Pasca
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ * http://www.apache.org/licenses/LICENSE-2.0
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
 
 #include "FGP_dTV_core.h"
 
@@ -45,18 +45,18 @@ limitations under the License.
 
 float dTV_FGP_CPU_main(float *Input, float *InputRef, float *Output, float *infovector, float lambdaPar, int iterationsNumb, float epsil, float eta, int methodTV, int nonneg, int dimX, int dimY, int dimZ)
 {
-	  int ll;
+    int ll;
     long j, DimTotal;
-	  float re, re1;
+    float re, re1;
     re = 0.0f; re1 = 0.0f;
-	  float tk = 1.0f;
+    float tk = 1.0f;
     float tkp1=1.0f;
     int count = 0;
-
-
+    
+    
     float *Output_prev=NULL, *P1=NULL, *P2=NULL, *P1_prev=NULL, *P2_prev=NULL, *R1=NULL, *R2=NULL, *InputRef_x=NULL, *InputRef_y=NULL;
     DimTotal = (long)(dimX*dimY*dimZ);
-
+    
     if (epsil != 0.0f) Output_prev = calloc(DimTotal, sizeof(float));
     P1 = calloc(DimTotal, sizeof(float));
     P2 = calloc(DimTotal, sizeof(float));
@@ -66,119 +66,119 @@ float dTV_FGP_CPU_main(float *Input, float *InputRef, float *Output, float *info
     R2 = calloc(DimTotal, sizeof(float));
     InputRef_x = calloc(DimTotal, sizeof(float));
     InputRef_y = calloc(DimTotal, sizeof(float));
-
-	if (dimZ <= 1) {
-    		/*2D case */
+    
+    if (dimZ <= 1) {
+        /*2D case */
         /* calculate gradient field (smoothed) for the reference image */
-		GradNorm_func2D(InputRef, InputRef_x, InputRef_y, eta, (long)(dimX), (long)(dimY));
-
-		/* begin iterations */
+        GradNorm_func2D(InputRef, InputRef_x, InputRef_y, eta, (long)(dimX), (long)(dimY));
+        
+        /* begin iterations */
         for(ll=0; ll<iterationsNumb; ll++) {
-
+            
             if ((epsil != 0.0f)  && (ll % 5 == 0)) copyIm(Output, Output_prev, (long)(dimX), (long)(dimY), 1l);
             /*projects a 2D vector field R-1,2 onto the orthogonal complement of another 2D vector field InputRef_xy*/
             ProjectVect_func2D(R1, R2, InputRef_x, InputRef_y, (long)(dimX), (long)(dimY));
-
+            
             /* computing the gradient of the objective function */
             Obj_dfunc2D(Input, Output, R1, R2, lambdaPar, (long)(dimX), (long)(dimY));
-
+            
             /* apply nonnegativity */
             if (nonneg == 1) for(j=0; j<DimTotal; j++) {if (Output[j] < 0.0f) Output[j] = 0.0f;}
-
+            
             /*Taking a step towards minus of the gradient*/
             Grad_dfunc2D(P1, P2, Output, R1, R2, InputRef_x, InputRef_y, lambdaPar, (long)(dimX), (long)(dimY));
-
+            
             /* projection step */
             Proj_dfunc2D(P1, P2, methodTV, DimTotal);
-
+            
             /*updating R and t*/
             tkp1 = (1.0f + sqrt(1.0f + 4.0f*tk*tk))*0.5f;
             Rupd_dfunc2D(P1, P1_prev, P2, P2_prev, R1, R2, tkp1, tk, DimTotal);
-
+            
             copyIm(P1, P1_prev, (long)(dimX), (long)(dimY), 1l);
             copyIm(P2, P2_prev, (long)(dimX), (long)(dimY), 1l);
             tk = tkp1;
-
+            
             /* check early stopping criteria */
             if ((epsil != 0.0f)  && (ll % 5 == 0)) {
-            re = 0.0f; re1 = 0.0f;
-              for(j=0; j<DimTotal; j++)
-              {
-                  re += powf(Output[j] - Output_prev[j],2);
-                  re1 += powf(Output[j],2);
-              }
-           re = sqrtf(re)/sqrtf(re1);
-           if (re < epsil)  count++;
-           if (count > 3) break;
+                re = 0.0f; re1 = 0.0f;
+                for(j=0; j<DimTotal; j++)
+                {
+                    re += powf(Output[j] - Output_prev[j],2);
+                    re1 += powf(Output[j],2);
+                }
+                re = sqrtf(re)/sqrtf(re1);
+                if (re < epsil)  count++;
+                if (count > 3) break;
             }
         }
-	}
-	else {
-		/*3D case*/
-		float *P3=NULL, *P3_prev=NULL, *R3=NULL, *InputRef_z=NULL;
-
+    }
+    else {
+        /*3D case*/
+        float *P3=NULL, *P3_prev=NULL, *R3=NULL, *InputRef_z=NULL;
+        
         P3 = calloc(DimTotal, sizeof(float));
         P3_prev = calloc(DimTotal, sizeof(float));
         R3 = calloc(DimTotal, sizeof(float));
         InputRef_z = calloc(DimTotal, sizeof(float));
-
-		/* calculate gradient field (smoothed) for the reference volume */
-		GradNorm_func3D(InputRef, InputRef_x, InputRef_y, InputRef_z, eta, (long)(dimX), (long)(dimY), (long)(dimZ));
-
-		/* begin iterations */
+        
+        /* calculate gradient field (smoothed) for the reference volume */
+        GradNorm_func3D(InputRef, InputRef_x, InputRef_y, InputRef_z, eta, (long)(dimX), (long)(dimY), (long)(dimZ));
+        
+        /* begin iterations */
         for(ll=0; ll<iterationsNumb; ll++) {
-
+            
             if ((epsil != 0.0f)  && (ll % 5 == 0)) copyIm(Output, Output_prev, (long)(dimX), (long)(dimY), (long)(dimZ));
-
-					 /*projects a 3D vector field R-1,2,3 onto the orthogonal complement of another 3D vector field InputRef_xyz*/
+            
+            /*projects a 3D vector field R-1,2,3 onto the orthogonal complement of another 3D vector field InputRef_xyz*/
             ProjectVect_func3D(R1, R2, R3, InputRef_x, InputRef_y, InputRef_z, (long)(dimX), (long)(dimY), (long)(dimZ));
-
+            
             /* computing the gradient of the objective function */
             Obj_dfunc3D(Input, Output, R1, R2, R3, lambdaPar, (long)(dimX), (long)(dimY), (long)(dimZ));
-
+            
             /* apply nonnegativity */
             if (nonneg == 1) for(j=0; j<DimTotal; j++) {if (Output[j] < 0.0f) Output[j] = 0.0f;}
-
+            
             /*Taking a step towards minus of the gradient*/
             Grad_dfunc3D(P1, P2, P3, Output, R1, R2, R3, InputRef_x, InputRef_y, InputRef_z, lambdaPar, (long)(dimX), (long)(dimY), (long)(dimZ));
-
+            
             /* projection step */
             Proj_dfunc3D(P1, P2, P3, methodTV, DimTotal);
-
+            
             /*updating R and t*/
             tkp1 = (1.0f + sqrt(1.0f + 4.0f*tk*tk))*0.5f;
             Rupd_dfunc3D(P1, P1_prev, P2, P2_prev, P3, P3_prev, R1, R2, R3, tkp1, tk, DimTotal);
-
-            /*storing old values*/            
+            
+            /*storing old values*/
             copyIm(P1, P1_prev, (long)(dimX), (long)(dimY), (long)(dimZ));
             copyIm(P2, P2_prev, (long)(dimX), (long)(dimY), (long)(dimZ));
             copyIm(P3, P3_prev, (long)(dimX), (long)(dimY), (long)(dimZ));
             tk = tkp1;
-
+            
             /* check early stopping criteria */
             if ((epsil != 0.0f)  && (ll % 5 == 0)) {
-            re = 0.0f; re1 = 0.0f;
-              for(j=0; j<DimTotal; j++)
-              {
-                  re += powf(Output[j] - Output_prev[j],2);
-                  re1 += powf(Output[j],2);
-              }
-           re = sqrtf(re)/sqrtf(re1);
-           if (re < epsil)  count++;
-           if (count > 3) break;
+                re = 0.0f; re1 = 0.0f;
+                for(j=0; j<DimTotal; j++)
+                {
+                    re += powf(Output[j] - Output_prev[j],2);
+                    re1 += powf(Output[j],2);
+                }
+                re = sqrtf(re)/sqrtf(re1);
+                if (re < epsil)  count++;
+                if (count > 3) break;
             }
         }
-
-		free(P3); free(P3_prev); free(R3); free(InputRef_z);
-	}
-  if (epsil != 0.0f) free(Output_prev);
-  free(P1); free(P2); free(P1_prev); free(P2_prev); free(R1); free(R2); free(InputRef_x); free(InputRef_y);
-
-   /*adding info into info_vector */
-   infovector[0] = (float)(ll);  /*iterations number (if stopped earlier based on tolerance)*/
-   infovector[1] = re;  /* reached tolerance */
-
-  return 0;
+        
+        free(P3); free(P3_prev); free(R3); free(InputRef_z);
+    }
+    if (epsil != 0.0f) free(Output_prev);
+    free(P1); free(P2); free(P1_prev); free(P2_prev); free(R1); free(R2); free(InputRef_x); free(InputRef_y);
+    
+    /*adding info into info_vector */
+    infovector[0] = (float)(ll);  /*iterations number (if stopped earlier based on tolerance)*/
+    infovector[1] = re;  /* reached tolerance */
+    
+    return 0;
 }
 
 
@@ -191,9 +191,9 @@ float GradNorm_func2D(float *B, float *B_x, float *B_y, float eta, long dimX, lo
     long i,j,index;
     float val1, val2, gradX, gradY, magn;
 #pragma omp parallel for shared(B, B_x, B_y) private(i,j,index,val1,val2,gradX,gradY,magn)
-    for(i=0; i<dimX; i++) {
-        for(j=0; j<dimY; j++) {
-			index = j*dimX+i;
+    for(j=0; j<dimY; j++) {
+        for(i=0; i<dimX; i++) {
+            index = j*dimX+i;
             /* zero boundary conditions */
             if (i == dimX-1) {val1 = 0.0f;} else {val1 = B[j*dimX + (i+1)];}
             if (j == dimY-1) {val2 = 0.0f;} else {val2 = B[(j+1)*dimX + i];}
@@ -212,9 +212,9 @@ float ProjectVect_func2D(float *R1, float *R2, float *B_x, float *B_y, long dimX
     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++) {
-        for(j=0; j<dimY; j++) {
-			index = j*dimX+i;
+    for(j=0; j<dimY; j++) {
+        for(i=0; i<dimX; i++) {
+            index = j*dimX+i;
             in_prod = R1[index]*B_x[index] + R2[index]*B_y[index];   /* calculate inner product */
             R1[index] = R1[index] - in_prod*B_x[index];
             R2[index] = R2[index] - in_prod*B_y[index];
@@ -227,9 +227,9 @@ float Obj_dfunc2D(float *A, float *D, float *R1, float *R2, float lambda, long d
     float val1, val2;
     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++) {
-			index = j*dimX+i;
+    for(j=0; j<dimY; j++) {
+        for(i=0; i<dimX; i++) {
+            index = j*dimX+i;
             /* boundary conditions  */
             if (i == 0) {val1 = 0.0f;} else {val1 = R1[j*dimX + (i-1)];}
             if (j == 0) {val2 = 0.0f;} else {val2 = R2[(j-1)*dimX + i];}
@@ -243,20 +243,20 @@ float Grad_dfunc2D(float *P1, float *P2, float *D, float *R1, float *R2, float *
     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++) {
-        for(j=0; j<dimY; j++) {
-			index = j*dimX+i;
+    for(j=0; j<dimY; j++) {
+        for(i=0; i<dimX; i++) {
+            index = j*dimX+i;
             /* boundary conditions */
             if (i == dimX-1) val1 = 0.0f; else val1 = D[index] - D[j*dimX + (i+1)];
             if (j == dimY-1) val2 = 0.0f; else val2 = D[index] - D[(j+1)*dimX + i];
-
+            
             in_prod = val1*B_x[index] + val2*B_y[index];   /* calculate inner product */
             val1 = val1 - in_prod*B_x[index];
             val2 = val2 - in_prod*B_y[index];
-
+            
             P1[index] = R1[index] + multip*val1;
             P2[index] = R2[index] + multip*val2;
-
+            
         }}
     return 1;
 }
@@ -268,25 +268,25 @@ float Proj_dfunc2D(float *P1, float *P2, int methTV, long DimTotal)
         /* isotropic TV*/
 #pragma omp parallel for shared(P1,P2) private(i,denom,sq_denom)
         for(i=0; i<DimTotal; i++) {
-                denom = powf(P1[i],2) +  powf(P2[i],2);
-                if (denom > 1.0f) {
-					sq_denom = 1.0f/sqrtf(denom);
-                    P1[i] = P1[i]*sq_denom;
-                    P2[i] = P2[i]*sq_denom;
-                }
+            denom = powf(P1[i],2) +  powf(P2[i],2);
+            if (denom > 1.0f) {
+                sq_denom = 1.0f/sqrtf(denom);
+                P1[i] = P1[i]*sq_denom;
+                P2[i] = P2[i]*sq_denom;
             }
+        }
     }
     else {
         /* anisotropic TV*/
 #pragma omp parallel for shared(P1,P2) private(i,val1,val2)
         for(i=0; i<DimTotal; i++) {
-                val1 = fabs(P1[i]);
-                val2 = fabs(P2[i]);
-                if (val1 < 1.0f) {val1 = 1.0f;}
-                if (val2 < 1.0f) {val2 = 1.0f;}
-                P1[i] = P1[i]/val1;
-                P2[i] = P2[i]/val2;
-            }
+            val1 = fabs(P1[i]);
+            val2 = fabs(P2[i]);
+            if (val1 < 1.0f) {val1 = 1.0f;}
+            if (val2 < 1.0f) {val2 = 1.0f;}
+            P1[i] = P1[i]/val1;
+            P2[i] = P2[i]/val2;
+        }
     }
     return 1;
 }
@@ -297,9 +297,9 @@ float Rupd_dfunc2D(float *P1, float *P1_old, float *P2, float *P2_old, float *R1
     multip = ((tk-1.0f)/tkp1);
 #pragma omp parallel for shared(P1,P2,P1_old,P2_old,R1,R2,multip) private(i)
     for(i=0; i<DimTotal; i++) {
-            R1[i] = P1[i] + multip*(P1[i] - P1_old[i]);
-            R2[i] = P2[i] + multip*(P2[i] - P2_old[i]);
-        }
+        R1[i] = P1[i] + multip*(P1[i] - P1_old[i]);
+        R2[i] = P2[i] + multip*(P2[i] - P2_old[i]);
+    }
     return 1;
 }
 
@@ -311,25 +311,26 @@ float GradNorm_func3D(float *B, float *B_x, float *B_y, float *B_z, float eta, l
     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++) {
+    for(k=0; k<dimZ; k++) {
         for(j=0; j<dimY; j++) {
-            for(k=0; k<dimZ; k++) {
-			index = (dimX*dimY)*k + j*dimX+i;
-
-            /* zero boundary conditions */
-            if (i == dimX-1) {val1 = 0.0f;} else {val1 = B[(dimX*dimY)*k + j*dimX+(i+1)];}
-            if (j == dimY-1) {val2 = 0.0f;} else {val2 = B[(dimX*dimY)*k + (j+1)*dimX+i];}
-            if (k == dimZ-1) {val3 = 0.0f;} else {val3 = B[(dimX*dimY)*(k+1) + (j)*dimX+i];}
-
-            gradX = val1 - B[index];
-            gradY = val2 - B[index];
-            gradZ = val3 - B[index];
-            magn = pow(gradX,2) + pow(gradY,2) + pow(gradZ,2);
-            magn = sqrt(magn + pow(eta,2)); /* the eta-smoothed gradients magnitude */
-            B_x[index] = gradX/magn;
-            B_y[index] = gradY/magn;
-            B_z[index] = gradZ/magn;
-        }}}
+            for(i=0; i<dimX; i++) {
+                
+                index = (dimX*dimY)*k + j*dimX+i;
+                
+                /* zero boundary conditions */
+                if (i == dimX-1) {val1 = 0.0f;} else {val1 = B[(dimX*dimY)*k + j*dimX+(i+1)];}
+                if (j == dimY-1) {val2 = 0.0f;} else {val2 = B[(dimX*dimY)*k + (j+1)*dimX+i];}
+                if (k == dimZ-1) {val3 = 0.0f;} else {val3 = B[(dimX*dimY)*(k+1) + (j)*dimX+i];}
+                
+                gradX = val1 - B[index];
+                gradY = val2 - B[index];
+                gradZ = val3 - B[index];
+                magn = pow(gradX,2) + pow(gradY,2) + pow(gradZ,2);
+                magn = sqrt(magn + pow(eta,2)); /* the eta-smoothed gradients magnitude */
+                B_x[index] = gradX/magn;
+                B_y[index] = gradY/magn;
+                B_z[index] = gradZ/magn;
+            }}}
     return 1;
 }
 
@@ -338,15 +339,15 @@ float ProjectVect_func3D(float *R1, float *R2, float *R3, float *B_x, float *B_y
     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++) {
+    for(k=0; k<dimZ; k++) {
         for(j=0; j<dimY; j++) {
-            for(k=0; k<dimZ; k++) {
-			index = (dimX*dimY)*k + j*dimX+i;
-            in_prod = R1[index]*B_x[index] + R2[index]*B_y[index] + R3[index]*B_z[index];   /* calculate inner product */
-            R1[index] = R1[index] - in_prod*B_x[index];
-            R2[index] = R2[index] - in_prod*B_y[index];
-            R3[index] = R3[index] - in_prod*B_z[index];
-        }}}
+            for(i=0; i<dimX; i++) {
+                index = (dimX*dimY)*k + j*dimX+i;
+                in_prod = R1[index]*B_x[index] + R2[index]*B_y[index] + R3[index]*B_z[index];   /* calculate inner product */
+                R1[index] = R1[index] - in_prod*B_x[index];
+                R2[index] = R2[index] - in_prod*B_y[index];
+                R3[index] = R3[index] - in_prod*B_z[index];
+            }}}
     return 1;
 }
 
@@ -355,10 +356,10 @@ float Obj_dfunc3D(float *A, float *D, float *R1, float *R2, float *R3, float lam
     float val1, val2, val3;
     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(k=0; k<dimZ; k++) {
         for(j=0; j<dimY; j++) {
-            for(k=0; k<dimZ; k++) {
-				index = (dimX*dimY)*k + j*dimX+i;
+            for(i=0; i<dimX; i++) {
+                index = (dimX*dimY)*k + j*dimX+i;
                 /* boundary conditions */
                 if (i == 0) {val1 = 0.0f;} else {val1 = R1[(dimX*dimY)*k + j*dimX + (i-1)];}
                 if (j == 0) {val2 = 0.0f;} else {val2 = R2[(dimX*dimY)*k + (j-1)*dimX + i];}
@@ -373,20 +374,20 @@ float Grad_dfunc3D(float *P1, float *P2, float *P3, float *D, float *R1, float *
     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++) {
+    for(k=0; k<dimZ; k++) {
         for(j=0; j<dimY; j++) {
-            for(k=0; k<dimZ; k++) {
-				index = (dimX*dimY)*k + j*dimX+i;
+            for(i=0; i<dimX; i++) {
+                index = (dimX*dimY)*k + j*dimX+i;
                 /* boundary conditions */
                 if (i == dimX-1) val1 = 0.0f; else val1 = D[index] - D[(dimX*dimY)*k + j*dimX + (i+1)];
                 if (j == dimY-1) val2 = 0.0f; else val2 = D[index] - D[(dimX*dimY)*k + (j+1)*dimX + i];
                 if (k == dimZ-1) val3 = 0.0f; else val3 = D[index] - D[(dimX*dimY)*(k+1) + j*dimX + i];
-
+                
                 in_prod = val1*B_x[index] + val2*B_y[index] + val3*B_z[index];   /* calculate inner product */
                 val1 = val1 - in_prod*B_x[index];
                 val2 = val2 - in_prod*B_y[index];
                 val3 = val3 - in_prod*B_z[index];
-
+                
                 P1[index] = R1[index] + multip*val1;
                 P2[index] = R2[index] + multip*val2;
                 P3[index] = R3[index] + multip*val3;
@@ -398,33 +399,33 @@ float Proj_dfunc3D(float *P1, float *P2, float *P3, int methTV, long DimTotal)
     float val1, val2, val3, denom, sq_denom;
     long i;
     if (methTV == 0) {
-	/* isotropic TV*/
-	#pragma omp parallel for shared(P1,P2,P3) private(i,val1,val2,val3,sq_denom)
-    for(i=0; i<DimTotal; i++) {
-				denom = powf(P1[i],2) + powf(P2[i],2) + powf(P3[i],2);
-                if (denom > 1.0f) {
-					sq_denom = 1.0f/sqrtf(denom);
-                    P1[i] = P1[i]*sq_denom;
-                    P2[i] = P2[i]*sq_denom;
-                    P3[i] = P3[i]*sq_denom;
-                }
-			}
-	}
+        /* isotropic TV*/
+#pragma omp parallel for shared(P1,P2,P3) private(i,val1,val2,val3,sq_denom)
+        for(i=0; i<DimTotal; i++) {
+            denom = powf(P1[i],2) + powf(P2[i],2) + powf(P3[i],2);
+            if (denom > 1.0f) {
+                sq_denom = 1.0f/sqrtf(denom);
+                P1[i] = P1[i]*sq_denom;
+                P2[i] = P2[i]*sq_denom;
+                P3[i] = P3[i]*sq_denom;
+            }
+        }
+    }
     else {
-    /* anisotropic TV*/
+        /* anisotropic TV*/
 #pragma omp parallel for shared(P1,P2,P3) private(i,val1,val2,val3)
-    for(i=0; i<DimTotal; i++) {
-                val1 = fabs(P1[i]);
-                val2 = fabs(P2[i]);
-                val3 = fabs(P3[i]);
-                if (val1 < 1.0f) {val1 = 1.0f;}
-                if (val2 < 1.0f) {val2 = 1.0f;}
-                if (val3 < 1.0f) {val3 = 1.0f;}
-                P1[i] = P1[i]/val1;
-                P2[i] = P2[i]/val2;
-                P3[i] = P3[i]/val3;
-            }
-		}
+        for(i=0; i<DimTotal; i++) {
+            val1 = fabs(P1[i]);
+            val2 = fabs(P2[i]);
+            val3 = fabs(P3[i]);
+            if (val1 < 1.0f) {val1 = 1.0f;}
+            if (val2 < 1.0f) {val2 = 1.0f;}
+            if (val3 < 1.0f) {val3 = 1.0f;}
+            P1[i] = P1[i]/val1;
+            P2[i] = P2[i]/val2;
+            P3[i] = P3[i]/val3;
+        }
+    }
     return 1;
 }
 float Rupd_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)
@@ -434,9 +435,9 @@ float Rupd_dfunc3D(float *P1, float *P1_old, float *P2, float *P2_old, float *P3
     multip = ((tk-1.0f)/tkp1);
 #pragma omp parallel for shared(P1,P2,P3,P1_old,P2_old,P3_old,R1,R2,R3,multip) private(i)
     for(i=0; i<DimTotal; i++) {
-                R1[i] = P1[i] + multip*(P1[i] - P1_old[i]);
-                R2[i] = P2[i] + multip*(P2[i] - P2_old[i]);
-                R3[i] = P3[i] + multip*(P3[i] - P3_old[i]);
-            }
+        R1[i] = P1[i] + multip*(P1[i] - P1_old[i]);
+        R2[i] = P2[i] + multip*(P2[i] - P2_old[i]);
+        R3[i] = P3[i] + multip*(P3[i] - P3_old[i]);
+    }
     return 1;
 }
diff --git a/src/Core/regularisers_CPU/LLT_ROF_core.c b/src/Core/regularisers_CPU/LLT_ROF_core.c
index 1064340..e097c90 100644
--- a/src/Core/regularisers_CPU/LLT_ROF_core.c
+++ b/src/Core/regularisers_CPU/LLT_ROF_core.c
@@ -1,21 +1,21 @@
 /*
-This work is part of the Core Imaging Library developed by
-Visual Analytics and Imaging System Group of the Science Technology
-Facilities Council, STFC
-
-Copyright 2017 Daniil Kazantsev
-Copyright 2017 Srikanth Nagella, Edoardo Pasca
-
-Licensed under the Apache License, Version 2.0 (the "License");
-you may not use this file except in compliance with the License.
-You may obtain a copy of the License at
-http://www.apache.org/licenses/LICENSE-2.0
-Unless required by applicable law or agreed to in writing, software
-distributed under the License is distributed on an "AS IS" BASIS,
-WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-See the License for the specific language governing permissions and
-limitations under the License.
-*/
+ * This work is part of the Core Imaging Library developed by
+ * Visual Analytics and Imaging System Group of the Science Technology
+ * Facilities Council, STFC
+ *
+ * Copyright 2017 Daniil Kazantsev
+ * Copyright 2017 Srikanth Nagella, Edoardo Pasca
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ * http://www.apache.org/licenses/LICENSE-2.0
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
 
 #include "LLT_ROF_core.h"
 #define EPS_LLT 1.0e-12
@@ -30,56 +30,56 @@ int signLLT(float x) {
 
 /* C-OMP implementation of Lysaker, Lundervold and Tai (LLT) model [1] combined with Rudin-Osher-Fatemi [2] TV regularisation penalty.
  *
-* This penalty can deliver visually pleasant piecewise-smooth recovery if regularisation parameters are selected well.
-* The rule of thumb for selection is to start with lambdaLLT = 0 (just the ROF-TV model) and then proceed to increase
-* lambdaLLT starting with smaller values.
-*
-* Input Parameters:
-* 1. U0 - original noise image/volume
-* 2. lambdaROF - ROF-related regularisation parameter
-* 3. lambdaLLT - LLT-related regularisation parameter
-* 4. tau - time-marching step
-* 5. iter - iterations number (for both models)
-* 6. eplsilon: tolerance constant
-*
-* Output:
-* [1] Filtered/regularized image/volume
-* [2] Information vector which contains [iteration no., reached tolerance]
-*
-* References:
-* [1] Lysaker, M., Lundervold, A. and Tai, X.C., 2003. Noise removal using fourth-order partial differential equation with applications to medical magnetic resonance images in space and time. IEEE Transactions on image processing, 12(12), pp.1579-1590.
-* [2] Rudin, Osher, Fatemi, "Nonlinear Total Variation based noise removal algorithms"
-*/
+ * This penalty can deliver visually pleasant piecewise-smooth recovery if regularisation parameters are selected well.
+ * The rule of thumb for selection is to start with lambdaLLT = 0 (just the ROF-TV model) and then proceed to increase
+ * lambdaLLT starting with smaller values.
+ *
+ * Input Parameters:
+ * 1. U0 - original noise image/volume
+ * 2. lambdaROF - ROF-related regularisation parameter
+ * 3. lambdaLLT - LLT-related regularisation parameter
+ * 4. tau - time-marching step
+ * 5. iter - iterations number (for both models)
+ * 6. eplsilon: tolerance constant
+ *
+ * Output:
+ * [1] Filtered/regularized image/volume
+ * [2] Information vector which contains [iteration no., reached tolerance]
+ *
+ * References:
+ * [1] Lysaker, M., Lundervold, A. and Tai, X.C., 2003. Noise removal using fourth-order partial differential equation with applications to medical magnetic resonance images in space and time. IEEE Transactions on image processing, 12(12), pp.1579-1590.
+ * [2] Rudin, Osher, Fatemi, "Nonlinear Total Variation based noise removal algorithms"
+ */
 
 float LLT_ROF_CPU_main(float *Input, float *Output, float *infovector, float lambdaROF, float lambdaLLT, int iterationsNumb, float tau, float epsil, int dimX, int dimY, int dimZ)
 {
-		long DimTotal;
+    long DimTotal;
     int ll, j;
     float re, re1;
     re = 0.0f; re1 = 0.0f;
     int count = 0;
-
-		float *D1_LLT=NULL, *D2_LLT=NULL, *D3_LLT=NULL, *D1_ROF=NULL, *D2_ROF=NULL, *D3_ROF=NULL, *Output_prev=NULL;
-				DimTotal = (long)(dimX*dimY*dimZ);
-
-        D1_ROF = calloc(DimTotal, sizeof(float));
-        D2_ROF = calloc(DimTotal, sizeof(float));
-        D3_ROF = calloc(DimTotal, sizeof(float));
-
-        D1_LLT = calloc(DimTotal, sizeof(float));
-        D2_LLT = calloc(DimTotal, sizeof(float));
-        D3_LLT = calloc(DimTotal, sizeof(float));
-
-        copyIm(Input, Output, (long)(dimX), (long)(dimY), (long)(dimZ)); /* initialize  */
-        if (epsil != 0.0f) Output_prev = calloc(DimTotal, sizeof(float));
-
-		for(ll = 0; ll < iterationsNumb; ll++) {
-            if ((epsil != 0.0f) && (ll % 5 == 0)) copyIm(Output, Output_prev, (long)(dimX), (long)(dimY), (long)(dimZ));
-
-            if (dimZ == 1) {
-			      /* 2D case */
-			      /****************ROF******************/
-			      /* calculate first-order differences */
+    
+    float *D1_LLT=NULL, *D2_LLT=NULL, *D3_LLT=NULL, *D1_ROF=NULL, *D2_ROF=NULL, *D3_ROF=NULL, *Output_prev=NULL;
+    DimTotal = (long)(dimX*dimY*dimZ);
+    
+    D1_ROF = calloc(DimTotal, sizeof(float));
+    D2_ROF = calloc(DimTotal, sizeof(float));
+    D3_ROF = calloc(DimTotal, sizeof(float));
+    
+    D1_LLT = calloc(DimTotal, sizeof(float));
+    D2_LLT = calloc(DimTotal, sizeof(float));
+    D3_LLT = calloc(DimTotal, sizeof(float));
+    
+    copyIm(Input, Output, (long)(dimX), (long)(dimY), (long)(dimZ)); /* initialize  */
+    if (epsil != 0.0f) Output_prev = calloc(DimTotal, sizeof(float));
+    
+    for(ll = 0; ll < iterationsNumb; ll++) {
+        if ((epsil != 0.0f) && (ll % 5 == 0)) copyIm(Output, Output_prev, (long)(dimX), (long)(dimY), (long)(dimZ));
+        
+        if (dimZ == 1) {
+            /* 2D case */
+            /****************ROF******************/
+            /* calculate first-order differences */
             D1_func_ROF(Output, D1_ROF, (long)(dimX), (long)(dimY), 1l);
             D2_func_ROF(Output, D2_ROF, (long)(dimX), (long)(dimY), 1l);
             /****************LLT******************/
@@ -87,10 +87,10 @@ float LLT_ROF_CPU_main(float *Input, float *Output, float *infovector, float lam
             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, (long)(dimX), (long)(dimY), 1l);
-            }
-            else {
-			      /* 3D case */
-			      /* calculate first-order differences */
+        }
+        else {
+            /* 3D case */
+            /* calculate first-order differences */
             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));
@@ -99,30 +99,30 @@ float LLT_ROF_CPU_main(float *Input, float *Output, float *infovector, float lam
             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, (long)(dimX), (long)(dimY), (long)(dimZ));
-			       }
-
-             /* check early stopping criteria */
-             if ((epsil != 0.0f) && (ll % 5 == 0)) {
-             re = 0.0f; re1 = 0.0f;
-              for(j=0; j<DimTotal; j++)
-              {
-                  re += powf(Output[j] - Output_prev[j],2);
-                  re1 += powf(Output[j],2);
-              }
-               re = sqrtf(re)/sqrtf(re1);
-               if (re < epsil)  count++;
-               if (count > 3) break;
-             }
-
+        }
+        
+        /* check early stopping criteria */
+        if ((epsil != 0.0f) && (ll % 5 == 0)) {
+            re = 0.0f; re1 = 0.0f;
+            for(j=0; j<DimTotal; j++)
+            {
+                re += powf(Output[j] - Output_prev[j],2);
+                re1 += powf(Output[j],2);
+            }
+            re = sqrtf(re)/sqrtf(re1);
+            if (re < epsil)  count++;
+            if (count > 3) break;
+        }
+        
     } /*end of iterations*/
     free(D1_LLT);free(D2_LLT);free(D3_LLT);
     free(D1_ROF);free(D2_ROF);free(D3_ROF);
     if (epsil != 0.0f) free(Output_prev);
-
+    
     /*adding info into info_vector */
     infovector[0] = (float)(ll);  /*iterations number (if stopped earlier based on tolerance)*/
     infovector[1] = re;  /* reached tolerance */
-	  return 0;
+    return 0;
 }
 
 /*************************************************************************/
@@ -130,65 +130,65 @@ float LLT_ROF_CPU_main(float *Input, float *Output, float *infovector, float lam
 /*************************************************************************/
 float der2D_LLT(float *U, float *D1, float *D2, long dimX, long dimY, long dimZ)
 {
-	long i, j, index, i_p, i_m, j_m, j_p;
-	float dxx, dyy, denom_xx, denom_yy;
+    long i, j, index, i_p, i_m, j_m, j_p;
+    float dxx, dyy, denom_xx, denom_yy;
 #pragma omp parallel for shared(U,D1,D2) private(i, j, index, i_p, i_m, j_m, j_p, denom_xx, denom_yy, dxx, dyy)
-	for (i = 0; i<dimX; i++) {
-		for (j = 0; j<dimY; j++) {
-			index = j*dimX+i;
-			/* symmetric boundary conditions (Neuman) */
-			i_p = i + 1; if (i_p == dimX) i_p = i - 1;
-			i_m = i - 1; if (i_m < 0) i_m = i + 1;
-			j_p = j + 1; if (j_p == dimY) j_p = j - 1;
-			j_m = j - 1; if (j_m < 0) j_m = j + 1;
-
-			dxx = U[j*dimX+i_p] - 2.0f*U[index] + U[j*dimX+i_m];
-			dyy = U[j_p*dimX+i] - 2.0f*U[index] + U[j_m*dimX+i];
-
-			denom_xx = fabs(dxx) + EPS_LLT;
-			denom_yy = fabs(dyy) + EPS_LLT;
-
-			D1[index] = dxx / denom_xx;
-			D2[index] = dyy / denom_yy;
-		}
-	}
-	return 1;
+    for (j = 0; j<dimY; j++) {
+        for (i = 0; i<dimX; i++) {
+            index = j*dimX+i;
+            /* symmetric boundary conditions (Neuman) */
+            i_p = i + 1; if (i_p == dimX) i_p = i - 1;
+            i_m = i - 1; if (i_m < 0) i_m = i + 1;
+            j_p = j + 1; if (j_p == dimY) j_p = j - 1;
+            j_m = j - 1; if (j_m < 0) j_m = j + 1;
+            
+            dxx = U[j*dimX+i_p] - 2.0f*U[index] + U[j*dimX+i_m];
+            dyy = U[j_p*dimX+i] - 2.0f*U[index] + U[j_m*dimX+i];
+            
+            denom_xx = fabs(dxx) + EPS_LLT;
+            denom_yy = fabs(dyy) + EPS_LLT;
+            
+            D1[index] = dxx / denom_xx;
+            D2[index] = dyy / denom_yy;
+        }
+    }
+    return 1;
 }
 
 float der3D_LLT(float *U, float *D1, float *D2, float *D3, long dimX, long dimY, long dimZ)
- {
- 	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++) {
- 		for (j = 0; j<dimY; j++) {
- 			for (k = 0; k<dimZ; k++) {
-				/* symmetric boundary conditions (Neuman) */
-				i_p = i + 1; if (i_p == dimX) i_p = i - 1;
-				i_m = i - 1; if (i_m < 0) i_m = i + 1;
-				j_p = j + 1; if (j_p == dimY) j_p = j - 1;
-				j_m = j - 1; if (j_m < 0) j_m = j + 1;
- 				k_p = k + 1; if (k_p == dimZ) k_p = k - 1;
- 				k_m = k - 1; if (k_m < 0) k_m = k + 1;
-
-				index = (dimX*dimY)*k + j*dimX+i;
-
- 				dxx = U[(dimX*dimY)*k + j*dimX+i_p] - 2.0f*U[index] + U[(dimX*dimY)*k + j*dimX+i_m];
- 				dyy = U[(dimX*dimY)*k + j_p*dimX+i] - 2.0f*U[index] + U[(dimX*dimY)*k + j_m*dimX+i];
- 				dzz = U[(dimX*dimY)*k_p + j*dimX+i] - 2.0f*U[index] + U[(dimX*dimY)*k_m + j*dimX+i];
-
- 				denom_xx = fabs(dxx) + EPS_LLT;
- 				denom_yy = fabs(dyy) + EPS_LLT;
- 				denom_zz = fabs(dzz) + EPS_LLT;
-
- 				D1[index] = dxx / denom_xx;
- 				D2[index] = dyy / denom_yy;
- 				D3[index] = dzz / denom_zz;
- 			}
- 		}
- 	}
- 	return 1;
- }
+{
+    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 (k = 0; k<dimZ; k++) {
+        for (j = 0; j<dimY; j++) {
+            for (i = 0; i<dimX; i++) {
+                /* symmetric boundary conditions (Neuman) */
+                i_p = i + 1; if (i_p == dimX) i_p = i - 1;
+                i_m = i - 1; if (i_m < 0) i_m = i + 1;
+                j_p = j + 1; if (j_p == dimY) j_p = j - 1;
+                j_m = j - 1; if (j_m < 0) j_m = j + 1;
+                k_p = k + 1; if (k_p == dimZ) k_p = k - 1;
+                k_m = k - 1; if (k_m < 0) k_m = k + 1;
+                
+                index = (dimX*dimY)*k + j*dimX+i;
+                
+                dxx = U[(dimX*dimY)*k + j*dimX+i_p] - 2.0f*U[index] + U[(dimX*dimY)*k + j*dimX+i_m];
+                dyy = U[(dimX*dimY)*k + j_p*dimX+i] - 2.0f*U[index] + U[(dimX*dimY)*k + j_m*dimX+i];
+                dzz = U[(dimX*dimY)*k_p + j*dimX+i] - 2.0f*U[index] + U[(dimX*dimY)*k_m + j*dimX+i];
+                
+                denom_xx = fabs(dxx) + EPS_LLT;
+                denom_yy = fabs(dyy) + EPS_LLT;
+                denom_zz = fabs(dzz) + EPS_LLT;
+                
+                D1[index] = dxx / denom_xx;
+                D2[index] = dyy / denom_yy;
+                D3[index] = dzz / denom_zz;
+            }
+        }
+    }
+    return 1;
+}
 
 /*************************************************************************/
 /**********************ROF-related functions *****************************/
@@ -199,13 +199,13 @@ 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;
     long i,j,k,i1,i2,k1,j1,j2,k2,index;
-
+    
     if (dimZ > 1) {
 #pragma omp parallel for shared (A, D1, dimX, dimY, dimZ) private(index, i, j, k, i1, j1, k1, i2, j2, k2, NOMx_1,NOMy_1,NOMy_0,NOMz_1,NOMz_0,denom1,denom2,denom3,T1)
-        for(j=0; j<dimY; j++) {
-            for(i=0; i<dimX; i++) {
-                for(k=0; k<dimZ; k++) {
-					index = (dimX*dimY)*k + j*dimX+i;
+        for (k = 0; k<dimZ; k++) {
+            for (j = 0; j<dimY; j++) {
+                for (i = 0; i<dimX; i++) {
+                    index = (dimX*dimY)*k + j*dimX+i;
                     /* symmetric boundary conditions (Neuman) */
                     i1 = i + 1; if (i1 >= dimX) i1 = i-1;
                     i2 = i - 1; if (i2 < 0) i2 = i+1;
@@ -213,17 +213,17 @@ float D1_func_ROF(float *A, float *D1, long dimX, long dimY, long dimZ)
                     j2 = j - 1; if (j2 < 0) j2 = j+1;
                     k1 = k + 1; if (k1 >= dimZ) k1 = k-1;
                     k2 = k - 1; if (k2 < 0) k2 = k+1;
-
+                    
                     /* Forward-backward differences */
                     NOMx_1 = A[(dimX*dimY)*k + j1*dimX + i] - A[index]; /* x+ */
                     NOMy_1 = A[(dimX*dimY)*k + j*dimX + i1] - A[index]; /* y+ */
                     /*NOMx_0 = (A[(i)*dimY + j] - A[(i2)*dimY + j]); */  /* x- */
                     NOMy_0 = A[index] - A[(dimX*dimY)*k + j*dimX + i2]; /* y- */
-
+                    
                     NOMz_1 = A[(dimX*dimY)*k1 + j*dimX + i] - A[index]; /* z+ */
                     NOMz_0 = A[index] - A[(dimX*dimY)*k2 + j*dimX + i]; /* z- */
-
-
+                    
+                    
                     denom1 = NOMx_1*NOMx_1;
                     denom2 = 0.5f*(signLLT(NOMy_1) + signLLT(NOMy_0))*(MIN(fabs(NOMy_1),fabs(NOMy_0)));
                     denom2 = denom2*denom2;
@@ -237,19 +237,19 @@ float D1_func_ROF(float *A, float *D1, long dimX, long dimY, long dimZ)
 #pragma omp parallel for shared (A, D1, dimX, dimY) private(i, j, i1, j1, i2, j2,NOMx_1,NOMy_1,NOMy_0,denom1,denom2,T1,index)
         for(j=0; j<dimY; j++) {
             for(i=0; i<dimX; i++) {
-				index = j*dimX+i;
+                index = j*dimX+i;
                 /* symmetric boundary conditions (Neuman) */
                 i1 = i + 1; if (i1 >= dimX) i1 = i-1;
                 i2 = i - 1; if (i2 < 0) i2 = i+1;
                 j1 = j + 1; if (j1 >= dimY) j1 = j-1;
                 j2 = j - 1; if (j2 < 0) j2 = j+1;
-
+                
                 /* Forward-backward differences */
                 NOMx_1 = A[j1*dimX + i] - A[index]; /* x+ */
                 NOMy_1 = A[j*dimX + i1] - A[index]; /* y+ */
                 /*NOMx_0 = (A[(i)*dimY + j] - A[(i2)*dimY + j]); */ /* x- */
                 NOMy_0 = A[index] - A[(j)*dimX + i2]; /* y- */
-
+                
                 denom1 = NOMx_1*NOMx_1;
                 denom2 = 0.5f*(signLLT(NOMy_1) + signLLT(NOMy_0))*(MIN(fabs(NOMy_1),fabs(NOMy_0)));
                 denom2 = denom2*denom2;
@@ -264,13 +264,13 @@ float D2_func_ROF(float *A, float *D2, long dimX, long dimY, long dimZ)
 {
     float NOMx_1, NOMy_1, NOMx_0, NOMz_1, NOMz_0, denom1, denom2, denom3, T2;
     long i,j,k,i1,i2,k1,j1,j2,k2,index;
-
+    
     if (dimZ > 1) {
 #pragma omp parallel for shared (A, D2, dimX, dimY, dimZ) private(index, i, j, k, i1, j1, k1, i2, j2, k2,  NOMx_1, NOMy_1, NOMx_0, NOMz_1, NOMz_0, denom1, denom2, denom3, T2)
-        for(j=0; j<dimY; j++) {
-            for(i=0; i<dimX; i++) {
-                for(k=0; k<dimZ; k++) {
-					index = (dimX*dimY)*k + j*dimX+i;
+        for (k = 0; k<dimZ; k++) {
+            for (j = 0; j<dimY; j++) {
+                for (i = 0; i<dimX; i++) {
+                    index = (dimX*dimY)*k + j*dimX+i;
                     /* symmetric boundary conditions (Neuman) */
                     i1 = i + 1; if (i1 >= dimX) i1 = i-1;
                     i2 = i - 1; if (i2 < 0) i2 = i+1;
@@ -278,16 +278,16 @@ float D2_func_ROF(float *A, float *D2, long dimX, long dimY, long dimZ)
                     j2 = j - 1; if (j2 < 0) j2 = j+1;
                     k1 = k + 1; if (k1 >= dimZ) k1 = k-1;
                     k2 = k - 1; if (k2 < 0) k2 = k+1;
-
-
+                    
+                    
                     /* Forward-backward differences */
                     NOMx_1 = A[(dimX*dimY)*k + (j1)*dimX + i] - A[index]; /* x+ */
                     NOMy_1 = A[(dimX*dimY)*k + (j)*dimX + i1] - A[index]; /* y+ */
                     NOMx_0 = A[index] - A[(dimX*dimY)*k + (j2)*dimX + i]; /* x- */
                     NOMz_1 = A[(dimX*dimY)*k1 + j*dimX + i] - A[index]; /* z+ */
                     NOMz_0 = A[index] - A[(dimX*dimY)*k2 + (j)*dimX + i]; /* z- */
-
-
+                    
+                    
                     denom1 = NOMy_1*NOMy_1;
                     denom2 = 0.5f*(signLLT(NOMx_1) + signLLT(NOMx_0))*(MIN(fabs(NOMx_1),fabs(NOMx_0)));
                     denom2 = denom2*denom2;
@@ -301,19 +301,19 @@ float D2_func_ROF(float *A, float *D2, long dimX, long dimY, long dimZ)
 #pragma omp parallel for shared (A, D2, dimX, dimY) private(i, j, i1, j1, i2, j2, NOMx_1,NOMy_1,NOMx_0,denom1,denom2,T2,index)
         for(j=0; j<dimY; j++) {
             for(i=0; i<dimX; i++) {
-				index = j*dimX+i;
+                index = j*dimX+i;
                 /* symmetric boundary conditions (Neuman) */
                 i1 = i + 1; if (i1 >= dimX) i1 = i-1;
                 i2 = i - 1; if (i2 < 0) i2 = i+1;
                 j1 = j + 1; if (j1 >= dimY) j1 = j-1;
                 j2 = j - 1; if (j2 < 0) j2 = j+1;
-
+                
                 /* Forward-backward differences */
                 NOMx_1 = A[j1*dimX + i] - A[index]; /* x+ */
                 NOMy_1 = A[j*dimX + i1] - A[index]; /* y+ */
                 NOMx_0 = A[index] - A[j2*dimX + i]; /* x- */
                 /*NOMy_0 = A[(i)*dimY + j] - A[(i)*dimY + j2]; */  /* y- */
-
+                
                 denom1 = NOMy_1*NOMy_1;
                 denom2 = 0.5f*(signLLT(NOMx_1) + signLLT(NOMx_0))*(MIN(fabs(NOMx_1),fabs(NOMx_0)));
                 denom2 = denom2*denom2;
@@ -329,12 +329,12 @@ float D3_func_ROF(float *A, float *D3, long dimX, long dimY, long dimZ)
 {
     float NOMx_1, NOMy_1, NOMx_0, NOMy_0, NOMz_1, denom1, denom2, denom3, T3;
     long index,i,j,k,i1,i2,k1,j1,j2,k2;
-
+    
 #pragma omp parallel for shared (A, D3, dimX, dimY, dimZ) private(index, i, j, k, i1, j1, k1, i2, j2, k2,  NOMx_1, NOMy_1, NOMy_0, NOMx_0, NOMz_1, denom1, denom2, denom3, T3)
-    for(j=0; j<dimY; j++) {
-        for(i=0; i<dimX; i++) {
-            for(k=0; k<dimZ; k++) {
-				index = (dimX*dimY)*k + j*dimX+i;
+    for (k = 0; k<dimZ; k++) {
+        for (j = 0; j<dimY; j++) {
+            for (i = 0; i<dimX; i++) {
+                index = (dimX*dimY)*k + j*dimX+i;
                 /* symmetric boundary conditions (Neuman) */
                 i1 = i + 1; if (i1 >= dimX) i1 = i-1;
                 i2 = i - 1; if (i2 < 0) i2 = i+1;
@@ -342,7 +342,7 @@ float D3_func_ROF(float *A, float *D3, long dimX, long dimY, long dimZ)
                 j2 = j - 1; if (j2 < 0) j2 = j+1;
                 k1 = k + 1; if (k1 >= dimZ) k1 = k-1;
                 k2 = k - 1; if (k2 < 0) k2 = k+1;
-
+                
                 /* Forward-backward differences */
                 NOMx_1 = A[(dimX*dimY)*k + (j1)*dimX + i] - A[index]; /* x+ */
                 NOMy_1 = A[(dimX*dimY)*k + (j)*dimX + i1] - A[index]; /* y+ */
@@ -350,7 +350,7 @@ float D3_func_ROF(float *A, float *D3, long dimX, long dimY, long dimZ)
                 NOMx_0 = A[index] - A[(dimX*dimY)*k + (j2)*dimX + i]; /* x- */
                 NOMz_1 = A[(dimX*dimY)*k1 + j*dimX + i] - A[index]; /* z+ */
                 /*NOMz_0 = A[(dimX*dimY)*k + (i)*dimY + j] - A[(dimX*dimY)*k2 + (i)*dimY + j]; */ /* z- */
-
+                
                 denom1 = NOMz_1*NOMz_1;
                 denom2 = 0.5f*(signLLT(NOMx_1) + signLLT(NOMx_0))*(MIN(fabs(NOMx_1),fabs(NOMx_0)));
                 denom2 = denom2*denom2;
@@ -368,69 +368,69 @@ float D3_func_ROF(float *A, float *D3, long dimX, long dimY, long dimZ)
 
 float Update2D_LLT_ROF(float *U0, float *U, float *D1_LLT, float *D2_LLT, float *D1_ROF, float *D2_ROF, float lambdaROF, float lambdaLLT, float tau, long dimX, long dimY, long dimZ)
 {
-	long i, j, index, i_p, i_m, j_m, j_p;
-	float div, laplc, dxx, dyy, dv1, dv2;
+    long i, j, index, i_p, i_m, j_m, j_p;
+    float div, laplc, dxx, dyy, dv1, dv2;
 #pragma omp parallel for shared(U,U0) private(i, j, index, i_p, i_m, j_m, j_p, laplc, div, dxx, dyy, dv1, dv2)
-	for (i = 0; i<dimX; i++) {
-		for (j = 0; j<dimY; j++) {
-			index = j*dimX+i;
-			/* symmetric boundary conditions (Neuman) */
-			i_p = i + 1; if (i_p == dimX) i_p = i - 1;
-			i_m = i - 1; if (i_m < 0) i_m = i + 1;
-			j_p = j + 1; if (j_p == dimY) j_p = j - 1;
-			j_m = j - 1; if (j_m < 0) j_m = j + 1;
-
-			/*LLT-related part*/
-			dxx = D1_LLT[j*dimX+i_p] - 2.0f*D1_LLT[index] + D1_LLT[j*dimX+i_m];
-			dyy = D2_LLT[j_p*dimX+i] - 2.0f*D2_LLT[index] + D2_LLT[j_m*dimX+i];
-			laplc = dxx + dyy; /*build Laplacian*/
-
-			/*ROF-related part*/
-			dv1 = D1_ROF[index] - D1_ROF[j_m*dimX + i];
+    for (j = 0; j<dimY; j++) {
+        for (i = 0; i<dimX; i++) {
+            index = j*dimX+i;
+            /* symmetric boundary conditions (Neuman) */
+            i_p = i + 1; if (i_p == dimX) i_p = i - 1;
+            i_m = i - 1; if (i_m < 0) i_m = i + 1;
+            j_p = j + 1; if (j_p == dimY) j_p = j - 1;
+            j_m = j - 1; if (j_m < 0) j_m = j + 1;
+            
+            /*LLT-related part*/
+            dxx = D1_LLT[j*dimX+i_p] - 2.0f*D1_LLT[index] + D1_LLT[j*dimX+i_m];
+            dyy = D2_LLT[j_p*dimX+i] - 2.0f*D2_LLT[index] + D2_LLT[j_m*dimX+i];
+            laplc = dxx + dyy; /*build Laplacian*/
+            
+            /*ROF-related part*/
+            dv1 = D1_ROF[index] - D1_ROF[j_m*dimX + i];
             dv2 = D2_ROF[index] - D2_ROF[j*dimX + i_m];
-			div = dv1 + dv2; /*build Divirgent*/
-
-			/*combine all into one cost function to minimise */
+            div = dv1 + dv2; /*build Divirgent*/
+            
+            /*combine all into one cost function to minimise */
             U[index] += tau*(lambdaROF*(div) - lambdaLLT*(laplc) - (U[index] - U0[index]));
-		}
-	}
-	return *U;
+        }
+    }
+    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, long dimX, long dimY, long dimZ)
 {
-	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;
+    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++) {
- 		for (j = 0; j<dimY; j++) {
- 			for (k = 0; k<dimZ; k++) {
-				/* symmetric boundary conditions (Neuman) */
-				i_p = i + 1; if (i_p == dimX) i_p = i - 1;
-				i_m = i - 1; if (i_m < 0) i_m = i + 1;
-				j_p = j + 1; if (j_p == dimY) j_p = j - 1;
-				j_m = j - 1; if (j_m < 0) j_m = j + 1;
- 				k_p = k + 1; if (k_p == dimZ) k_p = k - 1;
- 				k_m = k - 1; if (k_m < 0) k_m = k + 1;
-
-				index = (dimX*dimY)*k + j*dimX+i;
-
-				/*LLT-related part*/
-				dxx = D1_LLT[(dimX*dimY)*k + j*dimX+i_p] - 2.0f*D1_LLT[index] + D1_LLT[(dimX*dimY)*k + j*dimX+i_m];
-				dyy = D2_LLT[(dimX*dimY)*k + j_p*dimX+i] - 2.0f*D2_LLT[index] + D2_LLT[(dimX*dimY)*k + j_m*dimX+i];
-				dzz = D3_LLT[(dimX*dimY)*k_p + j*dimX+i] - 2.0f*D3_LLT[index] + D3_LLT[(dimX*dimY)*k_m + j*dimX+i];
-				laplc = dxx + dyy + dzz; /*build Laplacian*/
-
-				/*ROF-related part*/
-				dv1 = D1_ROF[index] - D1_ROF[(dimX*dimY)*k + j_m*dimX+i];
-				dv2 = D2_ROF[index] - D2_ROF[(dimX*dimY)*k + j*dimX+i_m];
-				dv3 = D3_ROF[index] - D3_ROF[(dimX*dimY)*k_m + j*dimX+i];
-				div = dv1 + dv2 + dv3; /*build Divirgent*/
-
-				/*combine all into one cost function to minimise */
-				U[index] += tau*(lambdaROF*(div) - lambdaLLT*(laplc) - (U[index] - U0[index]));
-			}
-		}
-	}
-	return *U;
+    for (k = 0; k<dimZ; k++) {
+        for (j = 0; j<dimY; j++) {
+            for (i = 0; i<dimX; i++) {
+                /* symmetric boundary conditions (Neuman) */
+                i_p = i + 1; if (i_p == dimX) i_p = i - 1;
+                i_m = i - 1; if (i_m < 0) i_m = i + 1;
+                j_p = j + 1; if (j_p == dimY) j_p = j - 1;
+                j_m = j - 1; if (j_m < 0) j_m = j + 1;
+                k_p = k + 1; if (k_p == dimZ) k_p = k - 1;
+                k_m = k - 1; if (k_m < 0) k_m = k + 1;
+                
+                index = (dimX*dimY)*k + j*dimX+i;
+                
+                /*LLT-related part*/
+                dxx = D1_LLT[(dimX*dimY)*k + j*dimX+i_p] - 2.0f*D1_LLT[index] + D1_LLT[(dimX*dimY)*k + j*dimX+i_m];
+                dyy = D2_LLT[(dimX*dimY)*k + j_p*dimX+i] - 2.0f*D2_LLT[index] + D2_LLT[(dimX*dimY)*k + j_m*dimX+i];
+                dzz = D3_LLT[(dimX*dimY)*k_p + j*dimX+i] - 2.0f*D3_LLT[index] + D3_LLT[(dimX*dimY)*k_m + j*dimX+i];
+                laplc = dxx + dyy + dzz; /*build Laplacian*/
+                
+                /*ROF-related part*/
+                dv1 = D1_ROF[index] - D1_ROF[(dimX*dimY)*k + j_m*dimX+i];
+                dv2 = D2_ROF[index] - D2_ROF[(dimX*dimY)*k + j*dimX+i_m];
+                dv3 = D3_ROF[index] - D3_ROF[(dimX*dimY)*k_m + j*dimX+i];
+                div = dv1 + dv2 + dv3; /*build Divirgent*/
+                
+                /*combine all into one cost function to minimise */
+                U[index] += tau*(lambdaROF*(div) - lambdaLLT*(laplc) - (U[index] - U0[index]));
+            }
+        }
+    }
+    return *U;
 }
diff --git a/src/Core/regularisers_CPU/Nonlocal_TV_core.c b/src/Core/regularisers_CPU/Nonlocal_TV_core.c
index de1c173..7b1368d 100644
--- a/src/Core/regularisers_CPU/Nonlocal_TV_core.c
+++ b/src/Core/regularisers_CPU/Nonlocal_TV_core.c
@@ -34,52 +34,52 @@
  * 5. Weights_ij(k) - associated weights
  * 6. regularisation parameter
  * 7. iterations number
-
+ *
  * Output:
  * 1. denoised image/volume
  * Elmoataz, Abderrahim, Olivier Lezoray, and Sébastien Bougleux. "Nonlocal discrete regularization on weighted graphs: a framework for image and manifold processing." IEEE Trans. Image Processing 17, no. 7 (2008): 1047-1060.
-
+ *
  */
 /*****************************************************************************/
 
 float Nonlocal_TV_CPU_main(float *A_orig, float *Output, unsigned short *H_i, unsigned short *H_j, unsigned short *H_k, float *Weights, int dimX, int dimY, int dimZ, int NumNeighb, float lambdaReg, int IterNumb, int switchM)
 {
-
+    
     long i, j, k;
     int iter;
     lambdaReg = 1.0f/lambdaReg;
-
+    
     /*****2D INPUT *****/
     if (dimZ == 0) {
-	  copyIm(A_orig, Output, (long)(dimX), (long)(dimY), 1l);
-    /* for each pixel store indeces of the most similar neighbours (patches) */
-     for(iter=0; iter<IterNumb; iter++) {
+        copyIm(A_orig, Output, (long)(dimX), (long)(dimY), 1l);
+        /* for each pixel store indeces of the most similar neighbours (patches) */
+        for(iter=0; iter<IterNumb; iter++) {
 #pragma omp parallel for shared (A_orig, Output, Weights, H_i, H_j, iter) private(i,j)
-   for(j=0; j<(long)(dimY); j++) {
-      for(i=0; i<(long)(dimX); i++) {
-             /*NLM_H1_2D(Output, A_orig, H_i, H_j, Weights, i, j, (long)(dimX), (long)(dimY), NumNeighb, lambdaReg);*/  /* NLM - H1 penalty */
-             if (switchM == 1) {
-               NLM_TV_2D(Output, A_orig, H_j, H_i, Weights, i, j, (long)(dimX), (long)(dimY), NumNeighb, lambdaReg);  /* NLM - TV penalty */
-             }
-             else {
-               NLM_TV_2D(Output, A_orig, H_i, H_j, Weights, i, j, (long)(dimX), (long)(dimY), NumNeighb, lambdaReg);  /* NLM - TV penalty */
-             }
-           }}
-          }
+            for(j=0; j<(long)(dimY); j++) {
+                for(i=0; i<(long)(dimX); i++) {
+                    /*NLM_H1_2D(Output, A_orig, H_i, H_j, Weights, i, j, (long)(dimX), (long)(dimY), NumNeighb, lambdaReg);*/  /* NLM - H1 penalty */
+                    if (switchM == 1) {
+                        NLM_TV_2D(Output, A_orig, H_j, H_i, Weights, i, j, (long)(dimX), (long)(dimY), NumNeighb, lambdaReg);  /* NLM - TV penalty */
+                    }
+                    else {
+                        NLM_TV_2D(Output, A_orig, H_i, H_j, Weights, i, j, (long)(dimX), (long)(dimY), NumNeighb, lambdaReg);  /* NLM - TV penalty */
+                    }
+                }}
+        }
     }
     else {
-     /*****3D INPUT *****/
+        /*****3D INPUT *****/
         copyIm(A_orig, Output, (long)(dimX), (long)(dimY), (long)(dimZ));
-    /* for each pixel store indeces of the most similar neighbours (patches) */
-     for(iter=0; iter<IterNumb; iter++) {
+        /* for each pixel store indeces of the most similar neighbours (patches) */
+        for(iter=0; iter<IterNumb; iter++) {
 #pragma omp parallel for shared (A_orig, Output, Weights, H_i, H_j, H_k, iter) private(i,j,k)
-      for(i=0; i<(long)(dimX); i++) {
-            for(j=0; j<(long)(dimY); j++) {
-               for(k=0; k<(long)(dimZ); k++) {
-            /* NLM_H1_3D(Output, A_orig, H_i, H_j, H_k, Weights, i, j, k, dimX, dimY, dimZ, NumNeighb, lambdaReg); */ /* NLM - H1 penalty */
-            NLM_TV_3D(Output, A_orig, H_i, H_j, H_k, Weights, i, j, k, (long)(dimX), (long)(dimY), (long)(dimZ), NumNeighb, lambdaReg);   /* NLM - TV penalty */
-           }}}
-          }
+            for(k=0; k<(long)(dimZ); k++) {
+                for(j=0; j<(long)(dimY); j++) {
+                    for(i=0; i<(long)(dimX); i++) {
+                        /* NLM_H1_3D(Output, A_orig, H_i, H_j, H_k, Weights, i, j, k, dimX, dimY, dimZ, NumNeighb, lambdaReg); */ /* NLM - H1 penalty */
+                        NLM_TV_3D(Output, A_orig, H_i, H_j, H_k, Weights, i, j, k, (long)(dimX), (long)(dimY), (long)(dimZ), NumNeighb, lambdaReg);   /* NLM - TV penalty */
+                    }}}
+        }
     }
     return *Output;
 }
@@ -87,61 +87,60 @@ float Nonlocal_TV_CPU_main(float *A_orig, float *Output, unsigned short *H_i, un
 /***********<<<<Main Function for NLM - H1 penalty>>>>**********/
 float NLM_H1_2D(float *A, float *A_orig, unsigned short *H_i, unsigned short *H_j, float *Weights, long i, long j, long dimX, long dimY, int NumNeighb, float lambdaReg)
 {
-	long x, i1, j1, index, index_m;
-	float value = 0.0f, normweight  = 0.0f;
-
-	index_m = j*dimX+i;
-	for(x=0; x < NumNeighb; x++) {
-	index =  (dimX*dimY*x) + j*dimX+i;
-		i1 = H_i[index];
-		j1 = H_j[index];
-		value += A[j1*dimX+i1]*Weights[index];
-		normweight += Weights[index];
-	}
-	 A[index_m] = (lambdaReg*A_orig[index_m] + value)/(lambdaReg + normweight);
+    long x, i1, j1, index, index_m;
+    float value = 0.0f, normweight  = 0.0f;
+    
+    index_m = j*dimX+i;
+    for(x=0; x < NumNeighb; x++) {
+        index =  (dimX*dimY*x) + j*dimX+i;
+        i1 = H_i[index];
+        j1 = H_j[index];
+        value += A[j1*dimX+i1]*Weights[index];
+        normweight += Weights[index];
+    }
+    A[index_m] = (lambdaReg*A_orig[index_m] + value)/(lambdaReg + normweight);
     return *A;
 }
 /*3D version*/
 float NLM_H1_3D(float *A, float *A_orig, unsigned short *H_i, unsigned short *H_j, unsigned short *H_k, float *Weights, long i, long j, long k, long dimX, long dimY, long dimZ, int NumNeighb, float lambdaReg)
 {
-	long x, i1, j1, k1, index;
-	float value = 0.0f, normweight  = 0.0f;
-
-	for(x=0; x < NumNeighb; x++) {
-	index = dimX*dimY*dimZ*x + (dimX*dimY*k) + j*dimX+i;
-		i1 = H_i[index];
-		j1 = H_j[index];
-		k1 = H_k[index];
-		value += A[(dimX*dimY*k1) + j1*dimX+i1]*Weights[index];
-		normweight += Weights[index];
-	}
+    long x, i1, j1, k1, index;
+    float value = 0.0f, normweight  = 0.0f;
+    
+    for(x=0; x < NumNeighb; x++) {
+        index = dimX*dimY*dimZ*x + (dimX*dimY*k) + j*dimX+i;
+        i1 = H_i[index];
+        j1 = H_j[index];
+        k1 = H_k[index];
+        value += A[(dimX*dimY*k1) + j1*dimX+i1]*Weights[index];
+        normweight += Weights[index];
+    }
     A[(dimX*dimY*k) + j*dimX+i] = (lambdaReg*A_orig[(dimX*dimY*k) + j*dimX+i] + value)/(lambdaReg + normweight);
     return *A;
 }
 
-
 /***********<<<<Main Function for NLM - TV penalty>>>>**********/
 float NLM_TV_2D(float *A, float *A_orig, unsigned short *H_i, unsigned short *H_j, float *Weights, long i, long j, long dimX, long dimY, int NumNeighb, float lambdaReg)
 {
-	long x, i1, j1, index, index_m;
-	float value = 0.0f, normweight  = 0.0f, NLgrad_magn = 0.0f, NLCoeff;
-
-	 index_m = j*dimX+i;
-
-	for(x=0; x < NumNeighb; x++) {
-		index =  (dimX*dimY*x) + j*dimX+i; /*c*/
-		i1 = H_i[index];
-		j1 = H_j[index];
-		NLgrad_magn += powf((A[j1*dimX+i1] - A[index_m]),2)*Weights[index];
-	}
-
+    long x, i1, j1, index, index_m;
+    float value = 0.0f, normweight  = 0.0f, NLgrad_magn = 0.0f, NLCoeff;
+    
+    index_m = j*dimX+i;
+    
+    for(x=0; x < NumNeighb; x++) {
+        index =  (dimX*dimY*x) + j*dimX+i; /*c*/
+        i1 = H_i[index];
+        j1 = H_j[index];
+        NLgrad_magn += powf((A[j1*dimX+i1] - A[index_m]),2)*Weights[index];
+    }
+    
     NLgrad_magn = sqrtf(NLgrad_magn); /*Non Local Gradients Magnitude */
     NLCoeff = 2.0f*(1.0f/(NLgrad_magn + EPS));
-
+    
     for(x=0; x < NumNeighb; x++) {
-	index =  (dimX*dimY*x) + j*dimX+i; /*c*/
-	i1 = H_i[index];
-	j1 = H_j[index];
+        index =  (dimX*dimY*x) + j*dimX+i; /*c*/
+        i1 = H_i[index];
+        j1 = H_j[index];
         value += A[j1*dimX+i1]*NLCoeff*Weights[index];
         normweight += Weights[index]*NLCoeff;
     }
@@ -151,25 +150,25 @@ float NLM_TV_2D(float *A, float *A_orig, unsigned short *H_i, unsigned short *H_
 /*3D version*/
 float NLM_TV_3D(float *A, float *A_orig, unsigned short *H_i, unsigned short *H_j, unsigned short *H_k, float *Weights, long i, long j, long k, long dimX, long dimY, long dimZ, int NumNeighb, float lambdaReg)
 {
-	long x, i1, j1, k1, index;
-	float value = 0.0f, normweight  = 0.0f, NLgrad_magn = 0.0f, NLCoeff;
-
-	for(x=0; x < NumNeighb; x++) {
-	index =  dimX*dimY*dimZ*x + (dimX*dimY*k) + j*dimX+i;
-		i1 = H_i[index];
-		j1 = H_j[index];
-		k1 = H_k[index];
-	        NLgrad_magn += powf((A[(dimX*dimY*k1) + j1*dimX+i1] - A[(dimX*dimY*k1) + j*dimX+i]),2)*Weights[index];
-	}
-
+    long x, i1, j1, k1, index;
+    float value = 0.0f, normweight  = 0.0f, NLgrad_magn = 0.0f, NLCoeff;
+    
+    for(x=0; x < NumNeighb; x++) {
+        index =  dimX*dimY*dimZ*x + (dimX*dimY*k) + j*dimX+i;
+        i1 = H_i[index];
+        j1 = H_j[index];
+        k1 = H_k[index];
+        NLgrad_magn += powf((A[(dimX*dimY*k1) + j1*dimX+i1] - A[(dimX*dimY*k1) + j*dimX+i]),2)*Weights[index];
+    }
+    
     NLgrad_magn = sqrtf(NLgrad_magn); /*Non Local Gradients Magnitude */
     NLCoeff = 2.0f*(1.0f/(NLgrad_magn + EPS));
-
+    
     for(x=0; x < NumNeighb; x++) {
-	index = dimX*dimY*dimZ*x + (dimX*dimY*k) + j*dimX+i;
-	i1 = H_i[index];
-	j1 = H_j[index];
-	k1 = H_k[index];
+        index = dimX*dimY*dimZ*x + (dimX*dimY*k) + j*dimX+i;
+        i1 = H_i[index];
+        j1 = H_j[index];
+        k1 = H_k[index];
         value += A[(dimX*dimY*k1) + j1*dimX+i1]*NLCoeff*Weights[index];
         normweight += Weights[index]*NLCoeff;
     }
diff --git a/src/Core/regularisers_CPU/PatchSelect_core.c b/src/Core/regularisers_CPU/PatchSelect_core.c
index 8581797..d80353b 100644
--- a/src/Core/regularisers_CPU/PatchSelect_core.c
+++ b/src/Core/regularisers_CPU/PatchSelect_core.c
@@ -77,13 +77,13 @@ float PatchSelect_CPU_main(float *A, unsigned short *H_i, unsigned short *H_j, u
             }} /*main neighb loop */
         /* for each pixel store indeces of the most similar neighbours (patches) */
 #pragma omp parallel for shared (A, Weights, H_i, H_j) private(i,j)
-    for(i=0; i<(long)(dimX); i++) {
-          for(j=0; j<(long)(dimY); j++) {
+        for(j=0; j<(long)(dimY); j++) {
+            for(i=0; i<(long)(dimX); i++) {
                 Indeces2D(A, H_i, H_j, Weights, i, j, (long)(dimX), (long)(dimY), Eucl_Vec, NumNeighb, SearchWindow, SimilarWin, h2);
             }}
     }
     else {
-    /****************3D INPUT ***************/
+        /****************3D INPUT ***************/
         /* generate a 3D Gaussian kernel for NLM procedure */
         Eucl_Vec = (float*) calloc ((2*SimilarWin+1)*(2*SimilarWin+1)*(2*SimilarWin+1),sizeof(float));
         counterG = 0;
@@ -93,12 +93,12 @@ float PatchSelect_CPU_main(float *A, unsigned short *H_i, unsigned short *H_j, u
                     Eucl_Vec[counterG] = (float)exp(-(pow(((float) i), 2) + pow(((float) j), 2) + pow(((float) k), 2))/(2*SimilarWin*SimilarWin*SimilarWin));
                     counterG++;
                 }}} /*main neighb loop */
-
+        
         /* for each voxel store indeces of the most similar neighbours (patches) */
 #pragma omp parallel for shared (A, Weights, H_i, H_j, H_k) private(i,j,k)
-        for(i=0; i<dimX; i++) {
+        for(k=0; k<dimZ; k++) {
             for(j=0; j<dimY; j++) {
-                for(k=0; k<dimZ; k++) {
+                for(i=0; i<dimX; i++) {
                     Indeces3D(A, H_i, H_j, H_k, Weights, i, j, k, (long)(dimX), (long)(dimY), (long)(dimZ), Eucl_Vec, NumNeighb, SearchWindow, SimilarWin, h2);
                 }}}
     }
@@ -111,13 +111,13 @@ float Indeces2D(float *Aorig, unsigned short *H_i, unsigned short *H_j, float *W
     long i1, j1, i_m, j_m, i_c, j_c, i2, j2, i3, j3, counter, x, y, index, sizeWin_tot, counterG;
     float *Weight_Vec, normsum;
     unsigned short *ind_i, *ind_j;
-
+    
     sizeWin_tot = (2*SearchWindow + 1)*(2*SearchWindow + 1);
-
+    
     Weight_Vec = (float*) calloc(sizeWin_tot, sizeof(float));
     ind_i = (unsigned short*) calloc(sizeWin_tot, sizeof(unsigned short));
     ind_j = (unsigned short*) calloc(sizeWin_tot, sizeof(unsigned short));
-
+    
     counter = 0;
     for(i_m=-SearchWindow; i_m<=SearchWindow; i_m++) {
         for(j_m=-SearchWindow; j_m<=SearchWindow; j_m++) {
@@ -149,15 +149,15 @@ float Indeces2D(float *Aorig, unsigned short *H_i, unsigned short *H_j, float *W
     /* do sorting to choose the most prominent weights [HIGH to LOW] */
     /* and re-arrange indeces accordingly */
     for (x = 0; x < counter-1; x++)  {
-       for (y = 0; y < counter-x-1; y++)  {
-           if (Weight_Vec[y] < Weight_Vec[y+1]) {
-            swap(&Weight_Vec[y], &Weight_Vec[y+1]);
-            swapUS(&ind_i[y], &ind_i[y+1]);
-            swapUS(&ind_j[y], &ind_j[y+1]);
+        for (y = 0; y < counter-x-1; y++)  {
+            if (Weight_Vec[y] < Weight_Vec[y+1]) {
+                swap(&Weight_Vec[y], &Weight_Vec[y+1]);
+                swapUS(&ind_i[y], &ind_i[y+1]);
+                swapUS(&ind_j[y], &ind_j[y+1]);
             }
-    	}
+        }
     }
-     /*sorting loop finished*/
+    /*sorting loop finished*/
     /*now select the NumNeighb more prominent weights and store into pre-allocated arrays */
     for(x=0; x < NumNeighb; x++) {
         index = (dimX*dimY*x) + j*dimX+i;
@@ -176,14 +176,14 @@ float Indeces3D(float *Aorig, unsigned short *H_i, unsigned short *H_j, unsigned
     long i1, j1, k1, i_m, j_m, k_m, i_c, j_c, k_c, i2, j2, k2, i3, j3, k3, counter, x, y, index, sizeWin_tot, counterG;
     float *Weight_Vec, normsum, temp;
     unsigned short *ind_i, *ind_j, *ind_k, temp_i, temp_j, temp_k;
-
+    
     sizeWin_tot = (2*SearchWindow + 1)*(2*SearchWindow + 1)*(2*SearchWindow + 1);
-
+    
     Weight_Vec = (float*) calloc(sizeWin_tot, sizeof(float));
     ind_i = (unsigned short*) calloc(sizeWin_tot, sizeof(unsigned short));
     ind_j = (unsigned short*) calloc(sizeWin_tot, sizeof(unsigned short));
     ind_k = (unsigned short*) calloc(sizeWin_tot, sizeof(unsigned short));
-
+    
     counter = 0l;
     for(i_m=-SearchWindow; i_m<=SearchWindow; i_m++) {
         for(j_m=-SearchWindow; j_m<=SearchWindow; j_m++) {
@@ -237,18 +237,18 @@ float Indeces3D(float *Aorig, unsigned short *H_i, unsigned short *H_j, unsigned
                 ind_k[y] = temp_k;
             }}}
     /*sorting loop finished*/
-
+    
     /*now select the NumNeighb more prominent weights and store into arrays */
     for(x=0; x < NumNeighb; x++) {
         index = dimX*dimY*dimZ*x + (dimX*dimY*k) + j*dimX+i;
-
+        
         H_i[index] = ind_i[x];
         H_j[index] = ind_j[x];
         H_k[index] = ind_k[x];
-
+        
         Weights[index] = Weight_Vec[x];
     }
-
+    
     free(ind_i);
     free(ind_j);
     free(ind_k);
diff --git a/src/Core/regularisers_CPU/ROF_TV_core.c b/src/Core/regularisers_CPU/ROF_TV_core.c
index 004a0f2..cf08652 100644
--- a/src/Core/regularisers_CPU/ROF_TV_core.c
+++ b/src/Core/regularisers_CPU/ROF_TV_core.c
@@ -56,46 +56,46 @@ float TV_ROF_CPU_main(float *Input, float *Output, float *infovector, float *lam
     int i;
     long DimTotal,j;
     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, (long)(dimX), (long)(dimY), (long)(dimZ));
     if (epsil != 0.0f) Output_prev = calloc(DimTotal, sizeof(float));
-
+    
     /* start TV iterations */
     for(i=0; i < iterationsNumb; i++) {
-            if ((epsil != 0.0f) && (i % 5 == 0)) copyIm(Output, Output_prev, (long)(dimX), (long)(dimY), (long)(dimZ));
-
-            /* calculate differences */
-            D1_func(Output, D1, (long)(dimX), (long)(dimY), (long)(dimZ));
-            D2_func(Output, D2, (long)(dimX), (long)(dimY), (long)(dimZ));
-            if (dimZ > 1) D3_func(Output, D3, (long)(dimX), (long)(dimY), (long)(dimZ));
-            TV_kernel(D1, D2, D3, Output, Input, lambdaPar, lambda_is_arr, tau, (long)(dimX), (long)(dimY), (long)(dimZ));
-
-            /* check early stopping criteria */
-            if ((epsil != 0.0f) && (i % 5 == 0)) {
+        if ((epsil != 0.0f) && (i % 5 == 0)) copyIm(Output, Output_prev, (long)(dimX), (long)(dimY), (long)(dimZ));
+        
+        /* calculate differences */
+        D1_func(Output, D1, (long)(dimX), (long)(dimY), (long)(dimZ));
+        D2_func(Output, D2, (long)(dimX), (long)(dimY), (long)(dimZ));
+        if (dimZ > 1) D3_func(Output, D3, (long)(dimX), (long)(dimY), (long)(dimZ));
+        TV_kernel(D1, D2, D3, Output, Input, lambdaPar, lambda_is_arr, tau, (long)(dimX), (long)(dimY), (long)(dimZ));
+        
+        /* check early stopping criteria */
+        if ((epsil != 0.0f) && (i % 5 == 0)) {
             re = 0.0f; re1 = 0.0f;
-	            for(j=0; j<DimTotal; j++)
-        	    {
-        	        re += powf(Output[j] - Output_prev[j],2);
-        	        re1 += powf(Output[j],2);
-        	    }
-              re = sqrtf(re)/sqrtf(re1);
-              if (re < epsil)  count++;
-              if (count > 3) break;
+            for(j=0; j<DimTotal; j++)
+            {
+                re += powf(Output[j] - Output_prev[j],2);
+                re1 += powf(Output[j],2);
             }
-		}
+            re = sqrtf(re)/sqrtf(re1);
+            if (re < epsil)  count++;
+            if (count > 3) break;
+        }
+    }
     free(D1);free(D2); free(D3);
     if (epsil != 0.0f) free(Output_prev);
-
+    
     /*adding info into info_vector */
     infovector[0] = (float)(i);  /*iterations number (if stopped earlier based on tolerance)*/
     infovector[1] = re;  /* reached tolerance */
-
-	return 0;
+    
+    return 0;
 }
 
 /* calculate differences 1 */
@@ -103,13 +103,13 @@ float D1_func(float *A, float *D1, long dimX, long dimY, long dimZ)
 {
     float NOMx_1, NOMy_1, NOMy_0, NOMz_1, NOMz_0, denom1, denom2,denom3, T1;
     long i,j,k,i1,i2,k1,j1,j2,k2,index;
-
+    
     if (dimZ > 1) {
 #pragma omp parallel for shared (A, D1, dimX, dimY, dimZ) private(index, i, j, k, i1, j1, k1, i2, j2, k2, NOMx_1,NOMy_1,NOMy_0,NOMz_1,NOMz_0,denom1,denom2,denom3,T1)
-        for(j=0; j<dimY; j++) {
-            for(i=0; i<dimX; i++) {
-                for(k=0; k<dimZ; k++) {
-					index = (dimX*dimY)*k + j*dimX+i;
+        for(k=0; k<dimZ; k++) {
+            for(j=0; j<dimY; j++) {
+                for(i=0; i<dimX; i++) {
+                    index = (dimX*dimY)*k + j*dimX+i;
                     /* symmetric boundary conditions (Neuman) */
                     i1 = i + 1; if (i1 >= dimX) i1 = i-1;
                     i2 = i - 1; if (i2 < 0) i2 = i+1;
@@ -117,17 +117,17 @@ float D1_func(float *A, float *D1, long dimX, long dimY, long dimZ)
                     j2 = j - 1; if (j2 < 0) j2 = j+1;
                     k1 = k + 1; if (k1 >= dimZ) k1 = k-1;
                     k2 = k - 1; if (k2 < 0) k2 = k+1;
-
+                    
                     /* Forward-backward differences */
                     NOMx_1 = A[(dimX*dimY)*k + j1*dimX + i] - A[index]; /* x+ */
                     NOMy_1 = A[(dimX*dimY)*k + j*dimX + i1] - A[index]; /* y+ */
                     /*NOMx_0 = (A[(i)*dimY + j] - A[(i2)*dimY + j]); */  /* x- */
                     NOMy_0 = A[index] - A[(dimX*dimY)*k + j*dimX + i2]; /* y- */
-
+                    
                     NOMz_1 = A[(dimX*dimY)*k1 + j*dimX + i] - A[index]; /* z+ */
                     NOMz_0 = A[index] - A[(dimX*dimY)*k2 + j*dimX + i]; /* z- */
-
-
+                    
+                    
                     denom1 = NOMx_1*NOMx_1;
                     denom2 = 0.5f*(sign(NOMy_1) + sign(NOMy_0))*(MIN(fabs(NOMy_1),fabs(NOMy_0)));
                     denom2 = denom2*denom2;
@@ -141,19 +141,19 @@ float D1_func(float *A, float *D1, long dimX, long dimY, long dimZ)
 #pragma omp parallel for shared (A, D1, dimX, dimY) private(i, j, i1, j1, i2, j2,NOMx_1,NOMy_1,NOMy_0,denom1,denom2,T1,index)
         for(j=0; j<dimY; j++) {
             for(i=0; i<dimX; i++) {
-				index = j*dimX+i;
+                index = j*dimX+i;
                 /* symmetric boundary conditions (Neuman) */
                 i1 = i + 1; if (i1 >= dimX) i1 = i-1;
                 i2 = i - 1; if (i2 < 0) i2 = i+1;
                 j1 = j + 1; if (j1 >= dimY) j1 = j-1;
                 j2 = j - 1; if (j2 < 0) j2 = j+1;
-
+                
                 /* Forward-backward differences */
                 NOMx_1 = A[j1*dimX + i] - A[index]; /* x+ */
                 NOMy_1 = A[j*dimX + i1] - A[index]; /* y+ */
                 /*NOMx_0 = (A[(i)*dimY + j] - A[(i2)*dimY + j]); */ /* x- */
                 NOMy_0 = A[index] - A[(j)*dimX + i2]; /* y- */
-
+                
                 denom1 = NOMx_1*NOMx_1;
                 denom2 = 0.5f*(sign(NOMy_1) + sign(NOMy_0))*(MIN(fabs(NOMy_1),fabs(NOMy_0)));
                 denom2 = denom2*denom2;
@@ -168,12 +168,12 @@ float D2_func(float *A, float *D2, long dimX, long dimY, long dimZ)
 {
     float NOMx_1, NOMy_1, NOMx_0, NOMz_1, NOMz_0, denom1, denom2, denom3, T2;
     long i,j,k,i1,i2,k1,j1,j2,k2,index;
-
+    
     if (dimZ > 1) {
 #pragma omp parallel for shared (A, D2, dimX, dimY, dimZ) private(index, i, j, k, i1, j1, k1, i2, j2, k2,  NOMx_1, NOMy_1, NOMx_0, NOMz_1, NOMz_0, denom1, denom2, denom3, T2)
-        for(j=0; j<dimY; j++) {
-            for(i=0; i<dimX; i++) {
-                for(k=0; k<dimZ; k++) {
+        for(k=0; k<dimZ; k++) {
+            for(j=0; j<dimY; j++) {
+                for(i=0; i<dimX; i++) {
                     index = (dimX*dimY)*k + j*dimX+i;
                     /* symmetric boundary conditions (Neuman) */
                     i1 = i + 1; if (i1 >= dimX) i1 = i-1;
@@ -182,15 +182,15 @@ float D2_func(float *A, float *D2, long dimX, long dimY, long dimZ)
                     j2 = j - 1; if (j2 < 0) j2 = j+1;
                     k1 = k + 1; if (k1 >= dimZ) k1 = k-1;
                     k2 = k - 1; if (k2 < 0) k2 = k+1;
-
+                    
                     /* Forward-backward differences */
                     NOMx_1 = A[(dimX*dimY)*k + (j1)*dimX + i] - A[index]; /* x+ */
                     NOMy_1 = A[(dimX*dimY)*k + (j)*dimX + i1] - A[index]; /* y+ */
                     NOMx_0 = A[index] - A[(dimX*dimY)*k + (j2)*dimX + i]; /* x- */
                     NOMz_1 = A[(dimX*dimY)*k1 + j*dimX + i] - A[index]; /* z+ */
                     NOMz_0 = A[index] - A[(dimX*dimY)*k2 + (j)*dimX + i]; /* z- */
-
-
+                    
+                    
                     denom1 = NOMy_1*NOMy_1;
                     denom2 = 0.5f*(sign(NOMx_1) + sign(NOMx_0))*(MIN(fabs(NOMx_1),fabs(NOMx_0)));
                     denom2 = denom2*denom2;
@@ -204,19 +204,19 @@ float D2_func(float *A, float *D2, long dimX, long dimY, long dimZ)
 #pragma omp parallel for shared (A, D2, dimX, dimY) private(i, j, i1, j1, i2, j2, NOMx_1,NOMy_1,NOMx_0,denom1,denom2,T2,index)
         for(j=0; j<dimY; j++) {
             for(i=0; i<dimX; i++) {
-		index = j*dimX+i;
+                index = j*dimX+i;
                 /* symmetric boundary conditions (Neuman) */
                 i1 = i + 1; if (i1 >= dimX) i1 = i-1;
                 i2 = i - 1; if (i2 < 0) i2 = i+1;
                 j1 = j + 1; if (j1 >= dimY) j1 = j-1;
                 j2 = j - 1; if (j2 < 0) j2 = j+1;
-
+                
                 /* Forward-backward differences */
                 NOMx_1 = A[j1*dimX + i] - A[index]; /* x+ */
                 NOMy_1 = A[j*dimX + i1] - A[index]; /* y+ */
                 NOMx_0 = A[index] - A[j2*dimX + i]; /* x- */
                 /*NOMy_0 = A[(i)*dimY + j] - A[(i)*dimY + j2]; */  /* y- */
-
+                
                 denom1 = NOMy_1*NOMy_1;
                 denom2 = 0.5f*(sign(NOMx_1) + sign(NOMx_0))*(MIN(fabs(NOMx_1),fabs(NOMx_0)));
                 denom2 = denom2*denom2;
@@ -232,12 +232,12 @@ float D3_func(float *A, float *D3, long dimX, long dimY, long dimZ)
 {
     float NOMx_1, NOMy_1, NOMx_0, NOMy_0, NOMz_1, denom1, denom2, denom3, T3;
     long index,i,j,k,i1,i2,k1,j1,j2,k2;
-
+    
 #pragma omp parallel for shared (A, D3, dimX, dimY, dimZ) private(index, i, j, k, i1, j1, k1, i2, j2, k2,  NOMx_1, NOMy_1, NOMy_0, NOMx_0, NOMz_1, denom1, denom2, denom3, T3)
-    for(j=0; j<dimY; j++) {
-        for(i=0; i<dimX; i++) {
-            for(k=0; k<dimZ; k++) {
-				index = (dimX*dimY)*k + j*dimX+i;
+    for(k=0; k<dimZ; k++) {
+        for(j=0; j<dimY; j++) {
+            for(i=0; i<dimX; i++) {
+                index = (dimX*dimY)*k + j*dimX+i;
                 /* symmetric boundary conditions (Neuman) */
                 i1 = i + 1; if (i1 >= dimX) i1 = i-1;
                 i2 = i - 1; if (i2 < 0) i2 = i+1;
@@ -245,7 +245,7 @@ float D3_func(float *A, float *D3, long dimX, long dimY, long dimZ)
                 j2 = j - 1; if (j2 < 0) j2 = j+1;
                 k1 = k + 1; if (k1 >= dimZ) k1 = k-1;
                 k2 = k - 1; if (k2 < 0) k2 = k+1;
-
+                
                 /* Forward-backward differences */
                 NOMx_1 = A[(dimX*dimY)*k + (j1)*dimX + i] - A[index]; /* x+ */
                 NOMy_1 = A[(dimX*dimY)*k + (j)*dimX + i1] - A[index]; /* y+ */
@@ -253,7 +253,7 @@ float D3_func(float *A, float *D3, long dimX, long dimY, long dimZ)
                 NOMx_0 = A[index] - A[(dimX*dimY)*k + (j2)*dimX + i]; /* x- */
                 NOMz_1 = A[(dimX*dimY)*k1 + j*dimX + i] - A[index]; /* z+ */
                 /*NOMz_0 = A[(dimX*dimY)*k + (i)*dimY + j] - A[(dimX*dimY)*k2 + (i)*dimY + j]; */ /* z- */
-
+                
                 denom1 = NOMz_1*NOMz_1;
                 denom2 = 0.5f*(sign(NOMx_1) + sign(NOMx_0))*(MIN(fabs(NOMx_1),fabs(NOMx_0)));
                 denom2 = denom2*denom2;
@@ -270,12 +270,12 @@ float TV_kernel(float *D1, float *D2, float *D3, float *B, float *A, float *lamb
 {
     float dv1, dv2, dv3, lambda_val;
     long index,i,j,k,i1,i2,k1,j1,j2,k2;
-
+    
     if (dimZ > 1) {
 #pragma omp parallel for shared (D1, D2, D3, B, dimX, dimY, dimZ) private(index, i, j, k, i1, j1, k1, i2, j2, k2, dv1,dv2,dv3,lambda_val)
-        for(j=0; j<dimY; j++) {
-            for(i=0; i<dimX; i++) {
-                for(k=0; k<dimZ; k++) {
+        for(k=0; k<dimZ; k++) {
+            for(j=0; j<dimY; j++) {
+                for(i=0; i<dimX; i++) {
                     index = (dimX*dimY)*k + j*dimX+i;
                     lambda_val = *(lambda + index* lambda_is_arr);
                     /* symmetric boundary conditions (Neuman) */
@@ -285,12 +285,12 @@ float TV_kernel(float *D1, float *D2, float *D3, float *B, float *A, float *lamb
                     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;
-
+                    
                     /*divergence components */
                     dv1 = D1[index] - D1[(dimX*dimY)*k + j2*dimX+i];
                     dv2 = D2[index] - D2[(dimX*dimY)*k + j*dimX+i2];
                     dv3 = D3[index] - D3[(dimX*dimY)*k2 + j*dimX+i];
-
+                    
                     B[index] += tau*(lambda_val*(dv1 + dv2 + dv3) - (B[index] - A[index]));
                 }}}
     }
@@ -305,11 +305,11 @@ float TV_kernel(float *D1, float *D2, float *D3, float *B, float *A, float *lamb
                 i2 = i - 1; if (i2 < 0) i2 = i+1;
                 j1 = j + 1; if (j1 >= dimY) j1 = j-1;
                 j2 = j - 1; if (j2 < 0) j2 = j+1;
-
+                
                 /* divergence components  */
                 dv1 = D1[index] - D1[j2*dimX + i];
                 dv2 = D2[index] - D2[j*dimX + i2];
-
+                
                 B[index] += tau*(lambda_val*(dv1 + dv2) - (B[index] - A[index]));
             }}
     }
diff --git a/src/Core/regularisers_CPU/SB_TV_core.c b/src/Core/regularisers_CPU/SB_TV_core.c
index 8d80787..1215baf 100755
--- a/src/Core/regularisers_CPU/SB_TV_core.c
+++ b/src/Core/regularisers_CPU/SB_TV_core.c
@@ -1,146 +1,146 @@
 /*
-This work is part of the Core Imaging Library developed by
-Visual Analytics and Imaging System Group of the Science Technology
-Facilities Council, STFC
-
-Copyright 2017 Daniil Kazantsev
-Copyright 2017 Srikanth Nagella, Edoardo Pasca
-
-Licensed under the Apache License, Version 2.0 (the "License");
-you may not use this file except in compliance with the License.
-You may obtain a copy of the License at
-http://www.apache.org/licenses/LICENSE-2.0
-Unless required by applicable law or agreed to in writing, software
-distributed under the License is distributed on an "AS IS" BASIS,
-WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-See the License for the specific language governing permissions and
-limitations under the License.
-*/
+ * This work is part of the Core Imaging Library developed by
+ * Visual Analytics and Imaging System Group of the Science Technology
+ * Facilities Council, STFC
+ *
+ * Copyright 2017 Daniil Kazantsev
+ * Copyright 2017 Srikanth Nagella, Edoardo Pasca
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ * http://www.apache.org/licenses/LICENSE-2.0
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
 
 #include "SB_TV_core.h"
 
 /* C-OMP implementation of Split Bregman - TV denoising-regularisation model (2D/3D) [1]
-*
-* Input Parameters:
-* 1. Noisy image/volume
-* 2. lambda - regularisation parameter
-* 3. Number of iterations [OPTIONAL parameter]
-* 4. eplsilon - tolerance constant [OPTIONAL parameter]
-* 5. TV-type: 'iso' or 'l1' [OPTIONAL parameter]
-*
-* Output:
-* [1] Filtered/regularized image/volume
-* [2] Information vector which contains [iteration no., reached tolerance]
-*
-* [1]. Goldstein, T. and Osher, S., 2009. The split Bregman method for L1-regularized problems. SIAM journal on imaging sciences, 2(2), pp.323-343.
-*/
+ *
+ * Input Parameters:
+ * 1. Noisy image/volume
+ * 2. lambda - regularisation parameter
+ * 3. Number of iterations [OPTIONAL parameter]
+ * 4. eplsilon - tolerance constant [OPTIONAL parameter]
+ * 5. TV-type: 'iso' or 'l1' [OPTIONAL parameter]
+ *
+ * Output:
+ * [1] Filtered/regularized image/volume
+ * [2] Information vector which contains [iteration no., reached tolerance]
+ *
+ * [1]. Goldstein, T. and Osher, S., 2009. The split Bregman method for L1-regularized problems. SIAM journal on imaging sciences, 2(2), pp.323-343.
+ */
 
 float SB_TV_CPU_main(float *Input, float *Output, float *infovector, float mu, int iter, float epsil, int methodTV, int dimX, int dimY, int dimZ)
 {
-	  int ll;
+    int ll;
     long j, DimTotal;
-	  float re, re1, lambda;
+    float re, re1, lambda;
     re = 0.0f; re1 = 0.0f;
     int count = 0;
     mu = 1.0f/mu;
     lambda = 2.0f*mu;
-
+    
     float *Output_prev=NULL, *Dx=NULL, *Dy=NULL, *Bx=NULL, *By=NULL;
     DimTotal = (long)(dimX*dimY*dimZ);
     Output_prev = calloc(DimTotal, sizeof(float));
-		Dx = calloc(DimTotal, sizeof(float));
-		Dy = calloc(DimTotal, sizeof(float));
-		Bx = calloc(DimTotal, sizeof(float));
-		By = calloc(DimTotal, sizeof(float));
-
-	if (dimZ == 1) {
-		    /* 2D case */
+    Dx = calloc(DimTotal, sizeof(float));
+    Dy = calloc(DimTotal, sizeof(float));
+    Bx = calloc(DimTotal, sizeof(float));
+    By = calloc(DimTotal, sizeof(float));
+    
+    if (dimZ == 1) {
+        /* 2D case */
         copyIm(Input, Output, (long)(dimX), (long)(dimY), 1l); /*initialize */
-
+        
         /* begin outer SB iterations */
         for(ll=0; ll<iter; ll++) {
-
+            
             /* storing old estimate */
             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, (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, (long)(dimX), (long)(dimY), lambda, mu);
-
+            
             /* TV-related step */
             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, (long)(dimX), (long)(dimY));
-
+            
             /* check early stopping criteria if epsilon not equal zero */
             if ((epsil != 0.0f)  && (ll % 5 == 0)) {
-            re = 0.0f; re1 = 0.0f;
-	            for(j=0; j<DimTotal; j++)
-        	    {
-        	        re += powf(Output[j] - Output_prev[j],2);
-        	        re1 += powf(Output[j],2);
-        	    }
-            re = sqrtf(re)/sqrtf(re1);
-            /* stop if the norm residual is less than the tolerance EPS */
-            if (re < epsil)  count++;
-            if (count > 3) break;
+                re = 0.0f; re1 = 0.0f;
+                for(j=0; j<DimTotal; j++)
+                {
+                    re += powf(Output[j] - Output_prev[j],2);
+                    re1 += powf(Output[j],2);
+                }
+                re = sqrtf(re)/sqrtf(re1);
+                /* stop if the norm residual is less than the tolerance EPS */
+                if (re < epsil)  count++;
+                if (count > 3) break;
             }
         }
-	}
-	else {
-		/* 3D case */
-		float *Dz=NULL, *Bz=NULL;
-
-		Dz = calloc(DimTotal, sizeof(float));
-		Bz = calloc(DimTotal, sizeof(float));
-
+    }
+    else {
+        /* 3D case */
+        float *Dz=NULL, *Bz=NULL;
+        
+        Dz = calloc(DimTotal, sizeof(float));
+        Bz = calloc(DimTotal, sizeof(float));
+        
         copyIm(Input, Output, (long)(dimX), (long)(dimY), (long)(dimZ)); /*initialize */
-
+        
         /* begin outer SB iterations */
         for(ll=0; ll<iter; ll++) {
-
+            
             /* storing old estimate */
             copyIm(Output, Output_prev, (long)(dimX), (long)(dimY), (long)(dimZ));
-
-             /* perform two GS iterations (normally 2 is enough for the convergence) */
+            
+            /* perform two GS iterations (normally 2 is enough for the convergence) */
             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, (long)(dimX), (long)(dimY), (long)(dimZ), lambda, mu);
-
+            
             /* TV-related step */
             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, (long)(dimX), (long)(dimY), (long)(dimZ));
-
+            
             /* check early stopping criteria if epsilon not equal zero */
             if ((epsil != 0.0f)  && (ll % 5 == 0)) {
-            re = 0.0f; re1 = 0.0f;
-              for(j=0; j<DimTotal; j++)
-              {
-                  re += powf(Output[j] - Output_prev[j],2);
-                  re1 += powf(Output[j],2);
-              }
-            re = sqrtf(re)/sqrtf(re1);
-            /* stop if the norm residual is less than the tolerance EPS */
-            if (re < epsil)  count++;
-            if (count > 3) break;
+                re = 0.0f; re1 = 0.0f;
+                for(j=0; j<DimTotal; j++)
+                {
+                    re += powf(Output[j] - Output_prev[j],2);
+                    re1 += powf(Output[j],2);
+                }
+                re = sqrtf(re)/sqrtf(re1);
+                /* stop if the norm residual is less than the tolerance EPS */
+                if (re < epsil)  count++;
+                if (count > 3) break;
             }
         }
-		free(Dz); free(Bz);
-	}
-
-  free(Output_prev); free(Dx); free(Dy); free(Bx); free(By);
-  /*adding info into info_vector */
-  infovector[0] = (float)(ll);  /*iterations number (if stopped earlier based on tolerance)*/
-  infovector[1] = re;  /* reached tolerance */
-	return 0;
+        free(Dz); free(Bz);
+    }
+    
+    free(Output_prev); free(Dx); free(Dy); free(Bx); free(By);
+    /*adding info into info_vector */
+    infovector[0] = (float)(ll);  /*iterations number (if stopped earlier based on tolerance)*/
+    infovector[1] = re;  /* reached tolerance */
+    return 0;
 }
 
 /********************************************************************/
@@ -151,18 +151,18 @@ float gauss_seidel2D(float *U, float *A, float *U_prev, float *Dx, float *Dy, fl
     float sum, normConst;
     long i,j,i1,i2,j1,j2,index;
     normConst = 1.0f/(mu + 4.0f*lambda);
-
+    
 #pragma omp parallel for shared(U) private(index,i,j,i1,i2,j1,j2,sum)
-    for(i=0; i<dimX; i++) {
+    for(j=0; j<dimY; j++) {
         /* symmetric boundary conditions (Neuman) */
-        i1 = i+1; if (i1 == dimX) i1 = i-1;
-        i2 = i-1; if (i2 < 0) i2 = i+1;
-        for(j=0; j<dimY; j++) {
+        j1 = j+1; if (j1 == dimY) j1 = j-1;
+        j2 = j-1; if (j2 < 0) j2 = j+1;
+        for(i=0; i<dimX; i++) {
             /* symmetric boundary conditions (Neuman) */
-            j1 = j+1; if (j1 == dimY) j1 = j-1;
-            j2 = j-1; if (j2 < 0) j2 = j+1;
+            i1 = i+1; if (i1 == dimX) i1 = i-1;
+            i2 = i-1; if (i2 < 0) i2 = i+1;
             index = j*dimX+i;
-
+            
             sum = Dx[j*dimX+i2] - Dx[index] + Dy[j2*dimX+i] - Dy[index] - Bx[j*dimX+i2] + Bx[index] - By[j2*dimX+i] + By[index];
             sum += U_prev[j*dimX+i1] + U_prev[j*dimX+i2] + U_prev[j1*dimX+i] + U_prev[j2*dimX+i];
             sum *= lambda;
@@ -178,22 +178,23 @@ float updDxDy_shrinkAniso2D(float *U, float *Dx, float *Dy, float *Bx, float *By
     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)
-    for(i=0; i<dimX; i++) {
-        for(j=0; j<dimY; j++) {
+    for(j=0; j<dimY; j++) {
+        j1 = j+1; if (j1 == dimY) j1 = j-1;
+        for(i=0; i<dimX; i++) {
             /* symmetric boundary conditions (Neuman) */
             i1 = i+1; if (i1 == dimX) i1 = i-1;
-            j1 = j+1; if (j1 == dimY) j1 = j-1;
+            
             index = j*dimX+i;
-
+            
             val1 = (U[j*dimX+i1] - U[index]) + Bx[index];
             val2 = (U[j1*dimX+i] - U[index]) + By[index];
-
+            
             val11 = fabs(val1) - denom_lam; if (val11 < 0) val11 = 0;
             val22 = fabs(val2) - denom_lam; if (val22 < 0) val22 = 0;
-
+            
             if (val1 !=0) Dx[index] = (val1/fabs(val1))*val11; else Dx[index] = 0;
             if (val2 !=0) Dy[index] = (val2/fabs(val2))*val22; else Dy[index] = 0;
-
+            
         }}
     return 1;
 }
@@ -202,22 +203,22 @@ float updDxDy_shrinkIso2D(float *U, float *Dx, float *Dy, float *Bx, float *By,
     long i,j,i1,j1,index;
     float val1, val11, val2, denom, denom_lam;
     denom_lam = 1.0f/lambda;
-
+    
 #pragma omp parallel for shared(U,denom_lam) private(index,i,j,i1,j1,val1,val11,val2,denom)
-    for(i=0; i<dimX; i++) {
-        for(j=0; j<dimY; j++) {
+    for(j=0; j<dimY; j++) {
+        j1 = j+1; if (j1 == dimY) j1 = j-1;
+        for(i=0; i<dimX; i++) {
             /* symmetric boundary conditions (Neuman) */
             i1 = i+1; if (i1 == dimX) i1 = i-1;
-            j1 = j+1; if (j1 == dimY) j1 = j-1;
             index = j*dimX+i;
-
+            
             val1 = (U[j*dimX+i1] - U[index]) + Bx[index];
             val2 = (U[j1*dimX+i] - U[index]) + By[index];
-
+            
             denom = sqrt(val1*val1 + val2*val2);
-
+            
             val11 = (denom - denom_lam); if (val11 < 0) val11 = 0.0f;
-
+            
             if (denom != 0.0f) {
                 Dx[index] = val11*(val1/denom);
                 Dy[index] = val11*(val2/denom);
@@ -233,13 +234,13 @@ float updBxBy2D(float *U, float *Dx, float *Dy, float *Bx, float *By, long dimX,
 {
     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++) {
+    for(j=0; j<dimY; j++) {
+        j1 = j+1; if (j1 == dimY) j1 = j-1;
+        for(i=0; i<dimX; i++) {
             /* symmetric boundary conditions (Neuman) */
             i1 = i+1; if (i1 == dimX) i1 = i-1;
-            j1 = j+1; if (j1 == dimY) j1 = j-1;
             index = j*dimX+i;
-
+            
             Bx[index] += (U[j*dimX+i1] - U[index]) - Dx[index];
             By[index] += (U[j1*dimX+i] - U[index]) - Dy[index];
         }}
@@ -256,18 +257,19 @@ float gauss_seidel3D(float *U, float *A, float *U_prev, float *Dx, float *Dy, fl
     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++) {
+    for(k=0; k<dimZ; k++) {
+        k1 = k+1; if (k1 == dimZ) k1 = k-1;
+        k2 = k-1; if (k2 < 0) k2 = k+1;
         for(j=0; j<dimY; j++) {
-            for(k=0; k<dimZ; k++) {
+            j1 = j+1; if (j1 == dimY) j1 = j-1;
+            j2 = j-1; if (j2 < 0) j2 = j+1;
+            for(i=0; i<dimX; i++) {
                 /* symmetric boundary conditions (Neuman) */
                 i1 = i+1; if (i1 == dimX) i1 = i-1;
                 i2 = i-1; if (i2 < 0) i2 = i+1;
-                j1 = j+1; if (j1 == dimY) j1 = j-1;
-                j2 = j-1; if (j2 < 0) j2 = j+1;
-                k1 = k+1; if (k1 == dimZ) k1 = k-1;
-                k2 = k-1; if (k2 < 0) k2 = k+1;
+                
                 index = (dimX*dimY)*k + j*dimX+i;
-
+                
                 d_val = Dx[(dimX*dimY)*k + j*dimX+i2] - Dx[index] + Dy[(dimX*dimY)*k + j2*dimX+i] - Dy[index] + Dz[(dimX*dimY)*k2 + j*dimX+i] - Dz[index];
                 b_val = -Bx[(dimX*dimY)*k + j*dimX+i2] + Bx[index] - By[(dimX*dimY)*k + j2*dimX+i] + By[index] - Bz[(dimX*dimY)*k2 + j*dimX+i] + Bz[index];
                 sum = d_val + b_val;
@@ -285,27 +287,28 @@ float updDxDyDz_shrinkAniso3D(float *U, float *Dx, float *Dy, float *Dz, float *
     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)
-    for(i=0; i<dimX; i++) {
+    for(k=0; k<dimZ; k++) {
+        k1 = k+1; if (k1 == dimZ) k1 = k-1;
         for(j=0; j<dimY; j++) {
-            for(k=0; k<dimZ; k++) {
+            j1 = j+1; if (j1 == dimY) j1 = j-1;
+            for(i=0; i<dimX; i++) {
+                
                 index = (dimX*dimY)*k + j*dimX+i;
                 /* symmetric boundary conditions (Neuman) */
                 i1 = i+1; if (i1 == dimX) i1 = i-1;
-                j1 = j+1; if (j1 == dimY) j1 = j-1;
-                k1 = k+1; if (k1 == dimZ) k1 = k-1;
-
+                
                 val1 = (U[(dimX*dimY)*k + j*dimX+i1] - U[index]) + Bx[index];
                 val2 = (U[(dimX*dimY)*k + j1*dimX+i] - U[index]) + By[index];
                 val3 = (U[(dimX*dimY)*k1 + j*dimX+i] - U[index]) + Bz[index];
-
+                
                 val11 = fabs(val1) - denom_lam; if (val11 < 0.0f) val11 = 0.0f;
                 val22 = fabs(val2) - denom_lam; if (val22 < 0.0f) val22 = 0.0f;
                 val33 = fabs(val3) - denom_lam; if (val33 < 0.0f) val33 = 0.0f;
-
+                
                 if (val1 !=0.0f) Dx[index] = (val1/fabs(val1))*val11; else Dx[index] = 0.0f;
                 if (val2 !=0.0f) Dy[index] = (val2/fabs(val2))*val22; else Dy[index] = 0.0f;
                 if (val3 !=0.0f) Dz[index] = (val3/fabs(val3))*val33; else Dz[index] = 0.0f;
-
+                
             }}}
     return 1;
 }
@@ -315,23 +318,25 @@ float updDxDyDz_shrinkIso3D(float *U, float *Dx, float *Dy, float *Dz, float *Bx
     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)
-    for(i=0; i<dimX; i++) {
+    for(k=0; k<dimZ; k++) {
+        k1 = k+1; if (k1 == dimZ) k1 = k-1;
         for(j=0; j<dimY; j++) {
-            for(k=0; k<dimZ; k++) {
-                index = (dimX*dimY)*k + j*dimX+i;
+            j1 = j+1; if (j1 == dimY) j1 = j-1;
+            for(i=0; i<dimX; i++) {          
+                
                 /* symmetric boundary conditions (Neuman) */
-                i1 = i+1; if (i1 == dimX) i1 = i-1;
-                j1 = j+1; if (j1 == dimY) j1 = j-1;
-                k1 = k+1; if (k1 == dimZ) k1 = k-1;
-
+                i1 = i+1; if (i1 == dimX) i1 = i-1;                
+                
+                index = (dimX*dimY)*k + j*dimX+i;
+                
                 val1 = (U[(dimX*dimY)*k + j*dimX+i1] - U[index]) + Bx[index];
                 val2 = (U[(dimX*dimY)*k + j1*dimX+i] - U[index]) + By[index];
                 val3 = (U[(dimX*dimY)*k1 + j*dimX+i] - U[index]) + Bz[index];
-
+                
                 denom = sqrt(val1*val1 + val2*val2 + val3*val3);
-
+                
                 val11 = (denom - denom_lam); if (val11 < 0) val11 = 0.0f;
-
+                
                 if (denom != 0.0f) {
                     Dx[index] = val11*(val1/denom);
                     Dy[index] = val11*(val2/denom);
@@ -349,15 +354,15 @@ float updBxByBz3D(float *U, float *Dx, float *Dy, float *Dz, float *Bx, float *B
 {
     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(k=0; k<dimZ; k++) {
+        k1 = k+1; if (k1 == dimZ) k1 = k-1;
         for(j=0; j<dimY; j++) {
-            for(k=0; k<dimZ; k++) {
-				index = (dimX*dimY)*k + j*dimX+i;
+            j1 = j+1; if (j1 == dimY) j1 = j-1;
+            for(i=0; i<dimX; i++) {
+                index = (dimX*dimY)*k + j*dimX+i;
                 /* symmetric boundary conditions (Neuman) */
                 i1 = i+1; if (i1 == dimX) i1 = i-1;
-                j1 = j+1; if (j1 == dimY) j1 = j-1;
-                k1 = k+1; if (k1 == dimZ) k1 = k-1;
-
+                
                 Bx[index] += (U[(dimX*dimY)*k + j*dimX+i1] - U[index]) - Dx[index];
                 By[index] += (U[(dimX*dimY)*k + j1*dimX+i] - U[index]) - Dy[index];
                 Bz[index] += (U[(dimX*dimY)*k1 + j*dimX+i] - U[index]) - Dz[index];
diff --git a/src/Core/regularisers_CPU/TGV_core.c b/src/Core/regularisers_CPU/TGV_core.c
index f43b56a..3f917de 100644
--- a/src/Core/regularisers_CPU/TGV_core.c
+++ b/src/Core/regularisers_CPU/TGV_core.c
@@ -48,148 +48,148 @@ float TGV_main(float *U0, float *U, float *infovector, float lambda, float alpha
     re = 0.0f; re1 = 0.0f;
     int count = 0;
     float *U_old, *P1, *P2, *Q1, *Q2, *Q3, *V1, *V1_old, *V2, *V2_old, tau, sigma;
-
+    
     DimTotal = (long)(dimX*dimY*dimZ);
     copyIm(U0, U, (long)(dimX), (long)(dimY), (long)(dimZ)); /* initialize */
     tau = pow(L2,-0.5);
     sigma = pow(L2,-0.5);
-
+    
     /* dual variables */
     P1 = calloc(DimTotal, sizeof(float));
     P2 = calloc(DimTotal, sizeof(float));
-
+    
     Q1 = calloc(DimTotal, sizeof(float));
     Q2 = calloc(DimTotal, sizeof(float));
     Q3 = calloc(DimTotal, sizeof(float));
-
+    
     U_old = calloc(DimTotal, sizeof(float));
-
+    
     V1 = calloc(DimTotal, sizeof(float));
     V1_old = calloc(DimTotal, sizeof(float));
     V2 = calloc(DimTotal, sizeof(float));
     V2_old = calloc(DimTotal, sizeof(float));
-
+    
     if (dimZ == 1) {
         /*2D case*/
-
+        
         /* Primal-dual iterations begin here */
         for(ll = 0; ll < iter; ll++) {
-
+            
             /* Calculate Dual Variable P */
             DualP_2D(U, V1, V2, P1, P2, (long)(dimX), (long)(dimY), sigma);
-
+            
             /*Projection onto convex set for P*/
             ProjP_2D(P1, P2, (long)(dimX), (long)(dimY), alpha1);
-
+            
             /* Calculate Dual Variable Q */
             DualQ_2D(V1, V2, Q1, Q2, Q3, (long)(dimX), (long)(dimY), sigma);
-
+            
             /*Projection onto convex set for Q*/
             ProjQ_2D(Q1, Q2, Q3, (long)(dimX), (long)(dimY), alpha0);
-
+            
             /*saving U into U_old*/
             copyIm(U, U_old, (long)(dimX), (long)(dimY), 1l);
-
+            
             /*adjoint operation  -> divergence and projection of P*/
             DivProjP_2D(U, U0, P1, P2, (long)(dimX), (long)(dimY), lambda, tau);
-
+            
             /*get updated solution U*/
             newU(U, U_old, (long)(dimX), (long)(dimY));
-
+            
             /*saving V into V_old*/
             copyIm(V1, V1_old, (long)(dimX), (long)(dimY), 1l);
             copyIm(V2, V2_old, (long)(dimX), (long)(dimY), 1l);
-
+            
             /* upd V*/
             UpdV_2D(V1, V2, P1, P2, Q1, Q2, Q3, (long)(dimX), (long)(dimY), tau);
-
+            
             /*get new V*/
             newU(V1, V1_old, (long)(dimX), (long)(dimY));
             newU(V2, V2_old, (long)(dimX), (long)(dimY));
-
+            
             /* check early stopping criteria */
             if ((epsil != 0.0f)  && (ll % 5 == 0)) {
-            re = 0.0f; re1 = 0.0f;
-	            for(j=0; j<DimTotal; j++)
-        	    {
-        	        re += powf(U[j] - U_old[j],2);
-        	        re1 += powf(U[j],2);
-        	    }
-           re = sqrtf(re)/sqrtf(re1);
-           if (re < epsil)  count++;
-           if (count > 3) break;
-           }
+                re = 0.0f; re1 = 0.0f;
+                for(j=0; j<DimTotal; j++)
+                {
+                    re += powf(U[j] - U_old[j],2);
+                    re1 += powf(U[j],2);
+                }
+                re = sqrtf(re)/sqrtf(re1);
+                if (re < epsil)  count++;
+                if (count > 3) break;
+            }
         } /*end of iterations*/
     }
     else {
         /*3D case*/
         float *P3, *Q4, *Q5, *Q6, *V3, *V3_old;
-
+        
         P3 = calloc(DimTotal, sizeof(float));
         Q4 = calloc(DimTotal, sizeof(float));
         Q5 = calloc(DimTotal, sizeof(float));
         Q6 = calloc(DimTotal, sizeof(float));
         V3 = calloc(DimTotal, sizeof(float));
         V3_old = calloc(DimTotal, sizeof(float));
-
+        
         /* Primal-dual iterations begin here */
         for(ll = 0; ll < iter; ll++) {
-
+            
             /* Calculate Dual Variable P */
             DualP_3D(U, V1, V2, V3, P1, P2, P3, (long)(dimX), (long)(dimY), (long)(dimZ), sigma);
-
+            
             /*Projection onto convex set for P*/
             ProjP_3D(P1, P2, P3, (long)(dimX), (long)(dimY), (long)(dimZ), alpha1);
-
+            
             /* Calculate Dual Variable Q */
             DualQ_3D(V1, V2, V3, Q1, Q2, Q3, Q4, Q5, Q6, (long)(dimX), (long)(dimY), (long)(dimZ), sigma);
-
+            
             /*Projection onto convex set for Q*/
             ProjQ_3D(Q1, Q2, Q3, Q4, Q5, Q6, (long)(dimX), (long)(dimY), (long)(dimZ), alpha0);
-
+            
             /*saving U into U_old*/
             copyIm(U, U_old, (long)(dimX), (long)(dimY), (long)(dimZ));
-
+            
             /*adjoint operation  -> divergence and projection of P*/
             DivProjP_3D(U, U0, P1, P2, P3, (long)(dimX), (long)(dimY), (long)(dimZ), lambda, tau);
-
+            
             /*get updated solution U*/
             newU3D(U, U_old, (long)(dimX), (long)(dimY), (long)(dimZ));
-
+            
             /*saving V into V_old*/
             copyIm_3Ar(V1, V2, V3, V1_old, V2_old, V3_old, (long)(dimX), (long)(dimY), (long)(dimZ));
-
+            
             /* upd V*/
             UpdV_3D(V1, V2, V3, P1, P2, P3, Q1, Q2, Q3, Q4, Q5, Q6, (long)(dimX), (long)(dimY), (long)(dimZ), tau);
-
+            
             /*get new V*/
             newU3D_3Ar(V1, V2, V3, V1_old, V2_old, V3_old, (long)(dimX), (long)(dimY), (long)(dimZ));
-
+            
             /* check early stopping criteria */
             if ((epsil != 0.0f)  && (ll % 5 == 0)) {
-            re = 0.0f; re1 = 0.0f;
-	            for(j=0; j<DimTotal; j++)
-        	    {
-        	        re += powf(U[j] - U_old[j],2);
-        	        re1 += powf(U[j],2);
-        	    }
-           re = sqrtf(re)/sqrtf(re1);
-           if (re < epsil)  count++;
-           if (count > 3) break;
-           }
-
+                re = 0.0f; re1 = 0.0f;
+                for(j=0; j<DimTotal; j++)
+                {
+                    re += powf(U[j] - U_old[j],2);
+                    re1 += powf(U[j],2);
+                }
+                re = sqrtf(re)/sqrtf(re1);
+                if (re < epsil)  count++;
+                if (count > 3) break;
+            }
+            
         } /*end of iterations*/
         free(P3);free(Q4);free(Q5);free(Q6);free(V3);free(V3_old);
     }
-
+    
     /*freeing*/
     free(P1);free(P2);free(Q1);free(Q2);free(Q3);free(U_old);
     free(V1);free(V2);free(V1_old);free(V2_old);
-
+    
     /*adding info into info_vector */
     infovector[0] = (float)(ll);  /*iterations number (if stopped earlier based on tolerance)*/
     infovector[1] = re;  /* reached tolerance */
-
+    
     return 0;
 }
 
@@ -201,15 +201,15 @@ float DualP_2D(float *U, float *V1, float *V2, float *P1, float *P2, long dimX,
 {
     long i,j, index;
 #pragma omp parallel for shared(U,V1,V2,P1,P2) private(i,j,index)
-    for(i=0; i<dimX; i++) {
-        for(j=0; j<dimY; j++) {
+    for(j=0; j<dimY; j++) {
+        for(i=0; i<dimX; i++) {
             index = j*dimX+i;
             /* symmetric boundary conditions (Neuman) */
             if (i == dimX-1) P1[index] += sigma*(-V1[index]);
             else P1[index] += sigma*((U[j*dimX+(i+1)] - U[index])  - V1[index]);
             if (j == dimY-1) P2[index] += sigma*(-V2[index]);
             else  P2[index] += sigma*((U[(j+1)*dimX+i] - U[index])  - V2[index]);
-
+            
         }}
     return 1;
 }
@@ -219,8 +219,8 @@ float ProjP_2D(float *P1, float *P2, long dimX, long dimY, float alpha1)
     float grad_magn;
     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++) {
+    for(j=0; j<dimY; j++) {
+        for(i=0; i<dimX; i++) {
             index = j*dimX+i;
             grad_magn = (sqrtf(pow(P1[index],2) + pow(P2[index],2)))/alpha1;
             if (grad_magn > 1.0f) {
@@ -236,8 +236,8 @@ float DualQ_2D(float *V1, float *V2, float *Q1, float *Q2, float *Q3, long dimX,
     long i,j,index;
     float q1, q2, q11, q22;
 #pragma omp parallel for shared(Q1,Q2,Q3,V1,V2) private(i,j,index,q1,q2,q11,q22)
-    for(i=0; i<dimX; i++) {
-        for(j=0; j<dimY; j++) {
+    for(j=0; j<dimY; j++) {
+        for(i=0; i<dimX; i++) {
             index = j*dimX+i;
             q1 = 0.0f; q11 = 0.0f; q2 = 0.0f; q22 = 0.0f;
             /* boundary conditions (Neuman) */
@@ -260,8 +260,8 @@ float ProjQ_2D(float *Q1, float *Q2, float *Q3, long dimX, long dimY, float alph
     float grad_magn;
     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++) {
+    for(j=0; j<dimY; j++) {
+        for(i=0; i<dimX; i++) {
             index = j*dimX+i;
             grad_magn = sqrtf(pow(Q1[index],2) + pow(Q2[index],2) + 2*pow(Q3[index],2));
             grad_magn = grad_magn/alpha0;
@@ -279,18 +279,18 @@ float DivProjP_2D(float *U, float *U0, float *P1, float *P2, long dimX, long dim
     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++) {
-        for(j=0; j<dimY; j++) {
+    for(j=0; j<dimY; j++) {
+        for(i=0; i<dimX; i++) {
             index = j*dimX+i;
-
+            
             if (i == 0) P_v1 = P1[index];
             else if (i == dimX-1) P_v1 = -P1[j*dimX+(i-1)];
             else P_v1 = P1[index] - P1[j*dimX+(i-1)];
-
+            
             if (j == 0) P_v2 = P2[index];
             else if (j == dimY-1) P_v2 = -P2[(j-1)*dimX+i];
             else P_v2 = P2[index] - P2[(j-1)*dimX+i];
-
+            
             div = P_v1 + P_v2;
             U[index] = (lambda*(U[index] + tau*div) + tau*U0[index])/(lambda + tau);
         }}
@@ -310,10 +310,10 @@ float UpdV_2D(float *V1, float *V2, float *P1, float *P2, float *Q1, float *Q2,
     long i, j, index;
     float q1, q3_x, q3_y, q2, div1, div2;
 #pragma omp parallel for shared(V1,V2,P1,P2,Q1,Q2,Q3) private(i, j, index, q1, q3_x, q3_y, q2, div1, div2)
-    for(i=0; i<dimX; i++) {
-        for(j=0; j<dimY; j++) {
+    for(j=0; j<dimY; j++) {
+        for(i=0; i<dimX; i++) {
             index = j*dimX+i;
-
+            
             /* boundary conditions (Neuman) */
             if (i == 0) {
                 q1 = Q1[index];
@@ -324,7 +324,7 @@ float UpdV_2D(float *V1, float *V2, float *P1, float *P2, float *Q1, float *Q2,
             else {
                 q1 = Q1[index] - Q1[j*dimX+(i-1)];
                 q3_x = Q3[index] - Q3[j*dimX+(i-1)];  }
-
+            
             if (j == 0) {
                 q2 = Q2[index];
                 q3_y = Q3[index]; }
@@ -334,8 +334,8 @@ float UpdV_2D(float *V1, float *V2, float *P1, float *P2, float *Q1, float *Q2,
             else {
                 q2 = Q2[index] - Q2[(j-1)*dimX+i];
                 q3_y = Q3[index] - Q3[(j-1)*dimX+i]; }
-
-
+            
+            
             div1 = q1 + q3_y;
             div2 = q3_x + q2;
             V1[index] += tau*(P1[index] + div1);
@@ -352,9 +352,9 @@ float DualP_3D(float *U, float *V1, float *V2, float *V3, float *P1, float *P2,
 {
     long i,j,k, index;
 #pragma omp parallel for shared(U,V1,V2,V3,P1,P2,P3) private(i,j,k,index)
-    for(i=0; i<dimX; i++) {
+    for(k=0; k<dimZ; k++) {
         for(j=0; j<dimY; j++) {
-            for(k=0; k<dimZ; k++) {
+            for(i=0; i<dimX; i++) {
                 index = (dimX*dimY)*k + j*dimX+i;
                 /* symmetric boundary conditions (Neuman) */
                 if (i == dimX-1) P1[index] += sigma*(-V1[index]);
@@ -372,9 +372,9 @@ float ProjP_3D(float *P1, float *P2, float *P3, long dimX, long dimY, long dimZ,
     float grad_magn;
     long i,j,k,index;
 #pragma omp parallel for shared(P1,P2,P3) private(i,j,k,index,grad_magn)
-    for(i=0; i<dimX; i++) {
+    for(k=0; k<dimZ; k++) {
         for(j=0; j<dimY; j++) {
-            for(k=0; k<dimZ; k++) {
+            for(i=0; i<dimX; i++) {
                 index = (dimX*dimY)*k + j*dimX+i;
                 grad_magn = (sqrtf(pow(P1[index],2) + pow(P2[index],2) + pow(P3[index],2)))/alpha1;
                 if (grad_magn > 1.0f) {
@@ -391,9 +391,9 @@ float DualQ_3D(float *V1, float *V2, float *V3, float *Q1, float *Q2, float *Q3,
     long i,j,k,index;
     float q1, q2, q3, q11, q22, q33, q44, q55, q66;
 #pragma omp parallel for shared(Q1,Q2,Q3,Q4,Q5,Q6,V1,V2,V3) private(i,j,k,index,q1,q2,q3,q11,q22,q33,q44,q55,q66)
-    for(i=0; i<dimX; i++) {
+    for(k=0; k<dimZ; k++) {
         for(j=0; j<dimY; j++) {
-            for(k=0; k<dimZ; k++) {
+            for(i=0; i<dimX; i++) {
                 index = (dimX*dimY)*k + j*dimX+i;
                 q1 = 0.0f; q11 = 0.0f; q33 = 0.0f; q2 = 0.0f; q22 = 0.0f; q55 = 0.0f; q3 = 0.0f; q44 = 0.0f; q66 = 0.0f;
                 /* symmetric boundary conditions (Neuman) */
@@ -412,7 +412,7 @@ float DualQ_3D(float *V1, float *V2, float *V3, float *Q1, float *Q2, float *Q3,
                     q44 = V1[(dimX*dimY)*(k+1) + j*dimX+i] - V1[index];
                     q66 = V2[(dimX*dimY)*(k+1) + j*dimX+i] - V2[index];
                 }
-
+                
                 Q1[index] += sigma*(q1); /*Q11*/
                 Q2[index] += sigma*(q2); /*Q22*/
                 Q3[index] += sigma*(q3); /*Q33*/
@@ -427,9 +427,9 @@ float ProjQ_3D(float *Q1, float *Q2, float *Q3, float *Q4, float *Q5, float *Q6,
     float grad_magn;
     long i,j,k,index;
 #pragma omp parallel for shared(Q1,Q2,Q3,Q4,Q5,Q6) private(i,j,k,index,grad_magn)
-    for(i=0; i<dimX; i++) {
+    for(k=0; k<dimZ; k++) {
         for(j=0; j<dimY; j++) {
-            for(k=0; k<dimZ; k++) {
+            for(i=0; i<dimX; i++) {
                 index = (dimX*dimY)*k + j*dimX+i;
                 grad_magn = sqrtf(pow(Q1[index],2) + pow(Q2[index],2) + pow(Q3[index],2) + 2.0f*pow(Q4[index],2) + 2.0f*pow(Q5[index],2) + 2.0f*pow(Q6[index],2));
                 grad_magn = grad_magn/alpha0;
@@ -450,11 +450,11 @@ float DivProjP_3D(float *U, float *U0, float *P1, float *P2, float *P3, long dim
     long i,j,k,index;
     float P_v1, P_v2, P_v3, div;
 #pragma omp parallel for shared(U,U0,P1,P2,P3) private(i,j,k,index,P_v1,P_v2,P_v3,div)
-    for(i=0; i<dimX; i++) {
+    for(k=0; k<dimZ; k++) {
         for(j=0; j<dimY; j++) {
-            for(k=0; k<dimZ; k++) {
+            for(i=0; i<dimX; i++) {
                 index = (dimX*dimY)*k + j*dimX+i;
-
+                
                 if (i == 0) P_v1 = P1[index];
                 else if (i == dimX-1)  P_v1 = -P1[(dimX*dimY)*k + j*dimX+(i-1)];
                 else P_v1 = P1[index] - P1[(dimX*dimY)*k + j*dimX+(i-1)];
@@ -464,7 +464,7 @@ float DivProjP_3D(float *U, float *U0, float *P1, float *P2, float *P3, long dim
                 if (k == 0) P_v3 = P3[index];
                 else if (k == dimZ-1) P_v3 = -P3[(dimX*dimY)*(k-1) + (j)*dimX+i];
                 else P_v3 = P3[index] - P3[(dimX*dimY)*(k-1) + (j)*dimX+i];
-
+                
                 div = P_v1 + P_v2 + P_v3;
                 U[index] = (lambda*(U[index] + tau*div) + tau*U0[index])/(lambda + tau);
             }}}
@@ -476,14 +476,14 @@ float UpdV_3D(float *V1, float *V2, float *V3, float *P1, float *P2, float *P3,
     long i,j,k,index;
     float q1, q4x, q5x, q2, q4y, q6y, q6z, q5z, q3, div1, div2, div3;
 #pragma omp parallel for shared(V1,V2,V3,P1,P2,P3,Q1,Q2,Q3,Q4,Q5,Q6) private(i,j,k,index,q1,q4x,q5x,q2,q4y,q6y,q6z,q5z,q3,div1,div2,div3)
-    for(i=0; i<dimX; i++) {
+    for(k=0; k<dimZ; k++) {
         for(j=0; j<dimY; j++) {
-            for(k=0; k<dimZ; k++) {
+            for(i=0; i<dimX; i++) {
                 index = (dimX*dimY)*k + j*dimX+i;
                 q1 = 0.0f; q4x= 0.0f; q5x= 0.0f; q2= 0.0f; q4y= 0.0f; q6y= 0.0f; q6z= 0.0f; q5z= 0.0f; q3= 0.0f;
                 /* Q1 - Q11, Q2 - Q22, Q3 -  Q33, Q4 - Q21/Q12, Q5 - Q31/Q13, Q6 - Q32/Q23*/
                 /* symmetric boundary conditions (Neuman) */
-
+                
                 if (i == 0) {
                     q1 = Q1[index];
                     q4x = Q4[index];
@@ -520,11 +520,11 @@ float UpdV_3D(float *V1, float *V2, float *V3, float *P1, float *P2, float *P3,
                     q6z = Q6[index] - Q6[(dimX*dimY)*(k-1) + (j)*dimX+i];
                     q5z = Q5[index] - Q5[(dimX*dimY)*(k-1) + (j)*dimX+i];
                     q3 = Q3[index] - Q3[(dimX*dimY)*(k-1) + (j)*dimX+i]; }
-
+                
                 div1 = q1 + q4y + q5z;
                 div2 = q4x + q2 + q6z;
                 div3 = q5x + q6y + q3;
-
+                
                 V1[index] += tau*(P1[index] + div1);
                 V2[index] += tau*(P2[index] + div2);
                 V3[index] += tau*(P3[index] + div3);