Merge pull request #53 from vais-ral/inpaint_Heat

Inpainting methods added to the toolkit

24 files changed, 1271 insertions, 95 deletions

diff --git a/Core/CMakeLists.txt b/Core/CMakeLists.txt
index 61986dc..86d69b2 100644
--- a/Core/CMakeLists.txt
+++ b/Core/CMakeLists.txt
@@ -24,15 +24,6 @@ if (PYTHONINTERP_FOUND)
   endif()
 endif()
 
-#find_package(Boost  REQUIRED
-#	COMPONENTS ${BOOST_PYTHON} ${BOOST_NUMPY})
-#
-#if (Boost_FOUND)
-#   message("Boost version " ${Boost_VERSION})
-#   message("Boost include dir " ${Boost_INCLUDE_DIRS})
-#   message("Boost library dir " ${Boost_LIBRARY_DIRS})
-#   message("Boost libraries " ${Boost_LIBRARIES})
-#endif()
 
 find_package(OpenMP)
 if (OPENMP_FOUND)
@@ -93,12 +84,15 @@ add_library(cilreg SHARED
         ${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/FGP_dTV_core.c
         ${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/TNV_core.c
 	    ${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/utils.c
+	    ${CMAKE_CURRENT_SOURCE_DIR}/inpainters_CPU/Diffusion_Inpaint_core.c
+	    ${CMAKE_CURRENT_SOURCE_DIR}/inpainters_CPU/NonlocalMarching_Inpaint_core.c
 	    )
 target_link_libraries(cilreg ${EXTRA_LIBRARIES} )
 include_directories(cilreg PUBLIC 
                       ${LIBRARY_INC}/include 
 					  ${CMAKE_CURRENT_SOURCE_DIR}
-		              ${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/  )
+		              ${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/
+		              ${CMAKE_CURRENT_SOURCE_DIR}/inpainters_CPU/  )
 #GENERATE_EXPORT_HEADER( cilreg
 #             BASE_NAME cilreg
 #             EXPORT_MACRO_NAME CCPiCore_EXPORTS
@@ -153,8 +147,3 @@ if (CUDA_FOUND)
 else()
   message("CUDA NOT FOUND")
 endif()
-
-
-#add_executable(regulariser_test ${CMAKE_CURRENT_SOURCE_DIR}/test/test_regulariser.cpp)
-
-#target_link_libraries (regulariser_test LINK_PUBLIC regularisers_lib)
diff --git a/Core/inpainters_CPU/Diffusion_Inpaint_core.c b/Core/inpainters_CPU/Diffusion_Inpaint_core.c
new file mode 100644
index 0000000..56e0421
--- /dev/null
+++ b/Core/inpainters_CPU/Diffusion_Inpaint_core.c
@@ -0,0 +1,322 @@
+/*
+ * 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 "Diffusion_Inpaint_core.h"
+#include "utils.h"
+
+/*sign function*/
+int signNDF_inc(float x) {
+    return (x > 0) - (x < 0);
+}
+
+/* C-OMP implementation of linear and nonlinear diffusion [1,2] for inpainting task (2D/3D case)
+ * The minimisation is performed using explicit scheme. 
+ *
+ * Input Parameters:
+ * 1. Image/volume to inpaint
+ * 2. Mask of the same size as (1) in 'unsigned char' format  (ones mark the region to inpaint, zeros belong to the data)
+ * 3. lambda - regularization parameter
+ * 4. Edge-preserving parameter (sigma), when sigma equals to zero nonlinear diffusion -> linear diffusion
+ * 5. Number of iterations, for explicit scheme >= 150 is recommended 
+ * 6. tau - time-marching step for explicit scheme
+ * 7. Penalty type: 1 - Huber, 2 - Perona-Malik, 3 - Tukey Biweight
+ *
+ * Output:
+ * [1] Inpainted image/volume 
+ *
+ * This function is based on the paper by
+ * [1] Perona, P. and Malik, J., 1990. Scale-space and edge detection using anisotropic diffusion. IEEE Transactions on pattern analysis and machine intelligence, 12(7), pp.629-639.
+ * [2] Black, M.J., Sapiro, G., Marimont, D.H. and Heeger, D., 1998. Robust anisotropic diffusion. IEEE Transactions on image processing, 7(3), pp.421-432.
+ */
+
+float Diffusion_Inpaint_CPU_main(float *Input, unsigned char *Mask, float *Output, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int penaltytype, int dimX, int dimY, int dimZ)
+{
+    int i,pointsone;
+    float sigmaPar2;
+    sigmaPar2 = sigmaPar/sqrt(2.0f);
+    
+    /* copy into output */
+    copyIm(Input, Output, dimX, dimY, dimZ);
+    
+    pointsone = 0;
+    for (i=0; i<dimY*dimX*dimZ; i++) if (Mask[i] == 1) pointsone++;
+        
+    if (pointsone == 0) printf("%s \n", "Nothing to inpaint, zero mask!");
+    else {
+    
+    if (dimZ == 1) {
+    /* running 2D diffusion iterations */
+    for(i=0; i < iterationsNumb; i++) {
+            if (sigmaPar == 0.0f) LinearDiff_Inp_2D(Input, Mask, Output, lambdaPar, tau, dimX, dimY); /* linear diffusion (heat equation) */
+            else NonLinearDiff_Inp_2D(Input, Mask, Output, lambdaPar, sigmaPar2, tau, penaltytype, dimX, dimY); /* nonlinear diffusion */
+		}
+	}
+	else {
+	/* running 3D diffusion iterations */
+    for(i=0; i < iterationsNumb; i++) {
+            if (sigmaPar == 0.0f) LinearDiff_Inp_3D(Input, Mask, Output, lambdaPar, tau, dimX, dimY, dimZ);
+            else NonLinearDiff_Inp_3D(Input, Mask, Output, lambdaPar, sigmaPar2, tau, penaltytype, dimX, dimY, dimZ);
+		}
+	}
+	}
+    return *Output;
+}
+/********************************************************************/
+/***************************2D Functions*****************************/
+/********************************************************************/
+/* linear diffusion (heat equation) */
+float LinearDiff_Inp_2D(float *Input, unsigned char *Mask, float *Output, float lambdaPar, float tau, int dimX, int dimY)
+{
+	int i,j,i1,i2,j1,j2,index;
+	float e,w,n,s,e1,w1,n1,s1;
+	
+#pragma omp parallel for shared(Input,Mask) private(index,i,j,i1,i2,j1,j2,e,w,n,s,e1,w1,n1,s1)
+    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(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 = j*dimX+i;
+             
+             if (Mask[index] > 0) {
+				/*inpainting process*/
+                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 NonLinearDiff_Inp_2D(float *Input, unsigned char *Mask, float *Output, float lambdaPar, float sigmaPar, float tau, int penaltytype, int dimX, int dimY)
+{
+	int i,j,i1,i2,j1,j2,index;
+	float e,w,n,s,e1,w1,n1,s1;
+	
+#pragma omp parallel for shared(Input,Mask) private(index,i,j,i1,i2,j1,j2,e,w,n,s,e1,w1,n1,s1)
+    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(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 = j*dimX+i;
+                
+        if (Mask[index] > 0) {
+		/*inpainting process*/
+                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 =  signNDF_inc(e1);
+            else e1 = e1/sigmaPar;
+            
+            if (fabs(w1) > sigmaPar) w1 =  signNDF_inc(w1);
+            else w1 = w1/sigmaPar;
+            
+            if (fabs(n1) > sigmaPar) n1 =  signNDF_inc(n1);
+            else n1 = n1/sigmaPar;
+            
+            if (fabs(s1) > sigmaPar) s1 =  signNDF_inc(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));
+            }
+            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;
+            }
+            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;
+}
+/********************************************************************/
+/***************************3D Functions*****************************/
+/********************************************************************/
+/* linear diffusion (heat equation) */
+float LinearDiff_Inp_3D(float *Input, unsigned char *Mask, float *Output, float lambdaPar, float tau, int dimX, int dimY, int dimZ)
+{
+	int i,j,k,i1,i2,j1,j2,k1,k2,index;
+	float e,w,n,s,u,d,e1,w1,n1,s1,u1,d1;
+	
+#pragma omp parallel for shared(Input,Mask) private(index,i,j,i1,i2,j1,j2,e,w,n,s,e1,w1,n1,s1,k,k1,k2,u1,d1,u,d)
+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(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;
+            
+            if (Mask[index] > 0) {
+			/*inpainting process*/
+            
+                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;
+}
+
+float NonLinearDiff_Inp_3D(float *Input, unsigned char *Mask, float *Output, float lambdaPar, float sigmaPar, float tau, int penaltytype, int dimX, int dimY, int dimZ)
+{
+	int i,j,k,i1,i2,j1,j2,k1,k2,index;
+	float e,w,n,s,u,d,e1,w1,n1,s1,u1,d1;
+	
+#pragma omp parallel for shared(Input,Mask) private(index,i,j,i1,i2,j1,j2,e,w,n,s,e1,w1,n1,s1,k,k1,k2,u1,d1,u,d)
+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(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;
+            
+        if (Mask[index] > 0) {
+			/*inpainting process*/
+                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 =  signNDF_inc(e1);
+            else e1 = e1/sigmaPar;
+            
+            if (fabs(w1) > sigmaPar) w1 =  signNDF_inc(w1);
+            else w1 = w1/sigmaPar;
+            
+            if (fabs(n1) > sigmaPar) n1 =  signNDF_inc(n1);
+            else n1 = n1/sigmaPar;
+            
+            if (fabs(s1) > sigmaPar) s1 =  signNDF_inc(s1);
+            else s1 = s1/sigmaPar;
+            
+            if (fabs(u1) > sigmaPar) u1 =  signNDF_inc(u1);
+            else u1 = u1/sigmaPar;
+            
+            if (fabs(d1) > sigmaPar) d1 =  signNDF_inc(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;
+}
diff --git a/Core/inpainters_CPU/Diffusion_Inpaint_core.h b/Core/inpainters_CPU/Diffusion_Inpaint_core.h
new file mode 100644
index 0000000..0fc2608
--- /dev/null
+++ b/Core/inpainters_CPU/Diffusion_Inpaint_core.h
@@ -0,0 +1,60 @@
+/*
+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 linear and nonlinear diffusion [1,2] for inpainting task (2D/3D case)
+ * The minimisation is performed using explicit scheme. 
+ *
+ * Input Parameters:
+ * 1. Image/volume to inpaint
+ * 2. Mask of the same size as (1) in 'unsigned char' format  (ones mark the region to inpaint, zeros belong to the data)
+ * 3. lambda - regularization parameter
+ * 4. Edge-preserving parameter (sigma), when sigma equals to zero nonlinear diffusion -> linear diffusion
+ * 5. Number of iterations, for explicit scheme >= 150 is recommended 
+ * 6. tau - time-marching step for explicit scheme
+ * 7. Penalty type: 1 - Huber, 2 - Perona-Malik, 3 - Tukey Biweight
+ *
+ * Output:
+ * [1] Inpainted image/volume 
+ *
+ * This function is based on the paper by
+ * [1] Perona, P. and Malik, J., 1990. Scale-space and edge detection using anisotropic diffusion. IEEE Transactions on pattern analysis and machine intelligence, 12(7), pp.629-639.
+ * [2] Black, M.J., Sapiro, G., Marimont, D.H. and Heeger, D., 1998. Robust anisotropic diffusion. IEEE Transactions on image processing, 7(3), pp.421-432.
+ */
+
+ 
+#ifdef __cplusplus
+extern "C" {
+#endif
+CCPI_EXPORT float Diffusion_Inpaint_CPU_main(float *Input, unsigned char *Mask, float *Output, float lambdaPar, float sigmaPar, int iterationsNumb,  float tau, int penaltytype, int dimX, int dimY, int dimZ);
+CCPI_EXPORT float LinearDiff_Inp_2D(float *Input, unsigned char *Mask, float *Output, float lambdaPar, float tau, int dimX, int dimY);
+CCPI_EXPORT float NonLinearDiff_Inp_2D(float *Input, unsigned char *Mask, float *Output, float lambdaPar, float sigmaPar, float tau, int penaltytype, int dimX, int dimY);
+CCPI_EXPORT float LinearDiff_Inp_3D(float *Input, unsigned char *Mask, float *Output, float lambdaPar, float tau, int dimX, int dimY, int dimZ);
+CCPI_EXPORT float NonLinearDiff_Inp_3D(float *Input, unsigned char *Mask, float *Output, float lambdaPar, float sigmaPar, float tau, int penaltytype, int dimX, int dimY, int dimZ);
+#ifdef __cplusplus
+}
+#endif
diff --git a/Core/inpainters_CPU/NonlocalMarching_Inpaint_core.c b/Core/inpainters_CPU/NonlocalMarching_Inpaint_core.c
new file mode 100644
index 0000000..66be15c
--- /dev/null
+++ b/Core/inpainters_CPU/NonlocalMarching_Inpaint_core.c
@@ -0,0 +1,188 @@
+/*
+ * 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 "NonlocalMarching_Inpaint_core.h"
+#include "utils.h"
+
+
+/* C-OMP implementation of Nonlocal Vertical Marching inpainting method (2D case)
+ * The method is heuristic but computationally efficent (especially for larger images).
+ * It developed specifically to smoothly inpaint horizontal or inclined missing data regions in sinograms
+ * The method WILL not work satisfactory if you have lengthy vertical stripes of missing data
+ *
+ * Input:
+ * 1. 2D image or sinogram with horizontal or inclined regions of missing data
+ * 2. Mask of the same size as A in 'unsigned char' format  (ones mark the region to inpaint, zeros belong to the data)
+ * 3. Linear increment to increase searching window size in iterations, values from 1-3 is a good choice
+ *
+ * Output:
+ * 1. Inpainted image or a sinogram
+ * 2. updated mask
+ *
+ * Reference: TBA
+ */
+
+float NonlocalMarching_Inpaint_main(float *Input, unsigned char *M, float *Output, unsigned char *M_upd, int SW_increment, int iterationsNumb, int trigger, int dimX, int dimY, int dimZ)
+{
+    int i, j, i_m, j_m, counter, iter, iterations_number, W_fullsize, switchmask, switchcurr, counterElements;
+    float *Gauss_weights;
+    
+    /* copying M to M_upd */
+    copyIm_unchar(M, M_upd, dimX, dimY, 1);
+    
+    /* Copying the image */
+    copyIm(Input, Output, dimX, dimY, 1);
+    
+    /* Find how many inpainting iterations (equal to the number of ones) required based on a mask  */
+    if (iterationsNumb == 0) {
+        iterations_number = 0;
+        for (i=0; i<dimY*dimX; i++) {
+            if (M[i] == 1) iterations_number++;
+        }
+        if ((int)(iterations_number/dimY) > dimX) iterations_number = dimX;
+    }
+    else iterations_number = iterationsNumb;
+    
+    if (iterations_number == 0) printf("%s \n", "Nothing to inpaint, zero mask!");
+    else {
+        
+        printf("%s %i \n", "Max iteration number equals to:", iterations_number);
+        
+        /* Inpainting iterations run here*/
+        int W_halfsize = 1;
+        for(iter=0; iter < iterations_number; iter++) {
+            
+            //if (mod (iter, 2) == 0) {W_halfsize += 1;}
+            // printf("%i \n", W_halfsize);
+            
+            /* pre-calculation of Gaussian distance weights  */
+            W_fullsize = (int)(2*W_halfsize + 1); /*full size of similarity window */
+            Gauss_weights = (float*)calloc(W_fullsize*W_fullsize,sizeof(float ));
+            counter = 0;
+            for(i_m=-W_halfsize; i_m<=W_halfsize; i_m++) {
+                for(j_m=-W_halfsize; j_m<=W_halfsize; j_m++) {
+                    Gauss_weights[counter] = exp(-(pow((i_m), 2) + pow((j_m), 2))/(2*W_halfsize*W_halfsize));
+                    counter++;
+                }
+            }
+            
+            if (trigger == 0) {
+                /*Matlab*/
+#pragma omp parallel for shared(Output, M_upd, Gauss_weights) private(i, j, switchmask, switchcurr)
+                for(j=0; j<dimY; j++) {
+                    switchmask = 0;
+                    for(i=0; i<dimX; i++) {
+                        switchcurr = 0;
+                        if ((M_upd[j*dimX + i] == 1) && (switchmask == 0)) {
+                            /* perform inpainting of the current pixel */
+                            inpaint_func(Output, M_upd, Gauss_weights, i, j, dimX, dimY, W_halfsize, W_fullsize);
+                            /* add value to the mask*/
+                            M_upd[j*dimX + i] = 0;
+                            switchmask = 1; switchcurr = 1;
+                        }
+                        if ((M_upd[j*dimX + i] == 0) && (switchmask == 1) && (switchcurr == 0)) {
+                            /* perform inpainting of the previous (i-1) pixel */
+                            inpaint_func(Output, M_upd, Gauss_weights, i-1, j, dimX, dimY, W_halfsize, W_fullsize);
+                            /* add value to the mask*/
+                            M_upd[(j)*dimX + i-1] = 0;
+                            switchmask = 0;
+                        }
+                    }
+                }
+            }
+            else {
+                /*Python*/
+                /* find a point in the mask to inpaint */
+#pragma omp parallel for shared(Output, M_upd, Gauss_weights) private(i, j, switchmask, switchcurr)
+                for(i=0; i<dimX; i++) {
+                    switchmask = 0;
+                    for(j=0; j<dimY; j++) {
+                        switchcurr = 0;
+                        if ((M_upd[j*dimX + i] == 1) && (switchmask == 0)) {
+                            /* perform inpainting of the current pixel */
+                            inpaint_func(Output, M_upd, Gauss_weights, i, j, dimX, dimY, W_halfsize, W_fullsize);
+                            /* add value to the mask*/
+                            M_upd[j*dimX + i] = 0;
+                            switchmask = 1; switchcurr = 1;
+                        }
+                        if ((M_upd[j*dimX + i] == 0) && (switchmask == 1) && (switchcurr == 0)) {
+                            /* perform inpainting of the previous (j-1) pixel */
+                            inpaint_func(Output, M_upd, Gauss_weights, i, j-1, dimX, dimY, W_halfsize, W_fullsize);
+                            /* add value to the mask*/
+                            M_upd[(j-1)*dimX + i] = 0;
+                            switchmask = 0;
+                        }
+                    }
+                }
+            }
+            free(Gauss_weights);
+            
+            /* check if possible to terminate iterations earlier */
+            counterElements = 0;
+            for(i=0; i<dimX*dimY; i++) if (M_upd[i] == 0) counterElements++;
+            
+            if (counterElements == dimX*dimY) {
+                printf("%s \n", "Padding completed!");
+                break;
+            }
+            W_halfsize += SW_increment;
+        }
+        printf("%s %i \n", "Iterations stopped at:", iter);
+    }
+    return *Output;
+}
+
+float inpaint_func(float *U, unsigned char *M_upd, float *Gauss_weights, int i, int j, int dimX, int dimY, int W_halfsize, int W_fullsize)
+{
+    int i1, j1, i_m, j_m, counter;
+    float sum_val, sumweight;
+    
+    /*method 1: inpainting based on Euclidian weights */
+    sumweight = 0.0f;
+    counter = 0; sum_val = 0.0f;
+    for(i_m=-W_halfsize; i_m<=W_halfsize; i_m++) {
+        i1 = i+i_m;
+        for(j_m=-W_halfsize; j_m<=W_halfsize; j_m++) {
+            j1 = j+j_m;
+            if (((i1 >= 0) && (i1 < dimX)) && ((j1 >= 0) && (j1 < dimY))) {
+                if (M_upd[j1*dimX + i1] == 0) {
+                    sumweight += Gauss_weights[counter];
+                }
+            }
+            counter++;
+        }
+    }
+    counter = 0; sum_val = 0.0f;
+    for(i_m=-W_halfsize; i_m<=W_halfsize; i_m++) {
+        i1 = i+i_m;
+        for(j_m=-W_halfsize; j_m<=W_halfsize; j_m++) {
+            j1 = j+j_m;
+            if (((i1 >= 0) && (i1 < dimX)) && ((j1 >= 0) && (j1 < dimY))) {
+                if ((M_upd[j1*dimX + i1] == 0) && (sumweight != 0.0f)) {
+                    /* we have data so add it with Euc weight */
+                    sum_val += (Gauss_weights[counter]/sumweight)*U[j1*dimX + i1];
+                }
+            }
+            counter++;
+        }
+    }
+    U[j*dimX + i] = sum_val;
+    return *U;
+}
+
diff --git a/Core/inpainters_CPU/NonlocalMarching_Inpaint_core.h b/Core/inpainters_CPU/NonlocalMarching_Inpaint_core.h
new file mode 100644
index 0000000..0f99ed4
--- /dev/null
+++ b/Core/inpainters_CPU/NonlocalMarching_Inpaint_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 Nonlocal Vertical Marching inpainting method (2D case)
+ * The method is heuristic but computationally efficent (especially for larger images).
+ * It developed specifically to smoothly inpaint horizontal or inclined missing data regions in sinograms
+ * The method WILL not work satisfactory if you have lengthy vertical stripes of missing data
+ *
+ * Inputs:
+ * 1. 2D image or sinogram with horizontal or inclined regions of missing data
+ * 2. Mask of the same size as A in 'unsigned char' format  (ones mark the region to inpaint, zeros belong to the data)
+ * 3. Linear increment to increase searching window size in iterations, values from 1-3 is a good choice
+
+ * Output:
+ * 1. Inpainted image or a sinogram
+ * 2. updated mask
+ *
+ * Reference: TBA
+ */
+
+ 
+#ifdef __cplusplus
+extern "C" {
+#endif
+CCPI_EXPORT float NonlocalMarching_Inpaint_main(float *Input, unsigned char *M, float *Output, unsigned char *M_upd, int SW_increment, int iterationsNumb, int trigger, int dimX, int dimY, int dimZ);
+CCPI_EXPORT float inpaint_func(float *U, unsigned char *M_upd, float *Gauss_weights, int i, int j, int dimX, int dimY, int W_halfsize, int W_fullsize);
+#ifdef __cplusplus
+}
+#endif
diff --git a/Core/regularisers_CPU/TNV_core.c b/Core/regularisers_CPU/TNV_core.c
index 8e61869..c51d6cd 100755
--- a/Core/regularisers_CPU/TNV_core.c
+++ b/Core/regularisers_CPU/TNV_core.c
@@ -251,18 +251,19 @@ float gradient(float *u_upd, float *gradx_upd, float *grady_upd, int dimX, int d
     // Compute discrete gradient using forward differences
 #pragma omp parallel for shared(gradx_upd,grady_upd,u_upd) private(i, j, k, l)
     for(k = 0; k < dimZ; k++)   {
-        for(j = 0; j < dimX; j++)   {
-            l = j * dimY;
-            for(i = 0; i < dimY; i++)   {
+        for(j = 0; j < dimY; j++)   {
+            l = j * dimX;           
+            for(i = 0; i < dimX; i++)   {
                 // Derivatives in the x-direction
-                if(i != dimY-1)
+                if(i != dimX-1)
                     gradx_upd[(dimX*dimY)*k + i+l] = u_upd[(dimX*dimY)*k + i+1+l] - u_upd[(dimX*dimY)*k + i+l];
                 else
                     gradx_upd[(dimX*dimY)*k + i+l] = 0.0f;
                 
                 // Derivatives in the y-direction
-                if(j != dimX-1)
-                    grady_upd[(dimX*dimY)*k + i+l] = u_upd[(dimX*dimY)*k + i+dimY+l] -u_upd[(dimX*dimY)*k + i+l];
+                if(j != dimY-1)
+                    //grady_upd[(dimX*dimY)*k + i+l] = u_upd[(dimX*dimY)*k + i+dimY+l] -u_upd[(dimX*dimY)*k + i+l];
+                    grady_upd[(dimX*dimY)*k + i+l] = u_upd[(dimX*dimY)*k + i+(j+1)*dimX] -u_upd[(dimX*dimY)*k + i+l];
                 else
                     grady_upd[(dimX*dimY)*k + i+l] = 0.0f;
             }}}
@@ -301,8 +302,8 @@ float proxF(float *gx, float *gy, float *vx, float *vy, float sigma, int p, int
             
             for(k = 0; k < dimZ; k++)
             {
-                valuex = vx[(dimX*dimY)*k + (i)*dimY + (j)];
-                valuey = vy[(dimX*dimY)*k + (i)*dimY + (j)];
+                valuex = vx[(dimX*dimY)*k + (j)*dimX + (i)];
+                valuey = vy[(dimX*dimY)*k + (j)*dimX + (i)];
                 
                 M1 += (valuex * valuex);
                 M2 += (valuex * valuey);
@@ -412,9 +413,9 @@ float proxF(float *gx, float *gy, float *vx, float *vy, float sigma, int p, int
             
             for(k = 0; k < dimZ; k++)
             {
-                gx[(dimX*dimY)*k + (i)*dimY + (j)] = vx[(dimX*dimY)*k + (i)*dimY + (j)] * t1 + vy[(dimX*dimY)*k + (i)*dimY + (j)] * t2;
-                gy[(dimX*dimY)*k + (i)*dimY + (j)] = vx[(dimX*dimY)*k + (i)*dimY + (j)] * t2 + vy[(dimX*dimY)*k + (i)*dimY + (j)] * t3;
-            }
+                gx[(dimX*dimY)*k + j*dimX + i] = vx[(dimX*dimY)*k + j*dimX + i] * t1 + vy[(dimX*dimY)*k + j*dimX + i] * t2;
+                gy[(dimX*dimY)*k + j*dimX + i] = vx[(dimX*dimY)*k + j*dimX + i] * t2 + vy[(dimX*dimY)*k + j*dimX + i] * t3;
+            }           
             
             // Delete allocated memory
             free(proj);
@@ -428,20 +429,21 @@ float divergence(float *qx_upd, float *qy_upd, float *div_upd, int dimX, int dim
     int i, j, k, l;
 #pragma omp parallel for shared(qx_upd,qy_upd,div_upd) private(i, j, k, l)
     for(k = 0; k < dimZ; k++)   {
-        for(j = 0; j < dimX; j++)   {
-            l = j * dimY;
-            for(i = 0; i < dimY; i++)   {
-                if(i != dimY-1)
+        for(j = 0; j < dimY; j++)   {
+            l = j * dimX;            
+            for(i = 0; i < dimX; i++)   {
+                if(i != dimX-1)
                 {
                     // ux[k][i+l] = u[k][i+1+l] - u[k][i+l]
                     div_upd[(dimX*dimY)*k + i+1+l] -= qx_upd[(dimX*dimY)*k + i+l];
                     div_upd[(dimX*dimY)*k + i+l] += qx_upd[(dimX*dimY)*k + i+l];
                 }
                 
-                if(j != dimX-1)
+                if(j != dimY-1)
                 {
                     // uy[k][i+l] = u[k][i+width+l] - u[k][i+l]
-                    div_upd[(dimX*dimY)*k + i+dimY+l] -= qy_upd[(dimX*dimY)*k + i+l];
+                    //div_upd[(dimX*dimY)*k + i+dimY+l] -= qy_upd[(dimX*dimY)*k + i+l];
+                    div_upd[(dimX*dimY)*k + i+(j+1)*dimX] -= qy_upd[(dimX*dimY)*k + i+l];                    
                     div_upd[(dimX*dimY)*k + i+l] += qy_upd[(dimX*dimY)*k + i+l];
                 }
             }
diff --git a/Core/regularisers_CPU/utils.c b/Core/regularisers_CPU/utils.c
index ca5c59a..3eebdaa 100644
--- a/Core/regularisers_CPU/utils.c
+++ b/Core/regularisers_CPU/utils.c
@@ -20,7 +20,7 @@ limitations under the License.
 #include "utils.h"
 #include <math.h>
 
-/* Copy Image */
+/* Copy Image (float) */
 float copyIm(float *A, float *U, int dimX, int dimY, int dimZ)
 {
 	int j;
@@ -29,12 +29,33 @@ float copyIm(float *A, float *U, int dimX, int dimY, int dimZ)
 	return *U;
 }
 
-/*static inline int8_t SIGN(int val) {
- if (val < 0) return -1;
- if (val==0) return 0;
- return 1;
+/* Copy Image */
+unsigned char copyIm_unchar(unsigned char *A, unsigned char *U, int dimX, int dimY, int dimZ)
+{
+	int j;
+#pragma omp parallel for shared(A, U) private(j)
+	for (j = 0; j<dimX*dimY*dimZ; j++)  U[j] = A[j];
+	return *U;
+}
+
+/*Roll image symmetrically from top to bottom*/
+float copyIm_roll(float *A, float *U, int dimX, int dimY, int roll_value, int switcher)
+{
+    int i, j;
+#pragma omp parallel for shared(U, A) private(i,j)
+    for (i=0; i<dimX; i++) {
+        for (j=0; j<dimY; j++) {
+            if (switcher == 0) {
+                if (j < (dimY - roll_value)) U[j*dimX + i] = A[(j+roll_value)*dimX + i];
+                else U[j*dimX + i] = A[(j - (dimY - roll_value))*dimX + i];
+            }
+            else {
+                if (j < roll_value) U[j*dimX + i] = A[(j+(dimY - roll_value))*dimX + i];
+                else U[j*dimX + i] = A[(j - roll_value)*dimX + i];
+            }
+        }}
+    return *U;
 }
-*/
 
 /* function that calculates TV energy
  * type - 1:  2*lambda*min||\nabla u|| + ||u -u0||^2
diff --git a/Core/regularisers_CPU/utils.h b/Core/regularisers_CPU/utils.h
index 866bc01..7ad83ab 100644
--- a/Core/regularisers_CPU/utils.h
+++ b/Core/regularisers_CPU/utils.h
@@ -17,17 +17,16 @@ See the License for the specific language governing permissions and
 limitations under the License.
 */
 
-//#include <matrix.h>
-//#include <math.h>
 #include <stdlib.h>
 #include <memory.h>
 #include "CCPiDefines.h"
-//#include <stdio.h>
 #include "omp.h"
 #ifdef __cplusplus
 extern "C" {
 #endif
 CCPI_EXPORT float copyIm(float *A, float *U, int dimX, int dimY, int dimZ);
+CCPI_EXPORT unsigned char copyIm_unchar(unsigned char *A, unsigned char *U, int dimX, int dimY, int dimZ);
+CCPI_EXPORT float copyIm_roll(float *A, float *U, int dimX, int dimY, int roll_value, int switcher);
 CCPI_EXPORT float TV_energy2D(float *U, float *U0, float *E_val, float lambda, int type, int dimX, int dimY);
 CCPI_EXPORT float TV_energy3D(float *U, float *U0, float *E_val, float lambda, int type, int dimX, int dimY, int dimZ);
 #ifdef __cplusplus
diff --git a/Readme.md b/Readme.md
index e73b4fb..8a67c60 100644
--- a/Readme.md
+++ b/Readme.md
@@ -16,17 +16,22 @@ the toolkit can be used independently to solve image denoising problems. The cor
  * C compilers
  * nvcc (CUDA SDK) compilers
 
-## Package modules (regularisers):
+## Package modules:
 
-### Single-channel
+### Single-channel (denoising):
 1. Rudin-Osher-Fatemi (ROF) Total Variation (explicit PDE minimisation scheme) **2D/3D CPU/GPU** (Ref. *1*)
 2. Fast-Gradient-Projection (FGP) Total Variation **2D/3D CPU/GPU** (Ref. *2*)
-3. Split-Bregman (SB) Total Variation **2D/3D CPU/GPU** (Ref. *4*)
-4. Linear and nonlinear diffusion (explicit PDE minimisation scheme) **2D/3D CPU/GPU** (Ref. *6*)
+3. Split-Bregman (SB) Total Variation **2D/3D CPU/GPU** (Ref. *5*)
+4. Linear and nonlinear diffusion (explicit PDE minimisation scheme) **2D/3D CPU/GPU** (Ref. *7*)
+
+### Multi-channel (denoising):
+1. Fast-Gradient-Projection (FGP) Directional Total Variation **2D/3D CPU/GPU** (Ref. *3,4,2*)
+2. Total Nuclear Variation (TNV) penalty **2D+channels CPU** (Ref. *6*)
+
+### Inpainting:
+1. Linear and nonlinear diffusion (explicit PDE minimisation scheme) **2D/3D CPU** (Ref. *7*)
+2. Iterative nonlocal vertical marching method  **2D CPU**
 
-### Multi-channel
-1. Fast-Gradient-Projection (FGP) Directional Total Variation **2D/3D CPU/GPU** (Ref. *3,2*)
-2. Total Nuclear Variation (TNV) penalty **2D+channels CPU** (Ref. *5*)
 
 ## Installation:
 
@@ -56,11 +61,13 @@ the toolkit can be used independently to solve image denoising problems. The cor
 
 *3. Ehrhardt, M.J. and Betcke, M.M., 2016. Multicontrast MRI reconstruction with structure-guided total variation. SIAM Journal on Imaging Sciences, 9(3), pp.1084-1106.*
 
-*4. Goldstein, T. and Osher, S., 2009. The split Bregman method for L1-regularized problems. SIAM journal on imaging sciences, 2(2), pp.323-343.*
+*4. Kazantsev, D., Jørgensen, J.S., Andersen, M., Lionheart, W.R., Lee, P.D. and Withers, P.J., 2018. Joint image reconstruction method with correlative multi-channel prior for X-ray spectral computed tomography. Inverse Problems*
+
+*5. Goldstein, T. and Osher, S., 2009. The split Bregman method for L1-regularized problems. SIAM journal on imaging sciences, 2(2), pp.323-343.*
 
-*5. Duran, J., Moeller, M., Sbert, C. and Cremers, D., 2016. Collaborative total variation: a general framework for vectorial TV models. SIAM Journal on Imaging Sciences, 9(1), pp.116-151.*
+*6. Duran, J., Moeller, M., Sbert, C. and Cremers, D., 2016. Collaborative total variation: a general framework for vectorial TV models. SIAM Journal on Imaging Sciences, 9(1), pp.116-151.*
 
-*6. Black, M.J., Sapiro, G., Marimont, D.H. and Heeger, D., 1998. Robust anisotropic diffusion. IEEE Transactions on image processing, 7(3), pp.421-432.*
+*7. Black, M.J., Sapiro, G., Marimont, D.H. and Heeger, D., 1998. Robust anisotropic diffusion. IEEE Transactions on image processing, 7(3), pp.421-432.*
 
 ### License:
 [Apache License, Version 2.0](http://www.apache.org/licenses/LICENSE-2.0)
diff --git a/Wrappers/Matlab/demos/demoMatlab_3Ddenoise.m b/Wrappers/Matlab/demos/demoMatlab_3Ddenoise.m
index 5a54d18..c087433 100644
--- a/Wrappers/Matlab/demos/demoMatlab_3Ddenoise.m
+++ b/Wrappers/Matlab/demos/demoMatlab_3Ddenoise.m
@@ -1,8 +1,9 @@
 % Volume (3D) denoising demo using CCPi-RGL
-clear
-close all
-addpath('../mex_compile/installed');
-addpath('../../../data/');
+clear; close all
+Path1 = sprintf(['..' filesep 'mex_compile' filesep 'installed'], 1i);
+Path2 = sprintf(['..' filesep '..' filesep '..' filesep 'data' filesep], 1i);
+addpath(Path1);
+addpath(Path2);
 
 N = 512; 
 slices = 30;
diff --git a/Wrappers/Matlab/demos/demoMatlab_denoise.m b/Wrappers/Matlab/demos/demoMatlab_denoise.m
index 151a604..d93f477 100644
--- a/Wrappers/Matlab/demos/demoMatlab_denoise.m
+++ b/Wrappers/Matlab/demos/demoMatlab_denoise.m
@@ -1,8 +1,9 @@
 % Image (2D) denoising demo using CCPi-RGL
-clear
-close all
-addpath('../mex_compile/installed');
-addpath('../../../data/');
+clear; close all
+Path1 = sprintf(['..' filesep 'mex_compile' filesep 'installed'], 1i);
+Path2 = sprintf(['..' filesep '..' filesep '..' filesep 'data' filesep], 1i);
+addpath(Path1);
+addpath(Path2);
 
 Im = double(imread('lena_gray_512.tif'))/255;  % loading image
 u0 = Im + .05*randn(size(Im)); u0(u0 < 0) = 0;
diff --git a/Wrappers/Matlab/demos/demoMatlab_inpaint.m b/Wrappers/Matlab/demos/demoMatlab_inpaint.m
new file mode 100644
index 0000000..66f9c15
--- /dev/null
+++ b/Wrappers/Matlab/demos/demoMatlab_inpaint.m
@@ -0,0 +1,35 @@
+% Image (2D) inpainting demo using CCPi-RGL
+clear; close all
+Path1 = sprintf(['..' filesep 'mex_compile' filesep 'installed'], 1i);
+Path2 = sprintf(['..' filesep '..' filesep '..' filesep 'data' filesep], 1i);
+addpath(Path1);
+addpath(Path2);
+
+load('SinoInpaint.mat');
+Sinogram = Sinogram./max(Sinogram(:));
+Sino_mask = Sinogram.*(1-single(Mask));
+figure; 
+subplot(1,2,1); imshow(Sino_mask, [0 1]); title('Missing data sinogram');
+subplot(1,2,2); imshow(Mask, [0 1]); title('Mask');
+%%
+fprintf('Inpaint using Linear-Diffusion model (CPU) \n');
+iter_diff = 5000; % number of diffusion iterations
+lambda_regDiff = 6000; % regularisation for the diffusivity 
+sigmaPar = 0.0; % edge-preserving parameter
+tau_param = 0.000075; % time-marching constant 
+tic; u_diff = NonlDiff_Inp(single(Sino_mask), Mask, lambda_regDiff, sigmaPar, iter_diff, tau_param); toc; 
+figure; imshow(u_diff, [0 1]); title('Linear-Diffusion inpainted sinogram (CPU)');
+%%
+fprintf('Inpaint using Nonlinear-Diffusion model (CPU) \n');
+iter_diff = 1500; % number of diffusion iterations
+lambda_regDiff = 80; % regularisation for the diffusivity 
+sigmaPar = 0.00009; % edge-preserving parameter
+tau_param = 0.000008; % time-marching constant 
+tic; u_diff = NonlDiff_Inp(single(Sino_mask), Mask, lambda_regDiff, sigmaPar, iter_diff, tau_param, 'Huber'); toc; 
+figure; imshow(u_diff, [0 1]); title('Non-Linear Diffusion inpainted sinogram (CPU)');
+%%
+fprintf('Inpaint using Nonlocal Vertical Marching model (CPU) \n');
+Increment = 1; % linear increment for the searching window
+tic; [u_nom,maskupd] = NonlocalMarching_Inpaint(single(Sino_mask), Mask, Increment); toc;
+figure; imshow(u_nom, [0 1]); title('NVM inpainted sinogram (CPU)');
+%%
\ No newline at end of file
diff --git a/Wrappers/Matlab/mex_compile/compileCPU_mex.m b/Wrappers/Matlab/mex_compile/compileCPU_mex.m
index ee0d99e..b232f33 100644
--- a/Wrappers/Matlab/mex_compile/compileCPU_mex.m
+++ b/Wrappers/Matlab/mex_compile/compileCPU_mex.m
@@ -2,9 +2,11 @@
 
 pathcopyFrom = sprintf(['..' filesep '..' filesep '..' filesep 'Core' filesep 'regularisers_CPU'], 1i);
 pathcopyFrom1 = sprintf(['..' filesep '..' filesep '..' filesep 'Core' filesep 'CCPiDefines.h'], 1i);
+pathcopyFrom2 = sprintf(['..' filesep '..' filesep '..' filesep 'Core' filesep 'inpainters_CPU'], 1i);
 
 copyfile(pathcopyFrom, 'regularisers_CPU');
 copyfile(pathcopyFrom1, 'regularisers_CPU');
+copyfile(pathcopyFrom2, 'regularisers_CPU');
 
 cd regularisers_CPU
 
@@ -32,9 +34,17 @@ movefile('NonlDiff.mex*',Pathmove);
 mex TV_energy.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
 movefile('TV_energy.mex*',Pathmove);
 
+%############Inpainters##############%
+mex NonlDiff_Inp.c Diffusion_Inpaint_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
+movefile('NonlDiff_Inp.mex*',Pathmove);
+
+mex NonlocalMarching_Inpaint.c NonlocalMarching_Inpaint_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
+movefile('NonlocalMarching_Inpaint.mex*',Pathmove);
+
 delete SB_TV_core* ROF_TV_core* FGP_TV_core* FGP_dTV_core* TNV_core* utils* Diffusion_core* CCPiDefines.h
-fprintf('%s \n', 'All successfully compiled!');
+delete Diffusion_Inpaint_core* NonlocalMarching_Inpaint_core*
+fprintf('%s \n', 'Regularisers successfully compiled!');
 
-pathA = sprintf(['..' filesep '..' filesep], 1i);
-cd(pathA);
+pathA2 = sprintf(['..' filesep '..' filesep], 1i);
+cd(pathA2);
 cd demos
diff --git a/Wrappers/Matlab/mex_compile/regularisers_CPU/NonlDiff_Inp.c b/Wrappers/Matlab/mex_compile/regularisers_CPU/NonlDiff_Inp.c
new file mode 100644
index 0000000..eaab4a7
--- /dev/null
+++ b/Wrappers/Matlab/mex_compile/regularisers_CPU/NonlDiff_Inp.c
@@ -0,0 +1,101 @@
+/*
+ * 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 "Diffusion_Inpaint_core.h"
+
+/* C-OMP implementation of linear and nonlinear diffusion [1,2] for inpainting task (2D/3D case)
+ * The minimisation is performed using explicit scheme. 
+ *
+ * Input Parameters:
+ * 1. Image/volume to inpaint
+ * 2. Inpainting Mask of the same size as (1) in 'unsigned char' format  (ones mark the region to inpaint, zeros belong to the data)
+ * 3. lambda - regularization parameter
+ * 4. Edge-preserving parameter (sigma), when sigma equals to zero nonlinear diffusion -> linear diffusion
+ * 5. Number of iterations, for explicit scheme >= 150 is recommended 
+ * 6. tau - time-marching step for explicit scheme
+ * 7. Penalty type: 1 - Huber, 2 - Perona-Malik, 3 - Tukey Biweight
+ *
+ * Output:
+ * [1] Inpainted image/volume 
+ *
+ * This function is based on the paper by
+ * [1] Perona, P. and Malik, J., 1990. Scale-space and edge detection using anisotropic diffusion. IEEE Transactions on pattern analysis and machine intelligence, 12(7), pp.629-639.
+ * [2] Black, M.J., Sapiro, G., Marimont, D.H. and Heeger, D., 1998. Robust anisotropic diffusion. IEEE Transactions on image processing, 7(3), pp.421-432.
+ */
+
+void mexFunction(
+        int nlhs, mxArray *plhs[],
+        int nrhs, const mxArray *prhs[])
+        
+{
+    int number_of_dims, iter_numb, dimX, dimY, dimZ, penaltytype, i, inpaint_elements;
+    const int  *dim_array;
+    const int  *dim_array2;
+    float *Input, *Output=NULL, lambda, tau, sigma;
+    unsigned char *Mask;
+    
+    dim_array = mxGetDimensions(prhs[0]);
+    dim_array2 = mxGetDimensions(prhs[1]);
+    number_of_dims = mxGetNumberOfDimensions(prhs[0]);
+    
+    /*Handling Matlab input data*/
+    Input  = (float *) mxGetData(prhs[0]);
+    Mask  = (unsigned char *) mxGetData(prhs[1]); /* MASK */
+    lambda =  (float) mxGetScalar(prhs[2]); /* regularization parameter */
+    sigma = (float) mxGetScalar(prhs[3]); /* Edge-preserving parameter */
+    iter_numb = 300; /* iterations number */
+    tau = 0.025; /* marching step parameter */
+    penaltytype = 1; /* Huber penalty by default */    
+  
+    if ((nrhs < 4) || (nrhs > 7)) mexErrMsgTxt("At least 4 parameters is required, all parameters are: Image(2D/3D), Mask(2D/3D), Regularisation parameter, Edge-preserving parameter, iterations number, time-marching constant, penalty type - Huber, PM or Tukey");
+    if ((nrhs == 5) || (nrhs == 6) || (nrhs == 7))  iter_numb = (int) mxGetScalar(prhs[4]); /* iterations number */
+    if ((nrhs == 6) || (nrhs == 7))  tau =  (float) mxGetScalar(prhs[5]); /* marching step parameter */
+    if (nrhs == 7)  {
+        char *penalty_type;
+        penalty_type = mxArrayToString(prhs[6]); /* Huber, PM or Tukey 'Huber' is the default */
+        if ((strcmp(penalty_type, "Huber") != 0) && (strcmp(penalty_type, "PM") != 0) && (strcmp(penalty_type, "Tukey") != 0)) mexErrMsgTxt("Choose penalty: 'Huber', 'PM' or 'Tukey',");
+        if (strcmp(penalty_type, "Huber") == 0)  penaltytype = 1;  /* enable 'Huber' penalty */
+        if (strcmp(penalty_type, "PM") == 0)  penaltytype = 2;  /* enable Perona-Malik penalty */
+        if (strcmp(penalty_type, "Tukey") == 0)  penaltytype = 3;  /* enable Tikey Biweight penalty */
+        mxFree(penalty_type);
+    }    
+    
+    if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); }
+    if (mxGetClassID(prhs[1]) != mxUINT8_CLASS) {mexErrMsgTxt("The mask must be in uint8 precision");}
+    
+    dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2];
+    
+    /* output arrays*/
+    if (number_of_dims == 2) {
+        dimZ = 1; /*2D case*/
+        /* output image/volume */
+        if ((dimX != dim_array2[0]) || (dimY != dim_array2[1])) mexErrMsgTxt("Input image and the provided mask are of different dimensions!");
+        Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
+    }    
+    if (number_of_dims == 3) {
+        if ((dimX != dim_array2[0]) || (dimY != dim_array2[1]) || (dimZ != dim_array2[2])) mexErrMsgTxt("Input image and the provided mask are of different dimensions!");
+        Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
+    }    
+    
+    inpaint_elements = 0;
+    for (i=0; i<dimY*dimX*dimZ; i++) if (Mask[i] == 1) inpaint_elements++;
+    if (inpaint_elements == 0) mexErrMsgTxt("The mask is full of zeros, nothing to inpaint");        
+    Diffusion_Inpaint_CPU_main(Input, Mask, Output, lambda, sigma, iter_numb, tau, penaltytype, dimX, dimY, dimZ);
+}
\ No newline at end of file
diff --git a/Wrappers/Matlab/mex_compile/regularisers_CPU/NonlocalMarching_Inpaint.c b/Wrappers/Matlab/mex_compile/regularisers_CPU/NonlocalMarching_Inpaint.c
new file mode 100644
index 0000000..36cf05c
--- /dev/null
+++ b/Wrappers/Matlab/mex_compile/regularisers_CPU/NonlocalMarching_Inpaint.c
@@ -0,0 +1,82 @@
+/*
+ * 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 "NonlocalMarching_Inpaint_core.h"
+
+/* C-OMP implementation of Nonlocal Vertical Marching inpainting method (2D case)
+ * The method is heuristic but computationally efficent (especially for larger images).
+ * It developed specifically to smoothly inpaint horizontal or inclined missing data regions in sinograms
+ * The method WILL not work satisfactory if you have lengthy vertical stripes of missing data
+ *
+ * Input:
+ * 1. 2D image or sinogram [REQUIRED]
+ * 2. Mask of the same size as A in 'unsigned char' format  (ones mark the region to inpaint, zeros belong to the data) [REQUIRED]
+ * 3. Linear increment to increase searching window size in iterations, values from 1-3 is a good choice [OPTIONAL, default 1]
+ * 4. Number of iterations [OPTIONAL, default - calculate based on the mask]
+ *
+ * Output:
+ * 1. Inpainted sinogram  
+ * 2. updated mask
+ * Reference: TBA
+ */
+
+void mexFunction(
+        int nlhs, mxArray *plhs[],
+        int nrhs, const mxArray *prhs[])
+        
+{
+    int number_of_dims, dimX, dimY, dimZ, iterations, SW_increment;
+    const int  *dim_array;
+    const int  *dim_array2;
+    float *Input, *Output=NULL;
+    unsigned char *Mask, *Mask_upd=NULL;
+    
+    dim_array = mxGetDimensions(prhs[0]);
+    dim_array2 = mxGetDimensions(prhs[1]);
+    number_of_dims = mxGetNumberOfDimensions(prhs[0]);
+    
+    /*Handling Matlab input data*/
+    Input  = (float *) mxGetData(prhs[0]);
+    Mask  = (unsigned char *) mxGetData(prhs[1]); /* MASK */    
+    SW_increment = 1;
+    iterations = 0;
+            
+    if ((nrhs < 2) || (nrhs > 4)) mexErrMsgTxt("At least 4 parameters is required, all parameters are: Image(2D/3D), Mask(2D/3D), Linear increment, Iterations number");
+    if ((nrhs == 3) || (nrhs == 4))  SW_increment =  (int) mxGetScalar(prhs[2]); /* linear increment */
+    if ((nrhs == 4))  iterations =  (int) mxGetScalar(prhs[3]); /* iterations number */
+       
+    if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); }
+    if (mxGetClassID(prhs[1]) != mxUINT8_CLASS) {mexErrMsgTxt("The mask must be in uint8 precision");}    
+    
+    dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2];
+    
+    /* output arrays*/
+    if (number_of_dims == 2) {
+        dimZ = 1; /*2D case*/
+        /* output image/volume */
+        if ((dimX != dim_array2[0]) || (dimY != dim_array2[1])) mexErrMsgTxt("Input image and the provided mask are of different dimensions!");
+        Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
+        Mask_upd = (unsigned char*)mxGetPr(plhs[1] = mxCreateNumericArray(2, dim_array, mxUINT8_CLASS, mxREAL));
+    }    
+    if (number_of_dims == 3) {
+        mexErrMsgTxt("Currently 2D supported only");        
+    }           
+    NonlocalMarching_Inpaint_main(Input, Mask, Output, Mask_upd, SW_increment, iterations, 0, dimX, dimY, dimZ);
+}
\ No newline at end of file
diff --git a/Wrappers/Python/CMakeLists.txt b/Wrappers/Python/CMakeLists.txt
index fb00706..7833b54 100644
--- a/Wrappers/Python/CMakeLists.txt
+++ b/Wrappers/Python/CMakeLists.txt
@@ -81,8 +81,3 @@ else()
 endif()
 
 configure_file("setup-regularisers.py.in" "setup-regularisers.py")
-
-
-#add_executable(regulariser_test ${CMAKE_CURRENT_SOURCE_DIR}/test/test_regulariser.cpp)
-
-#target_link_libraries (regulariser_test LINK_PUBLIC regularisers_lib)
diff --git a/Wrappers/Python/ccpi/filters/regularisers.py b/Wrappers/Python/ccpi/filters/regularisers.py
index eec8c4d..a07b39a 100644
--- a/Wrappers/Python/ccpi/filters/regularisers.py
+++ b/Wrappers/Python/ccpi/filters/regularisers.py
@@ -2,7 +2,7 @@
 script which assigns a proper device core function based on a flag ('cpu' or 'gpu')
 """
 
-from ccpi.filters.cpu_regularisers_cython import TV_ROF_CPU, TV_FGP_CPU, TV_SB_CPU, dTV_FGP_CPU, TNV_CPU, NDF_CPU
+from ccpi.filters.cpu_regularisers import TV_ROF_CPU, TV_FGP_CPU, TV_SB_CPU, dTV_FGP_CPU, TNV_CPU, NDF_CPU, NDF_INPAINT_CPU, NVM_INPAINT_CPU
 from ccpi.filters.gpu_regularisers import TV_ROF_GPU, TV_FGP_GPU, TV_SB_GPU, dTV_FGP_GPU, NDF_GPU
 
 def ROF_TV(inputData, regularisation_parameter, iterations,
@@ -110,3 +110,10 @@ def NDF(inputData, regularisation_parameter, edge_parameter, iterations,
     else:
         raise ValueError('Unknown device {0}. Expecting gpu or cpu'\
                          .format(device))
+def NDF_INP(inputData, maskData, regularisation_parameter, edge_parameter, iterations,
+                     time_marching_parameter, penalty_type):
+        return NDF_INPAINT_CPU(inputData, maskData, regularisation_parameter, 
+        edge_parameter, iterations, time_marching_parameter, penalty_type)
+        
+def NVM_INP(inputData, maskData, SW_increment, iterations):
+        return NVM_INPAINT_CPU(inputData, maskData, SW_increment, iterations)
diff --git a/Wrappers/Python/demos/demo_cpu_inpainters.py b/Wrappers/Python/demos/demo_cpu_inpainters.py
new file mode 100644
index 0000000..3b4191b
--- /dev/null
+++ b/Wrappers/Python/demos/demo_cpu_inpainters.py
@@ -0,0 +1,192 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Demonstration of CPU inpainters
+@authors: Daniil Kazantsev, Edoardo Pasca
+"""
+
+import matplotlib.pyplot as plt
+import numpy as np
+import os
+import timeit
+from scipy import io
+from ccpi.filters.regularisers import NDF_INP, NVM_INP
+from qualitymetrics import rmse
+###############################################################################
+def printParametersToString(pars):
+        txt = r''
+        for key, value in pars.items():
+            if key== 'algorithm' :
+                txt += "{0} = {1}".format(key, value.__name__)
+            elif key == 'input':
+                txt += "{0} = {1}".format(key, np.shape(value))
+            elif key == 'maskData':
+                txt += "{0} = {1}".format(key, np.shape(value))
+            else:
+                txt += "{0} = {1}".format(key, value)
+            txt += '\n'
+        return txt
+###############################################################################
+
+# read sinogram and the mask
+filename = os.path.join(".." , ".." , ".." , "data" ,"SinoInpaint.mat")
+sino = io.loadmat(filename)
+sino_full = sino.get('Sinogram')
+Mask = sino.get('Mask')
+[angles_dim,detectors_dim] = sino_full.shape
+sino_full = sino_full/np.max(sino_full)
+#apply mask to sinogram
+sino_cut = sino_full*(1-Mask)
+#sino_cut_new = np.zeros((angles_dim,detectors_dim),'float32')
+#sino_cut_new = sino_cut.copy(order='c')
+#sino_cut_new[:] = sino_cut[:]
+sino_cut_new = np.ascontiguousarray(sino_cut, dtype=np.float32);
+#mask = np.zeros((angles_dim,detectors_dim),'uint8')
+#mask =Mask.copy(order='c')
+#mask[:] = Mask[:]
+mask = np.ascontiguousarray(Mask, dtype=np.uint8);
+
+plt.figure(1)
+plt.subplot(121)
+plt.imshow(sino_cut_new,vmin=0.0, vmax=1)
+plt.title('Missing Data sinogram')
+plt.subplot(122)
+plt.imshow(mask)
+plt.title('Mask')
+plt.show()
+#%%
+print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
+print ("___Inpainting using linear diffusion (2D)__")
+print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
+
+## plot 
+fig = plt.figure(2)
+plt.suptitle('Performance of linear inpainting using the CPU')
+a=fig.add_subplot(1,2,1)
+a.set_title('Missing data sinogram')
+imgplot = plt.imshow(sino_cut_new,cmap="gray")
+
+# set parameters
+pars = {'algorithm' : NDF_INP, \
+        'input' : sino_cut_new,\
+        'maskData' : mask,\
+        'regularisation_parameter':5000,\
+        'edge_parameter':0,\
+        'number_of_iterations' :5000 ,\
+        'time_marching_parameter':0.000075,\
+        'penalty_type':0
+        }
+        
+start_time = timeit.default_timer()
+ndf_inp_linear = NDF_INP(pars['input'],
+              pars['maskData'],
+              pars['regularisation_parameter'],
+              pars['edge_parameter'], 
+              pars['number_of_iterations'],
+              pars['time_marching_parameter'], 
+              pars['penalty_type'])
+             
+rms = rmse(sino_full, ndf_inp_linear)
+pars['rmse'] = rms
+
+txtstr = printParametersToString(pars)
+txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
+print (txtstr)
+a=fig.add_subplot(1,2,2)
+
+# these are matplotlib.patch.Patch properties
+props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
+# place a text box in upper left in axes coords
+a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
+         verticalalignment='top', bbox=props)
+imgplot = plt.imshow(ndf_inp_linear, cmap="gray")
+plt.title('{}'.format('Linear diffusion inpainting results'))
+#%%
+print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
+print ("_Inpainting using nonlinear diffusion (2D)_")
+print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
+
+## plot 
+fig = plt.figure(3)
+plt.suptitle('Performance of nonlinear diffusion inpainting using the CPU')
+a=fig.add_subplot(1,2,1)
+a.set_title('Missing data sinogram')
+imgplot = plt.imshow(sino_cut_new,cmap="gray")
+
+# set parameters
+pars = {'algorithm' : NDF_INP, \
+        'input' : sino_cut_new,\
+        'maskData' : mask,\
+        'regularisation_parameter':80,\
+        'edge_parameter':0.00009,\
+        'number_of_iterations' :1500 ,\
+        'time_marching_parameter':0.000008,\
+        'penalty_type':1
+        }
+        
+start_time = timeit.default_timer()
+ndf_inp_nonlinear = NDF_INP(pars['input'],
+              pars['maskData'],
+              pars['regularisation_parameter'],
+              pars['edge_parameter'], 
+              pars['number_of_iterations'],
+              pars['time_marching_parameter'], 
+              pars['penalty_type'])
+             
+rms = rmse(sino_full, ndf_inp_nonlinear)
+pars['rmse'] = rms
+
+txtstr = printParametersToString(pars)
+txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
+print (txtstr)
+a=fig.add_subplot(1,2,2)
+
+# these are matplotlib.patch.Patch properties
+props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
+# place a text box in upper left in axes coords
+a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
+         verticalalignment='top', bbox=props)
+imgplot = plt.imshow(ndf_inp_nonlinear, cmap="gray")
+plt.title('{}'.format('Nonlinear diffusion inpainting results'))
+#%%
+print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
+print ("Inpainting using nonlocal vertical marching")
+print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
+
+## plot 
+fig = plt.figure(4)
+plt.suptitle('Performance of NVM inpainting using the CPU')
+a=fig.add_subplot(1,2,1)
+a.set_title('Missing data sinogram')
+imgplot = plt.imshow(sino_cut,cmap="gray")
+
+# set parameters
+pars = {'algorithm' : NVM_INP, \
+        'input' : sino_cut_new,\
+        'maskData' : mask,\
+        'SW_increment': 1,\
+        'number_of_iterations' : 150
+        }
+        
+start_time = timeit.default_timer()
+(nvm_inp, mask_upd) = NVM_INP(pars['input'],
+              pars['maskData'],
+              pars['SW_increment'],
+              pars['number_of_iterations'])
+             
+rms = rmse(sino_full, nvm_inp)
+pars['rmse'] = rms
+
+txtstr = printParametersToString(pars)
+txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
+print (txtstr)
+a=fig.add_subplot(1,2,2)
+
+# these are matplotlib.patch.Patch properties
+props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
+# place a text box in upper left in axes coords
+a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
+         verticalalignment='top', bbox=props)
+imgplot = plt.imshow(nvm_inp, cmap="gray")
+plt.title('{}'.format('Nonlocal Vertical Marching inpainting results'))
+#%%
diff --git a/Wrappers/Python/demos/demo_cpu_regularisers.py b/Wrappers/Python/demos/demo_cpu_regularisers.py
index 3567f91..986e3e9 100644
--- a/Wrappers/Python/demos/demo_cpu_regularisers.py
+++ b/Wrappers/Python/demos/demo_cpu_regularisers.py
@@ -44,13 +44,31 @@ u0 = Im + np.random.normal(loc = 0 ,
 u_ref = Im + np.random.normal(loc = 0 ,
                                   scale = 0.01 * Im , 
                                   size = np.shape(Im))
- 
+(N,M) = np.shape(u0)
 # map the u0 u0->u0>0
 # f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1)
 u0 = u0.astype('float32')
 u_ref = u_ref.astype('float32')
 
-
+# change dims to check that modules work with non-squared images
+"""
+M = M-100
+u_ref2 = np.zeros([N,M],dtype='float32')
+u_ref2[:,0:M] = u_ref[:,0:M]
+u_ref = u_ref2
+del u_ref2
+
+u02 = np.zeros([N,M],dtype='float32')
+u02[:,0:M] = u0[:,0:M]
+u0 = u02
+del u02
+
+Im2 = np.zeros([N,M],dtype='float32')
+Im2[:,0:M] = Im[:,0:M]
+Im = Im2
+del Im2
+"""
+#%%
 print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
 print ("_______________ROF-TV (2D)_________________")
 print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
@@ -288,7 +306,6 @@ a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
 imgplot = plt.imshow(fgp_dtv_cpu, cmap="gray")
 plt.title('{}'.format('CPU results'))
 
-
 print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
 print ("__________Total nuclear Variation__________")
 print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
@@ -301,9 +318,8 @@ a.set_title('Noisy Image')
 imgplot = plt.imshow(u0,cmap="gray")
 
 channelsNo = 5
-N = 512
-noisyVol = np.zeros((channelsNo,N,N),dtype='float32')
-idealVol = np.zeros((channelsNo,N,N),dtype='float32')
+noisyVol = np.zeros((channelsNo,N,M),dtype='float32')
+idealVol = np.zeros((channelsNo,N,M),dtype='float32')
 
 for i in range (channelsNo):
     noisyVol[i,:,:] = Im + np.random.normal(loc = 0 , scale = perc * Im , size = np.shape(Im))
@@ -344,25 +360,19 @@ plt.title('{}'.format('CPU results'))
 # Uncomment to test 3D regularisation performance 
 #%%
 """
-N = 512
 slices = 20
-
-filename = os.path.join(".." , ".." , ".." , "data" ,"lena_gray_512.tif")
-Im = plt.imread(filename)
-Im = np.asarray(Im, dtype='float32')
-
-Im = Im/255
 perc = 0.05
 
-noisyVol = np.zeros((slices,N,N),dtype='float32')
-noisyRef = np.zeros((slices,N,N),dtype='float32')
-idealVol = np.zeros((slices,N,N),dtype='float32')
+noisyVol = np.zeros((slices,N,M),dtype='float32')
+noisyRef = np.zeros((slices,N,M),dtype='float32')
+idealVol = np.zeros((slices,N,M),dtype='float32')
 
 for i in range (slices):
     noisyVol[i,:,:] = Im + np.random.normal(loc = 0 , scale = perc * Im , size = np.shape(Im))
     noisyRef[i,:,:] = Im + np.random.normal(loc = 0 , scale = 0.01 * Im , size = np.shape(Im))
     idealVol[i,:,:] = Im
 
+
 print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
 print ("_______________ROF-TV (3D)_________________")
 print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
@@ -403,6 +413,7 @@ a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
 imgplot = plt.imshow(rof_cpu3D[10,:,:], cmap="gray")
 plt.title('{}'.format('Recovered volume on the CPU using ROF-TV'))
 
+
 print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
 print ("_______________FGP-TV (3D)__________________")
 print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
diff --git a/Wrappers/Python/demos/demo_gpu_regularisers.py b/Wrappers/Python/demos/demo_gpu_regularisers.py
index b873700..f3ed50c 100644
--- a/Wrappers/Python/demos/demo_gpu_regularisers.py
+++ b/Wrappers/Python/demos/demo_gpu_regularisers.py
@@ -44,10 +44,28 @@ u0 = Im + np.random.normal(loc = 0 ,
 u_ref = Im + np.random.normal(loc = 0 ,
                                   scale = 0.01 * Im , 
                                   size = np.shape(Im))
+(N,M) = np.shape(u0)
 # map the u0 u0->u0>0
 # f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1)
 u0 = u0.astype('float32')
 u_ref = u_ref.astype('float32')
+"""
+M = M-100
+u_ref2 = np.zeros([N,M],dtype='float32')
+u_ref2[:,0:M] = u_ref[:,0:M]
+u_ref = u_ref2
+del u_ref2
+
+u02 = np.zeros([N,M],dtype='float32')
+u02[:,0:M] = u0[:,0:M]
+u0 = u02
+del u02
+
+Im2 = np.zeros([N,M],dtype='float32')
+Im2[:,0:M] = Im[:,0:M]
+Im = Im2
+del Im2
+"""
 
 print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
 print ("____________ROF-TV regulariser_____________")
diff --git a/Wrappers/Python/setup-regularisers.py.in b/Wrappers/Python/setup-regularisers.py.in
index b900efe..f55c6fe 100644
--- a/Wrappers/Python/setup-regularisers.py.in
+++ b/Wrappers/Python/setup-regularisers.py.in
@@ -34,6 +34,7 @@ extra_libraries = ['cilreg']
 
 extra_include_dirs += [os.path.join(".." , ".." , "Core"),
                        os.path.join(".." , ".." , "Core",  "regularisers_CPU"),
+                       os.path.join(".." , ".." , "Core",  "inpainters_CPU"),
                        os.path.join(".." , ".." , "Core",  "regularisers_GPU" , "TV_FGP" ) , 
                        os.path.join(".." , ".." , "Core",  "regularisers_GPU" , "TV_ROF" ) , 
                        os.path.join(".." , ".." , "Core",  "regularisers_GPU" , "TV_SB" ) ,
@@ -52,7 +53,7 @@ setup(
 	description='CCPi Core Imaging Library - Image regularisers',
 	version=cil_version,
     cmdclass = {'build_ext': build_ext},
-    ext_modules = [Extension("ccpi.filters.cpu_regularisers_cython",
+    ext_modules = [Extension("ccpi.filters.cpu_regularisers",
                              sources=[os.path.join("." , "src", "cpu_regularisers.pyx" ) ],
                              include_dirs=extra_include_dirs, 
 							 library_dirs=extra_library_dirs, 
diff --git a/Wrappers/Python/src/cpu_regularisers.pyx b/Wrappers/Python/src/cpu_regularisers.pyx
index 7ed8fa1..c934f1d 100644
--- a/Wrappers/Python/src/cpu_regularisers.pyx
+++ b/Wrappers/Python/src/cpu_regularisers.pyx
@@ -25,6 +25,8 @@ cdef extern float Diffusion_CPU_main(float *Input, float *Output, float lambdaPa
 cdef extern float TNV_CPU_main(float *Input, float *u, float lambdaPar, int maxIter, float tol, int dimX, int dimY, int dimZ);
 cdef extern float dTV_FGP_CPU_main(float *Input, float *InputRef, float *Output, float lambdaPar, int iterationsNumb, float epsil, float eta, int methodTV, int nonneg, int printM, int dimX, int dimY, int dimZ);
 
+cdef extern float Diffusion_Inpaint_CPU_main(float *Input, unsigned char *Mask, float *Output, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int penaltytype, int dimX, int dimY, int dimZ);
+cdef extern float NonlocalMarching_Inpaint_main(float *Input, unsigned char *M, float *Output, unsigned char *M_upd, int SW_increment, int iterationsNumb, int trigger, int dimX, int dimY, int dimZ);
 #****************************************************************#
 #********************** Total-variation ROF *********************#
 #****************************************************************#
@@ -46,7 +48,7 @@ def TV_ROF_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
             np.zeros([dims[0],dims[1]], dtype='float32')
                    
     # Run ROF iterations for 2D data 
-    TV_ROF_CPU_main(&inputData[0,0], &outputData[0,0], regularisation_parameter, iterationsNumb, marching_step_parameter, dims[0], dims[1], 1)
+    TV_ROF_CPU_main(&inputData[0,0], &outputData[0,0], regularisation_parameter, iterationsNumb, marching_step_parameter, dims[1], dims[0], 1)
     
     return outputData
             
@@ -99,7 +101,7 @@ def TV_FGP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
                        methodTV,
                        nonneg,
                        printM,
-                       dims[0], dims[1], 1)
+                       dims[1],dims[0],1)
     
     return outputData        
             
@@ -158,7 +160,7 @@ def TV_SB_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
                        tolerance_param,
                        methodTV,
                        printM,
-                       dims[0], dims[1], 1)
+                       dims[1],dims[0],1)
     
     return outputData        
             
@@ -219,7 +221,7 @@ def dTV_FGP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
                        methodTV,                       
                        nonneg,
                        printM,
-                       dims[0], dims[1], 1)
+                       dims[1], dims[0], 1)
     
     return outputData        
             
@@ -298,7 +300,7 @@ def NDF_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
             np.zeros([dims[0],dims[1]], dtype='float32')   
     
     # Run Nonlinear Diffusion iterations for 2D data 
-    Diffusion_CPU_main(&inputData[0,0], &outputData[0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type, dims[0], dims[1], 1)    
+    Diffusion_CPU_main(&inputData[0,0], &outputData[0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type, dims[1], dims[0], 1)
     return outputData
             
 def NDF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
@@ -319,3 +321,81 @@ def NDF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
     Diffusion_CPU_main(&inputData[0,0,0], &outputData[0,0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type, dims[2], dims[1], dims[0])
 
     return outputData
+
+#*********************Inpainting WITH****************************#
+#***************Nonlinear (Isotropic) Diffusion******************#
+#****************************************************************#
+def NDF_INPAINT_CPU(inputData, maskData, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type):
+    if inputData.ndim == 2:
+        return NDF_INP_2D(inputData, maskData, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type)
+    elif inputData.ndim == 3:
+        return NDF_INP_3D(inputData, maskData, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type)
+
+def NDF_INP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
+                     np.ndarray[np.uint8_t, ndim=2, mode="c"] maskData,
+                     float regularisation_parameter,
+                     float edge_parameter,
+                     int iterationsNumb,
+                     float time_marching_parameter,
+                     int penalty_type):
+
+    cdef long dims[2]
+    dims[0] = inputData.shape[0]
+    dims[1] = inputData.shape[1]
+
+
+    cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
+            np.zeros([dims[0],dims[1]], dtype='float32')
+    
+    # Run Inpaiting by Diffusion iterations for 2D data 
+    Diffusion_Inpaint_CPU_main(&inputData[0,0], &maskData[0,0], &outputData[0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type, dims[1], dims[0], 1)
+    return outputData
+            
+def NDF_INP_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
+                     np.ndarray[np.uint8_t, ndim=3, mode="c"] maskData,
+                     float regularisation_parameter,
+                     float edge_parameter,
+                     int iterationsNumb,
+                     float time_marching_parameter,
+                     int penalty_type):
+    cdef long dims[3]
+    dims[0] = inputData.shape[0]
+    dims[1] = inputData.shape[1]
+    dims[2] = inputData.shape[2]
+    
+    cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
+            np.zeros([dims[0],dims[1],dims[2]], dtype='float32')
+    
+    # Run Inpaiting by Diffusion iterations for 3D data 
+    Diffusion_Inpaint_CPU_main(&inputData[0,0,0], &maskData[0,0,0], &outputData[0,0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type, dims[2], dims[1], dims[0])
+
+    return outputData
+#*********************Inpainting WITH****************************#
+#***************Nonlocal Vertical Marching method****************#
+#****************************************************************#
+def NVM_INPAINT_CPU(inputData, maskData, SW_increment, iterationsNumb):
+    if inputData.ndim == 2:
+        return NVM_INP_2D(inputData, maskData, SW_increment, iterationsNumb)
+    elif inputData.ndim == 3:
+        return 
+
+def NVM_INP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
+               np.ndarray[np.uint8_t, ndim=2, mode="c"] maskData,
+                     int SW_increment,
+                     int iterationsNumb):
+    cdef long dims[2]
+    dims[0] = inputData.shape[0]
+    dims[1] = inputData.shape[1]
+    
+    cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
+            np.zeros([dims[0],dims[1]], dtype='float32')   
+    
+    cdef np.ndarray[np.uint8_t, ndim=2, mode="c"] maskData_upd = \
+            np.zeros([dims[0],dims[1]], dtype='uint8')
+    
+    # Run Inpaiting by Nonlocal vertical marching method for 2D data 
+    NonlocalMarching_Inpaint_main(&inputData[0,0], &maskData[0,0], &outputData[0,0], 
+                                  &maskData_upd[0,0],
+                                  SW_increment, iterationsNumb, 1, dims[1], dims[0], 1)
+    
+    return (outputData, maskData_upd)
diff --git a/Wrappers/Python/src/gpu_regularisers.pyx b/Wrappers/Python/src/gpu_regularisers.pyx
index b0775054..7eab5d5 100644
--- a/Wrappers/Python/src/gpu_regularisers.pyx
+++ b/Wrappers/Python/src/gpu_regularisers.pyx
@@ -157,7 +157,7 @@ def ROFTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
                        regularisation_parameter,
                        iterations , 
                        time_marching_parameter, 
-                       dims[0], dims[1], 1);   
+                       dims[1], dims[0], 1);   
      
     return outputData
     
@@ -210,7 +210,7 @@ def FGPTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
                        methodTV,
                        nonneg,
                        printM,
-                       dims[0], dims[1], 1);   
+                       dims[1], dims[0], 1);   
      
     return outputData
     
@@ -266,7 +266,7 @@ def SBTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
                        tolerance_param,
                        methodTV,
                        printM,
-                       dims[0], dims[1], 1);   
+                       dims[1], dims[0], 1);   
      
     return outputData
     
@@ -325,7 +325,7 @@ def FGPdTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
                        methodTV,
                        nonneg,
                        printM,
-                       dims[0], dims[1], 1);   
+                       dims[1], dims[0], 1);   
      
     return outputData
     
@@ -381,7 +381,7 @@ def NDF_GPU_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
     
     # Run Nonlinear Diffusion iterations for 2D data 
     # Running CUDA code here  
-    NonlDiff_GPU_main(&inputData[0,0], &outputData[0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type, dims[0], dims[1], 1)    
+    NonlDiff_GPU_main(&inputData[0,0], &outputData[0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type, dims[1], dims[0], 1)
     return outputData
             
 def NDF_GPU_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
diff --git a/data/SinoInpaint.mat b/data/SinoInpaint.mat
new file mode 100644
index 0000000..d748fb4
--- /dev/null
+++ b/data/SinoInpaint.mat