16 files changed, 182 insertions, 258 deletions
diff --git a/demos/DemoRD2.m b/demos/DemoRD2.m
index 0db3595..293c1cd 100644
--- a/demos/DemoRD2.m
+++ b/demos/DemoRD2.m
@@ -30,7 +30,7 @@ Weights3D = single(data_raw3D); % weights for PW model
 clear data_raw3D
 %%
 % set projection/reconstruction geometry here
-Z_slices = 3;
+Z_slices = 20;
 det_row_count = Z_slices;
 proj_geom = astra_create_proj_geom('parallel3d', 1, 1, det_row_count, size_det, angles_rad);
 vol_geom = astra_create_vol_geom(recon_size,recon_size,Z_slices);
@@ -44,7 +44,7 @@ clear params
 params.proj_geom = proj_geom; % pass geometry to the function
 params.vol_geom = vol_geom;
 params.sino = Sino3D;
-params.iterFISTA  = 35;
+params.iterFISTA  = 1;
 params.weights = Weights3D;
 params.show = 1;
 params.maxvalplot = 2.5; params.slice = 2;
diff --git a/main_func/FISTA_REC.m b/main_func/FISTA_REC.m
index 2823691..1e93719 100644
--- a/main_func/FISTA_REC.m
+++ b/main_func/FISTA_REC.m
@@ -99,7 +99,7 @@ else
         clear proj_geomT vol_geomT
     else
         % divergen beam geometry
-        fprintf('%s \n', 'Calculating Lipshitz constant for divergen beam geometry...');
+        fprintf('%s \n', 'Calculating Lipshitz constant for divergen beam geometry... will take some time!');
         niter = 8; % number of iteration for PM
         x1 = rand(N,N,SlicesZ);
         sqweight = sqrt(weights);
@@ -216,6 +216,39 @@ if (isfield(params,'initialize'))
 else
     X = zeros(N,N,SlicesZ, 'single'); % storage for the solution
 end
+if (isfield(params,'OS'))
+    % Ordered Subsets reorganisation of data and angles
+    OS = 1;
+    subsets = 8;
+    angles = proj_geom.ProjectionAngles;
+    binEdges = linspace(min(angles),max(angles),subsets+1);
+    
+    % assign values to bins
+    [binsDiscr,~] = histc(angles, [binEdges(1:end-1) Inf]);
+    
+    % rearrange angles into subsets
+    AnglesReorg = zeros(length(angles),1);
+    SinoReorg = zeros(Detectors, anglesNumb, SlicesZ, 'single');
+    
+    counterM = 0;
+    for ii = 1:max(binsDiscr(:))
+        counter = 0;
+        for jj = 1:subsets
+            curr_index = ii+jj-1 + counter;
+            if (binsDiscr(jj) >= ii)
+                counterM = counterM + 1;
+                AnglesReorg(counterM) = angles(curr_index);
+                SinoReorg(:,counterM,:) = squeeze(sino(:,curr_index,:));
+            end
+            counter = (counter + binsDiscr(jj)) - 1;
+        end
+    end
+    sino = SinoReorg;
+    clear SinoReorg;
+else 
+    OS = 0; % normal FISTA
+end
+
 
 %----------------Reconstruction part------------------------
 Resid_error = zeros(iterFISTA,1); % errors vector (if the ground truth is given)
@@ -241,7 +274,7 @@ for i = 1:iterFISTA
         % ring removal part (Group-Huber fidelity)
         for kkk = 1:anglesNumb
             residual(:,kkk,:) =  squeeze(weights(:,kkk,:)).*(squeeze(sino_updt(:,kkk,:)) - (squeeze(sino(:,kkk,:)) - alpha_ring.*r_x));
-        end        
+        end
         vec = sum(residual,2);
         if (SlicesZ > 1)
             vec = squeeze(vec(:,1,:));
@@ -249,7 +282,7 @@ for i = 1:iterFISTA
         r = r_x - (1./L_const).*vec;
     else
         % no ring removal
-        residual = weights.*(sino_updt - sino);        
+        residual = weights.*(sino_updt - sino);
     end
     
     objective(i) = (0.5*norm(residual(:))^2)/(Detectors*anglesNumb*SlicesZ); % for the objective function output
@@ -305,5 +338,6 @@ end
 output.Resid_error = Resid_error;
 output.objective = objective;
 output.L_const = L_const;
+output.sino = sino;
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 end
diff --git a/main_func/compile_mex.m b/main_func/compile_mex.m
index 7bfa8eb..66c05da 100644
--- a/main_func/compile_mex.m
+++ b/main_func/compile_mex.m
@@ -1,10 +1,11 @@
-% compile mex's once
+% compile mex's in Matlab once
 cd regularizers_CPU/
 
-mex LLT_model.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
-mex FGP_TV.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
-mex SplitBregman_TV.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
-mex TGV_PD.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
+mex LLT_model.c LLT_model_core.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
+mex FGP_TV.c FGP_TV_core.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
+mex SplitBregman_TV.c SplitBregman_TV_core.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
+mex TGV_PD.c TGV_PD_core.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
+mex PatchBased_Regul.c PatchBased_Regul_core.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
 
 cd ../../
 cd demos
diff --git a/main_func/regularizers_CPU/FGP_TV_core.c b/main_func/regularizers_CPU/FGP_TV_core.c
index c383a71..a55991c 100644
--- a/main_func/regularizers_CPU/FGP_TV_core.c
+++ b/main_func/regularizers_CPU/FGP_TV_core.c
@@ -19,6 +19,32 @@ limitations under the License.
 
 #include "FGP_TV_core.h"
 
+/* C-OMP implementation of FGP-TV [1] denoising/regularization model (2D/3D case)
+ *
+ * Input Parameters:
+ * 1. Noisy image/volume [REQUIRED]
+ * 2. lambda - regularization parameter [REQUIRED]
+ * 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
+ * [2] last function value 
+ *
+ * Example of image denoising:
+ * figure;
+ * Im = double(imread('lena_gray_256.tif'))/255;  % loading image
+ * u0 = Im + .05*randn(size(Im)); % adding noise
+ * u = FGP_TV(single(u0), 0.05, 100, 1e-04);
+ *
+ * This function is based on the Matlab's code and paper by
+ * [1] Amir Beck and Marc Teboulle, "Fast Gradient-Based Algorithms for Constrained Total Variation Image Denoising and Deblurring Problems"
+ *
+ * D. Kazantsev, 2016-17
+ *
+ */
+
 /* 2D-case related Functions */
 /*****************************************************************/
 float Obj_func2D(float *A, float *D, float *R1, float *R2, float lambda, int dimX, int dimY)
diff --git a/main_func/regularizers_CPU/FGP_TV_core.h b/main_func/regularizers_CPU/FGP_TV_core.h
index 8d6af3e..e5328fb 100644
--- a/main_func/regularizers_CPU/FGP_TV_core.h
+++ b/main_func/regularizers_CPU/FGP_TV_core.h
@@ -24,33 +24,6 @@ limitations under the License.
 #include <stdio.h>
 #include "omp.h"
 
-/* C-OMP implementation of FGP-TV [1] denoising/regularization model (2D/3D case)
-*
-* Input Parameters:
-* 1. Noisy image/volume [REQUIRED]
-* 2. lambda - regularization parameter [REQUIRED]
-* 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
-* [2] last function value
-*
-* Example of image denoising:
-* figure;
-* Im = double(imread('lena_gray_256.tif'))/255;  % loading image
-* u0 = Im + .05*randn(size(Im)); % adding noise
-* u = FGP_TV(single(u0), 0.05, 100, 1e-04);
-*
-* to compile with OMP support: mex FGP_TV.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
-* This function is based on the Matlab's code and paper by
-* [1] Amir Beck and Marc Teboulle, "Fast Gradient-Based Algorithms for Constrained Total Variation Image Denoising and Deblurring Problems"
-*
-* D. Kazantsev, 2016-17
-*
-*/
-
 float copyIm(float *A, float *B, int dimX, int dimY, int dimZ);
 float Obj_func2D(float *A, float *D, float *R1, float *R2, float lambda, int dimX, int dimY);
 float Grad_func2D(float *P1, float *P2, float *D, float *R1, float *R2, float lambda, int dimX, int dimY);
diff --git a/main_func/regularizers_CPU/LLT_model.c b/main_func/regularizers_CPU/LLT_model.c
index 6b50d33..47146a5 100644
--- a/main_func/regularizers_CPU/LLT_model.c
+++ b/main_func/regularizers_CPU/LLT_model.c
@@ -20,6 +20,32 @@ limitations under the License.
 #include "mex.h"
 #include "LLT_model_core.h"
 
+/* C-OMP implementation of Lysaker, Lundervold and Tai (LLT) model of higher order regularization penalty
+*
+* Input Parameters:
+* 1. U0 - origanal noise image/volume
+* 2. lambda - regularization parameter
+* 3. tau - time-step  for explicit scheme
+* 4. iter - iterations number
+* 5. epsil  - tolerance constant (to terminate earlier)
+* 6. switcher - default is 0, switch to (1) to restrictive smoothing in Z dimension (in test)
+*
+* Output:
+* Filtered/regularized image
+*
+* Example:
+* figure;
+* Im = double(imread('lena_gray_256.tif'))/255;  % loading image
+* u0 = Im + .03*randn(size(Im)); % adding noise
+* [Den] = LLT_model(single(u0), 10, 0.1, 1);
+*
+*
+* to compile with OMP support: mex LLT_model.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
+* References: Lysaker, Lundervold and Tai (LLT) 2003, IEEE
+*
+* 28.11.16/Harwell
+*/
+
 
 void mexFunction(
         int nlhs, mxArray *plhs[],
diff --git a/main_func/regularizers_CPU/LLT_model_core.c b/main_func/regularizers_CPU/LLT_model_core.c
index 3098269..7a1cdbe 100644
--- a/main_func/regularizers_CPU/LLT_model_core.c
+++ b/main_func/regularizers_CPU/LLT_model_core.c
@@ -19,6 +19,31 @@ limitations under the License.
 
 #include "LLT_model_core.h"
 
+/* C-OMP implementation of Lysaker, Lundervold and Tai (LLT) model of higher order regularization penalty
+*
+* Input Parameters:
+* 1. U0 - origanal noise image/volume
+* 2. lambda - regularization parameter
+* 3. tau - time-step  for explicit scheme
+* 4. iter - iterations number
+* 5. epsil  - tolerance constant (to terminate earlier)
+* 6. switcher - default is 0, switch to (1) to restrictive smoothing in Z dimension (in test)
+*
+* Output:
+* Filtered/regularized image
+*
+* Example:
+* figure;
+* Im = double(imread('lena_gray_256.tif'))/255;  % loading image
+* u0 = Im + .03*randn(size(Im)); % adding noise
+* [Den] = LLT_model(single(u0), 10, 0.1, 1);
+*
+* References: Lysaker, Lundervold and Tai (LLT) 2003, IEEE
+*
+* 28.11.16/Harwell
+*/
+
+
 float der2D(float *U, float *D1, float *D2, int dimX, int dimY, int dimZ)
 {
 	int i, j, i_p, i_m, j_m, j_p;
diff --git a/main_func/regularizers_CPU/LLT_model_core.h b/main_func/regularizers_CPU/LLT_model_core.h
index ef8803d..10f52fe 100644
--- a/main_func/regularizers_CPU/LLT_model_core.h
+++ b/main_func/regularizers_CPU/LLT_model_core.h
@@ -26,31 +26,6 @@ limitations under the License.
 
 #define EPS 0.01
 
-/* C-OMP implementation of Lysaker, Lundervold and Tai (LLT) model of higher order regularization penalty
-*
-* Input Parameters:
-* 1. U0 - origanal noise image/volume
-* 2. lambda - regularization parameter
-* 3. tau - time-step  for explicit scheme
-* 4. iter - iterations number
-* 5. epsil  - tolerance constant (to terminate earlier)
-* 6. switcher - default is 0, switch to (1) to restrictive smoothing in Z dimension (in test)
-*
-* Output:
-* Filtered/regularized image
-*
-* Example:
-* figure;
-* Im = double(imread('lena_gray_256.tif'))/255;  % loading image
-* u0 = Im + .03*randn(size(Im)); % adding noise
-* [Den] = LLT_model(single(u0), 10, 0.1, 1);
-*
-*
-* to compile with OMP support: mex LLT_model.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
-* References: Lysaker, Lundervold and Tai (LLT) 2003, IEEE
-*
-* 28.11.16/Harwell
-*/
 /* 2D functions */
 float der2D(float *U, float *D1, float *D2, int dimX, int dimY, int dimZ);
 float div_upd2D(float *U0, float *U, float *D1, float *D2, int dimX, int dimY, int dimZ, float lambda, float tau);
diff --git a/main_func/regularizers_CPU/PatchBased_Regul.c b/main_func/regularizers_CPU/PatchBased_Regul.c
index 24ee210..59eb3b3 100644
--- a/main_func/regularizers_CPU/PatchBased_Regul.c
+++ b/main_func/regularizers_CPU/PatchBased_Regul.c
@@ -27,27 +27,23 @@ limitations under the License.
  * References: 1. Yang Z. & Jacob M. "Nonlocal Regularization of Inverse Problems"
  *             2. Kazantsev D. et al. "4D-CT reconstruction with unified spatial-temporal patch-based regularization"
  *
- * Input Parameters (mandatory):
- * 1. Image (2D or 3D)
- * 2. ratio of the searching window (e.g. 3 = (2*3+1) = 7 pixels window)
- * 3. ratio of the similarity window (e.g. 1 = (2*1+1) = 3 pixels window)
- * 4. h - parameter for the PB penalty function
- * 5. lambda - regularization parameter 
+ * Input Parameters:
+ * 1. Image (2D or 3D) [required]
+ * 2. ratio of the searching window (e.g. 3 = (2*3+1) = 7 pixels window) [optional]
+ * 3. ratio of the similarity window (e.g. 1 = (2*1+1) = 3 pixels window) [optional]
+ * 4. h - parameter for the PB penalty function [optional]
+ * 5. lambda - regularization parameter  [optional]
 
  * Output:
  * 1. regularized (denoised) Image (N x N)/volume (N x N x N)
  *
- * Quick 2D denoising example in Matlab:   
+ * 2D denoising example in Matlab:   
    Im = double(imread('lena_gray_256.tif'))/255;  % loading image
    u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise
-   ImDen = PB_Regul_CPU(single(u0), 3, 1, 0.08, 0.05); 
- *
- * Please see more tests in a file: 
-   TestTemporalSmoothing.m
- 
+   ImDen = PatchBased_Regul(single(u0), 3, 1, 0.08, 0.05); 
  *
  * Matlab + C/mex compilers needed
- * to compile with OMP support: mex PB_Regul_CPU.c CFLAGS="\$CFLAGS -fopenmp -Wall" LDFLAGS="\$LDFLAGS -fopenmp"
+ * to compile with OMP support: mex PatchBased_Regul.c CFLAGS="\$CFLAGS -fopenmp -Wall" LDFLAGS="\$LDFLAGS -fopenmp"
  *
  * D. Kazantsev *
  * 02/07/2014
@@ -70,17 +66,23 @@ void mexFunction(
     M = dims[1];
     Z = dims[2];
     
-    if ((numdims < 2) || (numdims > 3)) {mexErrMsgTxt("The input should be 2D image or 3D volume");}
+    if ((numdims < 2) || (numdims > 3)) {mexErrMsgTxt("The input is 2D image or 3D volume");}
     if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input in single precision is required"); }
     
     if(nrhs != 5) mexErrMsgTxt("Five inputs reqired: Image(2D,3D), SearchW, SimilW, Threshold, Regularization parameter");
     
     /*Handling inputs*/
-    A  = (float *) mxGetData(prhs[0]);    /* the image to regularize/filter */
-    SearchW_real  = (int) mxGetScalar(prhs[1]); /* the searching window ratio */
-    SimilW =  (int) mxGetScalar(prhs[2]);  /* the similarity window ratio */
-    h =  (float) mxGetScalar(prhs[3]);  /* parameter for the PB filtering function */
-    lambda = (float) mxGetScalar(prhs[4]); /* regularization parameter */   
+    A  = (float *) mxGetData(prhs[0]);    /* the image/volume to regularize/filter */
+    SearchW_real = 3; /*default value*/
+    SimilW = 1; /*default value*/
+    h = 0.1; 
+    lambda = 0.1;
+    
+    if ((nrhs == 2) || (nrhs == 3) || (nrhs == 4) || (nrhs == 5))   SearchW_real  = (int) mxGetScalar(prhs[1]); /* the searching window ratio */
+    if ((nrhs == 3) || (nrhs == 4) || (nrhs == 5))   SimilW =  (int) mxGetScalar(prhs[2]);  /* the similarity window ratio */
+    if ((nrhs == 4) || (nrhs == 5))  h =  (float) mxGetScalar(prhs[3]);  /* parameter for the PB filtering function */
+    if ((nrhs == 5))  lambda = (float) mxGetScalar(prhs[4]); /* regularization parameter */   
+
 
     if (h <= 0) mexErrMsgTxt("Parmeter for the PB penalty function should be > 0");
     if (lambda <= 0) mexErrMsgTxt(" Regularization parmeter should be > 0");
@@ -89,7 +91,6 @@ void mexFunction(
     
     /* SearchW_full = 2*SearchW + 1; */ /* the full searching window  size */
     /* SimilW_full = 2*SimilW + 1;  */  /* the full similarity window  size */
-
     
     padXY = SearchW + 2*SimilW; /* padding sizes */
     newsizeX = N + 2*(padXY); /* the X size of the padded array */
@@ -135,4 +136,4 @@ void mexFunction(
         switchpad_crop = 1; /*cropping*/
         pad_crop(Bp, B, M, N, Z, newsizeY, newsizeX, newsizeZ, padXY, switchpad_crop);
     } /*end else ndims*/ 
-}    
-\ No newline at end of file
+}    
diff --git a/main_func/regularizers_CPU/PatchBased_Regul_core.c b/main_func/regularizers_CPU/PatchBased_Regul_core.c
index 6f0a48d..acfb464 100644
--- a/main_func/regularizers_CPU/PatchBased_Regul_core.c
+++ b/main_func/regularizers_CPU/PatchBased_Regul_core.c
@@ -19,40 +19,33 @@ limitations under the License.
 
 #include "PatchBased_Regul_core.h"
 
-/* C-OMP implementation of  patch-based (PB) regularization (2D and 3D cases).
-* This method finds self-similar patches in data and performs one fixed point iteration to mimimize the PB penalty function
-*
-* References: 1. Yang Z. & Jacob M. "Nonlocal Regularization of Inverse Problems"
-*             2. Kazantsev D. et al. "4D-CT reconstruction with unified spatial-temporal patch-based regularization"
-*
-* Input Parameters (mandatory):
-* 1. Image (2D or 3D)
-* 2. ratio of the searching window (e.g. 3 = (2*3+1) = 7 pixels window)
-* 3. ratio of the similarity window (e.g. 1 = (2*1+1) = 3 pixels window)
-* 4. h - parameter for the PB penalty function
-* 5. lambda - regularization parameter
+/* C-OMP implementation of  patch-based (PB) regularization (2D and 3D cases). 
+ * This method finds self-similar patches in data and performs one fixed point iteration to mimimize the PB penalty function
+ * 
+ * References: 1. Yang Z. & Jacob M. "Nonlocal Regularization of Inverse Problems"
+ *             2. Kazantsev D. et al. "4D-CT reconstruction with unified spatial-temporal patch-based regularization"
+ *
+ * Input Parameters:
+ * 1. Image (2D or 3D) [required]
+ * 2. ratio of the searching window (e.g. 3 = (2*3+1) = 7 pixels window) [optional]
+ * 3. ratio of the similarity window (e.g. 1 = (2*1+1) = 3 pixels window) [optional]
+ * 4. h - parameter for the PB penalty function [optional]
+ * 5. lambda - regularization parameter  [optional]
 
-* Output:
-* 1. regularized (denoised) Image (N x N)/volume (N x N x N)
-*
-* Quick 2D denoising example in Matlab:
-Im = double(imread('lena_gray_256.tif'))/255;  % loading image
-u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise
-ImDen = PB_Regul_CPU(single(u0), 3, 1, 0.08, 0.05);
-*
-* Please see more tests in a file:
-TestTemporalSmoothing.m
-
-*
-* Matlab + C/mex compilers needed
-* to compile with OMP support: mex PB_Regul_CPU.c CFLAGS="\$CFLAGS -fopenmp -Wall" LDFLAGS="\$LDFLAGS -fopenmp"
-*
-* D. Kazantsev *
-* 02/07/2014
-* Harwell, UK
-*/
+ * Output:
+ * 1. regularized (denoised) Image (N x N)/volume (N x N x N)
+ *
+ * 2D denoising example in Matlab:   
+   Im = double(imread('lena_gray_256.tif'))/255;  % loading image
+   u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise
+   ImDen = PatchBased_Regul(single(u0), 3, 1, 0.08, 0.05); 
+ 
+ * D. Kazantsev *
+ * 02/07/2014
+ * Harwell, UK
+ */
 
-/*2D version*/
+/*2D version function */
 float PB_FUNC2D(float *A, float *B, int dimX, int dimY, int padXY, int SearchW, int SimilW, float h, float lambda)
 {
     int i, j, i_n, j_n, i_m, j_m, i_p, j_p, i_l, j_l, i1, j1, i2, j2, i3, j3, i5,j5, count, SimilW_full;
diff --git a/main_func/regularizers_CPU/PatchBased_Regul_core.h b/main_func/regularizers_CPU/PatchBased_Regul_core.h
index b83cf10..5aa6415 100644
--- a/main_func/regularizers_CPU/PatchBased_Regul_core.h
+++ b/main_func/regularizers_CPU/PatchBased_Regul_core.h
@@ -26,39 +26,6 @@ limitations under the License.
 #include <stdio.h>
 #include "omp.h"
 
-/* C-OMP implementation of  patch-based (PB) regularization (2D and 3D cases).
-* This method finds self-similar patches in data and performs one fixed point iteration to mimimize the PB penalty function
-*
-* References: 1. Yang Z. & Jacob M. "Nonlocal Regularization of Inverse Problems"
-*             2. Kazantsev D. et al. "4D-CT reconstruction with unified spatial-temporal patch-based regularization"
-*
-* Input Parameters (mandatory):
-* 1. Image (2D or 3D)
-* 2. ratio of the searching window (e.g. 3 = (2*3+1) = 7 pixels window)
-* 3. ratio of the similarity window (e.g. 1 = (2*1+1) = 3 pixels window)
-* 4. h - parameter for the PB penalty function
-* 5. lambda - regularization parameter
-
-* Output:
-* 1. regularized (denoised) Image (N x N)/volume (N x N x N)
-*
-* Quick 2D denoising example in Matlab:
-Im = double(imread('lena_gray_256.tif'))/255;  % loading image
-u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise
-ImDen = PB_Regul_CPU(single(u0), 3, 1, 0.08, 0.05);
-*
-* Please see more tests in a file:
-TestTemporalSmoothing.m
-
-*
-* Matlab + C/mex compilers needed
-* to compile with OMP support: mex PB_Regul_CPU.c CFLAGS="\$CFLAGS -fopenmp -Wall" LDFLAGS="\$LDFLAGS -fopenmp"
-*
-* D. Kazantsev *
-* 02/07/2014
-* Harwell, UK
-*/
-
 float pad_crop(float *A, float *Ap, int OldSizeX, int OldSizeY, int OldSizeZ, int NewSizeX, int NewSizeY, int NewSizeZ, int padXY, int switchpad_crop);
 float PB_FUNC2D(float *A, float *B, int dimX, int dimY, int padXY, int SearchW, int SimilW, float h, float lambda);
 float PB_FUNC3D(float *A, float *B, int dimX, int dimY, int dimZ, int padXY, int SearchW, int SimilW, float h, float lambda);
diff --git a/main_func/regularizers_CPU/SplitBregman_TV_core.c b/main_func/regularizers_CPU/SplitBregman_TV_core.c
index 26ad5b1..283dd43 100644
--- a/main_func/regularizers_CPU/SplitBregman_TV_core.c
+++ b/main_func/regularizers_CPU/SplitBregman_TV_core.c
@@ -37,7 +37,6 @@ limitations under the License.
 * u0 = Im + .05*randn(size(Im)); u0(u0 < 0) = 0;
 * u = SplitBregman_TV(single(u0), 10, 30, 1e-04);
 *
-* to compile with OMP support: mex SplitBregman_TV.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
 * References:
 * The Split Bregman Method for L1 Regularized Problems, by Tom Goldstein and Stanley Osher.
 * D. Kazantsev, 2016*
diff --git a/main_func/regularizers_CPU/SplitBregman_TV_core.h b/main_func/regularizers_CPU/SplitBregman_TV_core.h
index ed9112c..a7aaabb 100644
--- a/main_func/regularizers_CPU/SplitBregman_TV_core.h
+++ b/main_func/regularizers_CPU/SplitBregman_TV_core.h
@@ -23,30 +23,6 @@ limitations under the License.
 #include <stdio.h>
 #include "omp.h"
 
-/* C-OMP implementation of Split Bregman - TV denoising-regularization model (2D/3D)
-*
-* Input Parameters:
-* 1. Noisy image/volume
-* 2. lambda - regularization parameter
-* 3. Number of iterations [OPTIONAL parameter]
-* 4. eplsilon - tolerance constant [OPTIONAL parameter]
-* 5. TV-type: 'iso' or 'l1' [OPTIONAL parameter]
-*
-* Output:
-* Filtered/regularized image
-*
-* Example:
-* figure;
-* Im = double(imread('lena_gray_256.tif'))/255;  % loading image
-* u0 = Im + .05*randn(size(Im)); u0(u0 < 0) = 0;
-* u = SplitBregman_TV(single(u0), 10, 30, 1e-04);
-*
-* to compile with OMP support: mex SplitBregman_TV.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
-* References:
-* The Split Bregman Method for L1 Regularized Problems, by Tom Goldstein and Stanley Osher.
-* D. Kazantsev, 2016*
-*/
-
 float copyIm(float *A, float *B, int dimX, int dimY, int dimZ);
 float gauss_seidel2D(float *U, float *A, float *Dx, float *Dy, float *Bx, float *By, int dimX, int dimY, float lambda, float mu);
 float updDxDy_shrinkAniso2D(float *U, float *Dx, float *Dy, float *Bx, float *By, int dimX, int dimY, float lambda);
diff --git a/main_func/regularizers_CPU/TGV_PD.c b/main_func/regularizers_CPU/TGV_PD.c
index 8bce18a..d4bb787 100644
--- a/main_func/regularizers_CPU/TGV_PD.c
+++ b/main_func/regularizers_CPU/TGV_PD.c
@@ -39,7 +39,7 @@ limitations under the License.
  * u0 = Im + .03*randn(size(Im)); % adding noise
  * tic; u = PrimalDual_TGV(single(u0), 0.02, 1.3, 1, 550); toc;
  *
- * to compile with OMP support: mex TGV_PD.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
+ * to compile with OMP support: mex TGV_PD.c  TGV_PD_core.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
  * References:
  * K. Bredies "Total Generalized Variation"
  *
@@ -53,7 +53,7 @@ void mexFunction(
 {
     int number_of_dims, iter, dimX, dimY, dimZ, ll;
     const int  *dim_array;
-    float *A, *U, *U_old, *P1, *P2, *P3, *Q1, *Q2, *Q3, *Q4, *Q5, *Q6, *Q7, *Q8, *Q9, *V1, *V1_old, *V2, *V2_old, *V3, *V3_old, lambda, L2, tau, sigma,  alpha1, alpha0;
+    float *A, *U, *U_old, *P1, *P2, *Q1, *Q2, *Q3, *V1, *V1_old, *V2, *V2_old, lambda, L2, tau, sigma,  alpha1, alpha0;
     
     number_of_dims = mxGetNumberOfDimensions(prhs[0]);
     dim_array = mxGetDimensions(prhs[0]);
@@ -88,40 +88,10 @@ void mexFunction(
         V1 = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
         V1_old = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
         V2 = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
-        V2_old = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));   }
-    else if (number_of_dims == 3) {
-        mexErrMsgTxt("The input data should be 2D");
-        /*3D case*/
-//         dimZ = dim_array[2];
-//         U = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
-//         
-//         P1 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
-//         P2 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
-//         P3 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
-//         
-//         Q1 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
-//         Q2 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
-//         Q3 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
-//         Q4 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
-//         Q5 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
-//         Q6 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
-//         Q7 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
-//         Q8 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
-//         Q9 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
-//         
-//         U_old = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
-//         
-//         V1 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
-//         V1_old = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
-//         V2 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
-//         V2_old = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
-//         V3 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
-//         V3_old = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));	  
-    }
-    else {mexErrMsgTxt("The input data should be 2D");}
-    
-    
-    /*printf("%i \n", i);*/
+        V2_old = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));   
+        
+        
+         /*printf("%i \n", i);*/
     L2 = 12.0; /*Lipshitz constant*/
     tau = 1.0/pow(L2,0.5);
     sigma = 1.0/pow(L2,0.5);
@@ -129,7 +99,6 @@ void mexFunction(
     /*Copy A to U*/
     copyIm(A, U, dimX, dimY, dimZ);
     
-    if (number_of_dims == 2) {
         /* Here primal-dual iterations begin for 2D */
         for(ll = 0; ll < iter; ll++) {
             
@@ -164,23 +133,12 @@ void mexFunction(
             /*get new V*/
             newU(V1, V1_old, dimX, dimY, dimZ);
             newU(V2, V2_old, dimX, dimY, dimZ);
-        } /*end of iterations*/
+        } /*end of iterations*/        
     }
-    
-//     /*3D version*/
-//     if (number_of_dims == 3) {
-//         /* Here primal-dual iterations begin for 3D */
-//         for(ll = 0; ll < iter; ll++) {
-//             
-//             /* Calculate Dual Variable P */
-//             DualP_3D(U, V1, V2, V3, P1, P2, P3, dimX, dimY, dimZ, sigma);
-//             
-//             /*Projection onto convex set for P*/
-//             ProjP_3D(P1, P2, P3, dimX, dimY, dimZ, alpha1);
-//             
-//             /* Calculate Dual Variable Q */
-//             DualQ_3D(V1, V2, V2, Q1, Q2, Q3, Q4, Q5, Q6, Q7, Q8, Q9, dimX, dimY, dimZ, sigma);
-//             
-//         } /*end of iterations*/
-//     }
+    else if (number_of_dims == 3) {
+        mexErrMsgTxt("The input data should be a 2D array");
+        /*3D case*/
+    }
+    else {mexErrMsgTxt("The input data should be a 2D array");}    
+   
 }
diff --git a/main_func/regularizers_CPU/TGV_PD_core.c b/main_func/regularizers_CPU/TGV_PD_core.c
index 4c0427c..1164b73 100644
--- a/main_func/regularizers_CPU/TGV_PD_core.c
+++ b/main_func/regularizers_CPU/TGV_PD_core.c
@@ -38,7 +38,6 @@ limitations under the License.
  * u0 = Im + .03*randn(size(Im)); % adding noise
  * tic; u = PrimalDual_TGV(single(u0), 0.02, 1.3, 1, 550); toc;
  *
- * to compile with OMP support: mex TGV_PD.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
  * References:
  * K. Bredies "Total Generalized Variation"
  *
diff --git a/main_func/regularizers_CPU/TGV_PD_core.h b/main_func/regularizers_CPU/TGV_PD_core.h
index e25ae68..04ba95c 100644
--- a/main_func/regularizers_CPU/TGV_PD_core.h
+++ b/main_func/regularizers_CPU/TGV_PD_core.h
@@ -24,32 +24,6 @@ limitations under the License.
 #include <stdio.h>
 #include "omp.h"
 
-/* C-OMP implementation of Primal-Dual denoising method for
-* Total Generilized Variation (TGV)-L2 model (2D case only)
-*
-* Input Parameters:
-* 1. Noisy image/volume (2D)
-* 2. lambda - regularization parameter
-* 3. parameter to control first-order term (alpha1)
-* 4. parameter to control the second-order term (alpha0)
-* 5. Number of CP iterations
-*
-* Output:
-* Filtered/regularized image
-*
-* Example:
-* figure;
-* Im = double(imread('lena_gray_256.tif'))/255;  % loading image
-* u0 = Im + .03*randn(size(Im)); % adding noise
-* tic; u = PrimalDual_TGV(single(u0), 0.02, 1.3, 1, 550); toc;
-*
-* to compile with OMP support: mex TGV_PD.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
-* References:
-* K. Bredies "Total Generalized Variation"
-*
-* 28.11.16/Harwell
-*/
-
 /* 2D functions */
 float DualP_2D(float *U, float *V1, float *V2, float *P1, float *P2, int dimX, int dimY, int dimZ, float sigma);
 float ProjP_2D(float *P1, float *P2, int dimX, int dimY, int dimZ, float alpha1);
@@ -57,8 +31,5 @@ float DualQ_2D(float *V1, float *V2, float *Q1, float *Q2, float *Q3, int dimX,
 float ProjQ_2D(float *Q1, float *Q2, float *Q3, int dimX, int dimY, int dimZ, float alpha0);
 float DivProjP_2D(float *U, float *A, float *P1, float *P2, int dimX, int dimY, int dimZ, float lambda, float tau);
 float UpdV_2D(float *V1, float *V2, float *P1, float *P2, float *Q1, float *Q2, float *Q3, int dimX, int dimY, int dimZ, float tau);
-/*3D functions*/
-float DualP_3D(float *U, float *V1, float *V2, float *V3, float *P1, float *P2, float *P3, int dimX, int dimY, int dimZ, float sigma);
-
 float newU(float *U, float *U_old, int dimX, int dimY, int dimZ);
 float copyIm(float *A, float *U, int dimX, int dimY, int dimZ);