ROF TV regularizer added #22

4 files changed, 421 insertions, 1 deletions

diff --git a/Core/regularizers_CPU/ROF_TV_core.c b/Core/regularizers_CPU/ROF_TV_core.c
new file mode 100644
index 0000000..0b24806
--- /dev/null
+++ b/Core/regularizers_CPU/ROF_TV_core.c
@@ -0,0 +1,259 @@
+/*
+ * 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 "ROF_TV_core.h"
+
+#define EPS 0.000001
+#define MAX(x, y) (((x) > (y)) ? (x) : (y))
+#define MIN(x, y) (((x) < (y)) ? (x) : (y))
+
+int sign(float x) {
+    return (x > 0) - (x < 0);
+}
+
+/* C-OMP implementation of ROF-TV denoising/regularization model [1] (2D/3D case)
+ *
+ * Input Parameters:
+ * 1. Noisy image/volume [REQUIRED]
+ * 2. lambda - regularization parameter [REQUIRED]
+ * 3. tau - marching step for explicit scheme, ~0.001 is recommended [REQUIRED]
+ * 4. Number of iterations, for explicit scheme >= 150 is recommended  [REQUIRED]
+ *
+ * Output:
+ * [1] Regularized image/volume
+ *
+ * This function is based on the paper by
+ * [1] Rudin, Osher, Fatemi, "Nonlinear Total Variation based noise removal algorithms"
+ *
+ * D. Kazantsev, 2016-18
+ */
+
+/* calculate differences 1 */
+float D1_func(float *A, float *D1, int dimY, int dimX, int dimZ)
+{
+    float NOMx_1, NOMy_1, NOMy_0, NOMz_1, NOMz_0, denom1, denom2,denom3, T1;
+    int i,j,k,i1,i2,k1,j1,j2,k2;
+    
+    if (dimZ == 0) {
+#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)
+        for(j=0; j<dimY; j++) {
+            for(i=0; i<dimX; i++) {
+                /* symmetric boundary conditions (Neuman) */
+                i1 = i + 1; if (i1 >= dimY) i1 = i-1;
+                i2 = i - 1; if (i2 < 0) i2 = i+1;
+                j1 = j + 1; if (j1 >= dimX) j1 = j-1;
+                j2 = j - 1; if (j2 < 0) j2 = j+1;
+                
+                /* Forward-backward differences */
+                NOMx_1 = A[i1*dimY + j] - A[(i)*dimY + j]; /* x+ */
+                NOMy_1 = A[i*dimY + j1] - A[(i)*dimY + j]; /* y+ */
+                /*NOMx_0 = (A[(i)*dimY + j] - A[(i2)*dimY + j]); */ /* x- */
+                NOMy_0 = A[(i)*dimY + j] - A[(i)*dimY + j2]; /* y- */
+                
+                denom1 = NOMx_1*NOMx_1;
+                denom2 = 0.5*(sign(NOMy_1) + sign(NOMy_0))*(MIN(fabs(NOMy_1),fabs(NOMy_0)));
+                denom2 = denom2*denom2;
+                T1 = sqrt(denom1 + denom2 + EPS);
+                D1[i*dimY+j] = NOMx_1/T1;
+            }}
+    }
+    else {
+#pragma omp parallel for shared (A, D1, dimX, dimY, dimZ) private(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++) {
+                    /* symmetric boundary conditions (Neuman) */
+                    i1 = i + 1; if (i1 >= dimY) i1 = i-1;
+                    i2 = i - 1; if (i2 < 0) i2 = i+1;
+                    j1 = j + 1; if (j1 >= dimX) 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;
+                    /*B[(dimX*dimY)*k + i*dimY+j] = 0.25*(A[(dimX*dimY)*k + (i1)*dimY + j] + A[(dimX*dimY)*k + (i2)*dimY + j] + A[(dimX*dimY)*k + (i)*dimY + j1] + A[(dimX*dimY)*k + (i)*dimY + j2]) -  A[(dimX*dimY)*k + i*dimY + j];*/
+                    
+                    /* Forward-backward differences */
+                    NOMx_1 = A[(dimX*dimY)*k + (i1)*dimY + j] - A[(dimX*dimY)*k + (i)*dimY + j]; /* x+ */
+                    NOMy_1 = A[(dimX*dimY)*k + (i)*dimY + j1] - A[(dimX*dimY)*k + (i)*dimY + j]; /* y+ */
+                    /*NOMx_0 = (A[(i)*dimY + j] - A[(i2)*dimY + j]); */  /* x- */
+                    NOMy_0 = A[(dimX*dimY)*k + (i)*dimY + j] - A[(dimX*dimY)*k + (i)*dimY + j2]; /* y- */
+                    
+                    NOMz_1 = A[(dimX*dimY)*k1 + (i)*dimY + j] - A[(dimX*dimY)*k + (i)*dimY + j]; /* z+ */
+                    NOMz_0 = A[(dimX*dimY)*k + (i)*dimY + j] - A[(dimX*dimY)*k2 + (i)*dimY + j]; /* z- */
+                    
+                    
+                    denom1 = NOMx_1*NOMx_1;
+                    denom2 = 0.5*(sign(NOMy_1) + sign(NOMy_0))*(MIN(fabs(NOMy_1),fabs(NOMy_0)));
+                    denom2 = denom2*denom2;
+                    denom3 = 0.5*(sign(NOMz_1) + sign(NOMz_0))*(MIN(fabs(NOMz_1),fabs(NOMz_0)));
+                    denom3 = denom3*denom3;
+                    T1 = sqrt(denom1 + denom2 + denom3 + EPS);
+                    D1[(dimX*dimY)*k + i*dimY+j] = NOMx_1/T1;
+                }}}
+    }
+    return *D1;
+}
+/* calculate differences 2 */
+float D2_func(float *A, float *D2, int dimY, int dimX, int dimZ)
+{
+    float NOMx_1, NOMy_1, NOMx_0, NOMz_1, NOMz_0, denom1, denom2, denom3, T2;
+    int i,j,k,i1,i2,k1,j1,j2,k2;
+    
+    if (dimZ == 0) {
+#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)
+        for(j=0; j<dimY; j++) {
+            for(i=0; i<dimX; i++) {
+                /* symmetric boundary conditions (Neuman) */
+                i1 = i + 1; if (i1 >= dimY) i1 = i-1;
+                i2 = i - 1; if (i2 < 0) i2 = i+1;
+                j1 = j + 1; if (j1 >= dimX) j1 = j-1;
+                j2 = j - 1; if (j2 < 0) j2 = j+1;
+                
+                /* Forward-backward differences */
+                NOMx_1 = A[(i1)*dimY + j] - A[(i)*dimY + j]; /* x+ */
+                NOMy_1 = A[i*dimY + j1] - A[(i)*dimY + j]; /* y+ */
+                NOMx_0 = A[(i)*dimY + j] - A[(i2)*dimY + j]; /* x- */
+                /*NOMy_0 = A[(i)*dimY + j] - A[(i)*dimY + j2]; */  /* y- */
+                
+                denom1 = NOMy_1*NOMy_1;
+                denom2 = 0.5*(sign(NOMx_1) + sign(NOMx_0))*(MIN(fabs(NOMx_1),fabs(NOMx_0)));
+                denom2 = denom2*denom2;
+                T2 = sqrt(denom1 + denom2 + EPS);
+                D2[i*dimY+j] = NOMy_1/T2;
+            }}
+    }
+    else {
+#pragma omp parallel for shared (A, D2, dimX, dimY, dimZ) private(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++) {
+                    /* symmetric boundary conditions (Neuman) */
+                    i1 = i + 1; if (i1 >= dimY) i1 = i-1;
+                    i2 = i - 1; if (i2 < 0) i2 = i+1;
+                    j1 = j + 1; if (j1 >= dimX) 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;
+                    
+                    
+                    /* Forward-backward differences */
+                    NOMx_1 = A[(dimX*dimY)*k + (i1)*dimY + j] - A[(dimX*dimY)*k + (i)*dimY + j]; /* x+ */
+                    NOMy_1 = A[(dimX*dimY)*k + (i)*dimY + j1] - A[(dimX*dimY)*k + (i)*dimY + j]; /* y+ */
+                    NOMx_0 = A[(dimX*dimY)*k + (i)*dimY + j] - A[(dimX*dimY)*k + (i2)*dimY + j]; /* x- */
+                    NOMz_1 = A[(dimX*dimY)*k1 + (i)*dimY + j] - A[(dimX*dimY)*k + (i)*dimY + j]; /* z+ */
+                    NOMz_0 = A[(dimX*dimY)*k + (i)*dimY + j] - A[(dimX*dimY)*k2 + (i)*dimY + j]; /* z- */
+                    
+                    
+                    denom1 = NOMy_1*NOMy_1;
+                    denom2 = 0.5*(sign(NOMx_1) + sign(NOMx_0))*(MIN(fabs(NOMx_1),fabs(NOMx_0)));
+                    denom2 = denom2*denom2;
+                    denom3 = 0.5*(sign(NOMz_1) + sign(NOMz_0))*(MIN(fabs(NOMz_1),fabs(NOMz_0)));
+                    denom3 = denom3*denom3;
+                    T2 = sqrt(denom1 + denom2 + denom3 + EPS);
+                    D2[(dimX*dimY)*k + i*dimY+j] = NOMy_1/T2;
+                }}}
+    }
+    return *D2;
+}
+
+/* calculate differences 3 */
+float D3_func(float *A, float *D3, int dimY, int dimX, int dimZ)
+{
+    float NOMx_1, NOMy_1, NOMx_0, NOMy_0, NOMz_1, denom1, denom2, denom3, T3;
+    int i,j,k,i1,i2,k1,j1,j2,k2;
+    
+#pragma omp parallel for shared (A, D3, dimX, dimY, dimZ) private(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++) {
+                /* symmetric boundary conditions (Neuman) */
+                i1 = i + 1; if (i1 >= dimY) i1 = i-1;
+                i2 = i - 1; if (i2 < 0) i2 = i+1;
+                j1 = j + 1; if (j1 >= dimX) 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;
+                
+                /* Forward-backward differences */
+                NOMx_1 = A[(dimX*dimY)*k + (i1)*dimY + j] - A[(dimX*dimY)*k + (i)*dimY + j]; /* x+ */
+                NOMy_1 = A[(dimX*dimY)*k + (i)*dimY + j1] - A[(dimX*dimY)*k + (i)*dimY + j]; /* y+ */
+                NOMy_0 = A[(dimX*dimY)*k + (i)*dimY + j] - A[(dimX*dimY)*k + (i)*dimY + j2]; /* y- */
+                NOMx_0 = A[(dimX*dimY)*k + (i)*dimY + j] - A[(dimX*dimY)*k + (i2)*dimY + j]; /* x- */
+                NOMz_1 = A[(dimX*dimY)*k1 + (i)*dimY + j] - A[(dimX*dimY)*k + (i)*dimY + j]; /* 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.5*(sign(NOMx_1) + sign(NOMx_0))*(MIN(fabs(NOMx_1),fabs(NOMx_0)));
+                denom2 = denom2*denom2;
+                denom3 = 0.5*(sign(NOMy_1) + sign(NOMy_0))*(MIN(fabs(NOMy_1),fabs(NOMy_0)));
+                denom3 = denom3*denom3;
+                T3 = sqrt(denom1 + denom2 + denom3 + EPS);
+                D3[(dimX*dimY)*k + i*dimY+j] = NOMz_1/T3;
+            }}}
+    return *D3;
+}
+
+/* calculate divergence */
+float TV_main(float *D1, float *D2, float *D3, float *B, float *A, float lambda, float tau, int dimY, int dimX, int dimZ)
+{
+    float dv1, dv2, dv3;
+    int index,i,j,k,i1,i2,k1,j1,j2,k2;
+    
+    if (dimZ == 0) {
+#pragma omp parallel for shared (D1, D2, B, dimX, dimY) private(index, i, j, i1, j1, i2, j2,dv1,dv2)
+        for(j=0; j<dimY; j++) {
+            for(i=0; i<dimX; i++) {
+                index = (i)*dimY + j;
+                /* symmetric boundary conditions (Neuman) */
+                i1 = i + 1; if (i1 >= dimY) i1 = i-1;
+                i2 = i - 1; if (i2 < 0) i2 = i+1;
+                j1 = j + 1; if (j1 >= dimX) j1 = j-1;
+                j2 = j - 1; if (j2 < 0) j2 = j+1;
+                
+                /* divergence components  */
+                dv1 = D1[index] - D1[(i2)*dimY + j];
+                dv2 = D2[index] - D2[(i)*dimY + j2];
+                
+                B[index] =  B[index] + tau*lambda*(dv1 + dv2) + tau*(A[index] - B[index]);
+                
+            }}
+    }
+    else {
+#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)
+        for(j=0; j<dimY; j++) {
+            for(i=0; i<dimX; i++) {
+                for(k=0; k<dimZ; k++) {
+                    index = (dimX*dimY)*k + i*dimY+j;
+                    /* symmetric boundary conditions (Neuman) */
+                    i1 = i + 1; if (i1 >= dimY) i1 = i-1;
+                    i2 = i - 1; if (i2 < 0) i2 = i+1;
+                    j1 = j + 1; if (j1 >= dimX) 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;
+                    
+                    /*divergence components */
+                    dv1 = D1[index] - D1[(dimX*dimY)*k + i2*dimY+j];
+                    dv2 = D2[index] - D2[(dimX*dimY)*k + i*dimY+j2];
+                    dv3 = D3[index] - D3[(dimX*dimY)*k2 + i*dimY+j];
+                    
+                    B[index] = B[index] + tau*lambda*(dv1 + dv2 + dv3) + tau*(A[index] - B[index]);
+                }}}
+    }
+    return *B;
+}
\ No newline at end of file
diff --git a/Core/regularizers_CPU/ROF_TV_core.h b/Core/regularizers_CPU/ROF_TV_core.h
new file mode 100644
index 0000000..6e4f961
--- /dev/null
+++ b/Core/regularizers_CPU/ROF_TV_core.h
@@ -0,0 +1,54 @@
+/*
+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 <math.h>
+#include <stdlib.h>
+#include <memory.h>
+#include <stdio.h>
+#include "omp.h"
+#include "utils.h"
+#include "CCPiDefines.h"
+
+/* C-OMP implementation of ROF-TV denoising/regularization model [1] (2D/3D case)
+*
+* Input Parameters:
+ * 1. Noisy image/volume [REQUIRED]
+ * 2. lambda - regularization parameter [REQUIRED]
+ * 3. tau - marching step for explicit scheme, ~0.001 is recommended [REQUIRED]
+ * 4. Number of iterations, for explicit scheme >= 150 is recommended [REQUIRED]
+*
+* Output:
+* [1] Regularized image/volume
+
+ * This function is based on the paper by
+* [1] Rudin, Osher, Fatemi, "Nonlinear Total Variation based noise removal algorithms"
+*
+* D. Kazantsev, 2016-18
+*/
+#ifdef __cplusplus
+extern "C" {
+#endif
+//float copyIm(float *A, float *B, int dimX, int dimY, int dimZ);
+CCPI_EXPORT float TV_main(float *D1, float *D2, float *D3, float *B, float *A, float lambda, float tau, int dimY, int dimX, int dimZ);
+CCPI_EXPORT float D1_func(float *A, float *D1, int dimY, int dimX, int dimZ);
+CCPI_EXPORT float D2_func(float *A, float *D2, int dimY, int dimX, int dimZ);
+CCPI_EXPORT float D3_func(float *A, float *D3, int dimY, int dimX, int dimZ);
+#ifdef __cplusplus
+}
+#endif
\ No newline at end of file
diff --git a/Wrappers/Matlab/mex_compile/compile_mex.m b/Wrappers/Matlab/mex_compile/compile_mex.m
index e1debf3..ee85b49 100644
--- a/Wrappers/Matlab/mex_compile/compile_mex.m
+++ b/Wrappers/Matlab/mex_compile/compile_mex.m
@@ -7,13 +7,14 @@ cd regularizers_CPU/
 
 % compile C regularizers
 
+mex ROF_TV.c ROF_TV_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
 mex LLT_model.c LLT_model_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
 mex FGP_TV.c FGP_TV_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
 mex SplitBregman_TV.c SplitBregman_TV_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
 mex TGV_PD.c TGV_PD_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
 mex PatchBased_Regul.c PatchBased_Regul_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
 
-delete LLT_model_core.c LLT_model_core.h FGP_TV_core.c FGP_TV_core.h SplitBregman_TV_core.c SplitBregman_TV_core.h  TGV_PD_core.c  TGV_PD_core.h PatchBased_Regul_core.c PatchBased_Regul_core.h utils.c utils.h CCPiDefines.h
+delete ROF_TV_core.c ROF_TV_core.h LLT_model_core.c LLT_model_core.h FGP_TV_core.c FGP_TV_core.h SplitBregman_TV_core.c SplitBregman_TV_core.h  TGV_PD_core.c  TGV_PD_core.h PatchBased_Regul_core.c PatchBased_Regul_core.h utils.c utils.h CCPiDefines.h
 
 % compile CUDA-based regularizers
 %cd regularizers_GPU/
diff --git a/Wrappers/Matlab/mex_compile/regularizers_CPU/ROF_TV.c b/Wrappers/Matlab/mex_compile/regularizers_CPU/ROF_TV.c
new file mode 100644
index 0000000..a800add
--- /dev/null
+++ b/Wrappers/Matlab/mex_compile/regularizers_CPU/ROF_TV.c
@@ -0,0 +1,106 @@
+/*
+ * 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 "matrix.h"
+#include "mex.h"
+#include "ROF_TV_core.h"
+
+/* C-OMP implementation of ROF-TV denoising/regularization model [1] (2D/3D case)
+ *
+ * Input Parameters:
+ * 1. Noisy image/volume [REQUIRED]
+ * 2. lambda - regularization parameter [REQUIRED]
+ * 3. tau - marching step for explicit scheme, ~0.001 is recommended [REQUIRED]
+ * 4. Number of iterations, for explicit scheme >= 150 is recommended  [REQUIRED]
+ *
+ * Output:
+ * [1] Regularized image/volume
+ *
+ * This function is based on the paper by
+ * [1] Rudin, Osher, Fatemi, "Nonlinear Total Variation based noise removal algorithms"
+ * compile: mex ROF_TV.c ROF_TV_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
+ * D. Kazantsev, 2016-18
+ */
+
+void mexFunction(
+        int nlhs, mxArray *plhs[],
+        int nrhs, const mxArray *prhs[])
+        
+{
+    int i, number_of_dims, iter_numb, dimX, dimY, dimZ;
+    const int  *dim_array;
+    float *A, *B, *D1, *D2, *D3, lambda, tau;
+    
+    dim_array = mxGetDimensions(prhs[0]);
+    number_of_dims = mxGetNumberOfDimensions(prhs[0]);
+    
+    /*Handling Matlab input data*/
+    A  = (float *) mxGetData(prhs[0]);
+    lambda =  (float) mxGetScalar(prhs[1]); /* regularization parameter */
+    tau =  (float) mxGetScalar(prhs[2]); /* marching step parameter */
+    iter_numb =  (int) mxGetScalar(prhs[3]); /* iterations number */
+    
+    if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); }
+    /*Handling Matlab output data*/
+    dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2];
+    
+    /* output arrays*/
+    if (number_of_dims == 2) {
+        dimZ = 0; /*2D case*/
+        B = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
+        D1 = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
+        D2 = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
+        
+        /* copy into B */
+        copyIm(A, B, dimX, dimY, 1);
+        
+        /* start TV iterations */
+        for(i=0; i < iter_numb; i++) {
+            
+            /* calculate differences */
+            D1_func(B, D1, dimX, dimY, dimZ);
+            D2_func(B, D2, dimX, dimY, dimZ);
+            
+            /* calculate divergence and image update*/
+            TV_main(D1, D2, D2, B, A, lambda, tau, dimX, dimY, dimZ);
+        }
+    }
+    
+    if (number_of_dims == 3) {
+        B = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
+        D1 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
+        D2 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
+        D3 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
+        
+        /* copy into B */
+        copyIm(A, B, dimX, dimY, dimZ);
+        
+        /* start TV iterations */
+        for(i=0; i < iter_numb; i++) {
+            
+        /* calculate differences */
+        D1_func(B, D1, dimX, dimY, dimZ);
+        D2_func(B, D2, dimX, dimY, dimZ);
+        D3_func(B, D3, dimX, dimY, dimZ);
+        
+        /* calculate divergence and image update*/
+        TV_main(D1, D2, D3, B, A, lambda, tau, dimX, dimY, dimZ);
+        }
+    }
+    
+}