From 09eb48ffbb4ad699e2eefd25678e10dc59d6a177 Mon Sep 17 00:00:00 2001
From: Daniil Kazantsev <dkazanc@hotmail.com>
Date: Tue, 1 May 2018 09:44:07 +0100
Subject: new inpainters

---
 Core/CMakeLists.txt                                |   4 +-
 Core/inpainters_CPU/Diffusion_Inpaint_core.c       | 317 +++++++++++++++++++++
 Core/inpainters_CPU/Diffusion_Inpaint_core.h       |  60 ++++
 .../inpainters_CPU/NonlocalMarching_Inpaint_core.c | 158 ++++++++++
 .../inpainters_CPU/NonlocalMarching_Inpaint_core.h |  54 ++++
 Core/regularisers_CPU/utils.c                      |  33 ++-
 Core/regularisers_CPU/utils.h                      |   5 +-
 Readme.md                                          |  27 +-
 Wrappers/Matlab/demos/demoMatlab_3Ddenoise.m       |   9 +-
 Wrappers/Matlab/demos/demoMatlab_denoise.m         |   9 +-
 Wrappers/Matlab/demos/demoMatlab_inpaint.m         |  35 +++
 Wrappers/Matlab/mex_compile/compileCPU_mex.m       |  16 +-
 .../mex_compile/regularisers_CPU/NonlDiff_Inp.c    | 101 +++++++
 .../regularisers_CPU/NonlocalMarching_Inpaint.c    |  82 ++++++
 Wrappers/Python/ccpi/filters/regularisers.py       |   7 +-
 Wrappers/Python/setup-regularisers.py.in           |   3 +-
 Wrappers/Python/src/cpu_regularisers.pyx           |  50 ++++
 data/SinoInpaint.mat                               | Bin 0 -> 3303762 bytes
 18 files changed, 937 insertions(+), 33 deletions(-)
 create mode 100644 Core/inpainters_CPU/Diffusion_Inpaint_core.c
 create mode 100644 Core/inpainters_CPU/Diffusion_Inpaint_core.h
 create mode 100644 Core/inpainters_CPU/NonlocalMarching_Inpaint_core.c
 create mode 100644 Core/inpainters_CPU/NonlocalMarching_Inpaint_core.h
 create mode 100644 Wrappers/Matlab/demos/demoMatlab_inpaint.m
 create mode 100644 Wrappers/Matlab/mex_compile/regularisers_CPU/NonlDiff_Inp.c
 create mode 100644 Wrappers/Matlab/mex_compile/regularisers_CPU/NonlocalMarching_Inpaint.c
 create mode 100644 data/SinoInpaint.mat

diff --git a/Core/CMakeLists.txt b/Core/CMakeLists.txt
index 61986dc..322732b 100644
--- a/Core/CMakeLists.txt
+++ b/Core/CMakeLists.txt
@@ -85,7 +85,7 @@ message("Adding regularisers as a shared library")
 add_library(cilreg SHARED
 	    ${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/FGP_TV_core.c
 	    ${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/SB_TV_core.c
-	    ${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/Diffusion_core.c
+	    ${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/Diffusion_core.c	    
 	    #${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/LLT_model_core.c
 	    #${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/PatchBased_Regul_core.c
 	    #${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/TGV_PD_core.c
@@ -93,6 +93,8 @@ 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 
diff --git a/Core/inpainters_CPU/Diffusion_Inpaint_core.c b/Core/inpainters_CPU/Diffusion_Inpaint_core.c
new file mode 100644
index 0000000..16e87de
--- /dev/null
+++ b/Core/inpainters_CPU/Diffusion_Inpaint_core.c
@@ -0,0 +1,317 @@
+/*
+ * 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 signNDFc(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;
+    float sigmaPar2;
+    sigmaPar2 = sigmaPar/sqrt(2.0f);
+    
+    /* copy into output */
+    copyIm(Input, Output, dimX, dimY, dimZ);
+    
+    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 =  signNDFc(e1);
+            else e1 = e1/sigmaPar;
+            
+            if (fabs(w1) > sigmaPar) w1 =  signNDFc(w1);
+            else w1 = w1/sigmaPar;
+            
+            if (fabs(n1) > sigmaPar) n1 =  signNDFc(n1);
+            else n1 = n1/sigmaPar;
+            
+            if (fabs(s1) > sigmaPar) s1 =  signNDFc(s1);
+            else s1 = s1/sigmaPar;
+            }
+            else if (penaltytype == 2) {
+            /* Perona-Malik */
+            e1 = (e1)/(1.0f + powf((e1/sigmaPar),2));
+            w1 = (w1)/(1.0f + powf((w1/sigmaPar),2));
+            n1 = (n1)/(1.0f + powf((n1/sigmaPar),2));
+            s1 = (s1)/(1.0f + powf((s1/sigmaPar),2));
+            }
+            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 =  signNDFc(e1);
+            else e1 = e1/sigmaPar;
+            
+            if (fabs(w1) > sigmaPar) w1 =  signNDFc(w1);
+            else w1 = w1/sigmaPar;
+            
+            if (fabs(n1) > sigmaPar) n1 =  signNDFc(n1);
+            else n1 = n1/sigmaPar;
+            
+            if (fabs(s1) > sigmaPar) s1 =  signNDFc(s1);
+            else s1 = s1/sigmaPar;
+            
+            if (fabs(u1) > sigmaPar) u1 =  signNDFc(u1);
+            else u1 = u1/sigmaPar;
+            
+            if (fabs(d1) > sigmaPar) d1 =  signNDFc(d1);
+            else d1 = d1/sigmaPar;            
+            }
+            else if (penaltytype == 2) {
+            /* Perona-Malik */
+            e1 = (e1)/(1.0f + powf((e1/sigmaPar),2));
+            w1 = (w1)/(1.0f + powf((w1/sigmaPar),2));
+            n1 = (n1)/(1.0f + powf((n1/sigmaPar),2));
+            s1 = (s1)/(1.0f + powf((s1/sigmaPar),2));
+            u1 = (u1)/(1.0f + powf((u1/sigmaPar),2));
+            d1 = (d1)/(1.0f + powf((d1/sigmaPar),2));
+            }
+            else if (penaltytype == 3) {
+            /* Tukey Biweight */
+            if (fabs(e1) <= sigmaPar) e1 =  e1*powf((1.0f - powf((e1/sigmaPar),2)), 2);
+            else e1 = 0.0f;
+            if (fabs(w1) <= sigmaPar) w1 =  w1*powf((1.0f - powf((w1/sigmaPar),2)), 2);
+            else w1 = 0.0f;
+            if (fabs(n1) <= sigmaPar) n1 =  n1*powf((1.0f - powf((n1/sigmaPar),2)), 2);
+            else n1 = 0.0f;
+            if (fabs(s1) <= sigmaPar) s1 =  s1*powf((1.0f - powf((s1/sigmaPar),2)), 2);
+            else s1 = 0.0f;
+            if (fabs(u1) <= sigmaPar) u1 =  u1*powf((1.0f - powf((u1/sigmaPar),2)), 2);
+            else u1 = 0.0f;
+            if (fabs(d1) <= sigmaPar) d1 =  d1*powf((1.0f - powf((d1/sigmaPar),2)), 2);
+            else d1 = 0.0f;
+            }
+            else {
+				printf("%s \n", "No penalty function selected! Use 1,2 or 3.");
+				break;
+				}
+
+                Output[index] += tau*(lambdaPar*(e1 + w1 + n1 + s1 + u1 + d1) - (Output[index] - Input[index]));  
+		}
+		}}}
+	return *Output;
+}
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..d6cdf93
--- /dev/null
+++ b/Core/inpainters_CPU/NonlocalMarching_Inpaint_core.c
@@ -0,0 +1,158 @@
+/*
+ * 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 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++;
+        }
+    }
+    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++;
+                }
+            }
+            /* find a point in the mask to inpaint */
+#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 (j-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;                        
+                    }
+                }
+            }
+            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..199820c
--- /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 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/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..5e4ab1f
--- /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, dimX, dimY, dimZ);
+}
\ No newline at end of file
diff --git a/Wrappers/Python/ccpi/filters/regularisers.py b/Wrappers/Python/ccpi/filters/regularisers.py
index eec8c4d..e62c020 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
 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,8 @@ 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, iterationsNumb, 
+        time_marching_parameter, penalty_type)
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..3625106 100644
--- a/Wrappers/Python/src/cpu_regularisers.pyx
+++ b/Wrappers/Python/src/cpu_regularisers.pyx
@@ -25,6 +25,9 @@ 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 dimX, int dimY, int dimZ);
+
 #****************************************************************#
 #********************** Total-variation ROF *********************#
 #****************************************************************#
@@ -319,3 +322,50 @@ 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[0], dims[1], 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
diff --git a/data/SinoInpaint.mat b/data/SinoInpaint.mat
new file mode 100644
index 0000000..d2fbe12
Binary files /dev/null and b/data/SinoInpaint.mat differ
-- 
cgit v1.2.3


From fa47bdc29ba4178254531174c02f790a9d10a187 Mon Sep 17 00:00:00 2001
From: Daniil Kazantsev <dkazanc@hotmail.com>
Date: Tue, 1 May 2018 10:03:16 +0100
Subject: make updates

---
 Core/CMakeLists.txt            | 19 +++----------------
 Wrappers/Python/CMakeLists.txt |  5 -----
 2 files changed, 3 insertions(+), 21 deletions(-)

diff --git a/Core/CMakeLists.txt b/Core/CMakeLists.txt
index 322732b..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)
@@ -85,7 +76,7 @@ message("Adding regularisers as a shared library")
 add_library(cilreg SHARED
 	    ${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/FGP_TV_core.c
 	    ${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/SB_TV_core.c
-	    ${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/Diffusion_core.c	    
+	    ${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/Diffusion_core.c
 	    #${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/LLT_model_core.c
 	    #${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/PatchBased_Regul_core.c
 	    #${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/TGV_PD_core.c
@@ -100,7 +91,8 @@ 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
@@ -155,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/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)
-- 
cgit v1.2.3


From c04d2d000abb5c4d98b11f60dbaadaff3e2b3ff8 Mon Sep 17 00:00:00 2001
From: Daniil Kazantsev <dkazanc@hotmail.com>
Date: Tue, 1 May 2018 11:05:12 +0100
Subject: bug in DiffCore fixed

---
 Core/inpainters_CPU/Diffusion_Inpaint_core.c | 26 ++++++++++++--------------
 1 file changed, 12 insertions(+), 14 deletions(-)

diff --git a/Core/inpainters_CPU/Diffusion_Inpaint_core.c b/Core/inpainters_CPU/Diffusion_Inpaint_core.c
index 16e87de..e6972ab 100644
--- a/Core/inpainters_CPU/Diffusion_Inpaint_core.c
+++ b/Core/inpainters_CPU/Diffusion_Inpaint_core.c
@@ -21,7 +21,7 @@
 #include "utils.h"
 
 /*sign function*/
-int signNDFc(float x) {
+int signNDF_inc(float x) {
     return (x > 0) - (x < 0);
 }
 
@@ -70,8 +70,6 @@ float Diffusion_Inpaint_CPU_main(float *Input, unsigned char *Mask, float *Outpu
 	}
     return *Output;
 }
-
-
 /********************************************************************/
 /***************************2D Functions*****************************/
 /********************************************************************/
@@ -141,16 +139,16 @@ float NonLinearDiff_Inp_2D(float *Input, unsigned char *Mask, float *Output, flo
 
             if (penaltytype == 1){
             /* Huber penalty */
-            if (fabs(e1) > sigmaPar) e1 =  signNDFc(e1);
+            if (fabs(e1) > sigmaPar) e1 =  signNDF_inc(e1);
             else e1 = e1/sigmaPar;
             
-            if (fabs(w1) > sigmaPar) w1 =  signNDFc(w1);
+            if (fabs(w1) > sigmaPar) w1 =  signNDF_inc(w1);
             else w1 = w1/sigmaPar;
             
-            if (fabs(n1) > sigmaPar) n1 =  signNDFc(n1);
+            if (fabs(n1) > sigmaPar) n1 =  signNDF_inc(n1);
             else n1 = n1/sigmaPar;
             
-            if (fabs(s1) > sigmaPar) s1 =  signNDFc(s1);
+            if (fabs(s1) > sigmaPar) s1 =  signNDF_inc(s1);
             else s1 = s1/sigmaPar;
             }
             else if (penaltytype == 2) {
@@ -263,23 +261,23 @@ for(k=0; k<dimZ; k++) {
                 
              if (penaltytype == 1){
             /* Huber penalty */
-            if (fabs(e1) > sigmaPar) e1 =  signNDFc(e1);
+            if (fabs(e1) > sigmaPar) e1 =  signNDF_inc(e1);
             else e1 = e1/sigmaPar;
             
-            if (fabs(w1) > sigmaPar) w1 =  signNDFc(w1);
+            if (fabs(w1) > sigmaPar) w1 =  signNDF_inc(w1);
             else w1 = w1/sigmaPar;
             
-            if (fabs(n1) > sigmaPar) n1 =  signNDFc(n1);
+            if (fabs(n1) > sigmaPar) n1 =  signNDF_inc(n1);
             else n1 = n1/sigmaPar;
             
-            if (fabs(s1) > sigmaPar) s1 =  signNDFc(s1);
+            if (fabs(s1) > sigmaPar) s1 =  signNDF_inc(s1);
             else s1 = s1/sigmaPar;
             
-            if (fabs(u1) > sigmaPar) u1 =  signNDFc(u1);
+            if (fabs(u1) > sigmaPar) u1 =  signNDF_inc(u1);
             else u1 = u1/sigmaPar;
             
-            if (fabs(d1) > sigmaPar) d1 =  signNDFc(d1);
-            else d1 = d1/sigmaPar;            
+            if (fabs(d1) > sigmaPar) d1 =  signNDF_inc(d1);
+            else d1 = d1/sigmaPar;
             }
             else if (penaltytype == 2) {
             /* Perona-Malik */
-- 
cgit v1.2.3


From 5ef87da22a31868fd88c7f0ab4c2201e816e92ed Mon Sep 17 00:00:00 2001
From: Daniil Kazantsev <dkazanc@hotmail.com>
Date: Tue, 1 May 2018 12:07:30 +0100
Subject: inpaint demo

---
 Wrappers/Python/demos/demo_cpu_inpainters.py | 143 +++++++++++++++++++++++++++
 1 file changed, 143 insertions(+)
 create mode 100644 Wrappers/Python/demos/demo_cpu_inpainters.py

diff --git a/Wrappers/Python/demos/demo_cpu_inpainters.py b/Wrappers/Python/demos/demo_cpu_inpainters.py
new file mode 100644
index 0000000..a022bc8
--- /dev/null
+++ b/Wrappers/Python/demos/demo_cpu_inpainters.py
@@ -0,0 +1,143 @@
+#!/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
+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 == 'refdata':
+                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_no_cut = sino.get('Sinogram')
+Mask = sino.get('Mask')
+[angles_dim,detectors_dim] = sino_no_cut.shape
+sinogram = sino_no_cut/np.max(sino_no_cut)
+#apply mask to sinogram
+sino_cut = sinogram*(1-Mask)
+
+plt.figure(1)
+plt.subplot(121)
+plt.imshow(sino_cut,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,cmap="gray")
+
+# set parameters
+pars = {'algorithm' : NDF_INP, \
+        'input' : sino_cut,\
+        'maskData' : Mask,\
+        'regularisation_parameter':6000,\
+        'edge_parameter':0.0,\
+        'number_of_iterations' :5000 ,\
+        'time_marching_parameter':0.000075,\
+        'penalty_type':1
+        }
+        
+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(sinogram, 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,cmap="gray")
+
+# set parameters
+pars = {'algorithm' : NDF_INP, \
+        'input' : sino_cut,\
+        '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(sinogram, 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'))
+#%%
\ No newline at end of file
-- 
cgit v1.2.3


From 033c2030a05c7aa4c832e7a5e9fd13346d05e33d Mon Sep 17 00:00:00 2001
From: algol <dkazanc@hotmail.com>
Date: Tue, 1 May 2018 14:48:45 +0100
Subject: some correction

---
 Core/inpainters_CPU/Diffusion_Inpaint_core.c |  9 ++++++-
 Wrappers/Python/ccpi/filters/regularisers.py |  3 +--
 Wrappers/Python/demos/demo_cpu_inpainters.py | 40 ++++++++++++++++------------
 3 files changed, 32 insertions(+), 20 deletions(-)

diff --git a/Core/inpainters_CPU/Diffusion_Inpaint_core.c b/Core/inpainters_CPU/Diffusion_Inpaint_core.c
index e6972ab..56e0421 100644
--- a/Core/inpainters_CPU/Diffusion_Inpaint_core.c
+++ b/Core/inpainters_CPU/Diffusion_Inpaint_core.c
@@ -47,13 +47,19 @@ int signNDF_inc(float x) {
 
 float Diffusion_Inpaint_CPU_main(float *Input, unsigned char *Mask, float *Output, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int penaltytype, int dimX, int dimY, int dimZ)
 {
-    int i;
+    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++) {
@@ -68,6 +74,7 @@ float Diffusion_Inpaint_CPU_main(float *Input, unsigned char *Mask, float *Outpu
             else NonLinearDiff_Inp_3D(Input, Mask, Output, lambdaPar, sigmaPar2, tau, penaltytype, dimX, dimY, dimZ);
 		}
 	}
+	}
     return *Output;
 }
 /********************************************************************/
diff --git a/Wrappers/Python/ccpi/filters/regularisers.py b/Wrappers/Python/ccpi/filters/regularisers.py
index e62c020..8120f72 100644
--- a/Wrappers/Python/ccpi/filters/regularisers.py
+++ b/Wrappers/Python/ccpi/filters/regularisers.py
@@ -113,5 +113,4 @@ def NDF(inputData, regularisation_parameter, edge_parameter, iterations,
 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, iterationsNumb, 
-        time_marching_parameter, penalty_type)
+        edge_parameter, iterations, time_marching_parameter, penalty_type)
diff --git a/Wrappers/Python/demos/demo_cpu_inpainters.py b/Wrappers/Python/demos/demo_cpu_inpainters.py
index a022bc8..b067b11 100644
--- a/Wrappers/Python/demos/demo_cpu_inpainters.py
+++ b/Wrappers/Python/demos/demo_cpu_inpainters.py
@@ -20,30 +20,36 @@ def printParametersToString(pars):
                 txt += "{0} = {1}".format(key, value.__name__)
             elif key == 'input':
                 txt += "{0} = {1}".format(key, np.shape(value))
-            elif key == 'refdata':
+            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_no_cut = sino.get('Sinogram')
+sino_full = sino.get('Sinogram')
 Mask = sino.get('Mask')
-[angles_dim,detectors_dim] = sino_no_cut.shape
-sinogram = sino_no_cut/np.max(sino_no_cut)
+[angles_dim,detectors_dim] = sino_full.shape
+sino_full = sino_full/np.max(sino_full)
 #apply mask to sinogram
-sino_cut = sinogram*(1-Mask)
+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[:]
+mask = np.zeros((angles_dim,detectors_dim),'uint8')
+#mask =Mask.copy(order='c')
+mask[:] = Mask[:]
 
 plt.figure(1)
 plt.subplot(121)
-plt.imshow(sino_cut,vmin=0.0, vmax=1)
+plt.imshow(sino_cut_new,vmin=0.0, vmax=1)
 plt.title('Missing Data sinogram')
 plt.subplot(122)
-plt.imshow(Mask)
+plt.imshow(mask)
 plt.title('Mask')
 plt.show()
 #%%
@@ -56,15 +62,15 @@ 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,cmap="gray")
+imgplot = plt.imshow(sino_cut_new,cmap="gray")
 
 # set parameters
 pars = {'algorithm' : NDF_INP, \
-        'input' : sino_cut,\
-        'maskData' : Mask,\
-        'regularisation_parameter':6000,\
+        'input' : sino_cut_new,\
+        'maskData' : mask,\
+        'regularisation_parameter':1000,\
         'edge_parameter':0.0,\
-        'number_of_iterations' :5000 ,\
+        'number_of_iterations' :1000 ,\
         'time_marching_parameter':0.000075,\
         'penalty_type':1
         }
@@ -78,7 +84,7 @@ ndf_inp_linear = NDF_INP(pars['input'],
               pars['time_marching_parameter'], 
               pars['penalty_type'])
              
-rms = rmse(sinogram, ndf_inp_linear)
+rms = rmse(sino_full, ndf_inp_linear)
 pars['rmse'] = rms
 
 txtstr = printParametersToString(pars)
@@ -107,8 +113,8 @@ imgplot = plt.imshow(sino_cut,cmap="gray")
 
 # set parameters
 pars = {'algorithm' : NDF_INP, \
-        'input' : sino_cut,\
-        'maskData' : Mask,\
+        'input' : sino_cut_new,\
+        'maskData' : mask,\
         'regularisation_parameter':80,\
         'edge_parameter':0.00009,\
         'number_of_iterations' :1500 ,\
@@ -125,7 +131,7 @@ ndf_inp_nonlinear = NDF_INP(pars['input'],
               pars['time_marching_parameter'], 
               pars['penalty_type'])
              
-rms = rmse(sinogram, ndf_inp_nonlinear)
+rms = rmse(sino_full, ndf_inp_nonlinear)
 pars['rmse'] = rms
 
 txtstr = printParametersToString(pars)
-- 
cgit v1.2.3


From 42a10faa06bd56bff3f0f1804ddcdf1a3e1283cd Mon Sep 17 00:00:00 2001
From: Daniil Kazantsev <dkazanc@hotmail.com>
Date: Tue, 1 May 2018 15:16:49 +0100
Subject: inpaint NVM added

---
 Wrappers/Python/ccpi/filters/regularisers.py |  5 +++-
 Wrappers/Python/demos/demo_cpu_inpainters.py | 45 ++++++++++++++++++++++++++--
 Wrappers/Python/src/cpu_regularisers.pyx     | 31 ++++++++++++++++++-
 3 files changed, 77 insertions(+), 4 deletions(-)

diff --git a/Wrappers/Python/ccpi/filters/regularisers.py b/Wrappers/Python/ccpi/filters/regularisers.py
index 8120f72..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 import TV_ROF_CPU, TV_FGP_CPU, TV_SB_CPU, dTV_FGP_CPU, TNV_CPU, NDF_CPU, NDF_INPAINT_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,
@@ -114,3 +114,6 @@ def NDF_INP(inputData, maskData, regularisation_parameter, edge_parameter, itera
                      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
index b067b11..ab7ed2f 100644
--- a/Wrappers/Python/demos/demo_cpu_inpainters.py
+++ b/Wrappers/Python/demos/demo_cpu_inpainters.py
@@ -10,7 +10,7 @@ import numpy as np
 import os
 import timeit
 from scipy import io
-from ccpi.filters.regularisers import NDF_INP
+from ccpi.filters.regularisers import NDF_INP, NVM_INP
 from qualitymetrics import rmse
 ###############################################################################
 def printParametersToString(pars):
@@ -146,4 +146,45 @@ 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'))
-#%%
\ No newline at end of file
+#%%
+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' :20
+        }
+        
+start_time = timeit.default_timer()
+nvm_inp = 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/src/cpu_regularisers.pyx b/Wrappers/Python/src/cpu_regularisers.pyx
index 3625106..19dd707 100644
--- a/Wrappers/Python/src/cpu_regularisers.pyx
+++ b/Wrappers/Python/src/cpu_regularisers.pyx
@@ -26,7 +26,7 @@ cdef extern float TNV_CPU_main(float *Input, float *u, float lambdaPar, int maxI
 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 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 dimX, int dimY, int dimZ);
 
 #****************************************************************#
 #********************** Total-variation ROF *********************#
@@ -368,4 +368,33 @@ def NDF_INP_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
     # 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, iterations):
+    if inputData.ndim == 2:
+        return NVM_INP_2D(inputData, maskData, SW_increment, iterations)
+    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,  
+    dims[0], dims[1], 1)    
+    
     return outputData
-- 
cgit v1.2.3


From c4f50db4f5b318aad785ae577908d37fe05f53d2 Mon Sep 17 00:00:00 2001
From: algol <dkazanc@hotmail.com>
Date: Tue, 1 May 2018 15:30:28 +0100
Subject: some updates in demo

---
 Wrappers/Python/demos/demo_cpu_inpainters.py | 2 +-
 Wrappers/Python/src/cpu_regularisers.pyx     | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/Wrappers/Python/demos/demo_cpu_inpainters.py b/Wrappers/Python/demos/demo_cpu_inpainters.py
index ab7ed2f..9a677c4 100644
--- a/Wrappers/Python/demos/demo_cpu_inpainters.py
+++ b/Wrappers/Python/demos/demo_cpu_inpainters.py
@@ -167,7 +167,7 @@ pars = {'algorithm' : NVM_INP, \
         }
         
 start_time = timeit.default_timer()
-nvm_inp = NVM_INP(pars['input'],
+(nvm_inp, mask_upd) = NVM_INP(pars['input'],
               pars['maskData'],
               pars['SW_increment'],
               pars['number_of_iterations'])
diff --git a/Wrappers/Python/src/cpu_regularisers.pyx b/Wrappers/Python/src/cpu_regularisers.pyx
index 19dd707..52befd7 100644
--- a/Wrappers/Python/src/cpu_regularisers.pyx
+++ b/Wrappers/Python/src/cpu_regularisers.pyx
@@ -397,4 +397,4 @@ def NVM_INP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
     SW_increment, iterationsNumb,  
     dims[0], dims[1], 1)    
     
-    return outputData
+    return (outputData, maskData_upd)
-- 
cgit v1.2.3


From 0fe584fdf3a3ce0b1c66bb7d25a27fb8f35daea6 Mon Sep 17 00:00:00 2001
From: Daniil Kazantsev <dkazanc@hotmail.com>
Date: Tue, 1 May 2018 16:14:50 +0100
Subject: data file change

---
 data/SinoInpaint.mat | Bin 3303762 -> 3584100 bytes
 1 file changed, 0 insertions(+), 0 deletions(-)

diff --git a/data/SinoInpaint.mat b/data/SinoInpaint.mat
index d2fbe12..58174a9 100644
Binary files a/data/SinoInpaint.mat and b/data/SinoInpaint.mat differ
-- 
cgit v1.2.3


From 73965b6b80c49a2867d54e4a42f3069fe35d9cc6 Mon Sep 17 00:00:00 2001
From: algol <dkazanc@hotmail.com>
Date: Wed, 2 May 2018 09:47:58 +0100
Subject: corrections to dimens issues

---
 Wrappers/Python/demos/demo_cpu_inpainters.py   | 18 ++++++++++--------
 Wrappers/Python/demos/demo_cpu_regularisers.py | 17 +++++++++++++++++
 Wrappers/Python/demos/demo_gpu_regularisers.py | 16 ++++++++++++++++
 Wrappers/Python/src/cpu_regularisers.pyx       | 10 ++++++----
 4 files changed, 49 insertions(+), 12 deletions(-)

diff --git a/Wrappers/Python/demos/demo_cpu_inpainters.py b/Wrappers/Python/demos/demo_cpu_inpainters.py
index 9a677c4..348d235 100644
--- a/Wrappers/Python/demos/demo_cpu_inpainters.py
+++ b/Wrappers/Python/demos/demo_cpu_inpainters.py
@@ -37,12 +37,14 @@ Mask = sino.get('Mask')
 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 = np.zeros((angles_dim,detectors_dim),'float32')
 #sino_cut_new = sino_cut.copy(order='c')
-sino_cut_new[:] = sino_cut[:]
-mask = np.zeros((angles_dim,detectors_dim),'uint8')
+#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[:] = Mask[:]
+mask = np.ascontiguousarray(Mask, dtype=np.uint8);
 
 plt.figure(1)
 plt.subplot(121)
@@ -68,11 +70,11 @@ imgplot = plt.imshow(sino_cut_new,cmap="gray")
 pars = {'algorithm' : NDF_INP, \
         'input' : sino_cut_new,\
         'maskData' : mask,\
-        'regularisation_parameter':1000,\
-        'edge_parameter':0.0,\
+        'regularisation_parameter':5000,\
+        'edge_parameter':0,\
         'number_of_iterations' :1000 ,\
         'time_marching_parameter':0.000075,\
-        'penalty_type':1
+        'penalty_type':0
         }
         
 start_time = timeit.default_timer()
@@ -163,7 +165,7 @@ pars = {'algorithm' : NVM_INP, \
         'input' : sino_cut_new,\
         'maskData' : mask,\
         'SW_increment': 1,\
-        'number_of_iterations' :20
+        'number_of_iterations' :0
         }
         
 start_time = timeit.default_timer()
diff --git a/Wrappers/Python/demos/demo_cpu_regularisers.py b/Wrappers/Python/demos/demo_cpu_regularisers.py
index 3567f91..f803870 100644
--- a/Wrappers/Python/demos/demo_cpu_regularisers.py
+++ b/Wrappers/Python/demos/demo_cpu_regularisers.py
@@ -50,7 +50,24 @@ u_ref = Im + np.random.normal(loc = 0 ,
 u0 = u0.astype('float32')
 u_ref = u_ref.astype('float32')
 
+# change dims to check that modules work with non-squared images
+(N,M) = np.shape(u0)
+u_ref2 = np.zeros([N,M-100],dtype='float32')
+u_ref2[:,0:M-100] = u_ref[:,0:M-100]
+u_ref = u_ref2
+del u_ref2
+
+u02 = np.zeros([N,M-100],dtype='float32')
+u02[:,0:M-100] = u0[:,0:M-100]
+u0 = u02
+del u02
+
+Im2 = np.zeros([N,M-100],dtype='float32')
+Im2[:,0:M-100] = Im[:,0:M-100]
+Im = Im2
+del Im2
 
+#%%
 print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
 print ("_______________ROF-TV (2D)_________________")
 print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
diff --git a/Wrappers/Python/demos/demo_gpu_regularisers.py b/Wrappers/Python/demos/demo_gpu_regularisers.py
index b873700..dfdceee 100644
--- a/Wrappers/Python/demos/demo_gpu_regularisers.py
+++ b/Wrappers/Python/demos/demo_gpu_regularisers.py
@@ -49,6 +49,22 @@ u_ref = Im + np.random.normal(loc = 0 ,
 u0 = u0.astype('float32')
 u_ref = u_ref.astype('float32')
 
+(N,M) = np.shape(u0)
+u_ref2 = np.zeros([N,M-100],dtype='float32')
+u_ref2[:,0:M-100] = u_ref[:,0:M-100]
+u_ref = u_ref2
+del u_ref2
+
+u02 = np.zeros([N,M-100],dtype='float32')
+u02[:,0:M-100] = u0[:,0:M-100]
+u0 = u02
+del u02
+
+Im2 = np.zeros([N,M-100],dtype='float32')
+Im2[:,0:M-100] = Im[:,0:M-100]
+Im = Im2
+del Im2
+
 print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
 print ("____________ROF-TV regulariser_____________")
 print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
diff --git a/Wrappers/Python/src/cpu_regularisers.pyx b/Wrappers/Python/src/cpu_regularisers.pyx
index 52befd7..21a1a00 100644
--- a/Wrappers/Python/src/cpu_regularisers.pyx
+++ b/Wrappers/Python/src/cpu_regularisers.pyx
@@ -49,7 +49,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
             
@@ -333,18 +333,20 @@ def NDF_INPAINT_CPU(inputData, maskData, regularisation_parameter, edge_paramete
         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,
+                     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')   
+            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[0], dims[1], 1)    
-- 
cgit v1.2.3


From a64fe4d083173cc67dd7585c3160a94ea24bca80 Mon Sep 17 00:00:00 2001
From: Daniil Kazantsev <dkazanc@hotmail.com>
Date: Wed, 2 May 2018 10:06:38 +0100
Subject: cyth corr

---
 Wrappers/Python/src/cpu_regularisers.pyx | 13 ++++++-------
 Wrappers/Python/src/gpu_regularisers.pyx | 10 +++++-----
 2 files changed, 11 insertions(+), 12 deletions(-)

diff --git a/Wrappers/Python/src/cpu_regularisers.pyx b/Wrappers/Python/src/cpu_regularisers.pyx
index 21a1a00..7c06c28 100644
--- a/Wrappers/Python/src/cpu_regularisers.pyx
+++ b/Wrappers/Python/src/cpu_regularisers.pyx
@@ -102,7 +102,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        
             
@@ -161,7 +161,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        
             
@@ -222,7 +222,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        
             
@@ -301,7 +301,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, 
@@ -349,7 +349,7 @@ def NDF_INP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
             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[0], dims[1], 1)    
+    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, 
@@ -396,7 +396,6 @@ def NVM_INP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
     
     # 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,  
-    dims[0], dims[1], 1)    
+    SW_increment, iterationsNumb,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, 
-- 
cgit v1.2.3


From 985fee04ac1abef2aaa69f282ae6c207e438b4af Mon Sep 17 00:00:00 2001
From: algol <dkazanc@hotmail.com>
Date: Wed, 2 May 2018 11:01:57 +0100
Subject: bugs in cython files

---
 Core/regularisers_CPU/TNV_core.c               |  36 +++++++++++-----------
 Wrappers/Python/demos/demo_cpu_inpainters.py   |   2 +-
 Wrappers/Python/demos/demo_cpu_regularisers.py |  40 +++++++++++--------------
 Wrappers/Python/demos/demo_gpu_regularisers.py |  18 ++++++-----
 data/SinoInpaint.mat                           | Bin 3584100 -> 3335061 bytes
 5 files changed, 47 insertions(+), 49 deletions(-)

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/Wrappers/Python/demos/demo_cpu_inpainters.py b/Wrappers/Python/demos/demo_cpu_inpainters.py
index 348d235..7f452c1 100644
--- a/Wrappers/Python/demos/demo_cpu_inpainters.py
+++ b/Wrappers/Python/demos/demo_cpu_inpainters.py
@@ -72,7 +72,7 @@ pars = {'algorithm' : NDF_INP, \
         'maskData' : mask,\
         'regularisation_parameter':5000,\
         'edge_parameter':0,\
-        'number_of_iterations' :1000 ,\
+        'number_of_iterations' :5000 ,\
         'time_marching_parameter':0.000075,\
         'penalty_type':0
         }
diff --git a/Wrappers/Python/demos/demo_cpu_regularisers.py b/Wrappers/Python/demos/demo_cpu_regularisers.py
index f803870..986e3e9 100644
--- a/Wrappers/Python/demos/demo_cpu_regularisers.py
+++ b/Wrappers/Python/demos/demo_cpu_regularisers.py
@@ -44,29 +44,30 @@ 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
-(N,M) = np.shape(u0)
-u_ref2 = np.zeros([N,M-100],dtype='float32')
-u_ref2[:,0:M-100] = u_ref[:,0:M-100]
+"""
+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-100],dtype='float32')
-u02[:,0:M-100] = u0[:,0:M-100]
+u02 = np.zeros([N,M],dtype='float32')
+u02[:,0:M] = u0[:,0:M]
 u0 = u02
 del u02
 
-Im2 = np.zeros([N,M-100],dtype='float32')
-Im2[:,0:M-100] = Im[:,0:M-100]
+Im2 = np.zeros([N,M],dtype='float32')
+Im2[:,0:M] = Im[:,0:M]
 Im = Im2
 del Im2
-
+"""
 #%%
 print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
 print ("_______________ROF-TV (2D)_________________")
@@ -305,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 ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
@@ -318,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))
@@ -361,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 ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
@@ -420,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 dfdceee..f3ed50c 100644
--- a/Wrappers/Python/demos/demo_gpu_regularisers.py
+++ b/Wrappers/Python/demos/demo_gpu_regularisers.py
@@ -44,26 +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')
-
-(N,M) = np.shape(u0)
-u_ref2 = np.zeros([N,M-100],dtype='float32')
-u_ref2[:,0:M-100] = u_ref[:,0:M-100]
+"""
+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-100],dtype='float32')
-u02[:,0:M-100] = u0[:,0:M-100]
+u02 = np.zeros([N,M],dtype='float32')
+u02[:,0:M] = u0[:,0:M]
 u0 = u02
 del u02
 
-Im2 = np.zeros([N,M-100],dtype='float32')
-Im2[:,0:M-100] = Im[:,0:M-100]
+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/data/SinoInpaint.mat b/data/SinoInpaint.mat
index 58174a9..d748fb4 100644
Binary files a/data/SinoInpaint.mat and b/data/SinoInpaint.mat differ
-- 
cgit v1.2.3


From 7066b1122a11bba244777052b24c103e612f7977 Mon Sep 17 00:00:00 2001
From: Daniil Kazantsev <dkazanc@hotmail.com>
Date: Wed, 2 May 2018 11:18:56 +0100
Subject: NVM

---
 Core/inpainters_CPU/NonlocalMarching_Inpaint_core.c | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/Core/inpainters_CPU/NonlocalMarching_Inpaint_core.c b/Core/inpainters_CPU/NonlocalMarching_Inpaint_core.c
index d6cdf93..3a846bb 100644
--- a/Core/inpainters_CPU/NonlocalMarching_Inpaint_core.c
+++ b/Core/inpainters_CPU/NonlocalMarching_Inpaint_core.c
@@ -56,7 +56,9 @@ float NonlocalMarching_Inpaint_main(float *Input, unsigned char *M, float *Outpu
             if (M[i] == 1) iterations_number++;
         }
     }
-    else iterations_number = iterationsNumb;
+    else iterations_number = (int)(iterationsNumb/dimX);
+    if (iterations_number > dimX) iterations_number = dimX;
+    
     if (iterations_number == 0) printf("%s \n", "Nothing to inpaint, zero mask!");
     else {
         
-- 
cgit v1.2.3


From 14edd18d07c871c0a355d70e68350a899014dbc7 Mon Sep 17 00:00:00 2001
From: algol <dkazanc@hotmail.com>
Date: Wed, 2 May 2018 13:11:50 +0100
Subject: bugs in NVM fixed

---
 .../inpainters_CPU/NonlocalMarching_Inpaint_core.c | 48 +++++++++++++++++-----
 Wrappers/Python/demos/demo_cpu_inpainters.py       |  6 ++-
 Wrappers/Python/src/cpu_regularisers.pyx           | 12 +++---
 3 files changed, 48 insertions(+), 18 deletions(-)

diff --git a/Core/inpainters_CPU/NonlocalMarching_Inpaint_core.c b/Core/inpainters_CPU/NonlocalMarching_Inpaint_core.c
index 3a846bb..a197375 100644
--- a/Core/inpainters_CPU/NonlocalMarching_Inpaint_core.c
+++ b/Core/inpainters_CPU/NonlocalMarching_Inpaint_core.c
@@ -42,7 +42,7 @@ float NonlocalMarching_Inpaint_main(float *Input, unsigned char *M, float *Outpu
 {
     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);
     
@@ -54,11 +54,11 @@ float NonlocalMarching_Inpaint_main(float *Input, unsigned char *M, float *Outpu
         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 = (int)(iterationsNumb/dimX);
-    if (iterations_number > dimX) iterations_number = dimX;
-    
+    else iterations_number = iterationsNumb;    
+  
     if (iterations_number == 0) printf("%s \n", "Nothing to inpaint, zero mask!");
     else {
         
@@ -81,11 +81,38 @@ float NonlocalMarching_Inpaint_main(float *Input, unsigned char *M, float *Outpu
                     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(j=0; j<dimY; j++) {
-                switchmask = 0;
-                for(i=0; i<dimX; i++) {                
+              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 */
@@ -96,13 +123,14 @@ float NonlocalMarching_Inpaint_main(float *Input, unsigned char *M, float *Outpu
                     }
                     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-1, j, dimX, dimY, W_halfsize, W_fullsize);
+                        inpaint_func(Output, M_upd, Gauss_weights, i, j-1, dimX, dimY, W_halfsize, W_fullsize);
                         /* add value to the mask*/
-                        M_upd[(j)*dimX + i-1] = 0;                 
+                        M_upd[(j-1)*dimX + i] = 0;                 
                         switchmask = 0;                        
                     }
                 }
             }
+		}
             free(Gauss_weights);
             
             /* check if possible to terminate iterations earlier */
diff --git a/Wrappers/Python/demos/demo_cpu_inpainters.py b/Wrappers/Python/demos/demo_cpu_inpainters.py
index 7f452c1..9197e91 100644
--- a/Wrappers/Python/demos/demo_cpu_inpainters.py
+++ b/Wrappers/Python/demos/demo_cpu_inpainters.py
@@ -55,6 +55,7 @@ plt.imshow(mask)
 plt.title('Mask')
 plt.show()
 #%%
+"""
 print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
 print ("___Inpainting using linear diffusion (2D)__")
 print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
@@ -111,7 +112,7 @@ 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,cmap="gray")
+imgplot = plt.imshow(sino_cut_new,cmap="gray")
 
 # set parameters
 pars = {'algorithm' : NDF_INP, \
@@ -148,6 +149,7 @@ 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")
@@ -165,7 +167,7 @@ pars = {'algorithm' : NVM_INP, \
         'input' : sino_cut_new,\
         'maskData' : mask,\
         'SW_increment': 1,\
-        'number_of_iterations' :0
+        'number_of_iterations' : 150
         }
         
 start_time = timeit.default_timer()
diff --git a/Wrappers/Python/src/cpu_regularisers.pyx b/Wrappers/Python/src/cpu_regularisers.pyx
index 7c06c28..732b4cb 100644
--- a/Wrappers/Python/src/cpu_regularisers.pyx
+++ b/Wrappers/Python/src/cpu_regularisers.pyx
@@ -27,7 +27,6 @@ cdef extern float dTV_FGP_CPU_main(float *Input, float *InputRef, float *Output,
 
 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 dimX, int dimY, int dimZ);
-
 #****************************************************************#
 #********************** Total-variation ROF *********************#
 #****************************************************************#
@@ -374,14 +373,14 @@ def NDF_INP_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
 #*********************Inpainting WITH****************************#
 #***************Nonlocal Vertical Marching method****************#
 #****************************************************************#
-def NVM_INPAINT_CPU(inputData, maskData, SW_increment, iterations):
+def NVM_INPAINT_CPU(inputData, maskData, SW_increment, iterationsNumb):
     if inputData.ndim == 2:
-        return NVM_INP_2D(inputData, maskData, SW_increment, iterations)
+        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,                     
+               np.ndarray[np.uint8_t, ndim=2, mode="c"] maskData,
                      int SW_increment,
                      int iterationsNumb):
     cdef long dims[2]
@@ -395,7 +394,8 @@ def NVM_INP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
             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,dims[1], dims[0], 1)
+    NonlocalMarching_Inpaint_main(&inputData[0,0], &maskData[0,0], &outputData[0,0], 
+                                  &maskData_upd[0,0],
+                                  SW_increment, iterationsNumb, dims[1], dims[0], 1)
     
     return (outputData, maskData_upd)
-- 
cgit v1.2.3


From 37ae2bdb0a15298f312e9f6545a465d4d20c57f1 Mon Sep 17 00:00:00 2001
From: Daniil Kazantsev <dkazanc@hotmail.com>
Date: Wed, 2 May 2018 15:47:19 +0100
Subject: bugs of NVM are fixed

---
 .../inpainters_CPU/NonlocalMarching_Inpaint_core.c | 100 ++++++++++-----------
 .../inpainters_CPU/NonlocalMarching_Inpaint_core.h |   2 +-
 .../regularisers_CPU/NonlocalMarching_Inpaint.c    |   2 +-
 Wrappers/Python/src/cpu_regularisers.pyx           |   4 +-
 4 files changed, 54 insertions(+), 54 deletions(-)

diff --git a/Core/inpainters_CPU/NonlocalMarching_Inpaint_core.c b/Core/inpainters_CPU/NonlocalMarching_Inpaint_core.c
index a197375..66be15c 100644
--- a/Core/inpainters_CPU/NonlocalMarching_Inpaint_core.c
+++ b/Core/inpainters_CPU/NonlocalMarching_Inpaint_core.c
@@ -30,7 +30,7 @@
  * 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
@@ -38,11 +38,11 @@
  * Reference: TBA
  */
 
-float NonlocalMarching_Inpaint_main(float *Input, unsigned char *M, float *Output, unsigned char *M_upd, int SW_increment, int iterationsNumb, int dimX, int dimY, int dimZ)
+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);
     
@@ -54,11 +54,11 @@ float NonlocalMarching_Inpaint_main(float *Input, unsigned char *M, float *Outpu
         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;    
+        }
+        if ((int)(iterations_number/dimY) > dimX) iterations_number = dimX;
     }
-    else iterations_number = iterationsNumb;    
-  
+    else iterations_number = iterationsNumb;
+    
     if (iterations_number == 0) printf("%s \n", "Nothing to inpaint, zero mask!");
     else {
         
@@ -83,54 +83,54 @@ float NonlocalMarching_Inpaint_main(float *Input, unsigned char *M, float *Outpu
             }
             
             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;                
+                /*Matlab*/
+#pragma omp parallel for shared(Output, M_upd, Gauss_weights) private(i, j, switchmask, switchcurr)
                 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;
+                    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;
+                        }
                     }
-                    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;                        
+                }
+            }
+            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 */
diff --git a/Core/inpainters_CPU/NonlocalMarching_Inpaint_core.h b/Core/inpainters_CPU/NonlocalMarching_Inpaint_core.h
index 199820c..0f99ed4 100644
--- a/Core/inpainters_CPU/NonlocalMarching_Inpaint_core.h
+++ b/Core/inpainters_CPU/NonlocalMarching_Inpaint_core.h
@@ -47,7 +47,7 @@ limitations under the License.
 #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 dimX, int dimY, int dimZ);
+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
 }
diff --git a/Wrappers/Matlab/mex_compile/regularisers_CPU/NonlocalMarching_Inpaint.c b/Wrappers/Matlab/mex_compile/regularisers_CPU/NonlocalMarching_Inpaint.c
index 5e4ab1f..36cf05c 100644
--- a/Wrappers/Matlab/mex_compile/regularisers_CPU/NonlocalMarching_Inpaint.c
+++ b/Wrappers/Matlab/mex_compile/regularisers_CPU/NonlocalMarching_Inpaint.c
@@ -78,5 +78,5 @@ void mexFunction(
     if (number_of_dims == 3) {
         mexErrMsgTxt("Currently 2D supported only");        
     }           
-    NonlocalMarching_Inpaint_main(Input, Mask, Output, Mask_upd, SW_increment, iterations, dimX, dimY, dimZ);
+    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/src/cpu_regularisers.pyx b/Wrappers/Python/src/cpu_regularisers.pyx
index 732b4cb..c934f1d 100644
--- a/Wrappers/Python/src/cpu_regularisers.pyx
+++ b/Wrappers/Python/src/cpu_regularisers.pyx
@@ -26,7 +26,7 @@ cdef extern float TNV_CPU_main(float *Input, float *u, float lambdaPar, int maxI
 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 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 *********************#
 #****************************************************************#
@@ -396,6 +396,6 @@ def NVM_INP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
     # 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, dims[1], dims[0], 1)
+                                  SW_increment, iterationsNumb, 1, dims[1], dims[0], 1)
     
     return (outputData, maskData_upd)
-- 
cgit v1.2.3


From 6e285c109938a43b5f8a84b7a48afaeb6b058c90 Mon Sep 17 00:00:00 2001
From: algol <dkazanc@hotmail.com>
Date: Wed, 2 May 2018 15:43:55 +0100
Subject: demo inpainters fixed for python

---
 Wrappers/Python/demos/demo_cpu_inpainters.py | 2 --
 1 file changed, 2 deletions(-)

diff --git a/Wrappers/Python/demos/demo_cpu_inpainters.py b/Wrappers/Python/demos/demo_cpu_inpainters.py
index 9197e91..3b4191b 100644
--- a/Wrappers/Python/demos/demo_cpu_inpainters.py
+++ b/Wrappers/Python/demos/demo_cpu_inpainters.py
@@ -55,7 +55,6 @@ plt.imshow(mask)
 plt.title('Mask')
 plt.show()
 #%%
-"""
 print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
 print ("___Inpainting using linear diffusion (2D)__")
 print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
@@ -149,7 +148,6 @@ 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")
-- 
cgit v1.2.3