From d6ee5585e696f855d1c687d34efa04328729e94c Mon Sep 17 00:00:00 2001
From: dkazanc <dkazanc@hotmail.com>
Date: Tue, 9 Apr 2019 16:44:39 +0100
Subject: [PATCH 1/7] 2D CPU version for constrained diffusion

---
 build/run.sh                                  |  18 +-
 src/Core/CMakeLists.txt                       |   1 +
 .../regularisers_CPU/DiffusionMASK_core.c     | 214 ++++++++++++++++++
 .../regularisers_CPU/DiffusionMASK_core.h     |  62 +++++
 src/Python/ccpi/filters/regularisers.py       |  27 ++-
 src/Python/src/cpu_regularisers.pyx           |  37 +++
 6 files changed, 349 insertions(+), 10 deletions(-)
 create mode 100644 src/Core/regularisers_CPU/DiffusionMASK_core.c
 create mode 100644 src/Core/regularisers_CPU/DiffusionMASK_core.h

diff --git a/build/run.sh b/build/run.sh
index e6f171e5..b40d222f 100755
--- a/build/run.sh
+++ b/build/run.sh
@@ -1,14 +1,14 @@
 #!/bin/bash
 echo "Building CCPi-regularisation Toolkit using CMake"
-rm -r build_proj
+rm -r ../build_proj
 # Requires Cython, install it first:
 # pip install cython
-mkdir build_proj
-cd build_proj/
+mkdir ../build_proj
+cd ../build_proj/
 #make clean
-export CIL_VERSION=19.03
+export CIL_VERSION=19.04
 # install Python modules without CUDA
-# cmake ../ -DBUILD_PYTHON_WRAPPER=ON -DBUILD_MATLAB_WRAPPER=OFF -DBUILD_CUDA=OFF -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=./install
+cmake ../ -DBUILD_PYTHON_WRAPPER=ON -DBUILD_MATLAB_WRAPPER=OFF -DBUILD_CUDA=OFF -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=./install
 # install Python modules with CUDA
 # cmake ../ -DBUILD_PYTHON_WRAPPER=ON -DBUILD_MATLAB_WRAPPER=OFF -DBUILD_CUDA=ON -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=./install
 # install Matlab modules without CUDA
@@ -18,12 +18,12 @@ export CIL_VERSION=19.03
 # cmake ../ -DBUILD_PYTHON_WRAPPER=OFF -DMatlab_ROOT_DIR=/home/algol/SOFT/MATLAB9/ -DBUILD_MATLAB_WRAPPER=ON -DBUILD_CUDA=ON -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=./install
 make install
 ############### Python(linux)###############
-#cp install/lib/libcilreg.so install/python/ccpi/filters
-# cd install/python
-# export LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:../lib
+cp install/lib/libcilreg.so install/python/ccpi/filters
+cd install/python
+export LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:../lib
 # spyder
 ############### Matlab(linux)###############
-export LD_PRELOAD=/home/algol/anaconda3/lib/libstdc++.so.6  # if there is libstdc error in matlab
+# export LD_PRELOAD=/home/algol/anaconda3/lib/libstdc++.so.6  # if there is libstdc error in matlab
 # PATH="/path/to/mex/:$PATH" LD_LIBRARY_PATH="/path/to/library:$LD_LIBRARY_PATH" matlab
 # PATH="/home/kjy41806/Documents/SOFT/CCPi-Regularisation-Toolkit/build_proj/install/matlab/:$PATH" LD_LIBRARY_PATH="/home/kjy41806/Documents/SOFT/CCPi-Regularisation-Toolkit/build_proj/install/lib:$LD_LIBRARY_PATH" matlab
 #PATH="/home/algol/Documents/DEV/CCPi-Regularisation-Toolkit/build_proj/install/matlab/:$PATH" LD_LIBRARY_PATH="/home/algol/Documents/DEV/CCPi-Regularisation-Toolkit/build_proj/install/lib:$LD_LIBRARY_PATH" /home/algol/SOFT/MATLAB9/bin/matlab
diff --git a/src/Core/CMakeLists.txt b/src/Core/CMakeLists.txt
index b3c0dfbe..6975a89f 100644
--- a/src/Core/CMakeLists.txt
+++ b/src/Core/CMakeLists.txt
@@ -68,6 +68,7 @@ add_library(cilreg SHARED
 	    ${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/SB_TV_core.c
 	    ${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/TGV_core.c
 	    ${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/Diffusion_core.c
+    	    ${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/DiffusionMASK_core.c
 	    ${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/Diffus4th_order_core.c
 	    ${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/LLT_ROF_core.c
         ${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/ROF_TV_core.c
diff --git a/src/Core/regularisers_CPU/DiffusionMASK_core.c b/src/Core/regularisers_CPU/DiffusionMASK_core.c
new file mode 100644
index 00000000..eef173d3
--- /dev/null
+++ b/src/Core/regularisers_CPU/DiffusionMASK_core.c
@@ -0,0 +1,214 @@
+/*
+ * 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 "DiffusionMASK_core.h"
+#include "utils.h"
+
+#define EPS 1.0e-5
+#define MAX(x, y) (((x) > (y)) ? (x) : (y))
+#define MIN(x, y) (((x) < (y)) ? (x) : (y))
+
+/*sign function*/
+int signNDF_m(float x) {
+    return (x > 0) - (x < 0);
+}
+
+/* C-OMP implementation of linear and nonlinear diffusion [1,2] which is constrained by the provided MASK.
+ * The minimisation is performed using explicit scheme.
+ * Implementation using the diffusivity window to increase the coverage area of the diffusivity
+ *
+ * Input Parameters:
+ * 1. Noisy image/volume
+ * 2. MASK (in unsigned char format)
+ * 3. Diffusivity window (half-size of the searching window, e.g. 3) 	 
+ * 4. lambda - regularization parameter
+ * 5. Edge-preserving parameter (sigma), when sigma equals to zero nonlinear diffusion -> linear diffusion
+ * 6. Number of iterations, for explicit scheme >= 150 is recommended
+ * 7. tau - time-marching step for explicit scheme
+ * 8. Penalty type: 1 - Huber, 2 - Perona-Malik, 3 - Tukey Biweight
+ * 9. eplsilon - tolerance constant
+
+ * Output:
+ * [1] Filtered/regularized image/volume
+ * [2] Information vector which contains [iteration no., reached tolerance]
+ *
+ * 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 DiffusionMASK_CPU_main(float *Input, unsigned char *MASK, float *Output, float *infovector, int DiffusWindow, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int penaltytype, float epsil, int dimX, int dimY, int dimZ)
+{
+    int i,j,k,counterG;
+    float sigmaPar2, *Output_prev=NULL, *Eucl_Vec;
+    int DiffusWindow_tot;
+    sigmaPar2 = sigmaPar/sqrt(2.0f);
+    long DimTotal;
+    float re, re1;
+    re = 0.0f; re1 = 0.0f;
+    int count = 0;
+    DimTotal = (long)(dimX*dimY*dimZ);
+
+    if (dimZ == 1) {
+	DiffusWindow_tot = (2*DiffusWindow + 1)*(2*DiffusWindow + 1);
+        /* generate a 2D Gaussian kernel for NLM procedure */
+        Eucl_Vec = (float*) calloc (DiffusWindow_tot,sizeof(float));
+        counterG = 0;
+        for(i=-DiffusWindow; i<=DiffusWindow; i++) {
+            for(j=-DiffusWindow; j<=DiffusWindow; j++) {
+                Eucl_Vec[counterG] = (float)expf(-(powf(((float) i), 2) + powf(((float) j), 2))/(2.0f*DiffusWindow*DiffusWindow));
+                counterG++;
+            }} /*main neighb loop */
+       }
+    else {
+	DiffusWindow_tot = (2*DiffusWindow + 1)*(2*DiffusWindow + 1)*(2*DiffusWindow + 1);
+	Eucl_Vec = (float*) calloc (DiffusWindow_tot,sizeof(float));
+        counterG = 0;
+        for(i=-DiffusWindow; i<=DiffusWindow; i++) {
+            for(j=-DiffusWindow; j<=DiffusWindow; j++) {
+                for(k=-DiffusWindow; k<=DiffusWindow; k++) {
+                    Eucl_Vec[counterG] = (float)expf(-(powf(((float) i), 2) + powf(((float) j), 2) + powf(((float) k), 2))/(2*DiffusWindow*DiffusWindow*DiffusWindow));
+                    counterG++;
+                }}} /*main neighb loop */            
+    }
+
+    if (epsil != 0.0f) Output_prev = calloc(DimTotal, sizeof(float));
+
+    /* copy into output */
+    copyIm(Input, Output, (long)(dimX), (long)(dimY), (long)(dimZ));
+
+    for(i=0; i < iterationsNumb; i++) {
+
+      if ((epsil != 0.0f)  && (i % 5 == 0)) copyIm(Output, Output_prev, (long)(dimX), (long)(dimY), (long)(dimZ));
+      if (dimZ == 1) {
+	    /* running 2D diffusion iterations */
+            if (sigmaPar == 0.0f) LinearDiff_MASK2D(Input, MASK, Output, Eucl_Vec, DiffusWindow, lambdaPar, tau, (long)(dimX), (long)(dimY)); /* constrained linear diffusion */
+            else NonLinearDiff_MASK2D(Input, MASK, Output, Eucl_Vec, DiffusWindow, lambdaPar, sigmaPar2, tau, penaltytype, (long)(dimX), (long)(dimY)); /* constrained nonlinear diffusion */
+          }
+      else {
+       	/* running 3D diffusion iterations */
+        //if (sigmaPar == 0.0f) LinearDiff3D(Input, Output, lambdaPar, tau, (long)(dimX), (long)(dimY), (long)(dimZ));
+//       else NonLinearDiff3D(Input, Output, lambdaPar, sigmaPar2, tau, penaltytype, (long)(dimX), (long)(dimY), (long)(dimZ));
+          }
+          /* check early stopping criteria if epsilon not equal zero */
+          if ((epsil != 0.0f)  && (i % 5 == 0)) {
+          re = 0.0f; re1 = 0.0f;
+            for(j=0; j<DimTotal; j++)
+            {
+                re += powf(Output[j] - Output_prev[j],2);
+                re1 += powf(Output[j],2);
+            }
+          re = sqrtf(re)/sqrtf(re1);
+          /* stop if the norm residual is less than the tolerance EPS */
+          if (re < epsil)  count++;
+          if (count > 3) break;
+          }
+		}
+
+    free(Output_prev);
+  /*adding info into info_vector */
+    infovector[0] = (float)(i);  /*iterations number (if stopped earlier based on tolerance)*/
+    infovector[1] = re;  /* reached tolerance */
+    return 0;
+}
+
+
+/********************************************************************/
+/***************************2D Functions*****************************/
+/********************************************************************/
+/* MASKED-constrained 2D linear diffusion (PDE heat equation) */
+float LinearDiff_MASK2D(float *Input, unsigned char *MASK, float *Output,  float *Eucl_Vec, int DiffusWindow, float lambdaPar, float tau, long dimX, long dimY)
+{
+
+long i,j,i1,j1,i_m,j_m,index,indexneighb,counter;
+unsigned char class_c, class_n;
+float diffVal;
+
+#pragma omp parallel for shared(Input) private(index,i,j,i1,j1,i_m,j_m,counter,diffVal,indexneighb,class_c,class_n)
+    for(i=0; i<dimX; i++) {
+        for(j=0; j<dimY; j++) {        
+        index = j*dimX+i; /* current pixel index */
+        counter = 0; diffVal = 0.0f;
+        for(i_m=-DiffusWindow; i_m<=DiffusWindow; i_m++) {
+            for(j_m=-DiffusWindow; j_m<=DiffusWindow; j_m++) {		
+            i1 = i+i_m;
+            j1 = j+j_m;
+            if (((i1 >= 0) && (i1 < dimX)) && ((j1 >= 0) && (j1 < dimY))) {                   
+            indexneighb = j1*dimX+i1; /* neighbour pixel index */
+	    class_c = MASK[index]; /* current class value */
+    	    class_n = MASK[indexneighb]; /* neighbour class value */
+	    
+	    /* perform diffusion only within the same class (given by MASK) */
+	    if (class_n == class_c) diffVal += Output[indexneighb] - Output[index];
+            }
+		counter++;
+	    }}
+            Output[index] += tau*(lambdaPar*(diffVal) - (Output[index] - Input[index]));
+        }}
+	return *Output;
+}
+
+/*  MASKED-constrained 2D nonlinear diffusion */
+float NonLinearDiff_MASK2D(float *Input, unsigned char *MASK, float *Output, float *Eucl_Vec, int DiffusWindow, float lambdaPar, float sigmaPar, float tau, int penaltytype, long dimX, long dimY)
+{
+	long i,j,i1,j1,i_m,j_m,index,indexneighb,counter;
+	unsigned char class_c, class_n;
+	float diffVal, funcVal;
+	
+#pragma omp parallel for shared(Input) private(index,i,j,i1,j1,i_m,j_m,counter,diffVal,funcVal,indexneighb,class_c,class_n)
+    for(i=0; i<dimX; i++) {
+        for(j=0; j<dimY; j++) {        
+        index = j*dimX+i; /* current pixel index */
+        counter = 0; diffVal = 0.0f; funcVal = 0.0f;
+        for(i_m=-DiffusWindow; i_m<=DiffusWindow; i_m++) {
+            for(j_m=-DiffusWindow; j_m<=DiffusWindow; j_m++) {		
+            i1 = i+i_m;
+            j1 = j+j_m;
+            if (((i1 >= 0) && (i1 < dimX)) && ((j1 >= 0) && (j1 < dimY))) {                   
+            indexneighb = j1*dimX+i1; /* neighbour pixel index */
+	    class_c = MASK[index]; /* current class value */
+    	    class_n = MASK[indexneighb]; /* neighbour class value */
+	    
+	    /* perform diffusion only within the same class (given by MASK) */
+	    if (class_n == class_c) {
+	    	diffVal = Output[indexneighb] - Output[index];
+	    	if (penaltytype == 1) {
+	        /* Huber penalty */
+                if (fabs(diffVal) > sigmaPar) funcVal += signNDF_m(diffVal);
+                else funcVal += diffVal/sigmaPar; }
+  		else if (penaltytype == 2) {
+  		/* Perona-Malik */
+  		funcVal += (diffVal)/(1.0f + powf((diffVal/sigmaPar),2)); }
+  		else if (penaltytype == 3) {
+  		/* Tukey Biweight */
+  		if (fabs(diffVal) <= sigmaPar) funcVal += diffVal*powf((1.0f - powf((diffVal/sigmaPar),2)), 2); }
+                else {
+                printf("%s \n", "No penalty function selected! Use Huber,2 or 3.");
+		break; }  		  
+            		}
+            	}
+		counter++;
+	    }}
+           Output[index] += tau*(lambdaPar*(funcVal) - (Output[index] - Input[index]));
+		}}
+	return *Output;
+}
+/********************************************************************/
+/***************************3D Functions*****************************/
+/********************************************************************/
diff --git a/src/Core/regularisers_CPU/DiffusionMASK_core.h b/src/Core/regularisers_CPU/DiffusionMASK_core.h
new file mode 100644
index 00000000..8890c736
--- /dev/null
+++ b/src/Core/regularisers_CPU/DiffusionMASK_core.h
@@ -0,0 +1,62 @@
+/*
+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] which is constrained by the provided MASK.
+ * The minimisation is performed using explicit scheme.
+ * Implementation using the Diffusivity window to increase the coverage area of the diffusivity
+ *
+ * Input Parameters:
+ * 1. Noisy image/volume
+ * 2. MASK (in unsigned short format)
+ * 3. Diffusivity window (half-size of the searching window, e.g. 3) 	 
+ * 4. lambda - regularization parameter
+ * 5. Edge-preserving parameter (sigma), when sigma equals to zero nonlinear diffusion -> linear diffusion
+ * 6. Number of iterations, for explicit scheme >= 150 is recommended
+ * 7. tau - time-marching step for explicit scheme
+ * 8. Penalty type: 1 - Huber, 2 - Perona-Malik, 3 - Tukey Biweight
+ * 9. eplsilon - tolerance constant
+
+ * Output:
+ * [1] Filtered/regularized image/volume
+ * [2] Information vector which contains [iteration no., reached tolerance]
+ *
+ * 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 DiffusionMASK_CPU_main(float *Input, unsigned char *MASK, float *Output, float *infovector, int DiffusWindow, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int penaltytype, float epsil, int dimX, int dimY, int dimZ);
+CCPI_EXPORT float LinearDiff_MASK2D(float *Input, unsigned char *MASK, float *Output,  float *Eucl_Vec, int DiffusWindow, float lambdaPar, float tau, long dimX, long dimY);
+CCPI_EXPORT float NonLinearDiff_MASK2D(float *Input, unsigned char *MASK, float *Output, float *Eucl_Vec, int DiffusWindow, float lambdaPar, float sigmaPar, float tau, int penaltytype, long dimX, long dimY);
+#ifdef __cplusplus
+}
+#endif
diff --git a/src/Python/ccpi/filters/regularisers.py b/src/Python/ccpi/filters/regularisers.py
index 398e11c2..1e427bf3 100644
--- a/src/Python/ccpi/filters/regularisers.py
+++ b/src/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, Diff4th_CPU, TGV_CPU, LLT_ROF_CPU, PATCHSEL_CPU, NLTV_CPU
+from ccpi.filters.cpu_regularisers import TV_ROF_CPU, TV_FGP_CPU, TV_SB_CPU, dTV_FGP_CPU, TNV_CPU, NDF_CPU, NDF_MASK_CPU, Diff4th_CPU, TGV_CPU, LLT_ROF_CPU, PATCHSEL_CPU, NLTV_CPU
 try:
     from ccpi.filters.gpu_regularisers import TV_ROF_GPU, TV_FGP_GPU, TV_SB_GPU, dTV_FGP_GPU, NDF_GPU, Diff4th_GPU, TGV_GPU, LLT_ROF_GPU, PATCHSEL_GPU
     gpu_enabled = True
@@ -127,6 +127,31 @@ def NDF(inputData, regularisation_parameter, edge_parameter, iterations,
     	    raise ValueError ('GPU is not available')
         raise ValueError('Unknown device {0}. Expecting gpu or cpu'\
                          .format(device))
+def NDF_MASK(inputData, diffuswindow, regularisation_parameter, edge_parameter, iterations,
+                     time_marching_parameter, penalty_type, tolerance_param, device='cpu'):
+    if device == 'cpu':
+        return NDF_MASK_CPU(inputData,
+                     diffuswindow,
+                     regularisation_parameter,
+                     edge_parameter,
+                     iterations,
+                     time_marching_parameter,
+                     penalty_type,
+                     tolerance_param)
+    elif device == 'gpu' and gpu_enabled:
+        return NDF_MASK_CPU(inputData,
+                     diffuswindow,
+                     regularisation_parameter,
+                     edge_parameter,
+                     iterations,
+                     time_marching_parameter,
+                     penalty_type,
+                     tolerance_param)
+    else:
+        if not gpu_enabled and device == 'gpu':
+    	    raise ValueError ('GPU is not available')
+        raise ValueError('Unknown device {0}. Expecting gpu or cpu'\
+                         .format(device))
 def Diff4th(inputData, regularisation_parameter, edge_parameter, iterations,
                      time_marching_parameter, tolerance_param, device='cpu'):
     if device == 'cpu':
diff --git a/src/Python/src/cpu_regularisers.pyx b/src/Python/src/cpu_regularisers.pyx
index add641b3..305ee1ff 100644
--- a/src/Python/src/cpu_regularisers.pyx
+++ b/src/Python/src/cpu_regularisers.pyx
@@ -24,6 +24,7 @@ cdef extern float SB_TV_CPU_main(float *Input, float *Output, float *infovector,
 cdef extern float LLT_ROF_CPU_main(float *Input, float *Output, float *infovector, float lambdaROF, float lambdaLLT, int iterationsNumb, float tau, float epsil, int dimX, int dimY, int dimZ);
 cdef extern float TGV_main(float *Input, float *Output, float *infovector, float lambdaPar, float alpha1, float alpha0, int iterationsNumb, float L2, float epsil, int dimX, int dimY, int dimZ);
 cdef extern float Diffusion_CPU_main(float *Input, float *Output, float *infovector, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int penaltytype, float epsil, int dimX, int dimY, int dimZ);
+cdef extern float DiffusionMASK_CPU_main(float *Input, unsigned char *MASK, float *Output, float *infovector, int DiffusWindow, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int penaltytype, float epsil, int dimX, int dimY, int dimZ);
 cdef extern float Diffus4th_CPU_main(float *Input, float *Output,  float *infovector, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, float epsil, int dimX, int dimY, int dimZ);
 cdef extern float dTV_FGP_CPU_main(float *Input, float *InputRef, float *Output, float *infovector, float lambdaPar, int iterationsNumb, float epsil, float eta, int methodTV, int nonneg, int dimX, int dimY, int dimZ);
 cdef extern float TNV_CPU_main(float *Input, float *u, float lambdaPar, int maxIter, float tol, int dimX, int dimY, int dimZ);
@@ -379,6 +380,42 @@ def NDF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
     tolerance_param,
     dims[2], dims[1], dims[0])
     return (outputData,infovec)
+
+#****************************************************************#
+#********Constrained Nonlinear(Isotropic) Diffusion**************#
+#****************************************************************#
+def NDF_MASK_CPU(inputData, maskData, diffuswindow, regularisation_parameter, edge_parameter, iterationsNumb,time_marching_parameter, penalty_type, tolerance_param):
+    if inputData.ndim == 2:
+        return NDF_MASK_2D(inputData, maskData, diffuswindow, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type, tolerance_param)
+    elif inputData.ndim == 3:
+        return 0
+
+def NDF_MASK_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
+                    np.ndarray[np.uint8_t, ndim=2, mode="c"] maskData,
+                     int diffuswindow,
+                     float regularisation_parameter,
+                     float edge_parameter,
+                     int iterationsNumb,
+                     float time_marching_parameter,
+                     int penalty_type,
+                     float tolerance_param):
+    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.float32_t, ndim=1, mode="c"] infovec = \
+                np.zeros([2], dtype='float32')
+
+    # Run constrained nonlinear diffusion iterations for 2D data
+    DiffusionMASK_CPU_main(&inputData[0,0], &maskData[0,0], &outputData[0,0], &infovec[0],
+    diffuswindow, regularisation_parameter, edge_parameter, iterationsNumb,
+    time_marching_parameter, penalty_type,
+    tolerance_param,
+    dims[1], dims[0], 1)
+    return (outputData,infovec)
+
 #****************************************************************#
 #*************Anisotropic Fourth-Order diffusion*****************#
 #****************************************************************#

From 2d436f37be6029d57d7f876d4c7c378ee712a11e Mon Sep 17 00:00:00 2001
From: Daniil Kazantsev <dkazanc@hotmail.com>
Date: Wed, 10 Apr 2019 22:39:11 +0100
Subject: [PATCH 2/7] progress with bresenham

---
 build/run.sh                                  |   2 +-
 .../regularisers_CPU/DiffusionMASK_core.c     | 210 ++++++++++++++++--
 .../regularisers_CPU/DiffusionMASK_core.h     |   7 +-
 src/Python/ccpi/filters/regularisers.py       |   4 +-
 src/Python/src/cpu_regularisers.pyx           |  10 +-
 5 files changed, 208 insertions(+), 25 deletions(-)

diff --git a/build/run.sh b/build/run.sh
index b40d222f..0a153aa4 100755
--- a/build/run.sh
+++ b/build/run.sh
@@ -21,7 +21,7 @@ make install
 cp install/lib/libcilreg.so install/python/ccpi/filters
 cd install/python
 export LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:../lib
-# spyder
+spyder
 ############### Matlab(linux)###############
 # export LD_PRELOAD=/home/algol/anaconda3/lib/libstdc++.so.6  # if there is libstdc error in matlab
 # PATH="/path/to/mex/:$PATH" LD_LIBRARY_PATH="/path/to/library:$LD_LIBRARY_PATH" matlab
diff --git a/src/Core/regularisers_CPU/DiffusionMASK_core.c b/src/Core/regularisers_CPU/DiffusionMASK_core.c
index eef173d3..2af3cb66 100644
--- a/src/Core/regularisers_CPU/DiffusionMASK_core.c
+++ b/src/Core/regularisers_CPU/DiffusionMASK_core.c
@@ -36,7 +36,7 @@ int signNDF_m(float x) {
  * Input Parameters:
  * 1. Noisy image/volume
  * 2. MASK (in unsigned char format)
- * 3. Diffusivity window (half-size of the searching window, e.g. 3) 	 
+ * 3. Diffusivity window (half-size of the searching window, e.g. 3)
  * 4. lambda - regularization parameter
  * 5. Edge-preserving parameter (sigma), when sigma equals to zero nonlinear diffusion -> linear diffusion
  * 6. Number of iterations, for explicit scheme >= 150 is recommended
@@ -53,9 +53,10 @@ int signNDF_m(float x) {
  * [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 DiffusionMASK_CPU_main(float *Input, unsigned char *MASK, float *Output, float *infovector, int DiffusWindow, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int penaltytype, float epsil, int dimX, int dimY, int dimZ)
+float DiffusionMASK_CPU_main(float *Input, unsigned char *MASK, unsigned char *MASK_upd, float *Output, float *infovector, int DiffusWindow, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int penaltytype, float epsil, int dimX, int dimY, int dimZ)
 {
-    int i,j,k,counterG;
+    long i,j,k;
+    int counterG;
     float sigmaPar2, *Output_prev=NULL, *Eucl_Vec;
     int DiffusWindow_tot;
     sigmaPar2 = sigmaPar/sqrt(2.0f);
@@ -63,6 +64,7 @@ float DiffusionMASK_CPU_main(float *Input, unsigned char *MASK, float *Output, f
     float re, re1;
     re = 0.0f; re1 = 0.0f;
     int count = 0;
+    unsigned char *MASK_temp;
     DimTotal = (long)(dimX*dimY*dimZ);
 
     if (dimZ == 1) {
@@ -85,21 +87,48 @@ float DiffusionMASK_CPU_main(float *Input, unsigned char *MASK, float *Output, f
                 for(k=-DiffusWindow; k<=DiffusWindow; k++) {
                     Eucl_Vec[counterG] = (float)expf(-(powf(((float) i), 2) + powf(((float) j), 2) + powf(((float) k), 2))/(2*DiffusWindow*DiffusWindow*DiffusWindow));
                     counterG++;
-                }}} /*main neighb loop */            
+                }}} /*main neighb loop */
     }
 
     if (epsil != 0.0f) Output_prev = calloc(DimTotal, sizeof(float));
 
-    /* copy into output */
+    MASK_temp = (unsigned char*) calloc (DimTotal,sizeof(unsigned char));
+
+    /* copy input into output */
     copyIm(Input, Output, (long)(dimX), (long)(dimY), (long)(dimZ));
+    /* copy given MASK  */
+    copyIm_unchar(MASK, MASK_upd, (long)(dimX), (long)(dimY), (long)(dimZ));
 
+    /********************** PERFORM MASK PROCESSING ************************/
+    if (dimZ == 1) {
+    #pragma omp parallel for shared(MASK,MASK_upd) private(i,j)
+    for(i=0; i<dimX; i++) {
+        for(j=0; j<dimY; j++) {
+    /* STEP1: in a smaller neighbourhood check that the current pixel is NOT an outlier */
+    OutiersRemoval2D(MASK, MASK_upd, i, j, (long)(dimX), (long)(dimY));
+    }}
+    /* copy the updated MASK (clean of outliers) */
+    copyIm_unchar(MASK_upd, MASK_temp, (long)(dimX), (long)(dimY), (long)(dimZ));
+
+//    #pragma omp parallel for shared(MASK_temp,MASK_upd) private(i,j)
+//  for(i=0; i<dimX; i++) {
+//        for(j=0; j<dimY; j++) {
+//      if (MASK_temp[j*dimX+i] == MASK[j*dimX+i]) {
+        /*the class of the central pixel has not changed, i.e. the central pixel is not an outlier -> continue */
+        i = 162; j = 258;
+        Mask_update2D(MASK_temp, MASK_upd, i, j, DiffusWindow, (long)(dimX), (long)(dimY));
+//       }
+      //}}
+    }
+
+    /* start diffusivity iterations */
     for(i=0; i < iterationsNumb; i++) {
 
       if ((epsil != 0.0f)  && (i % 5 == 0)) copyIm(Output, Output_prev, (long)(dimX), (long)(dimY), (long)(dimZ));
       if (dimZ == 1) {
-	    /* running 2D diffusion iterations */
-            if (sigmaPar == 0.0f) LinearDiff_MASK2D(Input, MASK, Output, Eucl_Vec, DiffusWindow, lambdaPar, tau, (long)(dimX), (long)(dimY)); /* constrained linear diffusion */
-            else NonLinearDiff_MASK2D(Input, MASK, Output, Eucl_Vec, DiffusWindow, lambdaPar, sigmaPar2, tau, penaltytype, (long)(dimX), (long)(dimY)); /* constrained nonlinear diffusion */
+	          /* running 2D diffusion iterations */
+            //if (sigmaPar == 0.0f) LinearDiff_MASK2D(Input, MASK, Output, Eucl_Vec, DiffusWindow, lambdaPar, tau, (long)(dimX), (long)(dimY)); /* constrained linear diffusion */
+            //else NonLinearDiff_MASK2D(Input, MASK, Output, Eucl_Vec, DiffusWindow, lambdaPar, sigmaPar2, tau, penaltytype, (long)(dimX), (long)(dimY)); /* constrained nonlinear diffusion */
           }
       else {
        	/* running 3D diffusion iterations */
@@ -122,6 +151,8 @@ float DiffusionMASK_CPU_main(float *Input, unsigned char *MASK, float *Output, f
 		}
 
     free(Output_prev);
+    free(Eucl_Vec);
+    free(MASK_temp);
   /*adding info into info_vector */
     infovector[0] = (float)(i);  /*iterations number (if stopped earlier based on tolerance)*/
     infovector[1] = re;  /* reached tolerance */
@@ -132,6 +163,62 @@ float DiffusionMASK_CPU_main(float *Input, unsigned char *MASK, float *Output, f
 /********************************************************************/
 /***************************2D Functions*****************************/
 /********************************************************************/
+float OutiersRemoval2D(unsigned char *MASK, unsigned char *MASK_upd, long i, long j, long dimX, long dimY)
+{
+  /*if the ROI pixel does not belong to the surrondings, turn it into the surronding*/
+  long i_m, j_m, i1, j1, counter;
+    counter = 0;
+    for(i_m=-1; i_m<=1; i_m++) {
+      for(j_m=-1; j_m<=1; j_m++) {
+        i1 = i+i_m;
+        j1 = j+j_m;
+        if (((i1 >= 0) && (i1 < dimX)) && ((j1 >= 0) && (j1 < dimY))) {
+          if (MASK[j*dimX+i] != MASK[j1*dimX+i1]) counter++;
+        }
+      }}
+      if (counter >= 8) MASK_upd[j*dimX+i] = MASK[j1*dimX+i1];
+      return *MASK_upd;
+}
+
+float Mask_update2D(unsigned char *MASK_temp, unsigned char *MASK_upd, long i, long j, int DiffusWindow, long dimX, long dimY)
+{
+  long i_m, j_m, i1, j1, CounterOtherClass;
+
+  /* STEP2: in a larger neighbourhood check that the other class is present  */
+  CounterOtherClass = 0;
+  for(i_m=-DiffusWindow; i_m<=DiffusWindow; i_m++) {
+      for(j_m=-DiffusWindow; j_m<=DiffusWindow; j_m++) {
+        i1 = i+i_m;
+        j1 = j+j_m;
+        if (((i1 >= 0) && (i1 < dimX)) && ((j1 >= 0) && (j1 < dimY))) {
+          if (MASK_temp[j*dimX+i] != MASK_temp[j1*dimX+i1]) CounterOtherClass++;
+        }
+      }}
+      if (CounterOtherClass > 0) {
+      /* the other class is present in the vicinity of DiffusWindow, continue to STEP 3 */
+      /*
+      STEP 3: Loop through all neighbours in DiffusWindow and check the spatial connection.
+      Meaning that we're instrested if there are any classes between points A and B that
+      does not belong to A and B (A,B \in C)
+      */
+      for(i_m=-DiffusWindow; i_m<=DiffusWindow; i_m++) {
+          for(j_m=-DiffusWindow; j_m<=DiffusWindow; j_m++) {
+            i1 = i+i_m;
+            j1 = j+j_m;
+            if (((i1 >= 0) && (i1 < dimX)) && ((j1 >= 0) && (j1 < dimY))) {
+              if (MASK_temp[j*dimX+i] == MASK_temp[j1*dimX+i1]) {
+                /*A and B points belong to the same class, do STEP 4*/
+                /* STEP 4: Run Bresenham line algorithm between A and B points
+                and convert all points on the way to the class of A/B.
+                */
+                bresenham2D(i, j, i1, j1, MASK_temp, MASK_upd, (long)(dimX), (long)(dimY));
+              }
+            }
+          }}
+      }
+  return *MASK_upd;
+}
+
 /* MASKED-constrained 2D linear diffusion (PDE heat equation) */
 float LinearDiff_MASK2D(float *Input, unsigned char *MASK, float *Output,  float *Eucl_Vec, int DiffusWindow, float lambdaPar, float tau, long dimX, long dimY)
 {
@@ -142,18 +229,18 @@ float diffVal;
 
 #pragma omp parallel for shared(Input) private(index,i,j,i1,j1,i_m,j_m,counter,diffVal,indexneighb,class_c,class_n)
     for(i=0; i<dimX; i++) {
-        for(j=0; j<dimY; j++) {        
+        for(j=0; j<dimY; j++) {
         index = j*dimX+i; /* current pixel index */
         counter = 0; diffVal = 0.0f;
         for(i_m=-DiffusWindow; i_m<=DiffusWindow; i_m++) {
-            for(j_m=-DiffusWindow; j_m<=DiffusWindow; j_m++) {		
+            for(j_m=-DiffusWindow; j_m<=DiffusWindow; j_m++) {
             i1 = i+i_m;
             j1 = j+j_m;
-            if (((i1 >= 0) && (i1 < dimX)) && ((j1 >= 0) && (j1 < dimY))) {                   
+            if (((i1 >= 0) && (i1 < dimX)) && ((j1 >= 0) && (j1 < dimY))) {
             indexneighb = j1*dimX+i1; /* neighbour pixel index */
 	    class_c = MASK[index]; /* current class value */
     	    class_n = MASK[indexneighb]; /* neighbour class value */
-	    
+
 	    /* perform diffusion only within the same class (given by MASK) */
 	    if (class_n == class_c) diffVal += Output[indexneighb] - Output[index];
             }
@@ -170,21 +257,21 @@ float NonLinearDiff_MASK2D(float *Input, unsigned char *MASK, float *Output, flo
 	long i,j,i1,j1,i_m,j_m,index,indexneighb,counter;
 	unsigned char class_c, class_n;
 	float diffVal, funcVal;
-	
+
 #pragma omp parallel for shared(Input) private(index,i,j,i1,j1,i_m,j_m,counter,diffVal,funcVal,indexneighb,class_c,class_n)
     for(i=0; i<dimX; i++) {
-        for(j=0; j<dimY; j++) {        
+        for(j=0; j<dimY; j++) {
         index = j*dimX+i; /* current pixel index */
         counter = 0; diffVal = 0.0f; funcVal = 0.0f;
         for(i_m=-DiffusWindow; i_m<=DiffusWindow; i_m++) {
-            for(j_m=-DiffusWindow; j_m<=DiffusWindow; j_m++) {		
+            for(j_m=-DiffusWindow; j_m<=DiffusWindow; j_m++) {
             i1 = i+i_m;
             j1 = j+j_m;
-            if (((i1 >= 0) && (i1 < dimX)) && ((j1 >= 0) && (j1 < dimY))) {                   
+            if (((i1 >= 0) && (i1 < dimX)) && ((j1 >= 0) && (j1 < dimY))) {
             indexneighb = j1*dimX+i1; /* neighbour pixel index */
 	    class_c = MASK[index]; /* current class value */
     	    class_n = MASK[indexneighb]; /* neighbour class value */
-	    
+
 	    /* perform diffusion only within the same class (given by MASK) */
 	    if (class_n == class_c) {
 	    	diffVal = Output[indexneighb] - Output[index];
@@ -200,7 +287,7 @@ float NonLinearDiff_MASK2D(float *Input, unsigned char *MASK, float *Output, flo
   		if (fabs(diffVal) <= sigmaPar) funcVal += diffVal*powf((1.0f - powf((diffVal/sigmaPar),2)), 2); }
                 else {
                 printf("%s \n", "No penalty function selected! Use Huber,2 or 3.");
-		break; }  		  
+		break; }
             		}
             	}
 		counter++;
@@ -212,3 +299,90 @@ float NonLinearDiff_MASK2D(float *Input, unsigned char *MASK, float *Output, flo
 /********************************************************************/
 /***************************3D Functions*****************************/
 /********************************************************************/
+
+
+int bresenham2D(int i, int j, int i1, int j1, unsigned char *MASK, unsigned char *MASK_upd, long dimX, long dimY) {
+
+                   int x[] = {i, i1};
+                   int y[] = {j, j1};
+                   int steep = (fabs(y[1]-y[0]) > fabs(x[1]-x[0]));
+                   int ystep = 0;
+
+                   //printf("[%i][%i][%i][%i]\n", x[1], y[1], steep, kk) ;
+
+                   //if (steep == 1) {swap(x[0],y[0]); swap(x[1],y[1]);}
+
+                   if (steep == 1) {
+
+                    ///swaping
+                   int a, b;
+
+                   a = x[0];
+                   b = y[0];
+                   x[0] = b;
+                   y[0] = a;
+
+                   a = x[1];
+                   b = y[1];
+                   x[1] = b;
+                   y[1] = a;
+                   }
+
+                   if (x[0] > x[1]) {
+                   int a, b;
+                   a = x[0];
+                   b = x[1];
+                   x[0] = b;
+                   x[1] = a;
+
+                   a = y[0];
+                   b = y[1];
+                   y[0] = b;
+                   y[1] = a;
+                   } //(x[0] > x[1])
+
+                  int delx = x[1]-x[0];
+                  int dely = fabs(y[1]-y[0]);
+                  int error = 0;
+                  int x_n = x[0];
+                  int y_n = y[0];
+                  if (y[0] < y[1]) {ystep = 1;}
+                  else {ystep = -1;}
+
+                  for(int n = 0; n<delx+1; n++) {
+                       if (steep == 1) {
+                        /*
+                        X_new[n] = x_n;
+                        Y_new[n] = y_n;
+                        */
+                        /*printf("[%i][%i][%u]\n", x_n, y_n, MASK[y_n*dimX+x_n]);*/
+
+                        MASK_upd[y_n*dimX+x_n] = 10;
+                        //if (MASK[j*dimX+i] != MASK[y_n*dimX+x_n]) {
+                        //  MASK_upd[y_n*dimX+x_n] = MASK[j*dimX+i];
+                        //}
+                       }
+                       else {
+                         /*
+                        X_new[n] = y_n;
+                        Y_new[n] = x_n;
+                        */
+                        // printf("[%i][%i][%u]\n", y_n, x_n, MASK[x_n*dimX+y_n]);
+                        // MASK_upd[x_n*dimX+y_n] = 1;
+                        //if (MASK[j*dimX+i] != MASK[x_n*dimX+y_n]) {
+                        //  MASK_upd[x_n*dimX+y_n] = MASK[j*dimX+i];
+                        //}
+                        MASK_upd[x_n*dimX+y_n] = 20;
+                       }
+                       x_n = x_n + 1;
+                       error = error + dely;
+
+                       if (2*error >= delx) {
+                          y_n = y_n + ystep;
+                         error = error - delx;
+                       } // (2*error >= delx)
+                       //printf("[%i][%i][%i]\n", X_new[n], Y_new[n], n) ;
+                  } // for(int n = 0; n<delx+1; n++)
+
+                  return 0;
+}
diff --git a/src/Core/regularisers_CPU/DiffusionMASK_core.h b/src/Core/regularisers_CPU/DiffusionMASK_core.h
index 8890c736..48142af1 100644
--- a/src/Core/regularisers_CPU/DiffusionMASK_core.h
+++ b/src/Core/regularisers_CPU/DiffusionMASK_core.h
@@ -33,7 +33,7 @@ limitations under the License.
  * Input Parameters:
  * 1. Noisy image/volume
  * 2. MASK (in unsigned short format)
- * 3. Diffusivity window (half-size of the searching window, e.g. 3) 	 
+ * 3. Diffusivity window (half-size of the searching window, e.g. 3)
  * 4. lambda - regularization parameter
  * 5. Edge-preserving parameter (sigma), when sigma equals to zero nonlinear diffusion -> linear diffusion
  * 6. Number of iterations, for explicit scheme >= 150 is recommended
@@ -54,9 +54,12 @@ limitations under the License.
 #ifdef __cplusplus
 extern "C" {
 #endif
-CCPI_EXPORT float DiffusionMASK_CPU_main(float *Input, unsigned char *MASK, float *Output, float *infovector, int DiffusWindow, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int penaltytype, float epsil, int dimX, int dimY, int dimZ);
+CCPI_EXPORT float DiffusionMASK_CPU_main(float *Input, unsigned char *MASK, unsigned char *MASK_upd, float *Output, float *infovector, int DiffusWindow, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int penaltytype, float epsil, int dimX, int dimY, int dimZ);
 CCPI_EXPORT float LinearDiff_MASK2D(float *Input, unsigned char *MASK, float *Output,  float *Eucl_Vec, int DiffusWindow, float lambdaPar, float tau, long dimX, long dimY);
 CCPI_EXPORT float NonLinearDiff_MASK2D(float *Input, unsigned char *MASK, float *Output, float *Eucl_Vec, int DiffusWindow, float lambdaPar, float sigmaPar, float tau, int penaltytype, long dimX, long dimY);
+CCPI_EXPORT float OutiersRemoval2D(unsigned char *MASK, unsigned char *MASK_upd, long i, long j, long dimX, long dimY);
+CCPI_EXPORT float Mask_update2D(unsigned char *MASK_temp, unsigned char *MASK_upd, long i, long j, int DiffusWindow, long dimX, long dimY);
+CCPI_EXPORT int bresenham2D(int i, int j, int i1, int j1, unsigned char *MASK, unsigned char *MASK_upd, long dimX, long dimY);
 #ifdef __cplusplus
 }
 #endif
diff --git a/src/Python/ccpi/filters/regularisers.py b/src/Python/ccpi/filters/regularisers.py
index 1e427bf3..00afcc20 100644
--- a/src/Python/ccpi/filters/regularisers.py
+++ b/src/Python/ccpi/filters/regularisers.py
@@ -127,10 +127,11 @@ def NDF(inputData, regularisation_parameter, edge_parameter, iterations,
     	    raise ValueError ('GPU is not available')
         raise ValueError('Unknown device {0}. Expecting gpu or cpu'\
                          .format(device))
-def NDF_MASK(inputData, diffuswindow, regularisation_parameter, edge_parameter, iterations,
+def NDF_MASK(inputData, maskdata, diffuswindow, regularisation_parameter, edge_parameter, iterations,
                      time_marching_parameter, penalty_type, tolerance_param, device='cpu'):
     if device == 'cpu':
         return NDF_MASK_CPU(inputData,
+                     maskdata,
                      diffuswindow,
                      regularisation_parameter,
                      edge_parameter,
@@ -140,6 +141,7 @@ def NDF_MASK(inputData, diffuswindow, regularisation_parameter, edge_parameter,
                      tolerance_param)
     elif device == 'gpu' and gpu_enabled:
         return NDF_MASK_CPU(inputData,
+                     maskdata,
                      diffuswindow,
                      regularisation_parameter,
                      edge_parameter,
diff --git a/src/Python/src/cpu_regularisers.pyx b/src/Python/src/cpu_regularisers.pyx
index 305ee1ff..ca402c19 100644
--- a/src/Python/src/cpu_regularisers.pyx
+++ b/src/Python/src/cpu_regularisers.pyx
@@ -24,7 +24,7 @@ cdef extern float SB_TV_CPU_main(float *Input, float *Output, float *infovector,
 cdef extern float LLT_ROF_CPU_main(float *Input, float *Output, float *infovector, float lambdaROF, float lambdaLLT, int iterationsNumb, float tau, float epsil, int dimX, int dimY, int dimZ);
 cdef extern float TGV_main(float *Input, float *Output, float *infovector, float lambdaPar, float alpha1, float alpha0, int iterationsNumb, float L2, float epsil, int dimX, int dimY, int dimZ);
 cdef extern float Diffusion_CPU_main(float *Input, float *Output, float *infovector, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int penaltytype, float epsil, int dimX, int dimY, int dimZ);
-cdef extern float DiffusionMASK_CPU_main(float *Input, unsigned char *MASK, float *Output, float *infovector, int DiffusWindow, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int penaltytype, float epsil, int dimX, int dimY, int dimZ);
+cdef extern float DiffusionMASK_CPU_main(float *Input, unsigned char *MASK, unsigned char *MASK_upd, float *Output, float *infovector, int DiffusWindow, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int penaltytype, float epsil, int dimX, int dimY, int dimZ);
 cdef extern float Diffus4th_CPU_main(float *Input, float *Output,  float *infovector, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, float epsil, int dimX, int dimY, int dimZ);
 cdef extern float dTV_FGP_CPU_main(float *Input, float *InputRef, float *Output, float *infovector, float lambdaPar, int iterationsNumb, float epsil, float eta, int methodTV, int nonneg, int dimX, int dimY, int dimZ);
 cdef extern float TNV_CPU_main(float *Input, float *u, float lambdaPar, int maxIter, float tol, int dimX, int dimY, int dimZ);
@@ -403,18 +403,22 @@ def NDF_MASK_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
     dims[0] = inputData.shape[0]
     dims[1] = inputData.shape[1]
 
+
+    cdef np.ndarray[np.uint8_t, ndim=2, mode="c"] mask_upd = \
+            np.zeros([dims[0],dims[1]], dtype='uint8')
     cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
             np.zeros([dims[0],dims[1]], dtype='float32')
     cdef np.ndarray[np.float32_t, ndim=1, mode="c"] infovec = \
                 np.zeros([2], dtype='float32')
 
+
     # Run constrained nonlinear diffusion iterations for 2D data
-    DiffusionMASK_CPU_main(&inputData[0,0], &maskData[0,0], &outputData[0,0], &infovec[0],
+    DiffusionMASK_CPU_main(&inputData[0,0], &maskData[0,0], &mask_upd[0,0], &outputData[0,0], &infovec[0],
     diffuswindow, regularisation_parameter, edge_parameter, iterationsNumb,
     time_marching_parameter, penalty_type,
     tolerance_param,
     dims[1], dims[0], 1)
-    return (outputData,infovec)
+    return (mask_upd,outputData,infovec)
 
 #****************************************************************#
 #*************Anisotropic Fourth-Order diffusion*****************#

From f736660128290df50aa9a431dfa40e96da600f0d Mon Sep 17 00:00:00 2001
From: dkazanc <dkazanc@hotmail.com>
Date: Thu, 11 Apr 2019 15:58:36 +0100
Subject: [PATCH 3/7] further progress

---
 .../regularisers_CPU/DiffusionMASK_core.c     | 106 ++++++++++++------
 .../regularisers_CPU/DiffusionMASK_core.h     |   2 +-
 src/Python/ccpi/filters/regularisers.py       |   6 +-
 src/Python/src/cpu_regularisers.pyx           |  15 ++-
 4 files changed, 86 insertions(+), 43 deletions(-)

diff --git a/src/Core/regularisers_CPU/DiffusionMASK_core.c b/src/Core/regularisers_CPU/DiffusionMASK_core.c
index 2af3cb66..9d128bc6 100644
--- a/src/Core/regularisers_CPU/DiffusionMASK_core.c
+++ b/src/Core/regularisers_CPU/DiffusionMASK_core.c
@@ -53,7 +53,14 @@ int signNDF_m(float x) {
  * [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 DiffusionMASK_CPU_main(float *Input, unsigned char *MASK, unsigned char *MASK_upd, float *Output, float *infovector, int DiffusWindow, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int penaltytype, float epsil, int dimX, int dimY, int dimZ)
+void swapVAL(unsigned char *xp, unsigned char *yp) 
+{ 
+    unsigned char temp = *xp; 
+    *xp = *yp; 
+    *yp = temp; 
+} 
+
+float DiffusionMASK_CPU_main(float *Input, unsigned char *MASK, unsigned char *MASK_upd, unsigned char *SelClassesList, int SelClassesList_length, float *Output, float *infovector, int classesNumb, int DiffusWindow, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int penaltytype, float epsil, int dimX, int dimY, int dimZ)
 {
     long i,j,k;
     int counterG;
@@ -64,9 +71,43 @@ float DiffusionMASK_CPU_main(float *Input, unsigned char *MASK, unsigned char *M
     float re, re1;
     re = 0.0f; re1 = 0.0f;
     int count = 0;
-    unsigned char *MASK_temp;
+    unsigned char *MASK_temp, *ClassesList, CurrClass, temp;
     DimTotal = (long)(dimX*dimY*dimZ);
 
+    /* defines the list for all classes in the mask */
+    ClassesList = (unsigned char*) calloc (classesNumb,sizeof(unsigned char));
+    /* indentify all classes in the MASK */
+    /* find the smaller class value first */
+    /*
+    unsigned char smallestClass;
+    smallestClass = MASK[0];
+    for(i=0; i<DimTotal; i++) {
+	if (MASK[i] < smallestClass) smallestClass = MASK[i];
+    }
+    */
+    MASK_temp = (unsigned char*) calloc (DimTotal,sizeof(unsigned char));   
+    copyIm_unchar(MASK, MASK_temp, (long)(dimX), (long)(dimY), (long)(dimZ));        
+     
+    int x, y;
+    for(x=0; x<classesNumb; x++)	{
+                for(y=0; y<classesNumb-1; y++) {
+                    if(MASK_temp[y] > MASK_temp[y+1]) {
+                        temp = MASK_temp[y+1];
+                        MASK_temp[y+1] = MASK_temp[y];
+                        MASK_temp[y] = temp;
+                    }}} 
+     
+     //printf("[%u]\n", MASK_temp[0]);
+    
+     CurrClass =  MASK_temp[0]; ClassesList[0]= MASK_temp[0]; counterG = 0;
+     for(i=0; i<DimTotal; i++) {
+        if (MASK_temp[i] > CurrClass) {
+        CurrClass = MASK_temp[i];
+        ClassesList[counterG] = MASK_temp[i];
+        printf("[%u]\n", ClassesList[counterG]);
+        counterG++; }    
+     }
+
     if (dimZ == 1) {
 	DiffusWindow_tot = (2*DiffusWindow + 1)*(2*DiffusWindow + 1);
         /* generate a 2D Gaussian kernel for NLM procedure */
@@ -92,7 +133,7 @@ float DiffusionMASK_CPU_main(float *Input, unsigned char *MASK, unsigned char *M
 
     if (epsil != 0.0f) Output_prev = calloc(DimTotal, sizeof(float));
 
-    MASK_temp = (unsigned char*) calloc (DimTotal,sizeof(unsigned char));
+    
 
     /* copy input into output */
     copyIm(Input, Output, (long)(dimX), (long)(dimY), (long)(dimZ));
@@ -110,25 +151,27 @@ float DiffusionMASK_CPU_main(float *Input, unsigned char *MASK, unsigned char *M
     /* copy the updated MASK (clean of outliers) */
     copyIm_unchar(MASK_upd, MASK_temp, (long)(dimX), (long)(dimY), (long)(dimZ));
 
-//    #pragma omp parallel for shared(MASK_temp,MASK_upd) private(i,j)
-//  for(i=0; i<dimX; i++) {
-//        for(j=0; j<dimY; j++) {
-//      if (MASK_temp[j*dimX+i] == MASK[j*dimX+i]) {
-        /*the class of the central pixel has not changed, i.e. the central pixel is not an outlier -> continue */
-        i = 162; j = 258;
+    #pragma omp parallel for shared(MASK_temp,MASK_upd) private(i,j)
+    for(i=0; i<dimX; i++) {
+        for(j=0; j<dimY; j++) {
+      if (MASK_temp[j*dimX+i] == MASK[j*dimX+i]) {
+	/* !One needs to work with a specific class to avoid overlaps! 
+	 hence it is crucial to establish relevant classes */
+       if (MASK_temp[j*dimX+i] == 149) {
+        /* The class of the central pixel has not changed, i.e. the central pixel is not an outlier -> continue */
+        /* i = 258; j = 165; */
         Mask_update2D(MASK_temp, MASK_upd, i, j, DiffusWindow, (long)(dimX), (long)(dimY));
-//       }
-      //}}
+       }}
+     }}
     }
 
-    /* start diffusivity iterations */
+    /* The mask has been processed, start diffusivity iterations */
     for(i=0; i < iterationsNumb; i++) {
-
       if ((epsil != 0.0f)  && (i % 5 == 0)) copyIm(Output, Output_prev, (long)(dimX), (long)(dimY), (long)(dimZ));
       if (dimZ == 1) {
-	          /* running 2D diffusion iterations */
-            //if (sigmaPar == 0.0f) LinearDiff_MASK2D(Input, MASK, Output, Eucl_Vec, DiffusWindow, lambdaPar, tau, (long)(dimX), (long)(dimY)); /* constrained linear diffusion */
-            //else NonLinearDiff_MASK2D(Input, MASK, Output, Eucl_Vec, DiffusWindow, lambdaPar, sigmaPar2, tau, penaltytype, (long)(dimX), (long)(dimY)); /* constrained nonlinear diffusion */
+             /* running 2D diffusion iterations */
+            if (sigmaPar == 0.0f) LinearDiff_MASK2D(Input, MASK_upd, Output, Eucl_Vec, DiffusWindow, lambdaPar, tau, (long)(dimX), (long)(dimY)); /* constrained linear diffusion */
+            else NonLinearDiff_MASK2D(Input, MASK_upd, Output, Eucl_Vec, DiffusWindow, lambdaPar, sigmaPar2, tau, penaltytype, (long)(dimX), (long)(dimY)); /* constrained nonlinear diffusion */
           }
       else {
        	/* running 3D diffusion iterations */
@@ -207,11 +250,10 @@ float Mask_update2D(unsigned char *MASK_temp, unsigned char *MASK_upd, long i, l
             j1 = j+j_m;
             if (((i1 >= 0) && (i1 < dimX)) && ((j1 >= 0) && (j1 < dimY))) {
               if (MASK_temp[j*dimX+i] == MASK_temp[j1*dimX+i1]) {
-                /*A and B points belong to the same class, do STEP 4*/
+                /* A and B points belong to the same class, do STEP 4*/
                 /* STEP 4: Run Bresenham line algorithm between A and B points
-                and convert all points on the way to the class of A/B.
-                */
-                bresenham2D(i, j, i1, j1, MASK_temp, MASK_upd, (long)(dimX), (long)(dimY));
+                and convert all points on the way to the class of A/B.  */
+               bresenham2D(i, j, i1, j1, MASK_temp, MASK_upd, (long)(dimX), (long)(dimY));
               }
             }
           }}
@@ -301,20 +343,20 @@ float NonLinearDiff_MASK2D(float *Input, unsigned char *MASK, float *Output, flo
 /********************************************************************/
 
 
-int bresenham2D(int i, int j, int i1, int j1, unsigned char *MASK, unsigned char *MASK_upd, long dimX, long dimY) {
-
+int bresenham2D(int i, int j, int i1, int j1, unsigned char *MASK, unsigned char *MASK_upd, long dimX, long dimY) 
+{               
+                   int n;
                    int x[] = {i, i1};
                    int y[] = {j, j1};
                    int steep = (fabs(y[1]-y[0]) > fabs(x[1]-x[0]));
                    int ystep = 0;
 
                    //printf("[%i][%i][%i][%i]\n", x[1], y[1], steep, kk) ;
-
                    //if (steep == 1) {swap(x[0],y[0]); swap(x[1],y[1]);}
 
                    if (steep == 1) {
 
-                    ///swaping
+                   // swaping
                    int a, b;
 
                    a = x[0];
@@ -349,18 +391,15 @@ int bresenham2D(int i, int j, int i1, int j1, unsigned char *MASK, unsigned char
                   if (y[0] < y[1]) {ystep = 1;}
                   else {ystep = -1;}
 
-                  for(int n = 0; n<delx+1; n++) {
+                  for(n = 0; n<delx+1; n++) {
                        if (steep == 1) {
                         /*
                         X_new[n] = x_n;
                         Y_new[n] = y_n;
                         */
                         /*printf("[%i][%i][%u]\n", x_n, y_n, MASK[y_n*dimX+x_n]);*/
-
-                        MASK_upd[y_n*dimX+x_n] = 10;
-                        //if (MASK[j*dimX+i] != MASK[y_n*dimX+x_n]) {
-                        //  MASK_upd[y_n*dimX+x_n] = MASK[j*dimX+i];
-                        //}
+                        // MASK_upd[x_n*dimX+y_n] = 10;                        
+                        if (MASK[j*dimX+i] != MASK[x_n*dimX+y_n]) MASK_upd[x_n*dimX+y_n] = MASK[j*dimX+i];
                        }
                        else {
                          /*
@@ -368,11 +407,8 @@ int bresenham2D(int i, int j, int i1, int j1, unsigned char *MASK, unsigned char
                         Y_new[n] = x_n;
                         */
                         // printf("[%i][%i][%u]\n", y_n, x_n, MASK[x_n*dimX+y_n]);
-                        // MASK_upd[x_n*dimX+y_n] = 1;
-                        //if (MASK[j*dimX+i] != MASK[x_n*dimX+y_n]) {
-                        //  MASK_upd[x_n*dimX+y_n] = MASK[j*dimX+i];
-                        //}
-                        MASK_upd[x_n*dimX+y_n] = 20;
+                        // MASK_upd[y_n*dimX+x_n] = 20;
+                        if (MASK[j*dimX+i] != MASK[y_n*dimX+x_n]) MASK_upd[y_n*dimX+x_n] = MASK[j*dimX+i];
                        }
                        x_n = x_n + 1;
                        error = error + dely;
diff --git a/src/Core/regularisers_CPU/DiffusionMASK_core.h b/src/Core/regularisers_CPU/DiffusionMASK_core.h
index 48142af1..a032905d 100644
--- a/src/Core/regularisers_CPU/DiffusionMASK_core.h
+++ b/src/Core/regularisers_CPU/DiffusionMASK_core.h
@@ -54,7 +54,7 @@ limitations under the License.
 #ifdef __cplusplus
 extern "C" {
 #endif
-CCPI_EXPORT float DiffusionMASK_CPU_main(float *Input, unsigned char *MASK, unsigned char *MASK_upd, float *Output, float *infovector, int DiffusWindow, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int penaltytype, float epsil, int dimX, int dimY, int dimZ);
+CCPI_EXPORT float DiffusionMASK_CPU_main(float *Input, unsigned char *MASK, unsigned char *MASK_upd, unsigned char *SelClassesList, int SelClassesList_length, float *Output, float *infovector, int classesNumb, int DiffusWindow, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int penaltytype, float epsil, int dimX, int dimY, int dimZ);
 CCPI_EXPORT float LinearDiff_MASK2D(float *Input, unsigned char *MASK, float *Output,  float *Eucl_Vec, int DiffusWindow, float lambdaPar, float tau, long dimX, long dimY);
 CCPI_EXPORT float NonLinearDiff_MASK2D(float *Input, unsigned char *MASK, float *Output, float *Eucl_Vec, int DiffusWindow, float lambdaPar, float sigmaPar, float tau, int penaltytype, long dimX, long dimY);
 CCPI_EXPORT float OutiersRemoval2D(unsigned char *MASK, unsigned char *MASK_upd, long i, long j, long dimX, long dimY);
diff --git a/src/Python/ccpi/filters/regularisers.py b/src/Python/ccpi/filters/regularisers.py
index 00afcc20..610907d3 100644
--- a/src/Python/ccpi/filters/regularisers.py
+++ b/src/Python/ccpi/filters/regularisers.py
@@ -127,11 +127,13 @@ def NDF(inputData, regularisation_parameter, edge_parameter, iterations,
     	    raise ValueError ('GPU is not available')
         raise ValueError('Unknown device {0}. Expecting gpu or cpu'\
                          .format(device))
-def NDF_MASK(inputData, maskdata, diffuswindow, regularisation_parameter, edge_parameter, iterations,
+def NDF_MASK(inputData, maskdata, select_classes, total_classesNum, diffuswindow, regularisation_parameter, edge_parameter, iterations,
                      time_marching_parameter, penalty_type, tolerance_param, device='cpu'):
     if device == 'cpu':
         return NDF_MASK_CPU(inputData,
                      maskdata,
+                     select_classes,
+                     total_classesNum,
                      diffuswindow,
                      regularisation_parameter,
                      edge_parameter,
@@ -142,6 +144,8 @@ def NDF_MASK(inputData, maskdata, diffuswindow, regularisation_parameter, edge_p
     elif device == 'gpu' and gpu_enabled:
         return NDF_MASK_CPU(inputData,
                      maskdata,
+                     select_classes,
+                     total_classesNum,
                      diffuswindow,
                      regularisation_parameter,
                      edge_parameter,
diff --git a/src/Python/src/cpu_regularisers.pyx b/src/Python/src/cpu_regularisers.pyx
index ca402c19..78c46e2f 100644
--- a/src/Python/src/cpu_regularisers.pyx
+++ b/src/Python/src/cpu_regularisers.pyx
@@ -24,7 +24,7 @@ cdef extern float SB_TV_CPU_main(float *Input, float *Output, float *infovector,
 cdef extern float LLT_ROF_CPU_main(float *Input, float *Output, float *infovector, float lambdaROF, float lambdaLLT, int iterationsNumb, float tau, float epsil, int dimX, int dimY, int dimZ);
 cdef extern float TGV_main(float *Input, float *Output, float *infovector, float lambdaPar, float alpha1, float alpha0, int iterationsNumb, float L2, float epsil, int dimX, int dimY, int dimZ);
 cdef extern float Diffusion_CPU_main(float *Input, float *Output, float *infovector, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int penaltytype, float epsil, int dimX, int dimY, int dimZ);
-cdef extern float DiffusionMASK_CPU_main(float *Input, unsigned char *MASK, unsigned char *MASK_upd, float *Output, float *infovector, int DiffusWindow, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int penaltytype, float epsil, int dimX, int dimY, int dimZ);
+cdef extern float DiffusionMASK_CPU_main(float *Input, unsigned char *MASK, unsigned char *MASK_upd, unsigned char *SelClassesList, int SelClassesList_length, float *Output, float *infovector, int classesNumb, int DiffusWindow, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int penaltytype, float epsil, int dimX, int dimY, int dimZ);
 cdef extern float Diffus4th_CPU_main(float *Input, float *Output,  float *infovector, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, float epsil, int dimX, int dimY, int dimZ);
 cdef extern float dTV_FGP_CPU_main(float *Input, float *InputRef, float *Output, float *infovector, float lambdaPar, int iterationsNumb, float epsil, float eta, int methodTV, int nonneg, int dimX, int dimY, int dimZ);
 cdef extern float TNV_CPU_main(float *Input, float *u, float lambdaPar, int maxIter, float tol, int dimX, int dimY, int dimZ);
@@ -384,14 +384,16 @@ def NDF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
 #****************************************************************#
 #********Constrained Nonlinear(Isotropic) Diffusion**************#
 #****************************************************************#
-def NDF_MASK_CPU(inputData, maskData, diffuswindow, regularisation_parameter, edge_parameter, iterationsNumb,time_marching_parameter, penalty_type, tolerance_param):
+def NDF_MASK_CPU(inputData, maskData, select_classes, total_classesNum, diffuswindow, regularisation_parameter, edge_parameter, iterationsNumb,time_marching_parameter, penalty_type, tolerance_param):
     if inputData.ndim == 2:
-        return NDF_MASK_2D(inputData, maskData, diffuswindow, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type, tolerance_param)
+        return NDF_MASK_2D(inputData, maskData, select_classes, total_classesNum, diffuswindow, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type, tolerance_param)
     elif inputData.ndim == 3:
         return 0
 
 def NDF_MASK_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=1, mode="c"] select_classes,
+                     int total_classesNum,
                      int diffuswindow,
                      float regularisation_parameter,
                      float edge_parameter,
@@ -403,7 +405,8 @@ def NDF_MASK_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
     dims[0] = inputData.shape[0]
     dims[1] = inputData.shape[1]
 
-
+    select_classes_length = select_classes.shape[0]
+    
     cdef np.ndarray[np.uint8_t, ndim=2, mode="c"] mask_upd = \
             np.zeros([dims[0],dims[1]], dtype='uint8')
     cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
@@ -413,8 +416,8 @@ def NDF_MASK_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
 
 
     # Run constrained nonlinear diffusion iterations for 2D data
-    DiffusionMASK_CPU_main(&inputData[0,0], &maskData[0,0], &mask_upd[0,0], &outputData[0,0], &infovec[0],
-    diffuswindow, regularisation_parameter, edge_parameter, iterationsNumb,
+    DiffusionMASK_CPU_main(&inputData[0,0], &maskData[0,0], &mask_upd[0,0], &select_classes[0], select_classes_length, &outputData[0,0], &infovec[0],
+    total_classesNum, diffuswindow, regularisation_parameter, edge_parameter, iterationsNumb,
     time_marching_parameter, penalty_type,
     tolerance_param,
     dims[1], dims[0], 1)

From 2e67930b6df118fb0dca129d3cdcf5515bf4827c Mon Sep 17 00:00:00 2001
From: Daniil Kazantsev <dkazanc@hotmail.com>
Date: Thu, 11 Apr 2019 22:05:35 +0100
Subject: [PATCH 4/7] sorting classes fixed

---
 .../regularisers_CPU/DiffusionMASK_core.c     | 88 +++++++++----------
 1 file changed, 44 insertions(+), 44 deletions(-)

diff --git a/src/Core/regularisers_CPU/DiffusionMASK_core.c b/src/Core/regularisers_CPU/DiffusionMASK_core.c
index 9d128bc6..dd4b1b48 100644
--- a/src/Core/regularisers_CPU/DiffusionMASK_core.c
+++ b/src/Core/regularisers_CPU/DiffusionMASK_core.c
@@ -53,17 +53,17 @@ int signNDF_m(float x) {
  * [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 swapVAL(unsigned char *xp, unsigned char *yp) 
-{ 
-    unsigned char temp = *xp; 
-    *xp = *yp; 
-    *yp = temp; 
-} 
+void swapVAL(unsigned char *xp, unsigned char *yp)
+{
+    unsigned char temp = *xp;
+    *xp = *yp;
+    *yp = temp;
+}
 
 float DiffusionMASK_CPU_main(float *Input, unsigned char *MASK, unsigned char *MASK_upd, unsigned char *SelClassesList, int SelClassesList_length, float *Output, float *infovector, int classesNumb, int DiffusWindow, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int penaltytype, float epsil, int dimX, int dimY, int dimZ)
 {
     long i,j,k;
-    int counterG;
+    int counterG, switcher;
     float sigmaPar2, *Output_prev=NULL, *Eucl_Vec;
     int DiffusWindow_tot;
     sigmaPar2 = sigmaPar/sqrt(2.0f);
@@ -76,38 +76,38 @@ float DiffusionMASK_CPU_main(float *Input, unsigned char *MASK, unsigned char *M
 
     /* defines the list for all classes in the mask */
     ClassesList = (unsigned char*) calloc (classesNumb,sizeof(unsigned char));
-    /* indentify all classes in the MASK */
-    /* find the smaller class value first */
-    /*
-    unsigned char smallestClass;
-    smallestClass = MASK[0];
-    for(i=0; i<DimTotal; i++) {
-	if (MASK[i] < smallestClass) smallestClass = MASK[i];
-    }
-    */
-    MASK_temp = (unsigned char*) calloc (DimTotal,sizeof(unsigned char));   
-    copyIm_unchar(MASK, MASK_temp, (long)(dimX), (long)(dimY), (long)(dimZ));        
-     
-    int x, y;
-    for(x=0; x<classesNumb; x++)	{
-                for(y=0; y<classesNumb-1; y++) {
-                    if(MASK_temp[y] > MASK_temp[y+1]) {
-                        temp = MASK_temp[y+1];
-                        MASK_temp[y+1] = MASK_temp[y];
-                        MASK_temp[y] = temp;
-                    }}} 
-     
-     //printf("[%u]\n", MASK_temp[0]);
-    
-     CurrClass =  MASK_temp[0]; ClassesList[0]= MASK_temp[0]; counterG = 0;
-     for(i=0; i<DimTotal; i++) {
-        if (MASK_temp[i] > CurrClass) {
-        CurrClass = MASK_temp[i];
-        ClassesList[counterG] = MASK_temp[i];
-        printf("[%u]\n", ClassesList[counterG]);
-        counterG++; }    
-     }
 
+     /* find which classes (values) are present in the segmented data */
+     CurrClass =  MASK[0]; ClassesList[0]= MASK[0]; counterG = 1;
+     for(i=0; i<DimTotal; i++) {
+       if (MASK[i] != CurrClass) {
+          switcher = 1;
+          for(j=0; j<counterG; j++) {
+            if (ClassesList[j] == MASK[i]) {
+              switcher = 0;
+              break;
+            }}
+            if (switcher == 1) {
+                CurrClass = MASK[i];
+                ClassesList[counterG] = MASK[i];
+                /*printf("[%u]\n", ClassesList[counterG]);*/
+                counterG++;
+              }
+        }
+        if (counterG == classesNumb) break;
+      }
+      /* sort the obtained values (classes) */
+      for(i=0; i<classesNumb; i++)	{
+                  for(j=0; j<classesNumb-1; j++) {
+                      if(ClassesList[j] > ClassesList[j+1]) {
+                          temp = ClassesList[j+1];
+                          ClassesList[j+1] = ClassesList[j];
+                          ClassesList[j] = temp;
+                      }}}
+
+    for(i=0; i<classesNumb; i++)	printf("[%u]\n", ClassesList[i]);
+
+    /*Euclidian weight for diffisuvuty window*/
     if (dimZ == 1) {
 	DiffusWindow_tot = (2*DiffusWindow + 1)*(2*DiffusWindow + 1);
         /* generate a 2D Gaussian kernel for NLM procedure */
@@ -133,11 +133,11 @@ float DiffusionMASK_CPU_main(float *Input, unsigned char *MASK, unsigned char *M
 
     if (epsil != 0.0f) Output_prev = calloc(DimTotal, sizeof(float));
 
-    
+    MASK_temp = (unsigned char*) calloc (DimTotal,sizeof(unsigned char));
 
     /* copy input into output */
     copyIm(Input, Output, (long)(dimX), (long)(dimY), (long)(dimZ));
-    /* copy given MASK  */
+    /* copy given MASK to MASK_upd*/
     copyIm_unchar(MASK, MASK_upd, (long)(dimX), (long)(dimY), (long)(dimZ));
 
     /********************** PERFORM MASK PROCESSING ************************/
@@ -155,7 +155,7 @@ float DiffusionMASK_CPU_main(float *Input, unsigned char *MASK, unsigned char *M
     for(i=0; i<dimX; i++) {
         for(j=0; j<dimY; j++) {
       if (MASK_temp[j*dimX+i] == MASK[j*dimX+i]) {
-	/* !One needs to work with a specific class to avoid overlaps! 
+	/* !One needs to work with a specific class to avoid overlaps!
 	 hence it is crucial to establish relevant classes */
        if (MASK_temp[j*dimX+i] == 149) {
         /* The class of the central pixel has not changed, i.e. the central pixel is not an outlier -> continue */
@@ -343,8 +343,8 @@ float NonLinearDiff_MASK2D(float *Input, unsigned char *MASK, float *Output, flo
 /********************************************************************/
 
 
-int bresenham2D(int i, int j, int i1, int j1, unsigned char *MASK, unsigned char *MASK_upd, long dimX, long dimY) 
-{               
+int bresenham2D(int i, int j, int i1, int j1, unsigned char *MASK, unsigned char *MASK_upd, long dimX, long dimY)
+{
                    int n;
                    int x[] = {i, i1};
                    int y[] = {j, j1};
@@ -398,7 +398,7 @@ int bresenham2D(int i, int j, int i1, int j1, unsigned char *MASK, unsigned char
                         Y_new[n] = y_n;
                         */
                         /*printf("[%i][%i][%u]\n", x_n, y_n, MASK[y_n*dimX+x_n]);*/
-                        // MASK_upd[x_n*dimX+y_n] = 10;                        
+                        // MASK_upd[x_n*dimX+y_n] = 10;
                         if (MASK[j*dimX+i] != MASK[x_n*dimX+y_n]) MASK_upd[x_n*dimX+y_n] = MASK[j*dimX+i];
                        }
                        else {

From 56a37d28b01078e43e742c47bba627e1a1a3ce86 Mon Sep 17 00:00:00 2001
From: Daniil Kazantsev <dkazanc@hotmail.com>
Date: Sun, 14 Apr 2019 18:40:07 +0100
Subject: [PATCH 5/7] progress2

---
 .../regularisers_CPU/DiffusionMASK_core.c     | 25 +++++++++++--------
 1 file changed, 14 insertions(+), 11 deletions(-)

diff --git a/src/Core/regularisers_CPU/DiffusionMASK_core.c b/src/Core/regularisers_CPU/DiffusionMASK_core.c
index dd4b1b48..a2110152 100644
--- a/src/Core/regularisers_CPU/DiffusionMASK_core.c
+++ b/src/Core/regularisers_CPU/DiffusionMASK_core.c
@@ -62,7 +62,7 @@ void swapVAL(unsigned char *xp, unsigned char *yp)
 
 float DiffusionMASK_CPU_main(float *Input, unsigned char *MASK, unsigned char *MASK_upd, unsigned char *SelClassesList, int SelClassesList_length, float *Output, float *infovector, int classesNumb, int DiffusWindow, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int penaltytype, float epsil, int dimX, int dimY, int dimZ)
 {
-    long i,j,k;
+    long i,j,k,l;
     int counterG, switcher;
     float sigmaPar2, *Output_prev=NULL, *Eucl_Vec;
     int DiffusWindow_tot;
@@ -96,7 +96,7 @@ float DiffusionMASK_CPU_main(float *Input, unsigned char *MASK, unsigned char *M
         }
         if (counterG == classesNumb) break;
       }
-      /* sort the obtained values (classes) */
+      /* sort from LOW->HIGH the obtained values (classes) */
       for(i=0; i<classesNumb; i++)	{
                   for(j=0; j<classesNumb-1; j++) {
                       if(ClassesList[j] > ClassesList[j+1]) {
@@ -105,8 +105,6 @@ float DiffusionMASK_CPU_main(float *Input, unsigned char *MASK, unsigned char *M
                           ClassesList[j] = temp;
                       }}}
 
-    for(i=0; i<classesNumb; i++)	printf("[%u]\n", ClassesList[i]);
-
     /*Euclidian weight for diffisuvuty window*/
     if (dimZ == 1) {
 	DiffusWindow_tot = (2*DiffusWindow + 1)*(2*DiffusWindow + 1);
@@ -151,18 +149,23 @@ float DiffusionMASK_CPU_main(float *Input, unsigned char *MASK, unsigned char *M
     /* copy the updated MASK (clean of outliers) */
     copyIm_unchar(MASK_upd, MASK_temp, (long)(dimX), (long)(dimY), (long)(dimZ));
 
-    #pragma omp parallel for shared(MASK_temp,MASK_upd) private(i,j)
+    for(l=0; l<SelClassesList_length; l++) {
+    /*printf("[%u]\n", ClassesList[SelClassesList[l]]);*/
+    #pragma omp parallel for shared(MASK_temp,MASK_upd,l) private(i,j)
     for(i=0; i<dimX; i++) {
         for(j=0; j<dimY; j++) {
+      /* The class of the central pixel has not changed, i.e. the central pixel is not an outlier -> continue */
       if (MASK_temp[j*dimX+i] == MASK[j*dimX+i]) {
-	/* !One needs to work with a specific class to avoid overlaps!
-	 hence it is crucial to establish relevant classes */
-       if (MASK_temp[j*dimX+i] == 149) {
-        /* The class of the central pixel has not changed, i.e. the central pixel is not an outlier -> continue */
+	    /* !One needs to work with a specific class to avoid overlaps! It is
+        crucial to establish relevant classes first (given as an input in SelClassesList) */
+       if (MASK_temp[j*dimX+i] == ClassesList[SelClassesList[l]]) {
         /* i = 258; j = 165; */
         Mask_update2D(MASK_temp, MASK_upd, i, j, DiffusWindow, (long)(dimX), (long)(dimY));
-       }}
-     }}
+        }}
+      }}
+      /* copy the updated mask */
+      copyIm_unchar(MASK_upd, MASK_temp, (long)(dimX), (long)(dimY), (long)(dimZ));
+      }
     }
 
     /* The mask has been processed, start diffusivity iterations */

From 7c79ab71d9d9613e03e9c822b9c8dd4de98d868d Mon Sep 17 00:00:00 2001
From: dkazanc <dkazanc@hotmail.com>
Date: Mon, 15 Apr 2019 14:46:28 +0100
Subject: [PATCH 6/7] created separate module for mask processing

---
 src/Core/CMakeLists.txt                     |   8 +-
 src/Core/regularisers_CPU/MASK_merge_core.c | 262 ++++++++++++++++++++
 src/Core/regularisers_CPU/MASK_merge_core.h |  56 +++++
 src/Python/ccpi/filters/regularisers.py     |  45 +---
 src/Python/src/cpu_regularisers.pyx         |  76 +++---
 5 files changed, 368 insertions(+), 79 deletions(-)
 create mode 100644 src/Core/regularisers_CPU/MASK_merge_core.c
 create mode 100644 src/Core/regularisers_CPU/MASK_merge_core.h

diff --git a/src/Core/CMakeLists.txt b/src/Core/CMakeLists.txt
index 6975a89f..8aec7495 100644
--- a/src/Core/CMakeLists.txt
+++ b/src/Core/CMakeLists.txt
@@ -68,12 +68,12 @@ add_library(cilreg SHARED
 	    ${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/SB_TV_core.c
 	    ${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/TGV_core.c
 	    ${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/Diffusion_core.c
-    	    ${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/DiffusionMASK_core.c
+    	    ${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/MASK_merge_core.c
 	    ${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/Diffus4th_order_core.c
 	    ${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/LLT_ROF_core.c
-        ${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/ROF_TV_core.c
-        ${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/ROF_TV_core.c
+            ${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/Nonlocal_TV_core.c
             ${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/PatchSelect_core.c
 	    ${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/utils.c
diff --git a/src/Core/regularisers_CPU/MASK_merge_core.c b/src/Core/regularisers_CPU/MASK_merge_core.c
new file mode 100644
index 00000000..d77c8127
--- /dev/null
+++ b/src/Core/regularisers_CPU/MASK_merge_core.c
@@ -0,0 +1,262 @@
+/*
+ * This work is part of the Core Imaging Library developed by
+ * Visual Analytics and Imaging System Group of the Science Technology
+ * Facilities Council, STFC
+ *
+ * Copyright 2019 Daniil Kazantsev
+ * Copyright 2019 Srikanth Nagella, Edoardo Pasca
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ * http://www.apache.org/licenses/LICENSE-2.0
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#include "MASK_merge_core.h"
+#include "utils.h"
+
+/* A method to ensure connectivity within regions of the segmented image/volume. Here we assume 
+ * that the MASK has been obtained using some classification/segmentation method such as k-means or gaussian
+ * mixture. Some pixels/voxels have been misclassified and we check the spatial dependences
+ * and correct the mask. We check the connectivity using the bresenham line algorithm within the non-local window
+ * surrounding the pixel of interest. 
+ *
+ * Input Parameters:
+ * 1. MASK [0:255], the result of some classification algorithm (information-based preferably)
+ * 2. The list of classes (e.g. [3,4]) to apply the method. The given order matters. 
+ * 3. The total number of classes in the MASK. 
+ * 4. The size of the Correction Window inside which the method works. 
+
+ * Output:
+ * 1. MASK_upd - the UPDATED MASK where some regions have been corrected (merged) or removed
+ * 2. CORRECTEDRegions - The array of the same size as MASK where all regions which were 
+ * changed are highlighted and the changes have been counted
+ */
+
+float Mask_merge_main(unsigned char *MASK, unsigned char *MASK_upd, unsigned char *CORRECTEDRegions, unsigned char *SelClassesList, int SelClassesList_length, int classesNumb, int CorrectionWindow, int dimX, int dimY, int dimZ)
+{
+    long i,j,l;
+    int counterG, switcher;
+    long DimTotal;
+    unsigned char *MASK_temp, *ClassesList, CurrClass, temp;
+    DimTotal = (long)(dimX*dimY*dimZ);
+
+    /* defines the list for all classes in the mask */
+    ClassesList = (unsigned char*) calloc (classesNumb,sizeof(unsigned char));
+
+     /* find which classes (values) are present in the segmented data */
+     CurrClass =  MASK[0]; ClassesList[0]= MASK[0]; counterG = 1;
+     for(i=0; i<DimTotal; i++) {
+       if (MASK[i] != CurrClass) {
+          switcher = 1;
+          for(j=0; j<counterG; j++) {
+            if (ClassesList[j] == MASK[i]) {
+              switcher = 0;
+              break;
+            }}
+            if (switcher == 1) {
+                CurrClass = MASK[i];
+                ClassesList[counterG] = MASK[i];
+                /*printf("[%u]\n", ClassesList[counterG]);*/
+                counterG++;
+              }
+        }
+        if (counterG == classesNumb) break;
+      }
+      /* sort from LOW->HIGH the obtained values (classes) */
+      for(i=0; i<classesNumb; i++)	{
+                  for(j=0; j<classesNumb-1; j++) {
+                      if(ClassesList[j] > ClassesList[j+1]) {
+                          temp = ClassesList[j+1];
+                          ClassesList[j+1] = ClassesList[j];
+                          ClassesList[j] = temp;
+                      }}}
+
+    MASK_temp = (unsigned char*) calloc (DimTotal,sizeof(unsigned char));
+
+    /* copy given MASK to MASK_upd*/
+    copyIm_unchar(MASK, MASK_upd, (long)(dimX), (long)(dimY), (long)(dimZ));
+
+    if (dimZ == 1) {
+    /********************** PERFORM 2D MASK PROCESSING ************************/
+    #pragma omp parallel for shared(MASK,MASK_upd) private(i,j)
+    for(i=0; i<dimX; i++) {
+        for(j=0; j<dimY; j++) {
+    /* STEP1: in a smaller neighbourhood check that the current pixel is NOT an outlier */
+    OutiersRemoval2D(MASK, MASK_upd, i, j, (long)(dimX), (long)(dimY));
+    }}
+    /* copy the updated MASK (clean of outliers) */
+    copyIm_unchar(MASK_upd, MASK_temp, (long)(dimX), (long)(dimY), (long)(dimZ));
+
+    for(l=0; l<SelClassesList_length; l++) {
+    /*printf("[%u]\n", ClassesList[SelClassesList[l]]);*/
+    #pragma omp parallel for shared(MASK_temp,MASK_upd,l) private(i,j)
+    for(i=0; i<dimX; i++) {
+        for(j=0; j<dimY; j++) {
+      /* The class of the central pixel has not changed, i.e. the central pixel is not an outlier -> continue */
+      if (MASK_temp[j*dimX+i] == MASK[j*dimX+i]) {
+	    /* !One needs to work with a specific class to avoid overlaps! It is
+        crucial to establish relevant classes first (given as an input in SelClassesList) */
+       if (MASK_temp[j*dimX+i] == ClassesList[SelClassesList[l]]) {
+        /* i = 258; j = 165; */
+        Mask_update2D(MASK_temp, MASK_upd, CORRECTEDRegions, i, j, CorrectionWindow, (long)(dimX), (long)(dimY));
+        }}
+      }}
+      /* copy the updated mask */
+      copyIm_unchar(MASK_upd, MASK_temp, (long)(dimX), (long)(dimY), (long)(dimZ));
+      }
+    }
+    else {
+    /********************** PERFORM 3D MASK PROCESSING ************************/
+    
+    
+    }
+
+    free(MASK_temp);   
+    return 0;
+}
+
+
+/********************************************************************/
+/***************************2D Functions*****************************/
+/********************************************************************/
+float OutiersRemoval2D(unsigned char *MASK, unsigned char *MASK_upd, long i, long j, long dimX, long dimY)
+{
+  /*if the ROI pixel does not belong to the surrondings, turn it into the surronding*/
+  long i_m, j_m, i1, j1, counter;
+    counter = 0;
+    for(i_m=-1; i_m<=1; i_m++) {
+      for(j_m=-1; j_m<=1; j_m++) {
+        i1 = i+i_m;
+        j1 = j+j_m;
+        if (((i1 >= 0) && (i1 < dimX)) && ((j1 >= 0) && (j1 < dimY))) {
+          if (MASK[j*dimX+i] != MASK[j1*dimX+i1]) counter++;
+        }
+      }}
+      if (counter >= 8) MASK_upd[j*dimX+i] = MASK[j1*dimX+i1];
+      return *MASK_upd;
+}
+
+float Mask_update2D(unsigned char *MASK_temp, unsigned char *MASK_upd, unsigned char *CORRECTEDRegions, long i, long j, int CorrectionWindow, long dimX, long dimY)
+{
+  long i_m, j_m, i1, j1, CounterOtherClass;
+
+  /* STEP2: in a larger neighbourhood check that the other class is present  */
+  CounterOtherClass = 0;
+  for(i_m=-CorrectionWindow; i_m<=CorrectionWindow; i_m++) {
+      for(j_m=-CorrectionWindow; j_m<=CorrectionWindow; j_m++) {
+        i1 = i+i_m;
+        j1 = j+j_m;
+        if (((i1 >= 0) && (i1 < dimX)) && ((j1 >= 0) && (j1 < dimY))) {
+          if (MASK_temp[j*dimX+i] != MASK_temp[j1*dimX+i1]) CounterOtherClass++;
+        }
+      }}
+      if (CounterOtherClass > 0) {
+      /* the other class is present in the vicinity of CorrectionWindow, continue to STEP 3 */
+      /*
+      STEP 3: Loop through all neighbours in CorrectionWindow and check the spatial connection.
+      Meaning that we're instrested if there are any classes between points A and B that
+      does not belong to A and B (A,B \in C)
+      */
+      for(i_m=-CorrectionWindow; i_m<=CorrectionWindow; i_m++) {
+          for(j_m=-CorrectionWindow; j_m<=CorrectionWindow; j_m++) {
+            i1 = i+i_m;
+            j1 = j+j_m;
+            if (((i1 >= 0) && (i1 < dimX)) && ((j1 >= 0) && (j1 < dimY))) {
+              if (MASK_temp[j*dimX+i] == MASK_temp[j1*dimX+i1]) {
+                /* A and B points belong to the same class, do STEP 4*/
+                /* STEP 4: Run Bresenham line algorithm between A and B points
+                and convert all points on the way to the class of A/B.  */
+               bresenham2D(i, j, i1, j1, MASK_temp, MASK_upd, CORRECTEDRegions, (long)(dimX), (long)(dimY));
+              }
+            }
+          }}
+      }
+  return 0;
+}
+int bresenham2D(int i, int j, int i1, int j1, unsigned char *MASK, unsigned char *MASK_upd, unsigned char *CORRECTEDRegions, long dimX, long dimY)
+{
+                   int n;
+                   int x[] = {i, i1};
+                   int y[] = {j, j1};
+                   int steep = (fabs(y[1]-y[0]) > fabs(x[1]-x[0]));
+                   int ystep = 0;
+
+                   //printf("[%i][%i][%i][%i]\n", x[1], y[1], steep, kk) ;
+                   //if (steep == 1) {swap(x[0],y[0]); swap(x[1],y[1]);}
+
+                   if (steep == 1) {
+                   // swaping
+                   int a, b;
+
+                   a = x[0];
+                   b = y[0];
+                   x[0] = b;
+                   y[0] = a;
+
+                   a = x[1];
+                   b = y[1];
+                   x[1] = b;
+                   y[1] = a;
+                   }
+
+                   if (x[0] > x[1]) {
+                   int a, b;
+                   a = x[0];
+                   b = x[1];
+                   x[0] = b;
+                   x[1] = a;
+
+                   a = y[0];
+                   b = y[1];
+                   y[0] = b;
+                   y[1] = a;
+                   } //(x[0] > x[1])
+
+                  int delx = x[1]-x[0];
+                  int dely = fabs(y[1]-y[0]);
+                  int error = 0;
+                  int x_n = x[0];
+                  int y_n = y[0];
+                  if (y[0] < y[1]) {ystep = 1;}
+                  else {ystep = -1;}
+
+                  for(n = 0; n<delx+1; n++) {
+                       if (steep == 1) {
+                        /*printf("[%i][%i][%u]\n", x_n, y_n, MASK[y_n*dimX+x_n]);*/
+                        // MASK_upd[x_n*dimX+y_n] = 10;
+                        if (MASK[j*dimX+i] != MASK[x_n*dimX+y_n]) {
+                        	MASK_upd[x_n*dimX+y_n] = MASK[j*dimX+i];
+                        	CORRECTEDRegions[x_n*dimX+y_n] += 1;
+                        }
+                       }
+                       else {
+                        // printf("[%i][%i][%u]\n", y_n, x_n, MASK[x_n*dimX+y_n]);
+                        // MASK_upd[y_n*dimX+x_n] = 20;
+                        if (MASK[j*dimX+i] != MASK[y_n*dimX+x_n]) {
+	                        MASK_upd[y_n*dimX+x_n] = MASK[j*dimX+i];
+                              	CORRECTEDRegions[y_n*dimX+x_n] += 1;
+                        }
+                       }
+                       x_n = x_n + 1;
+                       error = error + dely;
+
+                       if (2*error >= delx) {
+                          y_n = y_n + ystep;
+                         error = error - delx;
+                       } // (2*error >= delx)
+                       //printf("[%i][%i][%i]\n", X_new[n], Y_new[n], n) ;
+                  } // for(int n = 0; n<delx+1; n++)
+                  return 0;
+}
+/********************************************************************/
+/***************************3D Functions*****************************/
+/********************************************************************/
+
+
+
diff --git a/src/Core/regularisers_CPU/MASK_merge_core.h b/src/Core/regularisers_CPU/MASK_merge_core.h
new file mode 100644
index 00000000..287f8751
--- /dev/null
+++ b/src/Core/regularisers_CPU/MASK_merge_core.h
@@ -0,0 +1,56 @@
+/*
+This work is part of the Core Imaging Library developed by
+Visual Analytics and Imaging System Group of the Science Technology
+Facilities Council, STFC
+
+Copyright 2019 Daniil Kazantsev
+Copyright 2019 Srikanth Nagella, Edoardo Pasca
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+http://www.apache.org/licenses/LICENSE-2.0
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+*/
+
+#include <math.h>
+#include <stdlib.h>
+#include <memory.h>
+#include <stdio.h>
+#include "omp.h"
+#include "utils.h"
+#include "CCPiDefines.h"
+
+/* A method to ensure connectivity within regions of the segmented image/volume. Here we assume 
+ * that the MASK has been obtained using some classification/segmentation method such as k-means or gaussian
+ * mixture. Some pixels/voxels have been misclassified and we check the spatial dependences
+ * and correct the mask. We check the connectivity using the bresenham line algorithm within the non-local window
+ * surrounding the pixel of interest. 
+ *
+ * Input Parameters:
+ * 1. MASK [0:255], the result of some classification algorithm (information-based preferably)
+ * 2. The list of classes (e.g. [3,4]) to apply the method. The given order matters. 
+ * 3. The total number of classes in the MASK. 
+ * 4. The size of the Correction Window inside which the method works. 
+
+ * Output:
+ * 1. MASK_upd - the UPDATED MASK where some regions have been corrected (merged) or removed
+ * 2. CORRECTEDRegions - The array of the same size as MASK where all regions which were 
+ * changed are highlighted and the changes have been counted
+ */
+
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+CCPI_EXPORT float Mask_merge_main(unsigned char *MASK, unsigned char *MASK_upd, unsigned char *CORRECTEDRegions, unsigned char *SelClassesList, int SelClassesList_length, int classesNumb, int CorrectionWindow, int dimX, int dimY, int dimZ);
+CCPI_EXPORT float OutiersRemoval2D(unsigned char *MASK, unsigned char *MASK_upd, long i, long j, long dimX, long dimY);
+CCPI_EXPORT float Mask_update2D(unsigned char *MASK_temp, unsigned char *MASK_upd, unsigned char *CORRECTEDRegions, long i, long j, int CorrectionWindow, long dimX, long dimY);
+CCPI_EXPORT int bresenham2D(int i, int j, int i1, int j1, unsigned char *MASK, unsigned char *MASK_upd, unsigned char *CORRECTEDRegions, long dimX, long dimY);
+#ifdef __cplusplus
+}
+#endif
diff --git a/src/Python/ccpi/filters/regularisers.py b/src/Python/ccpi/filters/regularisers.py
index 610907d3..2fee8b3d 100644
--- a/src/Python/ccpi/filters/regularisers.py
+++ b/src/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_MASK_CPU, Diff4th_CPU, TGV_CPU, LLT_ROF_CPU, PATCHSEL_CPU, NLTV_CPU
+from ccpi.filters.cpu_regularisers import TV_ROF_CPU, TV_FGP_CPU, TV_SB_CPU, dTV_FGP_CPU, TNV_CPU, NDF_CPU, Diff4th_CPU, TGV_CPU, LLT_ROF_CPU, PATCHSEL_CPU, NLTV_CPU, MASK_CORR_CPU
 try:
     from ccpi.filters.gpu_regularisers import TV_ROF_GPU, TV_FGP_GPU, TV_SB_GPU, dTV_FGP_GPU, NDF_GPU, Diff4th_GPU, TGV_GPU, LLT_ROF_GPU, PATCHSEL_GPU
     gpu_enabled = True
@@ -127,37 +127,6 @@ def NDF(inputData, regularisation_parameter, edge_parameter, iterations,
     	    raise ValueError ('GPU is not available')
         raise ValueError('Unknown device {0}. Expecting gpu or cpu'\
                          .format(device))
-def NDF_MASK(inputData, maskdata, select_classes, total_classesNum, diffuswindow, regularisation_parameter, edge_parameter, iterations,
-                     time_marching_parameter, penalty_type, tolerance_param, device='cpu'):
-    if device == 'cpu':
-        return NDF_MASK_CPU(inputData,
-                     maskdata,
-                     select_classes,
-                     total_classesNum,
-                     diffuswindow,
-                     regularisation_parameter,
-                     edge_parameter,
-                     iterations,
-                     time_marching_parameter,
-                     penalty_type,
-                     tolerance_param)
-    elif device == 'gpu' and gpu_enabled:
-        return NDF_MASK_CPU(inputData,
-                     maskdata,
-                     select_classes,
-                     total_classesNum,
-                     diffuswindow,
-                     regularisation_parameter,
-                     edge_parameter,
-                     iterations,
-                     time_marching_parameter,
-                     penalty_type,
-                     tolerance_param)
-    else:
-        if not gpu_enabled and device == 'gpu':
-    	    raise ValueError ('GPU is not available')
-        raise ValueError('Unknown device {0}. Expecting gpu or cpu'\
-                         .format(device))
 def Diff4th(inputData, regularisation_parameter, edge_parameter, iterations,
                      time_marching_parameter, tolerance_param, device='cpu'):
     if device == 'cpu':
@@ -243,3 +212,15 @@ def NDF_INP(inputData, maskData, regularisation_parameter, edge_parameter, itera
 
 def NVM_INP(inputData, maskData, SW_increment, iterations):
         return NVM_INPAINT_CPU(inputData, maskData, SW_increment, iterations)
+    
+def MASK_CORR(maskdata, select_classes, total_classesNum, CorrectionWindow, device='cpu'):
+    if device == 'cpu':
+        return MASK_CORR_CPU(maskdata, select_classes, total_classesNum, CorrectionWindow)
+    elif device == 'gpu' and gpu_enabled:
+        return MASK_CORR_CPU(maskdata, select_classes, total_classesNum, CorrectionWindow)
+    else:
+        if not gpu_enabled and device == 'gpu':
+    	    raise ValueError ('GPU is not available')
+        raise ValueError('Unknown device {0}. Expecting gpu or cpu'\
+                         .format(device))
+
diff --git a/src/Python/src/cpu_regularisers.pyx b/src/Python/src/cpu_regularisers.pyx
index 78c46e2f..a63ecfaf 100644
--- a/src/Python/src/cpu_regularisers.pyx
+++ b/src/Python/src/cpu_regularisers.pyx
@@ -24,7 +24,7 @@ cdef extern float SB_TV_CPU_main(float *Input, float *Output, float *infovector,
 cdef extern float LLT_ROF_CPU_main(float *Input, float *Output, float *infovector, float lambdaROF, float lambdaLLT, int iterationsNumb, float tau, float epsil, int dimX, int dimY, int dimZ);
 cdef extern float TGV_main(float *Input, float *Output, float *infovector, float lambdaPar, float alpha1, float alpha0, int iterationsNumb, float L2, float epsil, int dimX, int dimY, int dimZ);
 cdef extern float Diffusion_CPU_main(float *Input, float *Output, float *infovector, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int penaltytype, float epsil, int dimX, int dimY, int dimZ);
-cdef extern float DiffusionMASK_CPU_main(float *Input, unsigned char *MASK, unsigned char *MASK_upd, unsigned char *SelClassesList, int SelClassesList_length, float *Output, float *infovector, int classesNumb, int DiffusWindow, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int penaltytype, float epsil, int dimX, int dimY, int dimZ);
+cdef extern float Mask_merge_main(unsigned char *MASK, unsigned char *MASK_upd, unsigned char *CORRECTEDRegions, unsigned char *SelClassesList, int SelClassesList_length, int classesNumb, int CorrectionWindow, int dimX, int dimY, int dimZ);
 cdef extern float Diffus4th_CPU_main(float *Input, float *Output,  float *infovector, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, float epsil, int dimX, int dimY, int dimZ);
 cdef extern float dTV_FGP_CPU_main(float *Input, float *InputRef, float *Output, float *infovector, float lambdaPar, int iterationsNumb, float epsil, float eta, int methodTV, int nonneg, int dimX, int dimY, int dimZ);
 cdef extern float TNV_CPU_main(float *Input, float *u, float lambdaPar, int maxIter, float tol, int dimX, int dimY, int dimZ);
@@ -381,48 +381,6 @@ def NDF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
     dims[2], dims[1], dims[0])
     return (outputData,infovec)
 
-#****************************************************************#
-#********Constrained Nonlinear(Isotropic) Diffusion**************#
-#****************************************************************#
-def NDF_MASK_CPU(inputData, maskData, select_classes, total_classesNum, diffuswindow, regularisation_parameter, edge_parameter, iterationsNumb,time_marching_parameter, penalty_type, tolerance_param):
-    if inputData.ndim == 2:
-        return NDF_MASK_2D(inputData, maskData, select_classes, total_classesNum, diffuswindow, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type, tolerance_param)
-    elif inputData.ndim == 3:
-        return 0
-
-def NDF_MASK_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=1, mode="c"] select_classes,
-                     int total_classesNum,
-                     int diffuswindow,
-                     float regularisation_parameter,
-                     float edge_parameter,
-                     int iterationsNumb,
-                     float time_marching_parameter,
-                     int penalty_type,
-                     float tolerance_param):
-    cdef long dims[2]
-    dims[0] = inputData.shape[0]
-    dims[1] = inputData.shape[1]
-
-    select_classes_length = select_classes.shape[0]
-    
-    cdef np.ndarray[np.uint8_t, ndim=2, mode="c"] mask_upd = \
-            np.zeros([dims[0],dims[1]], dtype='uint8')
-    cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
-            np.zeros([dims[0],dims[1]], dtype='float32')
-    cdef np.ndarray[np.float32_t, ndim=1, mode="c"] infovec = \
-                np.zeros([2], dtype='float32')
-
-
-    # Run constrained nonlinear diffusion iterations for 2D data
-    DiffusionMASK_CPU_main(&inputData[0,0], &maskData[0,0], &mask_upd[0,0], &select_classes[0], select_classes_length, &outputData[0,0], &infovec[0],
-    total_classesNum, diffuswindow, regularisation_parameter, edge_parameter, iterationsNumb,
-    time_marching_parameter, penalty_type,
-    tolerance_param,
-    dims[1], dims[0], 1)
-    return (mask_upd,outputData,infovec)
-
 #****************************************************************#
 #*************Anisotropic Fourth-Order diffusion*****************#
 #****************************************************************#
@@ -736,6 +694,38 @@ def NVM_INP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
 
     return (outputData, maskData_upd)
 
+##############################################################################
+#****************************************************************#
+#********Mask (segmented image) correction module **************#
+#****************************************************************#
+def MASK_CORR_CPU(maskData, select_classes, total_classesNum, CorrectionWindow):
+    if maskData.ndim == 2:
+        return MASK_CORR_CPU_2D(maskData, select_classes, total_classesNum, CorrectionWindow)
+    elif maskData.ndim == 3:
+        return 0
+
+def MASK_CORR_CPU_2D(np.ndarray[np.uint8_t, ndim=2, mode="c"] maskData,
+                    np.ndarray[np.uint8_t, ndim=1, mode="c"] select_classes,
+                     int total_classesNum,
+                     int CorrectionWindow):
+    
+    cdef long dims[2]
+    dims[0] = maskData.shape[0]
+    dims[1] = maskData.shape[1]
+
+    select_classes_length = select_classes.shape[0]
+    
+    cdef np.ndarray[np.uint8_t, ndim=2, mode="c"] mask_upd = \
+            np.zeros([dims[0],dims[1]], dtype='uint8')
+    cdef np.ndarray[np.uint8_t, ndim=2, mode="c"] corr_regions = \
+            np.zeros([dims[0],dims[1]], dtype='uint8')
+
+    # Run the function to process given MASK
+    Mask_merge_main(&maskData[0,0], &mask_upd[0,0], &corr_regions[0,0], &select_classes[0], select_classes_length, 
+    total_classesNum, CorrectionWindow, dims[1], dims[0], 1)
+    return (mask_upd,corr_regions)
+
+##############################################################################
 
 #****************************************************************#
 #***************Calculation of TV-energy functional**************#

From e6842ec7f2cdbd46a004758bc3a6543012c6a74a Mon Sep 17 00:00:00 2001
From: dkazanc <dkazanc@hotmail.com>
Date: Mon, 15 Apr 2019 14:47:36 +0100
Subject: [PATCH 7/7] old deleted

---
 .../regularisers_CPU/DiffusionMASK_core.c     | 427 ------------------
 .../regularisers_CPU/DiffusionMASK_core.h     |  65 ---
 2 files changed, 492 deletions(-)
 delete mode 100644 src/Core/regularisers_CPU/DiffusionMASK_core.c
 delete mode 100644 src/Core/regularisers_CPU/DiffusionMASK_core.h

diff --git a/src/Core/regularisers_CPU/DiffusionMASK_core.c b/src/Core/regularisers_CPU/DiffusionMASK_core.c
deleted file mode 100644
index a2110152..00000000
--- a/src/Core/regularisers_CPU/DiffusionMASK_core.c
+++ /dev/null
@@ -1,427 +0,0 @@
-/*
- * 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 "DiffusionMASK_core.h"
-#include "utils.h"
-
-#define EPS 1.0e-5
-#define MAX(x, y) (((x) > (y)) ? (x) : (y))
-#define MIN(x, y) (((x) < (y)) ? (x) : (y))
-
-/*sign function*/
-int signNDF_m(float x) {
-    return (x > 0) - (x < 0);
-}
-
-/* C-OMP implementation of linear and nonlinear diffusion [1,2] which is constrained by the provided MASK.
- * The minimisation is performed using explicit scheme.
- * Implementation using the diffusivity window to increase the coverage area of the diffusivity
- *
- * Input Parameters:
- * 1. Noisy image/volume
- * 2. MASK (in unsigned char format)
- * 3. Diffusivity window (half-size of the searching window, e.g. 3)
- * 4. lambda - regularization parameter
- * 5. Edge-preserving parameter (sigma), when sigma equals to zero nonlinear diffusion -> linear diffusion
- * 6. Number of iterations, for explicit scheme >= 150 is recommended
- * 7. tau - time-marching step for explicit scheme
- * 8. Penalty type: 1 - Huber, 2 - Perona-Malik, 3 - Tukey Biweight
- * 9. eplsilon - tolerance constant
-
- * Output:
- * [1] Filtered/regularized image/volume
- * [2] Information vector which contains [iteration no., reached tolerance]
- *
- * 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 swapVAL(unsigned char *xp, unsigned char *yp)
-{
-    unsigned char temp = *xp;
-    *xp = *yp;
-    *yp = temp;
-}
-
-float DiffusionMASK_CPU_main(float *Input, unsigned char *MASK, unsigned char *MASK_upd, unsigned char *SelClassesList, int SelClassesList_length, float *Output, float *infovector, int classesNumb, int DiffusWindow, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int penaltytype, float epsil, int dimX, int dimY, int dimZ)
-{
-    long i,j,k,l;
-    int counterG, switcher;
-    float sigmaPar2, *Output_prev=NULL, *Eucl_Vec;
-    int DiffusWindow_tot;
-    sigmaPar2 = sigmaPar/sqrt(2.0f);
-    long DimTotal;
-    float re, re1;
-    re = 0.0f; re1 = 0.0f;
-    int count = 0;
-    unsigned char *MASK_temp, *ClassesList, CurrClass, temp;
-    DimTotal = (long)(dimX*dimY*dimZ);
-
-    /* defines the list for all classes in the mask */
-    ClassesList = (unsigned char*) calloc (classesNumb,sizeof(unsigned char));
-
-     /* find which classes (values) are present in the segmented data */
-     CurrClass =  MASK[0]; ClassesList[0]= MASK[0]; counterG = 1;
-     for(i=0; i<DimTotal; i++) {
-       if (MASK[i] != CurrClass) {
-          switcher = 1;
-          for(j=0; j<counterG; j++) {
-            if (ClassesList[j] == MASK[i]) {
-              switcher = 0;
-              break;
-            }}
-            if (switcher == 1) {
-                CurrClass = MASK[i];
-                ClassesList[counterG] = MASK[i];
-                /*printf("[%u]\n", ClassesList[counterG]);*/
-                counterG++;
-              }
-        }
-        if (counterG == classesNumb) break;
-      }
-      /* sort from LOW->HIGH the obtained values (classes) */
-      for(i=0; i<classesNumb; i++)	{
-                  for(j=0; j<classesNumb-1; j++) {
-                      if(ClassesList[j] > ClassesList[j+1]) {
-                          temp = ClassesList[j+1];
-                          ClassesList[j+1] = ClassesList[j];
-                          ClassesList[j] = temp;
-                      }}}
-
-    /*Euclidian weight for diffisuvuty window*/
-    if (dimZ == 1) {
-	DiffusWindow_tot = (2*DiffusWindow + 1)*(2*DiffusWindow + 1);
-        /* generate a 2D Gaussian kernel for NLM procedure */
-        Eucl_Vec = (float*) calloc (DiffusWindow_tot,sizeof(float));
-        counterG = 0;
-        for(i=-DiffusWindow; i<=DiffusWindow; i++) {
-            for(j=-DiffusWindow; j<=DiffusWindow; j++) {
-                Eucl_Vec[counterG] = (float)expf(-(powf(((float) i), 2) + powf(((float) j), 2))/(2.0f*DiffusWindow*DiffusWindow));
-                counterG++;
-            }} /*main neighb loop */
-       }
-    else {
-	DiffusWindow_tot = (2*DiffusWindow + 1)*(2*DiffusWindow + 1)*(2*DiffusWindow + 1);
-	Eucl_Vec = (float*) calloc (DiffusWindow_tot,sizeof(float));
-        counterG = 0;
-        for(i=-DiffusWindow; i<=DiffusWindow; i++) {
-            for(j=-DiffusWindow; j<=DiffusWindow; j++) {
-                for(k=-DiffusWindow; k<=DiffusWindow; k++) {
-                    Eucl_Vec[counterG] = (float)expf(-(powf(((float) i), 2) + powf(((float) j), 2) + powf(((float) k), 2))/(2*DiffusWindow*DiffusWindow*DiffusWindow));
-                    counterG++;
-                }}} /*main neighb loop */
-    }
-
-    if (epsil != 0.0f) Output_prev = calloc(DimTotal, sizeof(float));
-
-    MASK_temp = (unsigned char*) calloc (DimTotal,sizeof(unsigned char));
-
-    /* copy input into output */
-    copyIm(Input, Output, (long)(dimX), (long)(dimY), (long)(dimZ));
-    /* copy given MASK to MASK_upd*/
-    copyIm_unchar(MASK, MASK_upd, (long)(dimX), (long)(dimY), (long)(dimZ));
-
-    /********************** PERFORM MASK PROCESSING ************************/
-    if (dimZ == 1) {
-    #pragma omp parallel for shared(MASK,MASK_upd) private(i,j)
-    for(i=0; i<dimX; i++) {
-        for(j=0; j<dimY; j++) {
-    /* STEP1: in a smaller neighbourhood check that the current pixel is NOT an outlier */
-    OutiersRemoval2D(MASK, MASK_upd, i, j, (long)(dimX), (long)(dimY));
-    }}
-    /* copy the updated MASK (clean of outliers) */
-    copyIm_unchar(MASK_upd, MASK_temp, (long)(dimX), (long)(dimY), (long)(dimZ));
-
-    for(l=0; l<SelClassesList_length; l++) {
-    /*printf("[%u]\n", ClassesList[SelClassesList[l]]);*/
-    #pragma omp parallel for shared(MASK_temp,MASK_upd,l) private(i,j)
-    for(i=0; i<dimX; i++) {
-        for(j=0; j<dimY; j++) {
-      /* The class of the central pixel has not changed, i.e. the central pixel is not an outlier -> continue */
-      if (MASK_temp[j*dimX+i] == MASK[j*dimX+i]) {
-	    /* !One needs to work with a specific class to avoid overlaps! It is
-        crucial to establish relevant classes first (given as an input in SelClassesList) */
-       if (MASK_temp[j*dimX+i] == ClassesList[SelClassesList[l]]) {
-        /* i = 258; j = 165; */
-        Mask_update2D(MASK_temp, MASK_upd, i, j, DiffusWindow, (long)(dimX), (long)(dimY));
-        }}
-      }}
-      /* copy the updated mask */
-      copyIm_unchar(MASK_upd, MASK_temp, (long)(dimX), (long)(dimY), (long)(dimZ));
-      }
-    }
-
-    /* The mask has been processed, start diffusivity iterations */
-    for(i=0; i < iterationsNumb; i++) {
-      if ((epsil != 0.0f)  && (i % 5 == 0)) copyIm(Output, Output_prev, (long)(dimX), (long)(dimY), (long)(dimZ));
-      if (dimZ == 1) {
-             /* running 2D diffusion iterations */
-            if (sigmaPar == 0.0f) LinearDiff_MASK2D(Input, MASK_upd, Output, Eucl_Vec, DiffusWindow, lambdaPar, tau, (long)(dimX), (long)(dimY)); /* constrained linear diffusion */
-            else NonLinearDiff_MASK2D(Input, MASK_upd, Output, Eucl_Vec, DiffusWindow, lambdaPar, sigmaPar2, tau, penaltytype, (long)(dimX), (long)(dimY)); /* constrained nonlinear diffusion */
-          }
-      else {
-       	/* running 3D diffusion iterations */
-        //if (sigmaPar == 0.0f) LinearDiff3D(Input, Output, lambdaPar, tau, (long)(dimX), (long)(dimY), (long)(dimZ));
-//       else NonLinearDiff3D(Input, Output, lambdaPar, sigmaPar2, tau, penaltytype, (long)(dimX), (long)(dimY), (long)(dimZ));
-          }
-          /* check early stopping criteria if epsilon not equal zero */
-          if ((epsil != 0.0f)  && (i % 5 == 0)) {
-          re = 0.0f; re1 = 0.0f;
-            for(j=0; j<DimTotal; j++)
-            {
-                re += powf(Output[j] - Output_prev[j],2);
-                re1 += powf(Output[j],2);
-            }
-          re = sqrtf(re)/sqrtf(re1);
-          /* stop if the norm residual is less than the tolerance EPS */
-          if (re < epsil)  count++;
-          if (count > 3) break;
-          }
-		}
-
-    free(Output_prev);
-    free(Eucl_Vec);
-    free(MASK_temp);
-  /*adding info into info_vector */
-    infovector[0] = (float)(i);  /*iterations number (if stopped earlier based on tolerance)*/
-    infovector[1] = re;  /* reached tolerance */
-    return 0;
-}
-
-
-/********************************************************************/
-/***************************2D Functions*****************************/
-/********************************************************************/
-float OutiersRemoval2D(unsigned char *MASK, unsigned char *MASK_upd, long i, long j, long dimX, long dimY)
-{
-  /*if the ROI pixel does not belong to the surrondings, turn it into the surronding*/
-  long i_m, j_m, i1, j1, counter;
-    counter = 0;
-    for(i_m=-1; i_m<=1; i_m++) {
-      for(j_m=-1; j_m<=1; j_m++) {
-        i1 = i+i_m;
-        j1 = j+j_m;
-        if (((i1 >= 0) && (i1 < dimX)) && ((j1 >= 0) && (j1 < dimY))) {
-          if (MASK[j*dimX+i] != MASK[j1*dimX+i1]) counter++;
-        }
-      }}
-      if (counter >= 8) MASK_upd[j*dimX+i] = MASK[j1*dimX+i1];
-      return *MASK_upd;
-}
-
-float Mask_update2D(unsigned char *MASK_temp, unsigned char *MASK_upd, long i, long j, int DiffusWindow, long dimX, long dimY)
-{
-  long i_m, j_m, i1, j1, CounterOtherClass;
-
-  /* STEP2: in a larger neighbourhood check that the other class is present  */
-  CounterOtherClass = 0;
-  for(i_m=-DiffusWindow; i_m<=DiffusWindow; i_m++) {
-      for(j_m=-DiffusWindow; j_m<=DiffusWindow; j_m++) {
-        i1 = i+i_m;
-        j1 = j+j_m;
-        if (((i1 >= 0) && (i1 < dimX)) && ((j1 >= 0) && (j1 < dimY))) {
-          if (MASK_temp[j*dimX+i] != MASK_temp[j1*dimX+i1]) CounterOtherClass++;
-        }
-      }}
-      if (CounterOtherClass > 0) {
-      /* the other class is present in the vicinity of DiffusWindow, continue to STEP 3 */
-      /*
-      STEP 3: Loop through all neighbours in DiffusWindow and check the spatial connection.
-      Meaning that we're instrested if there are any classes between points A and B that
-      does not belong to A and B (A,B \in C)
-      */
-      for(i_m=-DiffusWindow; i_m<=DiffusWindow; i_m++) {
-          for(j_m=-DiffusWindow; j_m<=DiffusWindow; j_m++) {
-            i1 = i+i_m;
-            j1 = j+j_m;
-            if (((i1 >= 0) && (i1 < dimX)) && ((j1 >= 0) && (j1 < dimY))) {
-              if (MASK_temp[j*dimX+i] == MASK_temp[j1*dimX+i1]) {
-                /* A and B points belong to the same class, do STEP 4*/
-                /* STEP 4: Run Bresenham line algorithm between A and B points
-                and convert all points on the way to the class of A/B.  */
-               bresenham2D(i, j, i1, j1, MASK_temp, MASK_upd, (long)(dimX), (long)(dimY));
-              }
-            }
-          }}
-      }
-  return *MASK_upd;
-}
-
-/* MASKED-constrained 2D linear diffusion (PDE heat equation) */
-float LinearDiff_MASK2D(float *Input, unsigned char *MASK, float *Output,  float *Eucl_Vec, int DiffusWindow, float lambdaPar, float tau, long dimX, long dimY)
-{
-
-long i,j,i1,j1,i_m,j_m,index,indexneighb,counter;
-unsigned char class_c, class_n;
-float diffVal;
-
-#pragma omp parallel for shared(Input) private(index,i,j,i1,j1,i_m,j_m,counter,diffVal,indexneighb,class_c,class_n)
-    for(i=0; i<dimX; i++) {
-        for(j=0; j<dimY; j++) {
-        index = j*dimX+i; /* current pixel index */
-        counter = 0; diffVal = 0.0f;
-        for(i_m=-DiffusWindow; i_m<=DiffusWindow; i_m++) {
-            for(j_m=-DiffusWindow; j_m<=DiffusWindow; j_m++) {
-            i1 = i+i_m;
-            j1 = j+j_m;
-            if (((i1 >= 0) && (i1 < dimX)) && ((j1 >= 0) && (j1 < dimY))) {
-            indexneighb = j1*dimX+i1; /* neighbour pixel index */
-	    class_c = MASK[index]; /* current class value */
-    	    class_n = MASK[indexneighb]; /* neighbour class value */
-
-	    /* perform diffusion only within the same class (given by MASK) */
-	    if (class_n == class_c) diffVal += Output[indexneighb] - Output[index];
-            }
-		counter++;
-	    }}
-            Output[index] += tau*(lambdaPar*(diffVal) - (Output[index] - Input[index]));
-        }}
-	return *Output;
-}
-
-/*  MASKED-constrained 2D nonlinear diffusion */
-float NonLinearDiff_MASK2D(float *Input, unsigned char *MASK, float *Output, float *Eucl_Vec, int DiffusWindow, float lambdaPar, float sigmaPar, float tau, int penaltytype, long dimX, long dimY)
-{
-	long i,j,i1,j1,i_m,j_m,index,indexneighb,counter;
-	unsigned char class_c, class_n;
-	float diffVal, funcVal;
-
-#pragma omp parallel for shared(Input) private(index,i,j,i1,j1,i_m,j_m,counter,diffVal,funcVal,indexneighb,class_c,class_n)
-    for(i=0; i<dimX; i++) {
-        for(j=0; j<dimY; j++) {
-        index = j*dimX+i; /* current pixel index */
-        counter = 0; diffVal = 0.0f; funcVal = 0.0f;
-        for(i_m=-DiffusWindow; i_m<=DiffusWindow; i_m++) {
-            for(j_m=-DiffusWindow; j_m<=DiffusWindow; j_m++) {
-            i1 = i+i_m;
-            j1 = j+j_m;
-            if (((i1 >= 0) && (i1 < dimX)) && ((j1 >= 0) && (j1 < dimY))) {
-            indexneighb = j1*dimX+i1; /* neighbour pixel index */
-	    class_c = MASK[index]; /* current class value */
-    	    class_n = MASK[indexneighb]; /* neighbour class value */
-
-	    /* perform diffusion only within the same class (given by MASK) */
-	    if (class_n == class_c) {
-	    	diffVal = Output[indexneighb] - Output[index];
-	    	if (penaltytype == 1) {
-	        /* Huber penalty */
-                if (fabs(diffVal) > sigmaPar) funcVal += signNDF_m(diffVal);
-                else funcVal += diffVal/sigmaPar; }
-  		else if (penaltytype == 2) {
-  		/* Perona-Malik */
-  		funcVal += (diffVal)/(1.0f + powf((diffVal/sigmaPar),2)); }
-  		else if (penaltytype == 3) {
-  		/* Tukey Biweight */
-  		if (fabs(diffVal) <= sigmaPar) funcVal += diffVal*powf((1.0f - powf((diffVal/sigmaPar),2)), 2); }
-                else {
-                printf("%s \n", "No penalty function selected! Use Huber,2 or 3.");
-		break; }
-            		}
-            	}
-		counter++;
-	    }}
-           Output[index] += tau*(lambdaPar*(funcVal) - (Output[index] - Input[index]));
-		}}
-	return *Output;
-}
-/********************************************************************/
-/***************************3D Functions*****************************/
-/********************************************************************/
-
-
-int bresenham2D(int i, int j, int i1, int j1, unsigned char *MASK, unsigned char *MASK_upd, long dimX, long dimY)
-{
-                   int n;
-                   int x[] = {i, i1};
-                   int y[] = {j, j1};
-                   int steep = (fabs(y[1]-y[0]) > fabs(x[1]-x[0]));
-                   int ystep = 0;
-
-                   //printf("[%i][%i][%i][%i]\n", x[1], y[1], steep, kk) ;
-                   //if (steep == 1) {swap(x[0],y[0]); swap(x[1],y[1]);}
-
-                   if (steep == 1) {
-
-                   // swaping
-                   int a, b;
-
-                   a = x[0];
-                   b = y[0];
-                   x[0] = b;
-                   y[0] = a;
-
-                   a = x[1];
-                   b = y[1];
-                   x[1] = b;
-                   y[1] = a;
-                   }
-
-                   if (x[0] > x[1]) {
-                   int a, b;
-                   a = x[0];
-                   b = x[1];
-                   x[0] = b;
-                   x[1] = a;
-
-                   a = y[0];
-                   b = y[1];
-                   y[0] = b;
-                   y[1] = a;
-                   } //(x[0] > x[1])
-
-                  int delx = x[1]-x[0];
-                  int dely = fabs(y[1]-y[0]);
-                  int error = 0;
-                  int x_n = x[0];
-                  int y_n = y[0];
-                  if (y[0] < y[1]) {ystep = 1;}
-                  else {ystep = -1;}
-
-                  for(n = 0; n<delx+1; n++) {
-                       if (steep == 1) {
-                        /*
-                        X_new[n] = x_n;
-                        Y_new[n] = y_n;
-                        */
-                        /*printf("[%i][%i][%u]\n", x_n, y_n, MASK[y_n*dimX+x_n]);*/
-                        // MASK_upd[x_n*dimX+y_n] = 10;
-                        if (MASK[j*dimX+i] != MASK[x_n*dimX+y_n]) MASK_upd[x_n*dimX+y_n] = MASK[j*dimX+i];
-                       }
-                       else {
-                         /*
-                        X_new[n] = y_n;
-                        Y_new[n] = x_n;
-                        */
-                        // printf("[%i][%i][%u]\n", y_n, x_n, MASK[x_n*dimX+y_n]);
-                        // MASK_upd[y_n*dimX+x_n] = 20;
-                        if (MASK[j*dimX+i] != MASK[y_n*dimX+x_n]) MASK_upd[y_n*dimX+x_n] = MASK[j*dimX+i];
-                       }
-                       x_n = x_n + 1;
-                       error = error + dely;
-
-                       if (2*error >= delx) {
-                          y_n = y_n + ystep;
-                         error = error - delx;
-                       } // (2*error >= delx)
-                       //printf("[%i][%i][%i]\n", X_new[n], Y_new[n], n) ;
-                  } // for(int n = 0; n<delx+1; n++)
-
-                  return 0;
-}
diff --git a/src/Core/regularisers_CPU/DiffusionMASK_core.h b/src/Core/regularisers_CPU/DiffusionMASK_core.h
deleted file mode 100644
index a032905d..00000000
--- a/src/Core/regularisers_CPU/DiffusionMASK_core.h
+++ /dev/null
@@ -1,65 +0,0 @@
-/*
-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] which is constrained by the provided MASK.
- * The minimisation is performed using explicit scheme.
- * Implementation using the Diffusivity window to increase the coverage area of the diffusivity
- *
- * Input Parameters:
- * 1. Noisy image/volume
- * 2. MASK (in unsigned short format)
- * 3. Diffusivity window (half-size of the searching window, e.g. 3)
- * 4. lambda - regularization parameter
- * 5. Edge-preserving parameter (sigma), when sigma equals to zero nonlinear diffusion -> linear diffusion
- * 6. Number of iterations, for explicit scheme >= 150 is recommended
- * 7. tau - time-marching step for explicit scheme
- * 8. Penalty type: 1 - Huber, 2 - Perona-Malik, 3 - Tukey Biweight
- * 9. eplsilon - tolerance constant
-
- * Output:
- * [1] Filtered/regularized image/volume
- * [2] Information vector which contains [iteration no., reached tolerance]
- *
- * 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 DiffusionMASK_CPU_main(float *Input, unsigned char *MASK, unsigned char *MASK_upd, unsigned char *SelClassesList, int SelClassesList_length, float *Output, float *infovector, int classesNumb, int DiffusWindow, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int penaltytype, float epsil, int dimX, int dimY, int dimZ);
-CCPI_EXPORT float LinearDiff_MASK2D(float *Input, unsigned char *MASK, float *Output,  float *Eucl_Vec, int DiffusWindow, float lambdaPar, float tau, long dimX, long dimY);
-CCPI_EXPORT float NonLinearDiff_MASK2D(float *Input, unsigned char *MASK, float *Output, float *Eucl_Vec, int DiffusWindow, float lambdaPar, float sigmaPar, float tau, int penaltytype, long dimX, long dimY);
-CCPI_EXPORT float OutiersRemoval2D(unsigned char *MASK, unsigned char *MASK_upd, long i, long j, long dimX, long dimY);
-CCPI_EXPORT float Mask_update2D(unsigned char *MASK_temp, unsigned char *MASK_upd, long i, long j, int DiffusWindow, long dimX, long dimY);
-CCPI_EXPORT int bresenham2D(int i, int j, int i1, int j1, unsigned char *MASK, unsigned char *MASK_upd, long dimX, long dimY);
-#ifdef __cplusplus
-}
-#endif