restructured sources

7 files changed, 1755 insertions, 0 deletions

diff --git a/src/Python/CMakeLists.txt b/src/Python/CMakeLists.txt
new file mode 100644
index 0000000..c2ef855
--- /dev/null
+++ b/src/Python/CMakeLists.txt
@@ -0,0 +1,141 @@
+#   Copyright 2018 Edoardo Pasca
+cmake_minimum_required (VERSION 3.0)
+
+project(regulariserPython)
+#https://stackoverflow.com/questions/13298504/using-cmake-with-setup-py
+
+# The version number.
+
+#set (CIL_VERSION $ENV{CIL_VERSION} CACHE INTERNAL "Core Imaging Library version" FORCE)
+
+# conda orchestrated build
+message("CIL_VERSION: ${CIL_VERSION}")
+#include (GenerateExportHeader)
+
+find_package(PythonInterp REQUIRED)
+if (PYTHONINTERP_FOUND)
+  message ("Current Python " ${PYTHON_VERSION_STRING} " found " ${PYTHON_EXECUTABLE})
+endif()
+
+	
+## Build the regularisers package as a library
+message("Creating Regularisers as shared library")
+
+message("CMAKE_CXX_FLAGS ${CMAKE_CXX_FLAGS}")
+
+set(CMAKE_BUILD_TYPE "Release")
+
+if(WIN32)
+  set (FLAGS "/DWIN32 /EHsc /openmp /DCCPiCore_EXPORTS")
+  set (CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} /NODEFAULTLIB:MSVCRT.lib")
+  
+  set (EXTRA_LIBRARIES)
+		
+  message("library lib: ${LIBRARY_LIB}")
+  
+elseif(UNIX)
+   set (FLAGS "-fopenmp -O2 -funsigned-char -Wall  -Wl,--no-undefined  -DCCPiReconstructionIterative_EXPORTS -std=c++0x")  
+   set (EXTRA_LIBRARIES 
+		"gomp"
+		)
+endif()
+
+# GPU regularisers
+if (BUILD_CUDA)
+    find_package(CUDA)
+    if (CUDA_FOUND)
+      message("CUDA FOUND")
+      set (SETUP_GPU_WRAPPERS "extra_libraries += ['cilregcuda']\n\
+setup( \n\
+        name='ccpi', \n\
+        description='CCPi Core Imaging Library - Image regularisers GPU',\n\
+        version=cil_version,\n\
+        cmdclass = {'build_ext': build_ext},\n\
+        ext_modules = [Extension('ccpi.filters.gpu_regularisers',\n\
+                                  sources=[ \n\
+                                          os.path.join('.' , 'src', 'gpu_regularisers.pyx' ),\n\
+                                            ],\n\
+                                 include_dirs=extra_include_dirs, \n\
+                                 library_dirs=extra_library_dirs, \n\
+                                 extra_compile_args=extra_compile_args, \n\
+                                 libraries=extra_libraries ), \n\
+        ],\n\
+        zip_safe = False,	\n\
+        packages = {'ccpi','ccpi.filters'},\n\
+    )")
+    else()
+      message("CUDA NOT FOUND")
+      set(SETUP_GPU_WRAPPERS "#CUDA NOT FOUND")
+    endif()
+endif()
+configure_file("${CMAKE_CURRENT_SOURCE_DIR}/setup-regularisers.py.in" "${CMAKE_CURRENT_BINARY_DIR}/setup-regularisers.py")
+
+
+find_package(PythonInterp)
+find_package(PythonLibs)
+if (PYTHONINTERP_FOUND)
+  message(STATUS "Found PYTHON_EXECUTABLE=${PYTHON_EXECUTABLE}")
+  message(STATUS "Python version ${PYTHON_VERSION_STRING}")
+endif()
+if (PYTHONLIBS_FOUND)
+  message(STATUS "Found PYTHON_INCLUDE_DIRS=${PYTHON_INCLUDE_DIRS}")
+  message(STATUS "Found PYTHON_LIBRARIES=${PYTHON_LIBRARIES}")
+endif()
+
+if (PYTHONINTERP_FOUND)
+    message("Python found " ${PYTHON_EXECUTABLE})
+    set(SETUP_PY_IN "${CMAKE_CURRENT_SOURCE_DIR}/setup-regularisers.py.in")
+    set(SETUP_PY    "${CMAKE_CURRENT_BINARY_DIR}/setup-regularisers.py")
+    #set(DEPS        "${CMAKE_CURRENT_SOURCE_DIR}/module/__init__.py")
+    set (DEPS       "${CMAKE_BINARY_DIR}/Core/")
+    set(OUTPUT      "${CMAKE_CURRENT_BINARY_DIR}/build/timestamp")
+
+    configure_file(${SETUP_PY_IN} ${SETUP_PY})
+
+    message("Core binary dir " ${CMAKE_BINARY_DIR}/Core/${CMAKE_BUILD_TYPE})
+    
+    if (CONDA_BUILD)
+      add_custom_command(OUTPUT ${OUTPUT}
+                       COMMAND ${CMAKE_COMMAND} -E copy_directory ${CMAKE_CURRENT_SOURCE_DIR}/src ${CMAKE_CURRENT_BINARY_DIR}/src
+                       COMMAND ${CMAKE_COMMAND} -E copy_directory ${CMAKE_CURRENT_SOURCE_DIR}/ccpi ${CMAKE_CURRENT_BINARY_DIR}/ccpi
+                       COMMAND ${CMAKE_COMMAND} -E env CIL_VERSION=${CIL_VERSION}
+                                                       PREFIX=${CMAKE_SOURCE_DIR}/Core 
+                                                       LIBRARY_INC=${CMAKE_SOURCE_DIR}/Core 
+                                                       LIBRARY_LIB=${CMAKE_BINARY_DIR}/Core
+                                                       ${PYTHON_EXECUTABLE} ${SETUP_PY} install
+                       COMMAND ${CMAKE_COMMAND} -E touch ${OUTPUT}
+                       DEPENDS cilreg)
+
+    else()
+      if (WIN32)
+        add_custom_command(OUTPUT ${OUTPUT}
+                       COMMAND ${CMAKE_COMMAND} -E copy_directory ${CMAKE_CURRENT_SOURCE_DIR}/src ${CMAKE_CURRENT_BINARY_DIR}/src
+                       COMMAND ${CMAKE_COMMAND} -E copy_directory ${CMAKE_CURRENT_SOURCE_DIR}/ccpi ${CMAKE_CURRENT_BINARY_DIR}/ccpi
+                       COMMAND ${CMAKE_COMMAND} -E env CIL_VERSION=${CIL_VERSION}
+                                                       PREFIX=${CMAKE_SOURCE_DIR}/Core 
+                                                       LIBRARY_INC=${CMAKE_SOURCE_DIR}/Core 
+                                                       LIBRARY_LIB=${CMAKE_BINARY_DIR}/Core/${CMAKE_BUILD_TYPE}
+                                                       ${PYTHON_EXECUTABLE} ${SETUP_PY} build_ext --inplace
+                       COMMAND ${CMAKE_COMMAND} -E touch ${OUTPUT}
+                       DEPENDS cilreg)
+      else()
+        add_custom_command(OUTPUT ${OUTPUT}
+                       COMMAND ${CMAKE_COMMAND} -E copy_directory ${CMAKE_CURRENT_SOURCE_DIR}/src ${CMAKE_CURRENT_BINARY_DIR}/src
+                       COMMAND ${CMAKE_COMMAND} -E copy_directory ${CMAKE_CURRENT_SOURCE_DIR}/ccpi ${CMAKE_CURRENT_BINARY_DIR}/ccpi
+                       COMMAND ${CMAKE_COMMAND} -E env CIL_VERSION=${CIL_VERSION}
+                                                       PREFIX=${CMAKE_SOURCE_DIR}/Core 
+                                                       LIBRARY_INC=${CMAKE_SOURCE_DIR}/Core 
+                                                       LIBRARY_LIB=${CMAKE_BINARY_DIR}/Core
+                                                       ${PYTHON_EXECUTABLE} ${SETUP_PY} build_ext --inplace
+                       COMMAND ${CMAKE_COMMAND} -E touch ${OUTPUT}
+                       DEPENDS cilreg)
+      endif()
+      install(DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/ccpi 
+              DESTINATION ${PYTHON_DEST})
+    endif()
+    
+    
+    add_custom_target(PythonWrapper ALL DEPENDS ${OUTPUT})
+
+    #install(CODE "execute_process(COMMAND ${PYTHON} ${SETUP_PY} install)")
+endif()
diff --git a/src/Python/ccpi/__init__.py b/src/Python/ccpi/__init__.py
new file mode 100644
index 0000000..e69de29
--- /dev/null
+++ b/src/Python/ccpi/__init__.py
diff --git a/src/Python/ccpi/filters/__init__.py b/src/Python/ccpi/filters/__init__.py
new file mode 100644
index 0000000..e69de29
--- /dev/null
+++ b/src/Python/ccpi/filters/__init__.py
diff --git a/src/Python/ccpi/filters/regularisers.py b/src/Python/ccpi/filters/regularisers.py
new file mode 100644
index 0000000..588ea32
--- /dev/null
+++ b/src/Python/ccpi/filters/regularisers.py
@@ -0,0 +1,214 @@
+"""
+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
+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
+except ImportError:
+    gpu_enabled = False    
+from ccpi.filters.cpu_regularisers import NDF_INPAINT_CPU, NVM_INPAINT_CPU
+
+def ROF_TV(inputData, regularisation_parameter, iterations,
+                     time_marching_parameter,device='cpu'):
+    if device == 'cpu':
+        return TV_ROF_CPU(inputData,
+                     regularisation_parameter,
+                     iterations, 
+                     time_marching_parameter)
+    elif device == 'gpu' and gpu_enabled:
+        return TV_ROF_GPU(inputData,
+                     regularisation_parameter,
+                     iterations, 
+                     time_marching_parameter)
+    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 FGP_TV(inputData, regularisation_parameter,iterations,
+                     tolerance_param, methodTV, nonneg, printM, device='cpu'):
+    if device == 'cpu':
+        return TV_FGP_CPU(inputData,
+                     regularisation_parameter,
+                     iterations, 
+                     tolerance_param,
+                     methodTV,
+                     nonneg,
+                     printM)
+    elif device == 'gpu' and gpu_enabled:
+        return TV_FGP_GPU(inputData,
+                     regularisation_parameter,
+                     iterations, 
+                     tolerance_param,
+                     methodTV,
+                     nonneg,
+                     printM)
+    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 SB_TV(inputData, regularisation_parameter, iterations,
+                     tolerance_param, methodTV, printM, device='cpu'):
+    if device == 'cpu':
+        return TV_SB_CPU(inputData,
+                     regularisation_parameter,
+                     iterations, 
+                     tolerance_param,
+                     methodTV,
+                     printM)
+    elif device == 'gpu' and gpu_enabled:
+        return TV_SB_GPU(inputData,
+                     regularisation_parameter,
+                     iterations, 
+                     tolerance_param,
+                     methodTV,
+                     printM)
+    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 FGP_dTV(inputData, refdata, regularisation_parameter, iterations,
+                     tolerance_param, eta_const, methodTV, nonneg, printM, device='cpu'):
+    if device == 'cpu':
+        return dTV_FGP_CPU(inputData,
+                     refdata,
+                     regularisation_parameter,
+                     iterations, 
+                     tolerance_param,
+                     eta_const,
+                     methodTV,
+                     nonneg,
+                     printM)
+    elif device == 'gpu' and gpu_enabled:
+        return dTV_FGP_GPU(inputData,
+                     refdata,
+                     regularisation_parameter,
+                     iterations, 
+                     tolerance_param,
+                     eta_const,
+                     methodTV,
+                     nonneg,
+                     printM)
+    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 TNV(inputData, regularisation_parameter, iterations, tolerance_param):
+        return TNV_CPU(inputData,
+                     regularisation_parameter,
+                     iterations, 
+                     tolerance_param)
+def NDF(inputData, regularisation_parameter, edge_parameter, iterations,
+                     time_marching_parameter, penalty_type, device='cpu'):
+    if device == 'cpu':
+        return NDF_CPU(inputData,
+                     regularisation_parameter,
+                     edge_parameter,
+                     iterations, 
+                     time_marching_parameter,
+                     penalty_type)
+    elif device == 'gpu' and gpu_enabled:
+        return NDF_GPU(inputData,
+                     regularisation_parameter,
+                     edge_parameter,
+                     iterations, 
+                     time_marching_parameter,
+                     penalty_type)
+    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, device='cpu'):
+    if device == 'cpu':
+        return Diff4th_CPU(inputData,
+                     regularisation_parameter,
+                     edge_parameter,
+                     iterations, 
+                     time_marching_parameter)
+    elif device == 'gpu' and gpu_enabled:
+        return Diff4th_GPU(inputData,
+                     regularisation_parameter,
+                     edge_parameter,
+                     iterations, 
+                     time_marching_parameter)
+    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 PatchSelect(inputData, searchwindow, patchwindow, neighbours, edge_parameter, device='cpu'):
+    if device == 'cpu':
+        return PATCHSEL_CPU(inputData,
+                     searchwindow,
+                     patchwindow,
+                     neighbours, 
+                     edge_parameter)
+    elif device == 'gpu' and gpu_enabled:
+        return PATCHSEL_GPU(inputData,
+                     searchwindow,
+                     patchwindow,
+                     neighbours, 
+                     edge_parameter)
+    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 NLTV(inputData, H_i, H_j, H_k, Weights, regularisation_parameter, iterations):
+    return NLTV_CPU(inputData,
+                     H_i,
+                     H_j,
+                     H_k, 
+                     Weights,
+                     regularisation_parameter,
+                     iterations)
+
+def TGV(inputData, regularisation_parameter, alpha1, alpha0, iterations,
+                     LipshitzConst, device='cpu'):
+    if device == 'cpu':
+        return TGV_CPU(inputData, 
+					regularisation_parameter, 
+					alpha1, 
+					alpha0, 
+					iterations,
+                    LipshitzConst)
+    elif device == 'gpu' and gpu_enabled:
+        return TGV_GPU(inputData, 
+					regularisation_parameter, 
+					alpha1, 
+					alpha0, 
+					iterations,
+                    LipshitzConst)
+    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 LLT_ROF(inputData, regularisation_parameterROF, regularisation_parameterLLT, iterations,
+                     time_marching_parameter, device='cpu'):
+    if device == 'cpu':
+        return LLT_ROF_CPU(inputData, regularisation_parameterROF, regularisation_parameterLLT, iterations, time_marching_parameter)
+    elif device == 'gpu' and gpu_enabled:
+        return LLT_ROF_GPU(inputData, regularisation_parameterROF, regularisation_parameterLLT, iterations, time_marching_parameter)
+    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 NDF_INP(inputData, maskData, regularisation_parameter, edge_parameter, iterations,
+                     time_marching_parameter, penalty_type):
+        return NDF_INPAINT_CPU(inputData, maskData, regularisation_parameter, 
+        edge_parameter, iterations, time_marching_parameter, penalty_type)
+        
+def NVM_INP(inputData, maskData, SW_increment, iterations):
+        return NVM_INPAINT_CPU(inputData, maskData, SW_increment, iterations)
diff --git a/src/Python/setup-regularisers.py.in b/src/Python/setup-regularisers.py.in
new file mode 100644
index 0000000..462edda
--- /dev/null
+++ b/src/Python/setup-regularisers.py.in
@@ -0,0 +1,75 @@
+#!/usr/bin/env python
+
+import setuptools
+from distutils.core import setup
+from distutils.extension import Extension
+from Cython.Distutils import build_ext
+
+import os
+import sys
+import numpy
+import platform	
+
+cil_version=os.environ['CIL_VERSION']
+if  cil_version == '':
+    print("Please set the environmental variable CIL_VERSION")
+    sys.exit(1)
+	
+library_include_path = ""
+library_lib_path = ""
+try:
+    library_include_path = os.environ['LIBRARY_INC']
+    library_lib_path = os.environ['LIBRARY_LIB']
+except:
+    library_include_path = os.environ['PREFIX']+'/include'
+    pass
+    
+extra_include_dirs = [numpy.get_include(), library_include_path]
+#extra_library_dirs = [os.path.join(library_include_path, "..", "lib")]
+extra_compile_args = []
+extra_library_dirs = [library_lib_path]
+extra_compile_args = []
+extra_link_args = []
+extra_libraries = ['cilreg']
+
+print ("extra_library_dirs " , extra_library_dirs)
+
+extra_include_dirs += [os.path.join(".." , ".." , "Core"),
+                       os.path.join(".." , ".." , "Core",  "regularisers_CPU"),
+                       os.path.join(".." , ".." , "Core",  "inpainters_CPU"),
+                       os.path.join(".." , ".." , "Core",  "regularisers_GPU" , "TV_FGP" ) , 
+                       os.path.join(".." , ".." , "Core",  "regularisers_GPU" , "TV_ROF" ) , 
+                       os.path.join(".." , ".." , "Core",  "regularisers_GPU" , "TV_SB" ) ,
+                       os.path.join(".." , ".." , "Core",  "regularisers_GPU" , "TGV" ) ,
+                       os.path.join(".." , ".." , "Core",  "regularisers_GPU" , "LLTROF" ) ,
+                       os.path.join(".." , ".." , "Core",  "regularisers_GPU" , "NDF" ) ,
+                       os.path.join(".." , ".." , "Core",  "regularisers_GPU" , "dTV_FGP" ) , 
+                       os.path.join(".." , ".." , "Core",  "regularisers_GPU" , "DIFF4th" ) , 
+                       os.path.join(".." , ".." , "Core",  "regularisers_GPU" , "PatchSelect" ) ,
+						   "."]
+
+if platform.system() == 'Windows':				   
+    extra_compile_args[0:] = ['/DWIN32','/EHsc','/DBOOST_ALL_NO_LIB' , '/openmp' ]   
+else:
+    extra_compile_args = ['-fopenmp','-O2', '-funsigned-char', '-Wall', '-std=c++0x']
+    extra_libraries += [@EXTRA_OMP_LIB@]
+    
+setup(
+    name='ccpi',
+	description='CCPi Core Imaging Library - Image regularisers',
+	version=cil_version,
+    cmdclass = {'build_ext': build_ext},
+    ext_modules = [Extension("ccpi.filters.cpu_regularisers",
+                             sources=[os.path.join("." , "src", "cpu_regularisers.pyx" ) ],
+                             include_dirs=extra_include_dirs, 
+							 library_dirs=extra_library_dirs, 
+							 extra_compile_args=extra_compile_args, 
+							 libraries=extra_libraries ), 
+    
+    ],
+	zip_safe = False,	
+	packages = {'ccpi','ccpi.filters'},
+)
+
+
+@SETUP_GPU_WRAPPERS@
diff --git a/src/Python/src/cpu_regularisers.pyx b/src/Python/src/cpu_regularisers.pyx
new file mode 100644
index 0000000..11a0617
--- /dev/null
+++ b/src/Python/src/cpu_regularisers.pyx
@@ -0,0 +1,685 @@
+# distutils: language=c++
+"""
+Copyright 2018 CCPi
+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.
+
+Author: Edoardo Pasca, Daniil Kazantsev
+"""
+
+import cython
+import numpy as np
+cimport numpy as np
+
+cdef extern float TV_ROF_CPU_main(float *Input, float *Output, float lambdaPar, int iterationsNumb, float tau, int dimX, int dimY, int dimZ);
+cdef extern float TV_FGP_CPU_main(float *Input, float *Output, float lambdaPar, int iterationsNumb, float epsil, int methodTV, int nonneg, int printM, int dimX, int dimY, int dimZ);
+cdef extern float SB_TV_CPU_main(float *Input, float *Output, float lambdaPar, int iterationsNumb, float epsil, int methodTV, int printM, int dimX, int dimY, int dimZ);
+cdef extern float LLT_ROF_CPU_main(float *Input, float *Output, float lambdaROF, float lambdaLLT, int iterationsNumb, float tau, int dimX, int dimY, int dimZ);
+cdef extern float TGV_main(float *Input, float *Output, float lambdaPar, float alpha1, float alpha0, int iterationsNumb, float L2, int dimX, int dimY, int dimZ);
+cdef extern float Diffusion_CPU_main(float *Input, float *Output, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int penaltytype, int dimX, int dimY, int dimZ);
+cdef extern float Diffus4th_CPU_main(float *Input, float *Output, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, 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);
+cdef extern float dTV_FGP_CPU_main(float *Input, float *InputRef, float *Output, float lambdaPar, int iterationsNumb, float epsil, float eta, int methodTV, int nonneg, int printM, int dimX, int dimY, int dimZ);
+cdef extern float PatchSelect_CPU_main(float *Input, unsigned short *H_i, unsigned short *H_j, unsigned short *H_k, float *Weights, int dimX, int dimY, int dimZ, int SearchWindow, int SimilarWin, int NumNeighb, float h, int switchM);
+cdef extern float Nonlocal_TV_CPU_main(float *A_orig, float *Output, unsigned short *H_i, unsigned short *H_j, unsigned short *H_k, float *Weights, int dimX, int dimY, int dimZ, int NumNeighb, float lambdaReg, int IterNumb);
+
+cdef extern float Diffusion_Inpaint_CPU_main(float *Input, unsigned char *Mask, float *Output, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int penaltytype, int dimX, int dimY, int dimZ);
+cdef extern float NonlocalMarching_Inpaint_main(float *Input, unsigned char *M, float *Output, unsigned char *M_upd, int SW_increment, int iterationsNumb, int trigger, int dimX, int dimY, int dimZ);
+cdef extern float TV_energy2D(float *U, float *U0, float *E_val, float lambdaPar, int type, int dimX, int dimY);
+cdef extern float TV_energy3D(float *U, float *U0, float *E_val, float lambdaPar, int type, int dimX, int dimY, int dimZ);
+#****************************************************************#
+#********************** Total-variation ROF *********************#
+#****************************************************************#
+def TV_ROF_CPU(inputData, regularisation_parameter, iterationsNumb, marching_step_parameter):
+    if inputData.ndim == 2:
+        return TV_ROF_2D(inputData, regularisation_parameter, iterationsNumb, marching_step_parameter)
+    elif inputData.ndim == 3:
+        return TV_ROF_3D(inputData, regularisation_parameter, iterationsNumb, marching_step_parameter)
+
+def TV_ROF_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
+                     float regularisation_parameter,
+                     int iterationsNumb,                     
+                     float marching_step_parameter):
+    cdef long dims[2]
+    dims[0] = inputData.shape[0]
+    dims[1] = inputData.shape[1]
+    
+    cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
+            np.zeros([dims[0],dims[1]], dtype='float32')
+                   
+    # Run ROF iterations for 2D data 
+    TV_ROF_CPU_main(&inputData[0,0], &outputData[0,0], regularisation_parameter, iterationsNumb, marching_step_parameter, dims[1], dims[0], 1)
+    
+    return outputData
+            
+def TV_ROF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
+                     float regularisation_parameter,
+                     int iterationsNumb,
+                     float marching_step_parameter):
+    cdef long dims[3]
+    dims[0] = inputData.shape[0]
+    dims[1] = inputData.shape[1]
+    dims[2] = inputData.shape[2]
+    
+    cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
+            np.zeros([dims[0],dims[1],dims[2]], dtype='float32')
+           
+    # Run ROF iterations for 3D data 
+    TV_ROF_CPU_main(&inputData[0,0,0], &outputData[0,0,0], regularisation_parameter, iterationsNumb, marching_step_parameter, dims[2], dims[1], dims[0])
+
+    return outputData
+
+#****************************************************************#
+#********************** Total-variation FGP *********************#
+#****************************************************************#
+#******** Total-variation Fast-Gradient-Projection (FGP)*********#
+def TV_FGP_CPU(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, printM):
+    if inputData.ndim == 2:
+        return TV_FGP_2D(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, printM)
+    elif inputData.ndim == 3:
+        return TV_FGP_3D(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, printM)
+
+def TV_FGP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
+                     float regularisation_parameter,
+                     int iterationsNumb, 
+                     float tolerance_param,
+                     int methodTV,
+                     int nonneg,
+                     int printM):
+                         
+    cdef long dims[2]
+    dims[0] = inputData.shape[0]
+    dims[1] = inputData.shape[1]
+    
+    cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
+            np.zeros([dims[0],dims[1]], dtype='float32')
+                   
+    #/* Run FGP-TV iterations for 2D data */
+    TV_FGP_CPU_main(&inputData[0,0], &outputData[0,0], regularisation_parameter, 
+                       iterationsNumb, 
+                       tolerance_param,
+                       methodTV,
+                       nonneg,
+                       printM,
+                       dims[1],dims[0],1)
+    
+    return outputData        
+            
+def TV_FGP_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
+                     float regularisation_parameter,
+                     int iterationsNumb, 
+                     float tolerance_param,
+                     int methodTV,
+                     int nonneg,
+                     int printM):
+    cdef long dims[3]
+    dims[0] = inputData.shape[0]
+    dims[1] = inputData.shape[1]
+    dims[2] = inputData.shape[2]
+    
+    cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
+            np.zeros([dims[0], dims[1], dims[2]], dtype='float32')
+           
+    #/* Run FGP-TV iterations for 3D data */
+    TV_FGP_CPU_main(&inputData[0,0,0], &outputData[0,0,0], regularisation_parameter,
+                       iterationsNumb, 
+                       tolerance_param,
+                       methodTV,
+                       nonneg,
+                       printM,
+                       dims[2], dims[1], dims[0])
+    return outputData 
+
+#***************************************************************#
+#********************** Total-variation SB *********************#
+#***************************************************************#
+#*************** Total-variation Split Bregman (SB)*************#
+def TV_SB_CPU(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, printM):
+    if inputData.ndim == 2:
+        return TV_SB_2D(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, printM)
+    elif inputData.ndim == 3:
+        return TV_SB_3D(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, printM)
+
+def TV_SB_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
+                     float regularisation_parameter,
+                     int iterationsNumb, 
+                     float tolerance_param,
+                     int methodTV,
+                     int printM):
+                         
+    cdef long dims[2]
+    dims[0] = inputData.shape[0]
+    dims[1] = inputData.shape[1]
+    
+    cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
+            np.zeros([dims[0],dims[1]], dtype='float32')
+                   
+    #/* Run SB-TV iterations for 2D data */
+    SB_TV_CPU_main(&inputData[0,0], &outputData[0,0], regularisation_parameter, 
+                       iterationsNumb, 
+                       tolerance_param,
+                       methodTV,
+                       printM,
+                       dims[1],dims[0],1)
+    
+    return outputData        
+            
+def TV_SB_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
+                     float regularisation_parameter,
+                     int iterationsNumb, 
+                     float tolerance_param,
+                     int methodTV,
+                     int printM):
+    cdef long dims[3]
+    dims[0] = inputData.shape[0]
+    dims[1] = inputData.shape[1]
+    dims[2] = inputData.shape[2]
+    
+    cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
+            np.zeros([dims[0], dims[1], dims[2]], dtype='float32')
+           
+    #/* Run SB-TV iterations for 3D data */
+    SB_TV_CPU_main(&inputData[0,0,0], &outputData[0,0,0], regularisation_parameter,
+                       iterationsNumb, 
+                       tolerance_param,
+                       methodTV,
+                       printM,
+                       dims[2], dims[1], dims[0])
+    return outputData 
+
+#***************************************************************#
+#***************** Total Generalised Variation *****************#
+#***************************************************************#
+def TGV_CPU(inputData, regularisation_parameter, alpha1, alpha0, iterations, LipshitzConst):
+    if inputData.ndim == 2:
+        return TGV_2D(inputData, regularisation_parameter, alpha1, alpha0, 
+                      iterations, LipshitzConst)
+    elif inputData.ndim == 3:
+        return TGV_3D(inputData, regularisation_parameter, alpha1, alpha0, 
+                      iterations, LipshitzConst)
+
+def TGV_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
+                     float regularisation_parameter,
+                     float alpha1,
+                     float alpha0,
+                     int iterationsNumb, 
+                     float LipshitzConst):
+                         
+    cdef long dims[2]
+    dims[0] = inputData.shape[0]
+    dims[1] = inputData.shape[1]
+    
+    cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
+            np.zeros([dims[0],dims[1]], dtype='float32')
+                   
+    #/* Run TGV iterations for 2D data */
+    TGV_main(&inputData[0,0], &outputData[0,0], regularisation_parameter, 
+                       alpha1,
+                       alpha0,
+                       iterationsNumb, 
+                       LipshitzConst,
+                       dims[1],dims[0],1)
+    return outputData
+def TGV_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
+                     float regularisation_parameter,
+                     float alpha1,
+                     float alpha0,
+                     int iterationsNumb, 
+                     float LipshitzConst):
+                         
+    cdef long dims[3]
+    dims[0] = inputData.shape[0]
+    dims[1] = inputData.shape[1]
+    dims[2] = inputData.shape[2]
+    
+    cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
+            np.zeros([dims[0], dims[1], dims[2]], dtype='float32')
+                   
+    #/* Run TGV iterations for 3D data */
+    TGV_main(&inputData[0,0,0], &outputData[0,0,0], regularisation_parameter, 
+                       alpha1,
+                       alpha0,
+                       iterationsNumb, 
+                       LipshitzConst,
+                       dims[2], dims[1], dims[0])
+    return outputData
+
+#***************************************************************#
+#******************* ROF - LLT regularisation ******************#
+#***************************************************************#
+def LLT_ROF_CPU(inputData, regularisation_parameterROF, regularisation_parameterLLT, iterations, time_marching_parameter):
+    if inputData.ndim == 2:
+        return LLT_ROF_2D(inputData, regularisation_parameterROF, regularisation_parameterLLT, iterations, time_marching_parameter)
+    elif inputData.ndim == 3:
+        return LLT_ROF_3D(inputData, regularisation_parameterROF, regularisation_parameterLLT, iterations, time_marching_parameter)
+
+def LLT_ROF_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
+                     float regularisation_parameterROF,
+                     float regularisation_parameterLLT,
+                     int iterations, 
+                     float time_marching_parameter):
+                         
+    cdef long dims[2]
+    dims[0] = inputData.shape[0]
+    dims[1] = inputData.shape[1]
+    
+    cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
+            np.zeros([dims[0],dims[1]], dtype='float32')
+                   
+    #/* Run ROF-LLT iterations for 2D data */
+    LLT_ROF_CPU_main(&inputData[0,0], &outputData[0,0], regularisation_parameterROF, regularisation_parameterLLT, iterations, time_marching_parameter, dims[1],dims[0],1)
+    return outputData
+
+def LLT_ROF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
+                     float regularisation_parameterROF,
+                     float regularisation_parameterLLT,
+                     int iterations, 
+                     float time_marching_parameter):
+						 
+    cdef long dims[3]
+    dims[0] = inputData.shape[0]
+    dims[1] = inputData.shape[1]
+    dims[2] = inputData.shape[2]
+    
+    cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
+            np.zeros([dims[0], dims[1], dims[2]], dtype='float32')
+           
+    #/* Run ROF-LLT iterations for 3D data */
+    LLT_ROF_CPU_main(&inputData[0,0,0], &outputData[0,0,0], regularisation_parameterROF, regularisation_parameterLLT, iterations, time_marching_parameter, dims[2], dims[1], dims[0])
+    return outputData 
+
+#****************************************************************#
+#**************Directional Total-variation FGP ******************#
+#****************************************************************#
+#******** Directional TV Fast-Gradient-Projection (FGP)*********#
+def dTV_FGP_CPU(inputData, refdata, regularisation_parameter, iterationsNumb, tolerance_param, eta_const, methodTV, nonneg, printM):
+    if inputData.ndim == 2:
+        return dTV_FGP_2D(inputData, refdata, regularisation_parameter, iterationsNumb, tolerance_param, eta_const, methodTV, nonneg, printM)
+    elif inputData.ndim == 3:
+        return dTV_FGP_3D(inputData, refdata, regularisation_parameter, iterationsNumb, tolerance_param, eta_const, methodTV, nonneg, printM)
+
+def dTV_FGP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
+               np.ndarray[np.float32_t, ndim=2, mode="c"] refdata,
+                     float regularisation_parameter,
+                     int iterationsNumb, 
+                     float tolerance_param,
+                     float eta_const,
+                     int methodTV,
+                     int nonneg,
+                     int printM):
+                         
+    cdef long dims[2]
+    dims[0] = inputData.shape[0]
+    dims[1] = inputData.shape[1]
+    
+    cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
+            np.zeros([dims[0],dims[1]], dtype='float32')
+                   
+    #/* Run FGP-dTV iterations for 2D data */
+    dTV_FGP_CPU_main(&inputData[0,0], &refdata[0,0], &outputData[0,0], regularisation_parameter, 
+                       iterationsNumb, 
+                       tolerance_param,
+                       eta_const,
+                       methodTV,                       
+                       nonneg,
+                       printM,
+                       dims[1], dims[0], 1)
+    
+    return outputData        
+            
+def dTV_FGP_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
+               np.ndarray[np.float32_t, ndim=3, mode="c"] refdata,
+                     float regularisation_parameter,
+                     int iterationsNumb, 
+                     float tolerance_param,
+                     float eta_const,
+                     int methodTV,
+                     int nonneg,
+                     int printM):
+    cdef long dims[3]
+    dims[0] = inputData.shape[0]
+    dims[1] = inputData.shape[1]
+    dims[2] = inputData.shape[2]
+    
+    cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
+            np.zeros([dims[0], dims[1], dims[2]], dtype='float32')
+           
+    #/* Run FGP-dTV iterations for 3D data */
+    dTV_FGP_CPU_main(&inputData[0,0,0], &refdata[0,0,0], &outputData[0,0,0], regularisation_parameter,
+                       iterationsNumb, 
+                       tolerance_param,
+                       eta_const,
+                       methodTV,
+                       nonneg,
+                       printM,
+                       dims[2], dims[1], dims[0])
+    return outputData
+    
+#****************************************************************#
+#*********************Total Nuclear Variation********************#
+#****************************************************************#
+def TNV_CPU(inputData, regularisation_parameter, iterationsNumb, tolerance_param):
+    if inputData.ndim == 2:
+        return 
+    elif inputData.ndim == 3:
+        return TNV_3D(inputData, regularisation_parameter, iterationsNumb, tolerance_param)
+
+def TNV_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
+                     float regularisation_parameter,
+                     int iterationsNumb,
+                     float tolerance_param):
+    cdef long dims[3]
+    dims[0] = inputData.shape[0]
+    dims[1] = inputData.shape[1]
+    dims[2] = inputData.shape[2]
+    
+    cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
+            np.zeros([dims[0],dims[1],dims[2]], dtype='float32')
+           
+    # Run TNV iterations for 3D (X,Y,Channels) data 
+    TNV_CPU_main(&inputData[0,0,0], &outputData[0,0,0], regularisation_parameter, iterationsNumb, tolerance_param, dims[2], dims[1], dims[0])
+    return outputData
+#****************************************************************#
+#***************Nonlinear (Isotropic) Diffusion******************#
+#****************************************************************#
+def NDF_CPU(inputData, regularisation_parameter, edge_parameter, iterationsNumb,time_marching_parameter, penalty_type):
+    if inputData.ndim == 2:
+        return NDF_2D(inputData, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type)
+    elif inputData.ndim == 3:
+        return NDF_3D(inputData, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type)
+
+def NDF_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
+                     float regularisation_parameter,
+                     float edge_parameter,
+                     int iterationsNumb,                     
+                     float time_marching_parameter,
+                     int penalty_type):
+    cdef long dims[2]
+    dims[0] = inputData.shape[0]
+    dims[1] = inputData.shape[1]
+    
+    cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
+            np.zeros([dims[0],dims[1]], dtype='float32')   
+    
+    # Run Nonlinear Diffusion iterations for 2D data 
+    Diffusion_CPU_main(&inputData[0,0], &outputData[0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type, dims[1], dims[0], 1)
+    return outputData
+            
+def NDF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
+                     float regularisation_parameter,
+                     float edge_parameter,
+                     int iterationsNumb,                     
+                     float time_marching_parameter,
+                     int penalty_type):
+    cdef long dims[3]
+    dims[0] = inputData.shape[0]
+    dims[1] = inputData.shape[1]
+    dims[2] = inputData.shape[2]
+    
+    cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
+            np.zeros([dims[0],dims[1],dims[2]], dtype='float32')
+    
+    # Run Nonlinear Diffusion iterations for  3D data 
+    Diffusion_CPU_main(&inputData[0,0,0], &outputData[0,0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type, dims[2], dims[1], dims[0])
+
+    return outputData
+
+#****************************************************************#
+#*************Anisotropic Fourth-Order diffusion*****************#
+#****************************************************************#
+def Diff4th_CPU(inputData, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter):
+    if inputData.ndim == 2:
+        return Diff4th_2D(inputData, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter)
+    elif inputData.ndim == 3:
+        return Diff4th_3D(inputData, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter)
+
+def Diff4th_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
+                     float regularisation_parameter,
+                     float edge_parameter,
+                     int iterationsNumb,                     
+                     float time_marching_parameter):
+    cdef long dims[2]
+    dims[0] = inputData.shape[0]
+    dims[1] = inputData.shape[1]
+    
+    cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
+            np.zeros([dims[0],dims[1]], dtype='float32')   
+    
+    # Run Anisotropic Fourth-Order diffusion for 2D data 
+    Diffus4th_CPU_main(&inputData[0,0], &outputData[0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, dims[1], dims[0], 1)
+    return outputData
+          
+def Diff4th_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
+                     float regularisation_parameter,
+                     float edge_parameter,
+                     int iterationsNumb,
+                     float time_marching_parameter):
+    cdef long dims[3]
+    dims[0] = inputData.shape[0]
+    dims[1] = inputData.shape[1]
+    dims[2] = inputData.shape[2]
+    
+    cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
+            np.zeros([dims[0],dims[1],dims[2]], dtype='float32')
+    
+    # Run Anisotropic Fourth-Order diffusion for  3D data 
+    Diffus4th_CPU_main(&inputData[0,0,0], &outputData[0,0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, dims[2], dims[1], dims[0])
+
+    return outputData
+
+#****************************************************************#
+#***************Patch-based weights calculation******************#
+#****************************************************************#
+def PATCHSEL_CPU(inputData, searchwindow, patchwindow, neighbours, edge_parameter):
+    if inputData.ndim == 2:
+        return PatchSel_2D(inputData, searchwindow, patchwindow, neighbours, edge_parameter)
+    elif inputData.ndim == 3:
+        return 1
+def PatchSel_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
+                     int searchwindow,
+                     int patchwindow,
+                     int neighbours,
+                     float edge_parameter):
+    cdef long dims[3]
+    dims[0] = neighbours
+    dims[1] = inputData.shape[0]
+    dims[2] = inputData.shape[1]
+    
+    
+    cdef np.ndarray[np.float32_t, ndim=3, mode="c"] Weights = \
+            np.zeros([dims[0], dims[1],dims[2]], dtype='float32')
+    
+    cdef np.ndarray[np.uint16_t, ndim=3, mode="c"] H_i = \
+            np.zeros([dims[0], dims[1],dims[2]], dtype='uint16')
+            
+    cdef np.ndarray[np.uint16_t, ndim=3, mode="c"] H_j = \
+            np.zeros([dims[0], dims[1],dims[2]], dtype='uint16')
+
+    # Run patch-based weight selection function
+    PatchSelect_CPU_main(&inputData[0,0], &H_j[0,0,0], &H_i[0,0,0], &H_i[0,0,0], &Weights[0,0,0], dims[2], dims[1], 0, searchwindow, patchwindow,  neighbours,  edge_parameter, 1)
+    return H_i, H_j, Weights
+"""
+def PatchSel_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
+                     int searchwindow,
+                     int patchwindow,
+                     int neighbours,
+                     float edge_parameter):
+    cdef long dims[4]
+    dims[0] = inputData.shape[0]
+    dims[1] = inputData.shape[1]
+    dims[2] = inputData.shape[2]
+    dims[3] = neighbours
+    
+    cdef np.ndarray[np.float32_t, ndim=4, mode="c"] Weights = \
+            np.zeros([dims[3],dims[0],dims[1],dims[2]], dtype='float32')
+    
+    cdef np.ndarray[np.uint16_t, ndim=4, mode="c"] H_i = \
+            np.zeros([dims[3],dims[0],dims[1],dims[2]], dtype='uint16')
+            
+    cdef np.ndarray[np.uint16_t, ndim=4, mode="c"] H_j = \
+            np.zeros([dims[3],dims[0],dims[1],dims[2]], dtype='uint16')
+            
+    cdef np.ndarray[np.uint16_t, ndim=4, mode="c"] H_k = \
+            np.zeros([dims[3],dims[0],dims[1],dims[2]], dtype='uint16')
+
+    # Run patch-based weight selection function
+    PatchSelect_CPU_main(&inputData[0,0,0], &H_i[0,0,0,0], &H_j[0,0,0,0], &H_k[0,0,0,0], &Weights[0,0,0,0], dims[2], dims[1], dims[0], searchwindow, patchwindow,  neighbours, edge_parameter, 1)
+    return H_i, H_j, H_k, Weights
+"""
+
+#****************************************************************#
+#***************Non-local Total Variation******************#
+#****************************************************************#
+def NLTV_CPU(inputData, H_i, H_j, H_k, Weights, regularisation_parameter, iterations):
+    if inputData.ndim == 2:
+        return NLTV_2D(inputData, H_i, H_j, Weights, regularisation_parameter, iterations)
+    elif inputData.ndim == 3:
+        return 1
+def NLTV_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
+                     np.ndarray[np.uint16_t, ndim=3, mode="c"] H_i,
+                     np.ndarray[np.uint16_t, ndim=3, mode="c"] H_j,
+                     np.ndarray[np.float32_t, ndim=3, mode="c"] Weights,
+                     float regularisation_parameter,
+                     int iterations):
+
+    cdef long dims[2]
+    dims[0] = inputData.shape[0]
+    dims[1] = inputData.shape[1]
+    neighbours = H_i.shape[0]
+    
+    cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
+            np.zeros([dims[0],dims[1]], dtype='float32')
+    
+    # Run nonlocal TV regularisation
+    Nonlocal_TV_CPU_main(&inputData[0,0], &outputData[0,0], &H_i[0,0,0], &H_j[0,0,0], &H_i[0,0,0], &Weights[0,0,0], dims[1], dims[0], 0, neighbours, regularisation_parameter, iterations)
+    return outputData
+
+#*********************Inpainting WITH****************************#
+#***************Nonlinear (Isotropic) Diffusion******************#
+#****************************************************************#
+def NDF_INPAINT_CPU(inputData, maskData, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type):
+    if inputData.ndim == 2:
+        return NDF_INP_2D(inputData, maskData, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type)
+    elif inputData.ndim == 3:
+        return NDF_INP_3D(inputData, maskData, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type)
+
+def NDF_INP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
+                     np.ndarray[np.uint8_t, ndim=2, mode="c"] maskData,
+                     float regularisation_parameter,
+                     float edge_parameter,
+                     int iterationsNumb,
+                     float time_marching_parameter,
+                     int penalty_type):
+
+    cdef long dims[2]
+    dims[0] = inputData.shape[0]
+    dims[1] = inputData.shape[1]
+
+
+    cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
+            np.zeros([dims[0],dims[1]], dtype='float32')
+    
+    # Run Inpaiting by Diffusion iterations for 2D data 
+    Diffusion_Inpaint_CPU_main(&inputData[0,0], &maskData[0,0], &outputData[0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type, dims[1], dims[0], 1)
+    return outputData
+            
+def NDF_INP_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
+                     np.ndarray[np.uint8_t, ndim=3, mode="c"] maskData,
+                     float regularisation_parameter,
+                     float edge_parameter,
+                     int iterationsNumb,
+                     float time_marching_parameter,
+                     int penalty_type):
+    cdef long dims[3]
+    dims[0] = inputData.shape[0]
+    dims[1] = inputData.shape[1]
+    dims[2] = inputData.shape[2]
+    
+    cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
+            np.zeros([dims[0],dims[1],dims[2]], dtype='float32')
+    
+    # Run Inpaiting by Diffusion iterations for 3D data 
+    Diffusion_Inpaint_CPU_main(&inputData[0,0,0], &maskData[0,0,0], &outputData[0,0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type, dims[2], dims[1], dims[0])
+
+    return outputData
+#*********************Inpainting WITH****************************#
+#***************Nonlocal Vertical Marching method****************#
+#****************************************************************#
+def NVM_INPAINT_CPU(inputData, maskData, SW_increment, iterationsNumb):
+    if inputData.ndim == 2:
+        return NVM_INP_2D(inputData, maskData, SW_increment, iterationsNumb)
+    elif inputData.ndim == 3:
+        return 
+
+def NVM_INP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
+               np.ndarray[np.uint8_t, ndim=2, mode="c"] maskData,
+                     int SW_increment,
+                     int iterationsNumb):
+    cdef long dims[2]
+    dims[0] = inputData.shape[0]
+    dims[1] = inputData.shape[1]
+    
+    cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
+            np.zeros([dims[0],dims[1]], dtype='float32')   
+    
+    cdef np.ndarray[np.uint8_t, ndim=2, mode="c"] maskData_upd = \
+            np.zeros([dims[0],dims[1]], dtype='uint8')
+    
+    # Run Inpaiting by Nonlocal vertical marching method for 2D data 
+    NonlocalMarching_Inpaint_main(&inputData[0,0], &maskData[0,0], &outputData[0,0], 
+                                  &maskData_upd[0,0],
+                                  SW_increment, iterationsNumb, 1, dims[1], dims[0], 1)
+    
+    return (outputData, maskData_upd)
+
+
+#****************************************************************#
+#***************Calculation of TV-energy functional**************#
+#****************************************************************#
+def TV_ENERGY(inputData, inputData0, regularisation_parameter, typeFunctional):
+    if inputData.ndim == 2:
+        return TV_ENERGY_2D(inputData, inputData0, regularisation_parameter, typeFunctional)
+    elif inputData.ndim == 3:
+        return TV_ENERGY_3D(inputData, inputData0, regularisation_parameter, typeFunctional)
+
+def TV_ENERGY_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
+                 np.ndarray[np.float32_t, ndim=2, mode="c"] inputData0, 
+                     float regularisation_parameter,
+                     int typeFunctional):
+    
+    cdef long dims[2]
+    dims[0] = inputData.shape[0]
+    dims[1] = inputData.shape[1]
+    
+    cdef np.ndarray[np.float32_t, ndim=1, mode="c"] outputData = \
+            np.zeros([1], dtype='float32')
+                   
+    # run function    
+    TV_energy2D(&inputData[0,0], &inputData0[0,0], &outputData[0], regularisation_parameter, typeFunctional, dims[1], dims[0])
+    
+    return outputData
+            
+def TV_ENERGY_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
+                 np.ndarray[np.float32_t, ndim=3, mode="c"] inputData0, 
+                     float regularisation_parameter,
+                     int typeFunctional):
+						 
+    cdef long dims[3]
+    dims[0] = inputData.shape[0]
+    dims[1] = inputData.shape[1]
+    dims[2] = inputData.shape[2]
+    
+    cdef np.ndarray[np.float32_t, ndim=1, mode="c"] outputData = \
+            np.zeros([1], dtype='float32')
+           
+    # Run function
+    TV_energy3D(&inputData[0,0,0], &inputData0[0,0,0], &outputData[0], regularisation_parameter, typeFunctional, dims[2], dims[1], dims[0])
+
+    return outputData
diff --git a/src/Python/src/gpu_regularisers.pyx b/src/Python/src/gpu_regularisers.pyx
new file mode 100644
index 0000000..b52f669
--- /dev/null
+++ b/src/Python/src/gpu_regularisers.pyx
@@ -0,0 +1,640 @@
+# distutils: language=c++
+"""
+Copyright 2018 CCPi
+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.
+
+Author: Edoardo Pasca, Daniil Kazantsev
+"""
+
+import cython
+import numpy as np
+cimport numpy as np
+
+CUDAErrorMessage = 'CUDA error'
+
+cdef extern int TV_ROF_GPU_main(float* Input, float* Output, float lambdaPar, int iter, float tau, int N, int M, int Z);
+cdef extern int TV_FGP_GPU_main(float *Input, float *Output, float lambdaPar, int iter, float epsil, int methodTV, int nonneg, int printM, int N, int M, int Z);
+cdef extern int TV_SB_GPU_main(float *Input, float *Output, float lambdaPar, int iter, float epsil, int methodTV, int printM, int N, int M, int Z);
+cdef extern int TGV_GPU_main(float *Input, float *Output, float lambdaPar, float alpha1, float alpha0, int iterationsNumb, float L2, int dimX, int dimY, int dimZ);
+cdef extern int LLT_ROF_GPU_main(float *Input, float *Output, float lambdaROF, float lambdaLLT, int iterationsNumb, float tau, int N, int M, int Z);
+cdef extern int NonlDiff_GPU_main(float *Input, float *Output, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int penaltytype, int N, int M, int Z);
+cdef extern int dTV_FGP_GPU_main(float *Input, float *InputRef, float *Output, float lambdaPar, int iterationsNumb, float epsil, float eta, int methodTV, int nonneg, int printM, int N, int M, int Z);
+cdef extern int Diffus4th_GPU_main(float *Input, float *Output, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int N, int M, int Z);
+cdef extern int PatchSelect_GPU_main(float *Input, unsigned short *H_i, unsigned short *H_j, float *Weights, int N, int M, int SearchWindow, int SimilarWin, int NumNeighb, float h);
+
+# Total-variation Rudin-Osher-Fatemi (ROF)
+def TV_ROF_GPU(inputData,
+                     regularisation_parameter,
+                     iterations, 
+                     time_marching_parameter):
+    if inputData.ndim == 2:
+        return ROFTV2D(inputData, 
+                     regularisation_parameter,
+                     iterations,
+                     time_marching_parameter)
+    elif inputData.ndim == 3:
+        return ROFTV3D(inputData, 
+                     regularisation_parameter,
+                     iterations, 
+                     time_marching_parameter)
+                     
+# Total-variation Fast-Gradient-Projection (FGP)
+def TV_FGP_GPU(inputData,
+                     regularisation_parameter,
+                     iterations, 
+                     tolerance_param,
+                     methodTV,
+                     nonneg,
+                     printM):
+    if inputData.ndim == 2:
+        return FGPTV2D(inputData,
+                     regularisation_parameter,
+                     iterations, 
+                     tolerance_param,
+                     methodTV,
+                     nonneg,
+                     printM)
+    elif inputData.ndim == 3:
+        return FGPTV3D(inputData,
+                     regularisation_parameter,
+                     iterations, 
+                     tolerance_param,
+                     methodTV,
+                     nonneg,
+                     printM)
+# Total-variation Split Bregman (SB)
+def TV_SB_GPU(inputData,
+                     regularisation_parameter,
+                     iterations, 
+                     tolerance_param,
+                     methodTV,
+                     printM):
+    if inputData.ndim == 2:
+        return SBTV2D(inputData,
+                     regularisation_parameter,
+                     iterations, 
+                     tolerance_param,
+                     methodTV,
+                     printM)
+    elif inputData.ndim == 3:
+        return SBTV3D(inputData,
+                     regularisation_parameter,
+                     iterations, 
+                     tolerance_param,
+                     methodTV,
+                     printM)
+# LLT-ROF model
+def LLT_ROF_GPU(inputData, regularisation_parameterROF, regularisation_parameterLLT, iterations, time_marching_parameter):
+    if inputData.ndim == 2:
+        return LLT_ROF_GPU2D(inputData, regularisation_parameterROF, regularisation_parameterLLT, iterations, time_marching_parameter)
+    elif inputData.ndim == 3:
+        return LLT_ROF_GPU3D(inputData, regularisation_parameterROF, regularisation_parameterLLT, iterations, time_marching_parameter)
+# Total Generilised Variation (TGV)
+def TGV_GPU(inputData, regularisation_parameter, alpha1, alpha0, iterations, LipshitzConst):
+    if inputData.ndim == 2:
+        return TGV2D(inputData, regularisation_parameter, alpha1, alpha0, iterations, LipshitzConst)
+    elif inputData.ndim == 3:
+        return TGV3D(inputData, regularisation_parameter, alpha1, alpha0, iterations, LipshitzConst)
+# Directional Total-variation Fast-Gradient-Projection (FGP)
+def dTV_FGP_GPU(inputData,
+                     refdata,
+                     regularisation_parameter,
+                     iterations, 
+                     tolerance_param,
+                     eta_const,
+                     methodTV,
+                     nonneg,
+                     printM):
+    if inputData.ndim == 2:
+        return FGPdTV2D(inputData,
+                     refdata,
+                     regularisation_parameter,
+                     iterations, 
+                     tolerance_param,
+                     eta_const,
+                     methodTV,
+                     nonneg,
+                     printM)
+    elif inputData.ndim == 3:
+        return FGPdTV3D(inputData,
+                     refdata,
+                     regularisation_parameter,
+                     iterations, 
+                     tolerance_param,
+                     eta_const,
+                     methodTV,
+                     nonneg,
+                     printM)
+# Nonlocal Isotropic Diffusion (NDF)
+def NDF_GPU(inputData,
+                     regularisation_parameter,
+                     edge_parameter,
+                     iterations, 
+                     time_marching_parameter,
+                     penalty_type):
+    if inputData.ndim == 2:
+        return NDF_GPU_2D(inputData,
+                     regularisation_parameter,
+                     edge_parameter,
+                     iterations, 
+                     time_marching_parameter,
+                     penalty_type)
+    elif inputData.ndim == 3:
+        return NDF_GPU_3D(inputData,
+                     regularisation_parameter,
+                     edge_parameter,
+                     iterations, 
+                     time_marching_parameter,
+                     penalty_type)
+# Anisotropic Fourth-Order diffusion
+def Diff4th_GPU(inputData,
+                     regularisation_parameter,
+                     edge_parameter,
+                     iterations, 
+                     time_marching_parameter):
+    if inputData.ndim == 2:
+        return Diff4th_2D(inputData,
+                     regularisation_parameter,
+                     edge_parameter,
+                     iterations, 
+                     time_marching_parameter)
+    elif inputData.ndim == 3:
+        return Diff4th_3D(inputData,
+                     regularisation_parameter,
+                     edge_parameter,
+                     iterations, 
+                     time_marching_parameter)
+                     
+#****************************************************************#
+#********************** Total-variation ROF *********************#
+#****************************************************************#
+def ROFTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
+                     float regularisation_parameter,
+                     int iterations, 
+                     float time_marching_parameter):
+    
+    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')
+          
+    # Running CUDA code here
+    if (TV_ROF_GPU_main(
+            &inputData[0,0], &outputData[0,0], 
+                       regularisation_parameter,
+                       iterations , 
+                       time_marching_parameter, 
+                       dims[1], dims[0], 1)==0):
+        return outputData;
+    else:
+        raise ValueError(CUDAErrorMessage);
+    
+def ROFTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
+                     float regularisation_parameter,
+                     int iterations, 
+                     float time_marching_parameter):
+    
+    cdef long dims[3]
+    dims[0] = inputData.shape[0]
+    dims[1] = inputData.shape[1]
+    dims[2] = inputData.shape[2]
+
+    cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
+		    np.zeros([dims[0],dims[1],dims[2]], dtype='float32')
+          
+    # Running CUDA code here    
+    if (TV_ROF_GPU_main(
+            &inputData[0,0,0], &outputData[0,0,0], 
+                       regularisation_parameter,
+                       iterations , 
+                       time_marching_parameter, 
+                       dims[2], dims[1], dims[0])==0):
+        return outputData;
+    else:
+        raise ValueError(CUDAErrorMessage);
+#****************************************************************#
+#********************** Total-variation FGP *********************#
+#****************************************************************#
+#******** Total-variation Fast-Gradient-Projection (FGP)*********#
+def FGPTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
+                     float regularisation_parameter,
+                     int iterations, 
+                     float tolerance_param,
+                     int methodTV,
+                     int nonneg,
+                     int printM):
+    
+    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')
+          
+    # Running CUDA code here    
+    if (TV_FGP_GPU_main(&inputData[0,0], &outputData[0,0],
+                       regularisation_parameter, 
+                       iterations, 
+                       tolerance_param,
+                       methodTV,
+                       nonneg,
+                       printM,
+                       dims[1], dims[0], 1)==0):
+        return outputData;
+    else:
+        raise ValueError(CUDAErrorMessage);
+
+    
+def FGPTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
+                     float regularisation_parameter,
+                     int iterations, 
+                     float tolerance_param,
+                     int methodTV,
+                     int nonneg,
+                     int printM):
+    
+    cdef long dims[3]
+    dims[0] = inputData.shape[0]
+    dims[1] = inputData.shape[1]
+    dims[2] = inputData.shape[2]
+
+    cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
+		    np.zeros([dims[0],dims[1],dims[2]], dtype='float32')
+          
+    # Running CUDA code here    
+    if (TV_FGP_GPU_main(&inputData[0,0,0], &outputData[0,0,0],
+                       regularisation_parameter , 
+                       iterations, 
+                       tolerance_param,
+                       methodTV,
+                       nonneg,
+                       printM,
+                       dims[2], dims[1], dims[0])==0):
+        return outputData;
+    else:
+        raise ValueError(CUDAErrorMessage);
+
+#***************************************************************#
+#********************** Total-variation SB *********************#
+#***************************************************************#
+#*************** Total-variation Split Bregman (SB)*************#
+def SBTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
+                     float regularisation_parameter,
+                     int iterations, 
+                     float tolerance_param,
+                     int methodTV,
+                     int printM):
+    
+    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')
+          
+    # Running CUDA code here    
+    if (TV_SB_GPU_main(&inputData[0,0], &outputData[0,0],
+                       regularisation_parameter, 
+                       iterations, 
+                       tolerance_param,
+                       methodTV,
+                       printM,
+                       dims[1], dims[0], 1)==0):
+        return outputData;
+    else:
+        raise ValueError(CUDAErrorMessage);
+
+    
+def SBTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
+                     float regularisation_parameter,
+                     int iterations, 
+                     float tolerance_param,
+                     int methodTV,
+                     int printM):
+    
+    cdef long dims[3]
+    dims[0] = inputData.shape[0]
+    dims[1] = inputData.shape[1]
+    dims[2] = inputData.shape[2]
+
+    cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
+		    np.zeros([dims[0],dims[1],dims[2]], dtype='float32')
+          
+    # Running CUDA code here    
+    if (TV_SB_GPU_main(&inputData[0,0,0], &outputData[0,0,0],
+                       regularisation_parameter , 
+                       iterations, 
+                       tolerance_param,
+                       methodTV,
+                       printM,
+                       dims[2], dims[1], dims[0])==0):
+        return outputData;
+    else:
+        raise ValueError(CUDAErrorMessage);
+
+
+#***************************************************************#
+#************************ LLT-ROF model ************************#
+#***************************************************************#
+#************Joint LLT-ROF model for higher order **************#
+def LLT_ROF_GPU2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
+                     float regularisation_parameterROF,
+                     float regularisation_parameterLLT,
+                     int iterations, 
+                     float time_marching_parameter):
+    
+    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')
+          
+    # Running CUDA code here    
+    if (LLT_ROF_GPU_main(&inputData[0,0], &outputData[0,0],regularisation_parameterROF, regularisation_parameterLLT, iterations, time_marching_parameter, dims[1],dims[0],1)==0):
+        return outputData;
+    else:
+        raise ValueError(CUDAErrorMessage);
+
+    
+def LLT_ROF_GPU3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
+                     float regularisation_parameterROF,
+                     float regularisation_parameterLLT,
+                     int iterations, 
+                     float time_marching_parameter):
+    
+    cdef long dims[3]
+    dims[0] = inputData.shape[0]
+    dims[1] = inputData.shape[1]
+    dims[2] = inputData.shape[2]
+
+    cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
+		    np.zeros([dims[0],dims[1],dims[2]], dtype='float32')
+          
+    # Running CUDA code here    
+    if (LLT_ROF_GPU_main(&inputData[0,0,0], &outputData[0,0,0], regularisation_parameterROF, regularisation_parameterLLT, iterations, time_marching_parameter, dims[2], dims[1], dims[0])==0):
+        return outputData;
+    else:
+        raise ValueError(CUDAErrorMessage);
+
+
+#***************************************************************#
+#***************** Total Generalised Variation *****************#
+#***************************************************************#
+def TGV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
+                     float regularisation_parameter,
+                     float alpha1,
+                     float alpha0,
+                     int iterationsNumb, 
+                     float LipshitzConst):
+                         
+    cdef long dims[2]
+    dims[0] = inputData.shape[0]
+    dims[1] = inputData.shape[1]
+    
+    cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
+            np.zeros([dims[0],dims[1]], dtype='float32')
+                   
+    #/* Run TGV iterations for 2D data */
+    if (TGV_GPU_main(&inputData[0,0], &outputData[0,0], regularisation_parameter,
+                       alpha1,
+                       alpha0,
+                       iterationsNumb, 
+                       LipshitzConst,
+                       dims[1],dims[0], 1)==0):
+        return outputData
+    else:
+        raise ValueError(CUDAErrorMessage);
+
+def TGV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
+                     float regularisation_parameter,
+                     float alpha1,
+                     float alpha0,
+                     int iterationsNumb, 
+                     float LipshitzConst):
+    
+    cdef long dims[3]
+    dims[0] = inputData.shape[0]
+    dims[1] = inputData.shape[1]
+    dims[2] = inputData.shape[2]
+
+    cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
+		    np.zeros([dims[0],dims[1],dims[2]], dtype='float32')
+          
+    # Running CUDA code here    
+    if (TGV_GPU_main(
+            &inputData[0,0,0], &outputData[0,0,0], regularisation_parameter,
+                       alpha1,
+                       alpha0,
+                       iterationsNumb, 
+                       LipshitzConst,
+                       dims[2], dims[1], dims[0])==0):
+        return outputData;
+    else:
+        raise ValueError(CUDAErrorMessage);
+
+
+#****************************************************************#
+#**************Directional Total-variation FGP ******************#
+#****************************************************************#
+#******** Directional TV Fast-Gradient-Projection (FGP)*********#
+def FGPdTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
+             np.ndarray[np.float32_t, ndim=2, mode="c"] refdata,
+                     float regularisation_parameter,
+                     int iterations, 
+                     float tolerance_param,
+                     float eta_const,
+                     int methodTV,
+                     int nonneg,
+                     int printM):
+    
+    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')
+          
+    # Running CUDA code here    
+    if (dTV_FGP_GPU_main(&inputData[0,0], &refdata[0,0], &outputData[0,0],
+                       regularisation_parameter, 
+                       iterations, 
+                       tolerance_param,
+                       eta_const,
+                       methodTV,
+                       nonneg,
+                       printM,
+                       dims[1], dims[0], 1)==0):
+        return outputData
+    else:
+        raise ValueError(CUDAErrorMessage);
+
+    
+def FGPdTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
+             np.ndarray[np.float32_t, ndim=3, mode="c"] refdata, 
+                     float regularisation_parameter,
+                     int iterations, 
+                     float tolerance_param,
+                     float eta_const,
+                     int methodTV,
+                     int nonneg,
+                     int printM):
+    
+    cdef long dims[3]
+    dims[0] = inputData.shape[0]
+    dims[1] = inputData.shape[1]
+    dims[2] = inputData.shape[2]
+
+    cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
+		    np.zeros([dims[0],dims[1],dims[2]], dtype='float32')
+          
+    # Running CUDA code here    
+    if (dTV_FGP_GPU_main(&inputData[0,0,0], &refdata[0,0,0], &outputData[0,0,0],
+                       regularisation_parameter , 
+                       iterations, 
+                       tolerance_param,
+                       eta_const,
+                       methodTV,
+                       nonneg,
+                       printM,
+                       dims[2], dims[1], dims[0])==0):
+        return outputData;
+    else:
+        raise ValueError(CUDAErrorMessage);
+
+
+#****************************************************************#
+#***************Nonlinear (Isotropic) Diffusion******************#
+#****************************************************************#
+def NDF_GPU_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
+                     float regularisation_parameter,
+                     float edge_parameter,
+                     int iterationsNumb,                     
+                     float time_marching_parameter,
+                     int penalty_type):
+    cdef long dims[2]
+    dims[0] = inputData.shape[0]
+    dims[1] = inputData.shape[1]
+    
+    cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
+            np.zeros([dims[0],dims[1]], dtype='float32')
+    
+    #rangecheck = penalty_type < 1 and penalty_type > 3
+    #if not rangecheck:
+#        raise ValueError('Choose penalty type as 1 for Huber, 2 - Perona-Malik, 3 - Tukey Biweight')
+    
+    # Run Nonlinear Diffusion iterations for 2D data 
+    # Running CUDA code here  
+    if (NonlDiff_GPU_main(&inputData[0,0], &outputData[0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type, dims[1], dims[0], 1)==0):
+        return outputData;
+    else:
+        raise ValueError(CUDAErrorMessage);
+
+            
+def NDF_GPU_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
+                     float regularisation_parameter,
+                     float edge_parameter,
+                     int iterationsNumb,                     
+                     float time_marching_parameter,
+                     int penalty_type):
+    cdef long dims[3]
+    dims[0] = inputData.shape[0]
+    dims[1] = inputData.shape[1]
+    dims[2] = inputData.shape[2]
+    
+    cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
+            np.zeros([dims[0],dims[1],dims[2]], dtype='float32')    
+       
+    # Run Nonlinear Diffusion iterations for  3D data 
+    # Running CUDA code here  
+    if (NonlDiff_GPU_main(&inputData[0,0,0], &outputData[0,0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type, dims[2], dims[1], dims[0])==0):
+        return outputData;
+    else:
+        raise ValueError(CUDAErrorMessage);
+
+#****************************************************************#
+#************Anisotropic Fourth-Order diffusion******************#
+#****************************************************************#
+def Diff4th_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
+                     float regularisation_parameter,
+                     float edge_parameter,
+                     int iterationsNumb,
+                     float time_marching_parameter):
+    cdef long dims[2]
+    dims[0] = inputData.shape[0]
+    dims[1] = inputData.shape[1]
+    
+    cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
+            np.zeros([dims[0],dims[1]], dtype='float32')
+    
+    # Run Anisotropic Fourth-Order diffusion for 2D data 
+    # Running CUDA code here  
+    if (Diffus4th_GPU_main(&inputData[0,0], &outputData[0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, dims[1], dims[0], 1)==0):
+        return outputData
+    else:
+        raise ValueError(CUDAErrorMessage);
+
+            
+def Diff4th_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
+                     float regularisation_parameter,
+                     float edge_parameter,
+                     int iterationsNumb,
+                     float time_marching_parameter):
+    cdef long dims[3]
+    dims[0] = inputData.shape[0]
+    dims[1] = inputData.shape[1]
+    dims[2] = inputData.shape[2]
+    
+    cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
+            np.zeros([dims[0],dims[1],dims[2]], dtype='float32')    
+       
+    # Run Anisotropic Fourth-Order diffusion for  3D data 
+    # Running CUDA code here  
+    if (Diffus4th_GPU_main(&inputData[0,0,0], &outputData[0,0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, dims[2], dims[1], dims[0])==0):
+        return outputData;
+    else:
+        raise ValueError(CUDAErrorMessage);
+
+#****************************************************************#
+#************Patch-based weights pre-selection******************#
+#****************************************************************#
+def PATCHSEL_GPU(inputData, searchwindow, patchwindow, neighbours, edge_parameter):
+    if inputData.ndim == 2:
+        return PatchSel_2D(inputData, searchwindow, patchwindow, neighbours, edge_parameter)
+    elif inputData.ndim == 3:
+        return 1
+def PatchSel_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
+                     int searchwindow,
+                     int patchwindow,
+                     int neighbours,
+                     float edge_parameter):
+    cdef long dims[3]
+    dims[0] = neighbours
+    dims[1] = inputData.shape[0]
+    dims[2] = inputData.shape[1]    
+    
+    cdef np.ndarray[np.float32_t, ndim=3, mode="c"] Weights = \
+            np.zeros([dims[0], dims[1],dims[2]], dtype='float32')
+    
+    cdef np.ndarray[np.uint16_t, ndim=3, mode="c"] H_i = \
+            np.zeros([dims[0], dims[1],dims[2]], dtype='uint16')
+            
+    cdef np.ndarray[np.uint16_t, ndim=3, mode="c"] H_j = \
+            np.zeros([dims[0], dims[1],dims[2]], dtype='uint16')
+
+    # Run patch-based weight selection function
+    if (PatchSelect_GPU_main(&inputData[0,0], &H_j[0,0,0], &H_i[0,0,0], &Weights[0,0,0], dims[2], dims[1], searchwindow, patchwindow,  neighbours,  edge_parameter)==0):
+        return H_i, H_j, Weights;
+    else:
+        raise ValueError(CUDAErrorMessage);
+