restructured sources

29 files changed, 2373 insertions, 0 deletions

diff --git a/src/Matlab/CMakeLists.txt b/src/Matlab/CMakeLists.txt
new file mode 100755
index 0000000..b97f845
--- /dev/null
+++ b/src/Matlab/CMakeLists.txt
@@ -0,0 +1,147 @@
+project(regulariserMatlab)

+

+

+find_package(Matlab REQUIRED COMPONENTS MAIN_PROGRAM MX_LIBRARY ENG_LIBRARY )

+

+

+

+#C:\Users\ofn77899\Documents\Projects\CCPi\GitHub\CCPi-FISTA_Reconstruction\Core\regularisers_CPU

+# matlab_add_mex(

+    # NAME CPU_ROF

+    # SRC 

+      # ${CMAKE_SOURCE_DIR}/Matlab/mex_compile/regularisers_CPU/ROF_TV.c

+    # LINK_TO cilreg ${Matlab_LIBRARIES}

+    # )

+    

+# target_include_directories(CPU_ROF 

+   # PUBLIC ${CMAKE_SOURCE_DIR}/Core/regularisers_CPU

+   # ${CMAKE_SOURCE_DIR}/Core/regularisers_GPU

+   # ${CMAKE_SOURCE_DIR}/Core/inpainters_CPU

+   # ${CMAKE_SOURCE_DIR}/Core/

+   # ${MATLAB_INCLUDE_DIR})

+   

+   # matlab_add_mex(

+    # NAME CPU_TNV

+    # SRC 

+      # ${CMAKE_SOURCE_DIR}/Matlab/mex_compile/regularisers_CPU/TNV.c 

+    # LINK_TO cilreg ${Matlab_LIBRARIES}

+    # )

+    

+# target_include_directories(CPU_TNV 

+   # PUBLIC ${CMAKE_SOURCE_DIR}/Core/regularisers_CPU

+   # ${CMAKE_SOURCE_DIR}/Core/regularisers_GPU

+   # ${CMAKE_SOURCE_DIR}/Core/inpainters_CPU

+   # ${CMAKE_SOURCE_DIR}/Core/

+   # ${MATLAB_INCLUDE_DIR})

+   

+#set (CPU_MEX_FILES "regularisers_CPU/TNV.c;regularisers_CPU/ROF_TV.c")

+#set (MEX_TARGETS "CPU_TNV;CPU_ROF")

+#list(APPEND MEX_TARGETS "CPU_TNV")

+#list(APPEND MEX_TARGETS "CPU_ROF")

+

+file(GLOB CPU_MEX_FILES

+    "${CMAKE_SOURCE_DIR}/Matlab/mex_compile/regularisers_CPU/*.c"

+    #"${CMAKE_SOURCE_DIR}/Matlab/mex_compile/regularisers_GPU/*.c"

+)

+

+#message("CPU_MEX_FILES " ${CPU_MEX_FILES})

+

+list(LENGTH CPU_MEX_FILES num)

+

+

+MATH(EXPR num "${num}-1")

+#set(num "-1")

+message("found ${num} files")

+

+foreach(tgt RANGE 0 ${num})

+  message("number " ${tgt})

+  list(LENGTH CPU_MEX_FILES num2)

+  message("the list is ${num2}")

+  #list(GET CPU_TARGETS ${tgt} current_target)

+  list(GET CPU_MEX_FILES ${tgt} current_file_name)

+  get_filename_component(current_file ${current_file_name} NAME)

+  string(REGEX MATCH "(.+).c" match ${current_file})

+  if (NOT ${match} EQUAL "" )

+  set (current_target ${CMAKE_MATCH_1})

+  endif()

+  message("matlab_add_mex target " ${current_file} " and " ${current_target})

+  matlab_add_mex(

+    NAME ${current_target}

+    SRC 

+      ${current_file_name} 

+            #${CMAKE_SOURCE_DIR}/Core/regularisers_CPU/FGP_TV_core.c

+	    #${CMAKE_SOURCE_DIR}/Core/regularisers_CPU/SB_TV_core.c

+	    #${CMAKE_SOURCE_DIR}/Core/regularisers_CPU/TGV_core.c

+	    #${CMAKE_SOURCE_DIR}/Core/regularisers_CPU/Diffusion_core.c

+	    #${CMAKE_SOURCE_DIR}/Core/regularisers_CPU/Diffus4th_order_core.c

+	    #${CMAKE_SOURCE_DIR}/Core/regularisers_CPU/LLT_ROF_core.c

+            #${CMAKE_SOURCE_DIR}/Core/regularisers_CPU/ROF_TV_core.c

+            #${CMAKE_SOURCE_DIR}/Core/regularisers_CPU/FGP_dTV_core.c

+            #${CMAKE_SOURCE_DIR}/Core/regularisers_CPU/TNV_core.c

+            #${CMAKE_SOURCE_DIR}/Core/regularisers_CPU/utils.c

+	    #${CMAKE_SOURCE_DIR}/Core/inpainters_CPU/Diffusion_Inpaint_core.c

+	    #${CMAKE_SOURCE_DIR}/Core/inpainters_CPU/NonlocalMarching_Inpaint_core.c

+    LINK_TO cilreg ${Matlab_LIBRARIES}

+    )

+    

+target_include_directories(${current_target}

+   PUBLIC ${CMAKE_SOURCE_DIR}/Core/regularisers_CPU

+   ${CMAKE_SOURCE_DIR}/Core/regularisers_GPU

+   ${CMAKE_SOURCE_DIR}/Core/inpainters_CPU

+   ${CMAKE_SOURCE_DIR}/Core/

+   ${MATLAB_INCLUDE_DIR})

+   set_property(TARGET ${current_target} PROPERTY C_STANDARD 99)

+   list(APPEND CPU_MEX_TARGETS ${current_target})

+   INSTALL(TARGETS ${current_target} DESTINATION "${MATLAB_DEST}")

+endforeach()

+   

+add_custom_target(MatlabWrapper DEPENDS ${CPU_MEX_TARGETS})

+

+if (BUILD_CUDA)

+    find_package(CUDA)

+    if (CUDA_FOUND)

+      file(GLOB GPU_MEX_FILES

+        "${CMAKE_SOURCE_DIR}/Matlab/mex_compile/regularisers_GPU/*.cpp"

+      )

+

+      list(LENGTH GPU_MEX_FILES num)

+message("number of GPU files  " ${num})

+

+      MATH(EXPR num "${num}-1")

+    #set(num "-1")

+

+      foreach(tgt RANGE ${num})

+        message("number " ${tgt})

+  list(LENGTH GPU_MEX_FILES num2)

+  message("the list is ${num2}")

+  #list(GET CPU_TARGETS ${tgt} current_target)

+  list(GET GPU_MEX_FILES ${tgt} current_file_name)

+  get_filename_component(current_file ${current_file_name} NAME)

+  string(REGEX MATCH "(.+).c" match ${current_file})

+  if (NOT ${match} EQUAL "" )

+  set (current_target ${CMAKE_MATCH_1})

+  endif()

+  message("matlab_add_mex target " ${current_file} " and " ${current_target})

+        message("matlab_add_mex " ${current_target})

+        matlab_add_mex(

+          NAME ${current_target}

+          SRC 

+            ${current_file_name} 

+          LINK_TO cilregcuda ${Matlab_LIBRARIES}

+          )

+        

+        target_include_directories(${current_target}

+        PUBLIC ${CMAKE_SOURCE_DIR}/Core/regularisers_CPU

+               ${CMAKE_SOURCE_DIR}/Core/regularisers_GPU

+               ${CMAKE_SOURCE_DIR}/Core/inpainters_CPU

+               ${CMAKE_SOURCE_DIR}/Core/

+               ${MATLAB_INCLUDE_DIR})

+       

+        list(APPEND GPU_MEX_TARGETS ${current_target})

+        INSTALL(TARGETS ${current_target} DESTINATION "${MATLAB_DEST}")

+      endforeach()

+       

+      add_custom_target(MatlabWrapperGPU DEPENDS ${GPU_MEX_TARGETS})

+      

+    endif()

+endif()

diff --git a/src/Matlab/mex_compile/compileCPU_mex_Linux.m b/src/Matlab/mex_compile/compileCPU_mex_Linux.m
new file mode 100644
index 0000000..72a828e
--- /dev/null
+++ b/src/Matlab/mex_compile/compileCPU_mex_Linux.m
@@ -0,0 +1,81 @@
+% execute this mex file on Linux in Matlab once
+
+fsep = '/';
+
+pathcopyFrom = sprintf(['..' fsep '..' fsep '..' fsep 'Core' fsep 'regularisers_CPU'], 1i);
+pathcopyFrom1 = sprintf(['..' fsep '..' fsep '..' fsep 'Core' fsep 'CCPiDefines.h'], 1i);
+pathcopyFrom2 = sprintf(['..' fsep '..' fsep '..' fsep 'Core' fsep 'inpainters_CPU'], 1i);
+
+copyfile(pathcopyFrom, 'regularisers_CPU');
+copyfile(pathcopyFrom1, 'regularisers_CPU');
+copyfile(pathcopyFrom2, 'regularisers_CPU');
+
+cd regularisers_CPU
+
+Pathmove = sprintf(['..' fsep 'installed' fsep], 1i);
+
+fprintf('%s \n', '<<<<<<<<<<<Compiling CPU regularisers>>>>>>>>>>>>>');
+
+fprintf('%s \n', 'Compiling ROF-TV...');
+mex ROF_TV.c ROF_TV_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
+movefile('ROF_TV.mex*',Pathmove);
+
+fprintf('%s \n', 'Compiling FGP-TV...');
+mex FGP_TV.c FGP_TV_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
+movefile('FGP_TV.mex*',Pathmove);
+
+fprintf('%s \n', 'Compiling SB-TV...');
+mex SB_TV.c SB_TV_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
+movefile('SB_TV.mex*',Pathmove);
+
+fprintf('%s \n', 'Compiling dFGP-TV...');
+mex FGP_dTV.c FGP_dTV_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
+movefile('FGP_dTV.mex*',Pathmove);
+
+fprintf('%s \n', 'Compiling TNV...');
+mex TNV.c TNV_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
+movefile('TNV.mex*',Pathmove);
+
+fprintf('%s \n', 'Compiling NonLinear Diffusion...');
+mex NonlDiff.c Diffusion_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
+movefile('NonlDiff.mex*',Pathmove);
+
+fprintf('%s \n', 'Compiling Anisotropic diffusion of higher order...');
+mex Diffusion_4thO.c Diffus4th_order_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
+movefile('Diffusion_4thO.mex*',Pathmove);
+
+fprintf('%s \n', 'Compiling TGV...');
+mex TGV.c TGV_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
+movefile('TGV.mex*',Pathmove);
+
+fprintf('%s \n', 'Compiling ROF-LLT...');
+mex LLT_ROF.c LLT_ROF_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
+movefile('LLT_ROF.mex*',Pathmove);
+
+fprintf('%s \n', 'Compiling NonLocal-TV...');
+mex PatchSelect.c PatchSelect_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
+mex Nonlocal_TV.c Nonlocal_TV_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
+movefile('Nonlocal_TV.mex*',Pathmove);
+movefile('PatchSelect.mex*',Pathmove);
+
+fprintf('%s \n', 'Compiling additional tools...');
+mex TV_energy.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
+movefile('TV_energy.mex*',Pathmove);
+
+%############Inpainters##############%
+fprintf('%s \n', 'Compiling Nonlinear/Linear diffusion inpainting...');
+mex NonlDiff_Inp.c Diffusion_Inpaint_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
+movefile('NonlDiff_Inp.mex*',Pathmove);
+
+fprintf('%s \n', 'Compiling Nonlocal marching method for inpainting...');
+mex NonlocalMarching_Inpaint.c NonlocalMarching_Inpaint_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp"
+movefile('NonlocalMarching_Inpaint.mex*',Pathmove);
+
+delete SB_TV_core* ROF_TV_core* FGP_TV_core* FGP_dTV_core* TNV_core* utils* Diffusion_core* Diffus4th_order_core* TGV_core* LLT_ROF_core* CCPiDefines.h
+delete PatchSelect_core* Nonlocal_TV_core*
+delete Diffusion_Inpaint_core* NonlocalMarching_Inpaint_core*
+fprintf('%s \n', '<<<<<<< Regularisers successfully compiled! >>>>>>>');
+
+pathA2 = sprintf(['..' fsep '..' fsep], 1i);
+cd(pathA2);
+cd demos
diff --git a/src/Matlab/mex_compile/compileCPU_mex_WINDOWS.m b/src/Matlab/mex_compile/compileCPU_mex_WINDOWS.m
new file mode 100644
index 0000000..6f7541c
--- /dev/null
+++ b/src/Matlab/mex_compile/compileCPU_mex_WINDOWS.m
@@ -0,0 +1,135 @@
+% execute this mex file on Windows in Matlab once
+
+% >>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+% I've been able to compile on Windows 7 with MinGW and Matlab 2016b, however, 
+% not sure if openmp is enabled after the compilation. 
+
+% Here I present two ways how software can be compiled, if you have some
+% other suggestions/remarks please contact me at dkazanc@hotmail.com 
+% >>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+
+fsep = '/';
+
+pathcopyFrom = sprintf(['..' fsep '..' fsep '..' fsep 'Core' fsep 'regularisers_CPU'], 1i);
+pathcopyFrom1 = sprintf(['..' fsep '..' fsep '..' fsep 'Core' fsep 'CCPiDefines.h'], 1i);
+pathcopyFrom2 = sprintf(['..' fsep '..' fsep '..' fsep 'Core' fsep 'inpainters_CPU'], 1i);
+
+copyfile(pathcopyFrom, 'regularisers_CPU');
+copyfile(pathcopyFrom1, 'regularisers_CPU');
+copyfile(pathcopyFrom2, 'regularisers_CPU');
+
+cd regularisers_CPU
+
+Pathmove = sprintf(['..' fsep 'installed' fsep], 1i);
+
+fprintf('%s \n', '<<<<<<<<<<<Compiling CPU regularisers>>>>>>>>>>>>>');
+
+fprintf('%s \n', 'Compiling ROF-TV...');
+mex ROF_TV.c ROF_TV_core.c utils.c COMPFLAGS="\$COMPFLAGS -fopenmp -Wall -std=c99"
+movefile('ROF_TV.mex*',Pathmove);
+
+fprintf('%s \n', 'Compiling FGP-TV...');
+mex FGP_TV.c FGP_TV_core.c utils.c COMPFLAGS="\$COMPFLAGS -fopenmp -Wall -std=c99"
+movefile('FGP_TV.mex*',Pathmove);
+
+fprintf('%s \n', 'Compiling SB-TV...');
+mex SB_TV.c SB_TV_core.c utils.c COMPFLAGS="\$COMPFLAGS -fopenmp -Wall -std=c99"
+movefile('SB_TV.mex*',Pathmove);
+
+fprintf('%s \n', 'Compiling dFGP-TV...');
+mex FGP_dTV.c FGP_dTV_core.c utils.c COMPFLAGS="\$COMPFLAGS -fopenmp -Wall -std=c99"
+movefile('FGP_dTV.mex*',Pathmove);
+
+fprintf('%s \n', 'Compiling TNV...');
+mex TNV.c TNV_core.c utils.c COMPFLAGS="\$COMPFLAGS -fopenmp -Wall -std=c99"
+movefile('TNV.mex*',Pathmove);
+
+fprintf('%s \n', 'Compiling NonLinear Diffusion...');
+mex NonlDiff.c Diffusion_core.c utils.c COMPFLAGS="\$COMPFLAGS -fopenmp -Wall -std=c99"
+movefile('NonlDiff.mex*',Pathmove);
+
+fprintf('%s \n', 'Compiling Anisotropic diffusion of higher order...');
+mex Diffusion_4thO.c Diffus4th_order_core.c utils.c COMPFLAGS="\$COMPFLAGS -fopenmp -Wall -std=c99"
+movefile('Diffusion_4thO.mex*',Pathmove);
+
+fprintf('%s \n', 'Compiling TGV...');
+mex TGV.c TGV_core.c utils.c COMPFLAGS="\$COMPFLAGS -fopenmp -Wall -std=c99"
+movefile('TGV.mex*',Pathmove);
+
+fprintf('%s \n', 'Compiling ROF-LLT...');
+mex LLT_ROF.c LLT_ROF_core.c utils.c COMPFLAGS="\$COMPFLAGS -fopenmp -Wall -std=c99"
+movefile('LLT_ROF.mex*',Pathmove);
+
+fprintf('%s \n', 'Compiling NonLocal-TV...');
+mex PatchSelect.c PatchSelect_core.c utils.c COMPFLAGS="\$COMPFLAGS -fopenmp -Wall -std=c99"
+mex Nonlocal_TV.c Nonlocal_TV_core.c utils.c COMPFLAGS="\$COMPFLAGS -fopenmp -Wall -std=c99"
+movefile('Nonlocal_TV.mex*',Pathmove);
+movefile('PatchSelect.mex*',Pathmove);
+
+fprintf('%s \n', 'Compiling additional tools...');
+mex TV_energy.c utils.c COMPFLAGS="\$COMPFLAGS -fopenmp -Wall -std=c99"
+movefile('TV_energy.mex*',Pathmove);
+
+%############Inpainters##############%
+fprintf('%s \n', 'Compiling Nonlinear/Linear diffusion inpainting...');
+mex NonlDiff_Inp.c Diffusion_Inpaint_core.c utils.c COMPFLAGS="\$COMPFLAGS -fopenmp -Wall -std=c99"
+movefile('NonlDiff_Inp.mex*',Pathmove);
+
+fprintf('%s \n', 'Compiling Nonlocal marching method for inpaiting...');
+mex NonlocalMarching_Inpaint.c NonlocalMarching_Inpaint_core.c utils.c COMPFLAGS="\$COMPFLAGS -fopenmp -Wall -std=c99"
+movefile('NonlocalMarching_Inpaint.mex*',Pathmove);
+
+
+%%
+%%% The second approach to compile using TDM-GCC which follows this
+%%% discussion:
+%%% https://uk.mathworks.com/matlabcentral/answers/279171-using-mingw-compiler-and-open-mp#comment_359122
+%%% 1. Install TDM-GCC independently from http://tdm-gcc.tdragon.net/ (I installed 5.1.0)
+%%% Install openmp version: http://sourceforge.net/projects/tdm-gcc/files/TDM-GCC%205%20series/5.1.0-tdm64-1/gcc-5.1.0-tdm64-1-openmp.zip/download
+%%% 2. Link til libgomp.a in that installation when compilling your mex file.
+
+%%% assuming you unzipped TDM GCC (OpenMp) in folder TDMGCC on C drive, uncomment
+%%% bellow
+% fprintf('%s \n', 'Compiling CPU regularisers...');
+% mex C:\TDMGCC\lib\gcc\x86_64-w64-mingw32\5.1.0\libgomp.a CXXFLAGS="$CXXFLAGS -std=c++11 -fopenmp" ROF_TV.c ROF_TV_core.c utils.c
+% movefile('ROF_TV.mex*',Pathmove);
+% mex C:\TDMGCC\lib\gcc\x86_64-w64-mingw32\5.1.0\libgomp.a CXXFLAGS="$CXXFLAGS -std=c++11 -fopenmp" FGP_TV.c FGP_TV_core.c utils.c
+% movefile('FGP_TV.mex*',Pathmove);
+% mex C:\TDMGCC\lib\gcc\x86_64-w64-mingw32\5.1.0\libgomp.a CXXFLAGS="$CXXFLAGS -std=c++11 -fopenmp" SB_TV.c SB_TV_core.c utils.c
+% movefile('SB_TV.mex*',Pathmove);
+% mex C:\TDMGCC\lib\gcc\x86_64-w64-mingw32\5.1.0\libgomp.a CXXFLAGS="$CXXFLAGS -std=c++11 -fopenmp" FGP_dTV.c FGP_dTV_core.c utils.c
+% movefile('FGP_dTV.mex*',Pathmove);
+% mex C:\TDMGCC\lib\gcc\x86_64-w64-mingw32\5.1.0\libgomp.a CXXFLAGS="$CXXFLAGS -std=c++11 -fopenmp" TNV.c TNV_core.c utils.c
+% movefile('TNV.mex*',Pathmove);
+% mex C:\TDMGCC\lib\gcc\x86_64-w64-mingw32\5.1.0\libgomp.a CXXFLAGS="$CXXFLAGS -std=c++11 -fopenmp" NonlDiff.c Diffusion_core.c utils.c
+% movefile('NonlDiff.mex*',Pathmove);
+% mex C:\TDMGCC\lib\gcc\x86_64-w64-mingw32\5.1.0\libgomp.a CXXFLAGS="$CXXFLAGS -std=c++11 -fopenmp" Diffusion_4thO.c Diffus4th_order_core.c utils.c
+% movefile('Diffusion_4thO.mex*',Pathmove);
+% mex C:\TDMGCC\lib\gcc\x86_64-w64-mingw32\5.1.0\libgomp.a CXXFLAGS="$CXXFLAGS -std=c++11 -fopenmp" TGV.c TGV_core.c utils.c
+% movefile('TGV.mex*',Pathmove);
+% mex C:\TDMGCC\lib\gcc\x86_64-w64-mingw32\5.1.0\libgomp.a CXXFLAGS="$CXXFLAGS -std=c++11 -fopenmp" LLT_ROF.c LLT_ROF_core.c utils.c
+% movefile('LLT_ROF.mex*',Pathmove);
+% mex C:\TDMGCC\lib\gcc\x86_64-w64-mingw32\5.1.0\libgomp.a CXXFLAGS="$CXXFLAGS -std=c++11 -fopenmp" PatchSelect.c PatchSelect_core.c utils.c
+% mex C:\TDMGCC\lib\gcc\x86_64-w64-mingw32\5.1.0\libgomp.a CXXFLAGS="$CXXFLAGS -std=c++11 -fopenmp" Nonlocal_TV.c Nonlocal_TV_core.c utils.c
+% movefile('Nonlocal_TV.mex*',Pathmove);
+% movefile('PatchSelect.mex*',Pathmove);
+% mex C:\TDMGCC\lib\gcc\x86_64-w64-mingw32\5.1.0\libgomp.a CXXFLAGS="$CXXFLAGS -std=c++11 -fopenmp" TV_energy.c utils.c
+% movefile('TV_energy.mex*',Pathmove);
+% mex C:\TDMGCC\lib\gcc\x86_64-w64-mingw32\5.1.0\libgomp.a CXXFLAGS="$CXXFLAGS -std=c++11 -fopenmp" NonlDiff_Inp.c Diffusion_Inpaint_core.c utils.c
+% movefile('NonlDiff_Inp.mex*',Pathmove);
+% mex C:\TDMGCC\lib\gcc\x86_64-w64-mingw32\5.1.0\libgomp.a CXXFLAGS="$CXXFLAGS -std=c++11 -fopenmp" NonlocalMarching_Inpaint.c NonlocalMarching_Inpaint_core.c utils.c
+% movefile('NonlocalMarching_Inpaint.mex*',Pathmove);
+
+
+delete SB_TV_core* ROF_TV_core* FGP_TV_core* FGP_dTV_core* TNV_core* utils* Diffusion_core* Diffus4th_order_core* TGV_core* CCPiDefines.h
+delete PatchSelect_core* Nonlocal_TV_core*
+delete Diffusion_Inpaint_core* NonlocalMarching_Inpaint_core*
+fprintf('%s \n', 'Regularisers successfully compiled!');
+
+
+%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%pathA2 = sprintf(['..' fsep '..' fsep], 1i);
+%cd(pathA2);
+%cd demos
diff --git a/src/Matlab/mex_compile/compileGPU_mex.m b/src/Matlab/mex_compile/compileGPU_mex.m
new file mode 100644
index 0000000..dd1475c
--- /dev/null
+++ b/src/Matlab/mex_compile/compileGPU_mex.m
@@ -0,0 +1,74 @@
+% execute this mex file in Matlab once
+
+%>>>>>>>>>>>>>>>>>Important<<<<<<<<<<<<<<<<<<<
+% In order to compile CUDA modules one needs to have nvcc-compiler
+% installed (see CUDA SDK), check it under MATLAB with !nvcc --version
+
+% In the code bellow we provide a full explicit path to nvcc compiler 
+% ! paths to matlab and CUDA sdk can be different, modify accordingly !
+
+% Tested on Ubuntu 18.04/MATLAB 2016b/cuda10.0/gcc7.3
+
+% Installation HAS NOT been tested on Windows, please you Cmake build or
+% modify the code bellow accordingly
+fsep = '/';
+
+pathcopyFrom = sprintf(['..' fsep '..' fsep '..' fsep 'Core' fsep 'regularisers_GPU'], 1i);
+pathcopyFrom1 = sprintf(['..' fsep '..' fsep '..' fsep 'Core' fsep 'CCPiDefines.h'], 1i);
+
+copyfile(pathcopyFrom, 'regularisers_GPU');
+copyfile(pathcopyFrom1, 'regularisers_GPU');
+
+cd regularisers_GPU
+
+Pathmove = sprintf(['..' fsep 'installed' fsep], 1i);
+
+fprintf('%s \n', '<<<<<<<<<<<Compiling GPU regularisers (CUDA)>>>>>>>>>>>>>');
+
+fprintf('%s \n', 'Compiling ROF-TV...');
+!/usr/local/cuda/bin/nvcc -O0 -c TV_ROF_GPU_core.cu -Xcompiler -fPIC -I~/SOFT/MATLAB9/extern/include/
+mex -g -I/usr/local/cuda-10.0/include -L/usr/local/cuda-10.0/lib64 -lcudart -lcufft -lmwgpu ROF_TV_GPU.cpp TV_ROF_GPU_core.o
+movefile('ROF_TV_GPU.mex*',Pathmove);
+
+fprintf('%s \n', 'Compiling FGP-TV...');
+!/usr/local/cuda/bin/nvcc -O0 -c TV_FGP_GPU_core.cu -Xcompiler -fPIC -I~/SOFT/MATLAB9/extern/include/
+mex -g -I/usr/local/cuda-10.0/include -L/usr/local/cuda-10.0/lib64 -lcudart -lcufft -lmwgpu FGP_TV_GPU.cpp TV_FGP_GPU_core.o
+movefile('FGP_TV_GPU.mex*',Pathmove);
+
+fprintf('%s \n', 'Compiling SB-TV...');
+!/usr/local/cuda/bin/nvcc -O0 -c TV_SB_GPU_core.cu -Xcompiler -fPIC -I~/SOFT/MATLAB9/extern/include/
+mex -g -I/usr/local/cuda-10.0/include -L/usr/local/cuda-10.0/lib64 -lcudart -lcufft -lmwgpu SB_TV_GPU.cpp TV_SB_GPU_core.o
+movefile('SB_TV_GPU.mex*',Pathmove);
+
+fprintf('%s \n', 'Compiling TGV...');
+!/usr/local/cuda/bin/nvcc -O0 -c TGV_GPU_core.cu -Xcompiler -fPIC -I~/SOFT/MATLAB9/extern/include/
+mex -g -I/usr/local/cuda-10.0/include -L/usr/local/cuda-10.0/lib64 -lcudart -lcufft -lmwgpu TGV_GPU.cpp TGV_GPU_core.o
+movefile('TGV_GPU.mex*',Pathmove);
+
+fprintf('%s \n', 'Compiling dFGP-TV...');
+!/usr/local/cuda/bin/nvcc -O0 -c dTV_FGP_GPU_core.cu -Xcompiler -fPIC -I~/SOFT/MATLAB9/extern/include/
+mex -g -I/usr/local/cuda-10.0/include -L/usr/local/cuda-10.0/lib64 -lcudart -lcufft -lmwgpu FGP_dTV_GPU.cpp dTV_FGP_GPU_core.o
+movefile('FGP_dTV_GPU.mex*',Pathmove);
+
+fprintf('%s \n', 'Compiling NonLinear Diffusion...');
+!/usr/local/cuda/bin/nvcc -O0 -c NonlDiff_GPU_core.cu -Xcompiler -fPIC -I~/SOFT/MATLAB9/extern/include/
+mex -g -I/usr/local/cuda-10.0/include -L/usr/local/cuda-10.0/lib64 -lcudart -lcufft -lmwgpu NonlDiff_GPU.cpp NonlDiff_GPU_core.o
+movefile('NonlDiff_GPU.mex*',Pathmove);
+
+fprintf('%s \n', 'Compiling Anisotropic diffusion of higher order...');
+!/usr/local/cuda/bin/nvcc -O0 -c Diffus_4thO_GPU_core.cu -Xcompiler -fPIC -I~/SOFT/MATLAB9/extern/include/
+mex -g -I/usr/local/cuda-10.0/include -L/usr/local/cuda-10.0/lib64 -lcudart -lcufft -lmwgpu Diffusion_4thO_GPU.cpp Diffus_4thO_GPU_core.o
+movefile('Diffusion_4thO_GPU.mex*',Pathmove);
+
+fprintf('%s \n', 'Compiling ROF-LLT...');
+!/usr/local/cuda/bin/nvcc -O0 -c LLT_ROF_GPU_core.cu -Xcompiler -fPIC -I~/SOFT/MATLAB9/extern/include/
+mex -g -I/usr/local/cuda-10.0/include -L/usr/local/cuda-10.0/lib64 -lcudart -lcufft -lmwgpu LLT_ROF_GPU.cpp LLT_ROF_GPU_core.o
+movefile('LLT_ROF_GPU.mex*',Pathmove);
+
+
+delete TV_ROF_GPU_core* TV_FGP_GPU_core* TV_SB_GPU_core* dTV_FGP_GPU_core* NonlDiff_GPU_core* Diffus_4thO_GPU_core* TGV_GPU_core* LLT_ROF_GPU_core* CCPiDefines.h
+fprintf('%s \n', 'All successfully compiled!');
+
+pathA2 = sprintf(['..' fsep '..' fsep], 1i);
+cd(pathA2);
+cd demos
\ No newline at end of file
diff --git a/src/Matlab/mex_compile/installed/MEXed_files_location.txt b/src/Matlab/mex_compile/installed/MEXed_files_location.txt
new file mode 100644
index 0000000..e69de29
--- /dev/null
+++ b/src/Matlab/mex_compile/installed/MEXed_files_location.txt
diff --git a/src/Matlab/mex_compile/regularisers_CPU/Diffusion_4thO.c b/src/Matlab/mex_compile/regularisers_CPU/Diffusion_4thO.c
new file mode 100644
index 0000000..66ea9be
--- /dev/null
+++ b/src/Matlab/mex_compile/regularisers_CPU/Diffusion_4thO.c
@@ -0,0 +1,77 @@
+/*
+ * This work is part of the Core Imaging Library developed by
+ * Visual Analytics and Imaging System Group of the Science Technology
+ * Facilities Council, STFC
+ *
+ * Copyright 2017 Daniil Kazantsev
+ * Copyright 2017 Srikanth Nagella, Edoardo Pasca
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ * http://www.apache.org/licenses/LICENSE-2.0
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+#include "matrix.h"
+#include "mex.h"
+#include "Diffus4th_order_core.h"
+
+/* C-OMP implementation of fourth-order diffusion scheme [1] for piecewise-smooth recovery (2D/3D case)
+ * The minimisation is performed using explicit scheme. 
+ *
+ * Input Parameters:
+ * 1. Noisy image/volume [REQUIRED]
+ * 2. lambda - regularization parameter [REQUIRED]
+ * 3. Edge-preserving parameter (sigma) [REQUIRED]
+ * 4. Number of iterations, for explicit scheme >= 150 is recommended [OPTIONAL, default 300]
+ * 5. tau - time-marching step for the explicit scheme [OPTIONAL, default 0.015]
+ *
+ * Output:
+ * [1] Regularized image/volume 
+ *
+ * This function is based on the paper by
+ * [1] Hajiaboli, M.R., 2011. An anisotropic fourth-order diffusion filter for image noise removal. International Journal of Computer Vision, 92(2), pp.177-191.
+ */
+
+void mexFunction(
+        int nlhs, mxArray *plhs[],
+        int nrhs, const mxArray *prhs[])
+        
+{
+    int number_of_dims, iter_numb;
+    mwSize dimX, dimY, dimZ;
+    const mwSize *dim_array;
+    float *Input, *Output=NULL, lambda, tau, sigma;
+    
+    dim_array = mxGetDimensions(prhs[0]);
+    number_of_dims = mxGetNumberOfDimensions(prhs[0]);
+    
+    /*Handling Matlab input data*/
+    Input  = (float *) mxGetData(prhs[0]);
+    lambda =  (float) mxGetScalar(prhs[1]); /* regularization parameter */
+    sigma = (float) mxGetScalar(prhs[2]); /* Edge-preserving parameter */
+    iter_numb = 300; /* iterations number */
+    tau = 0.01; /* marching step parameter */
+    
+    if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); }
+    if ((nrhs < 3) || (nrhs > 5)) mexErrMsgTxt("At least 3 parameters is required, all parameters are: Image(2D/3D), Regularisation parameter, Edge-preserving parameter, iterations number, time-marching constant");
+    if ((nrhs == 4) || (nrhs == 5))  iter_numb = (int) mxGetScalar(prhs[3]); /* iterations number */
+    if (nrhs == 5)  tau =  (float) mxGetScalar(prhs[4]); /* marching step parameter */
+    
+    /*Handling Matlab output data*/
+    dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2];
+    
+    /* output arrays*/
+    if (number_of_dims == 2) {
+        dimZ = 1; /*2D case*/
+        /* output image/volume */
+        Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
+    }
+    if (number_of_dims == 3) Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
+    
+    Diffus4th_CPU_main(Input, Output, lambda, sigma, iter_numb, tau, dimX, dimY, dimZ);
+}
\ No newline at end of file
diff --git a/src/Matlab/mex_compile/regularisers_CPU/FGP_TV.c b/src/Matlab/mex_compile/regularisers_CPU/FGP_TV.c
new file mode 100644
index 0000000..642362f
--- /dev/null
+++ b/src/Matlab/mex_compile/regularisers_CPU/FGP_TV.c
@@ -0,0 +1,97 @@
+/*
+ * This work is part of the Core Imaging Library developed by
+ * Visual Analytics and Imaging System Group of the Science Technology
+ * Facilities Council, STFC
+ *
+ * Copyright 2017 Daniil Kazantsev
+ * Copyright 2017 Srikanth Nagella, Edoardo Pasca
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ * http://www.apache.org/licenses/LICENSE-2.0
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+#include "matrix.h"
+#include "mex.h"
+#include "FGP_TV_core.h"
+
+/* C-OMP implementation of FGP-TV [1] denoising/regularization model (2D/3D case)
+ *
+ * Input Parameters:
+ * 1. Noisy image/volume
+ * 2. lambdaPar - regularization parameter
+ * 3. Number of iterations
+ * 4. eplsilon: tolerance constant
+ * 5. TV-type: methodTV - 'iso' (0) or 'l1' (1)
+ * 6. nonneg: 'nonnegativity (0 is OFF by default)
+ * 7. print information: 0 (off) or 1 (on)
+ *
+ * Output:
+ * [1] Filtered/regularized image
+ *
+ * This function is based on the Matlab's code and paper by
+ * [1] Amir Beck and Marc Teboulle, "Fast Gradient-Based Algorithms for Constrained Total Variation Image Denoising and Deblurring Problems"
+ */
+
+
+void mexFunction(
+        int nlhs, mxArray *plhs[],
+        int nrhs, const mxArray *prhs[])
+        
+{
+    int number_of_dims, iter, methTV, printswitch, nonneg;
+    mwSize dimX, dimY, dimZ;
+    const mwSize *dim_array;
+    float *Input, *Output=NULL, lambda, epsil;
+    
+    number_of_dims = mxGetNumberOfDimensions(prhs[0]);
+    dim_array = mxGetDimensions(prhs[0]);
+    
+    /*Handling Matlab input data*/
+    if ((nrhs < 2) || (nrhs > 7)) mexErrMsgTxt("At least 2 parameters is required, all parameters are: Image(2D/3D), Regularization parameter, Regularization parameter, iterations number, tolerance, penalty type ('iso' or 'l1'), nonnegativity switch, print switch");
+    
+    Input  = (float *) mxGetData(prhs[0]); /*noisy image (2D/3D) */
+    lambda =  (float) mxGetScalar(prhs[1]); /* regularization parameter */
+    iter = 300; /* default iterations number */
+    epsil = 0.0001; /* default tolerance constant */
+    methTV = 0;  /* default isotropic TV penalty */
+    nonneg = 0; /* default nonnegativity switch, off - 0 */
+    printswitch = 0; /*default print is switched, off - 0 */
+    
+    if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); }
+    
+    if ((nrhs == 3) || (nrhs == 4) || (nrhs == 5) || (nrhs == 6) || (nrhs == 7))  iter = (int) mxGetScalar(prhs[2]); /* iterations number */
+    if ((nrhs == 4) || (nrhs == 5) || (nrhs == 6) || (nrhs == 7))  epsil =  (float) mxGetScalar(prhs[3]); /* tolerance constant */
+    if ((nrhs == 5) || (nrhs == 6) || (nrhs == 7))  {
+        char *penalty_type;
+        penalty_type = mxArrayToString(prhs[4]); /* choosing TV penalty: 'iso' or 'l1', 'iso' is the default */
+        if ((strcmp(penalty_type, "l1") != 0) && (strcmp(penalty_type, "iso") != 0)) mexErrMsgTxt("Choose TV type: 'iso' or 'l1',");
+        if (strcmp(penalty_type, "l1") == 0)  methTV = 1;  /* enable 'l1' penalty */
+        mxFree(penalty_type);
+    }
+    if ((nrhs == 6) || (nrhs == 7))  {
+        nonneg = (int) mxGetScalar(prhs[5]);
+        if ((nonneg != 0) && (nonneg != 1)) mexErrMsgTxt("Nonnegativity constraint can be enabled by choosing 1 or off - 0");
+    }
+    if (nrhs == 7)  {
+        printswitch = (int) mxGetScalar(prhs[6]);
+        if ((printswitch != 0) && (printswitch != 1)) mexErrMsgTxt("Print can be enabled by choosing 1 or off - 0");
+    }
+    
+    /*Handling Matlab output data*/
+    dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2];
+    
+    if (number_of_dims == 2) {
+        dimZ = 1; /*2D case*/
+        Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
+    }
+    if (number_of_dims == 3) Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
+    
+    /* running the function */
+    TV_FGP_CPU_main(Input, Output, lambda, iter, epsil, methTV, nonneg, printswitch, dimX, dimY, dimZ);
+}
\ No newline at end of file
diff --git a/src/Matlab/mex_compile/regularisers_CPU/FGP_dTV.c b/src/Matlab/mex_compile/regularisers_CPU/FGP_dTV.c
new file mode 100644
index 0000000..1a0c070
--- /dev/null
+++ b/src/Matlab/mex_compile/regularisers_CPU/FGP_dTV.c
@@ -0,0 +1,114 @@
+/*
+ * This work is part of the Core Imaging Library developed by
+ * Visual Analytics and Imaging System Group of the Science Technology
+ * Facilities Council, STFC
+ *
+ * Copyright 2017 Daniil Kazantsev
+ * Copyright 2017 Srikanth Nagella, Edoardo Pasca
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ * http://www.apache.org/licenses/LICENSE-2.0
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+#include "matrix.h"
+#include "mex.h"
+#include "FGP_dTV_core.h"
+
+/* C-OMP implementation of FGP-dTV [1,2] denoising/regularization model (2D/3D case)
+ * which employs structural similarity of the level sets of two images/volumes, see [1,2]
+ * The current implementation updates image 1 while image 2 is being fixed.
+ *
+ * Input Parameters:
+ * 1. Noisy image/volume [REQUIRED]
+ * 2. Additional reference image/volume of the same dimensions as (1) [REQUIRED]
+ * 3. lambdaPar - regularization parameter [REQUIRED]
+ * 4. Number of iterations [OPTIONAL]
+ * 5. eplsilon: tolerance constant [OPTIONAL]
+ * 6. eta: smoothing constant to calculate gradient of the reference [OPTIONAL] * 
+ * 7. TV-type: methodTV - 'iso' (0) or 'l1' (1) [OPTIONAL]
+ * 8. nonneg: 'nonnegativity (0 is OFF by default) [OPTIONAL]
+ * 9. print information: 0 (off) or 1 (on) [OPTIONAL]
+ *
+ * Output:
+ * [1] Filtered/regularized image/volume
+ *
+ * This function is based on the Matlab's codes and papers by
+ * [1] Amir Beck and Marc Teboulle, "Fast Gradient-Based Algorithms for Constrained Total Variation Image Denoising and Deblurring Problems"
+ * [2] M. J. Ehrhardt and M. M. Betcke, Multi-Contrast MRI Reconstruction with Structure-Guided Total Variation, SIAM Journal on Imaging Sciences 9(3), pp. 1084–1106
+ */
+
+
+void mexFunction(
+        int nlhs, mxArray *plhs[],
+        int nrhs, const mxArray *prhs[])
+        
+{
+    int number_of_dims, iter, methTV, printswitch, nonneg;
+    mwSize dimX, dimY, dimZ;
+    const mwSize *dim_array;
+    const mwSize *dim_array2;    
+    float *Input, *InputRef, *Output=NULL, lambda, epsil, eta;
+    
+    number_of_dims = mxGetNumberOfDimensions(prhs[0]);
+    dim_array = mxGetDimensions(prhs[0]);
+    dim_array2 = mxGetDimensions(prhs[1]);
+    
+    /*Handling Matlab input data*/
+    if ((nrhs < 3) || (nrhs > 9)) mexErrMsgTxt("At least 3 parameters is required, all parameters are: Image(2D/3D), Reference(2D/3D), Regularization parameter, iterations number, tolerance, smoothing constant, penalty type ('iso' or 'l1'), nonnegativity switch, print switch");
+    
+    Input  = (float *) mxGetData(prhs[0]); /*noisy image (2D/3D) */
+    InputRef  = (float *) mxGetData(prhs[1]); /* reference image (2D/3D) */
+    lambda =  (float) mxGetScalar(prhs[2]); /* regularization parameter */
+    iter = 300; /* default iterations number */
+    epsil = 0.0001; /* default tolerance constant */
+    eta = 0.01; /* default smoothing constant */
+    methTV = 0;  /* default isotropic TV penalty */
+    nonneg = 0; /* default nonnegativity switch, off - 0 */
+    printswitch = 0; /*default print is switched, off - 0 */
+    
+        
+    if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); }
+    if (mxGetClassID(prhs[1]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); }
+    
+    /*Handling Matlab output data*/
+    dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2];
+    if (number_of_dims == 2) { if ((dimX != dim_array2[0]) || (dimY != dim_array2[1])) mexErrMsgTxt("The input images have different dimensionalities");}
+    if (number_of_dims == 3) { if ((dimX != dim_array2[0]) || (dimY != dim_array2[1]) || (dimZ != dim_array2[2])) mexErrMsgTxt("The input volumes have different dimensionalities");}   
+    
+    
+    if ((nrhs == 4) || (nrhs == 5) || (nrhs == 6) || (nrhs == 7) || (nrhs == 8) || (nrhs == 9))  iter = (int) mxGetScalar(prhs[3]); /* iterations number */
+    if ((nrhs == 5) || (nrhs == 6) || (nrhs == 7) || (nrhs == 8) || (nrhs == 9))  epsil =  (float) mxGetScalar(prhs[4]); /* tolerance constant */
+    if ((nrhs == 6) || (nrhs == 7) || (nrhs == 8) || (nrhs == 9))  {
+    eta =  (float) mxGetScalar(prhs[5]); /* smoothing constant for the gradient of InputRef */
+    }
+    if ((nrhs == 7) || (nrhs == 8) || (nrhs == 9))  {        
+        char *penalty_type;
+        penalty_type = mxArrayToString(prhs[6]); /* choosing TV penalty: 'iso' or 'l1', 'iso' is the default */
+        if ((strcmp(penalty_type, "l1") != 0) && (strcmp(penalty_type, "iso") != 0)) mexErrMsgTxt("Choose TV type: 'iso' or 'l1',");
+        if (strcmp(penalty_type, "l1") == 0)  methTV = 1;  /* enable 'l1' penalty */
+        mxFree(penalty_type);
+    }    
+    if ((nrhs == 8) || (nrhs == 9))  {
+        nonneg = (int) mxGetScalar(prhs[7]);
+        if ((nonneg != 0) && (nonneg != 1)) mexErrMsgTxt("Nonnegativity constraint can be enabled by choosing 1 or off - 0");
+    }
+    if (nrhs == 9)  {
+        printswitch = (int) mxGetScalar(prhs[8]);
+        if ((printswitch != 0) && (printswitch != 1)) mexErrMsgTxt("Print can be enabled by choosing 1 or off - 0");
+    }    
+   
+    if (number_of_dims == 2) {
+        dimZ = 1; /*2D case*/
+        Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
+    }
+    if (number_of_dims == 3) Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
+    
+    /* running the function */
+    dTV_FGP_CPU_main(Input, InputRef, Output, lambda, iter, epsil, eta, methTV, nonneg, printswitch, dimX, dimY, dimZ);
+}
\ No newline at end of file
diff --git a/src/Matlab/mex_compile/regularisers_CPU/LLT_ROF.c b/src/Matlab/mex_compile/regularisers_CPU/LLT_ROF.c
new file mode 100644
index 0000000..ab45446
--- /dev/null
+++ b/src/Matlab/mex_compile/regularisers_CPU/LLT_ROF.c
@@ -0,0 +1,82 @@
+/*
+ * This work is part of the Core Imaging Library developed by
+ * Visual Analytics and Imaging System Group of the Science Technology
+ * Facilities Council, STFC
+ *
+ * Copyright 2017 Daniil Kazantsev
+ * Copyright 2017 Srikanth Nagella, Edoardo Pasca
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ * http://www.apache.org/licenses/LICENSE-2.0
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+#include "matrix.h"
+#include "mex.h"
+#include "LLT_ROF_core.h"
+
+/* C-OMP implementation of Lysaker, Lundervold and Tai (LLT) model [1] combined with Rudin-Osher-Fatemi [2] TV regularisation penalty.
+* 
+* This penalty can deliver visually pleasant piecewise-smooth recovery if regularisation parameters are selected well. 
+* The rule of thumb for selection is to start with lambdaLLT = 0 (just the ROF-TV model) and then proceed to increase 
+* lambdaLLT starting with smaller values. 
+*
+* Input Parameters:
+* 1. U0 - original noise image/volume
+* 2. lambdaROF - ROF-related regularisation parameter
+* 3. lambdaLLT - LLT-related regularisation parameter
+* 4. tau - time-marching step 
+* 5. iter - iterations number (for both models)
+*
+* Output:
+* Filtered/regularised image
+*
+* References: 
+* [1] Lysaker, M., Lundervold, A. and Tai, X.C., 2003. Noise removal using fourth-order partial differential equation with applications to medical magnetic resonance images in space and time. IEEE Transactions on image processing, 12(12), pp.1579-1590.
+* [2] Rudin, Osher, Fatemi, "Nonlinear Total Variation based noise removal algorithms"
+*/
+
+void mexFunction(
+        int nlhs, mxArray *plhs[],
+        int nrhs, const mxArray *prhs[])
+        
+{
+    int number_of_dims, iterationsNumb;
+    mwSize dimX, dimY, dimZ;
+    const mwSize *dim_array;    
+    float *Input, *Output=NULL, lambdaROF, lambdaLLT, tau;
+    
+    dim_array = mxGetDimensions(prhs[0]);
+    number_of_dims = mxGetNumberOfDimensions(prhs[0]);
+    
+    if ((nrhs < 3) || (nrhs > 5)) mexErrMsgTxt("At least 3 parameters is required, all parameters are: Image(2D/3D), Regularisation parameter (ROF), Regularisation parameter (LTT), iterations number, time-marching parameter");
+    
+    /*Handling Matlab input data*/
+    Input  = (float *) mxGetData(prhs[0]);
+    lambdaROF =  (float) mxGetScalar(prhs[1]); /* ROF regularization parameter */
+    lambdaLLT =  (float) mxGetScalar(prhs[2]); /* ROF regularization parameter */    
+    iterationsNumb = 250;
+    tau =  0.0025;
+    
+    if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); }   
+    if ((nrhs == 4) || (nrhs == 5)) iterationsNumb =  (int) mxGetScalar(prhs[3]); /* iterations number */    
+    if (nrhs == 5) tau =  (float) mxGetScalar(prhs[4]); /* marching step parameter */  
+        
+    /*Handling Matlab output data*/
+    dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2];
+    
+    /* output arrays*/
+    if (number_of_dims == 2) {
+        dimZ = 1; /*2D case*/
+        /* output image/volume */
+        Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));                        
+    }    
+    if (number_of_dims == 3) Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));   
+  
+    LLT_ROF_CPU_main(Input, Output, lambdaROF, lambdaLLT, iterationsNumb, tau, dimX, dimY, dimZ);    
+}
\ No newline at end of file
diff --git a/src/Matlab/mex_compile/regularisers_CPU/NonlDiff.c b/src/Matlab/mex_compile/regularisers_CPU/NonlDiff.c
new file mode 100644
index 0000000..ec35b8b
--- /dev/null
+++ b/src/Matlab/mex_compile/regularisers_CPU/NonlDiff.c
@@ -0,0 +1,89 @@
+/*
+ * This work is part of the Core Imaging Library developed by
+ * Visual Analytics and Imaging System Group of the Science Technology
+ * Facilities Council, STFC
+ *
+ * Copyright 2017 Daniil Kazantsev
+ * Copyright 2017 Srikanth Nagella, Edoardo Pasca
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ * http://www.apache.org/licenses/LICENSE-2.0
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+#include "matrix.h"
+#include "mex.h"
+#include "Diffusion_core.h"
+
+/* C-OMP implementation of linear and nonlinear diffusion with the regularisation model [1] (2D/3D case)
+ * The minimisation is performed using explicit scheme.
+ *
+ * Input Parameters:
+ * 1. Noisy image/volume
+ * 2. lambda - regularization parameter
+ * 3. Edge-preserving parameter (sigma), when sigma equals to zero nonlinear diffusion -> linear diffusion
+ * 4. Number of iterations, for explicit scheme >= 150 is recommended  [OPTIONAL parameter]
+ * 5. tau - time-marching step for explicit scheme [OPTIONAL parameter]
+ * 6. Penalty type: 1 - Huber, 2 - Perona-Malik, 3 - Tukey Biweight [OPTIONAL parameter]
+ *
+ * Output:
+ * [1] Regularized image/volume
+ *
+ * This function is based on the paper by
+ * [1] Perona, P. and Malik, J., 1990. Scale-space and edge detection using anisotropic diffusion. IEEE Transactions on pattern analysis and machine intelligence, 12(7), pp.629-639.
+ */
+
+void mexFunction(
+        int nlhs, mxArray *plhs[],
+        int nrhs, const mxArray *prhs[])
+        
+{
+    int number_of_dims, iter_numb, penaltytype;
+    mwSize dimX, dimY, dimZ;
+    const mwSize *dim_array;   
+    
+    float *Input, *Output=NULL, lambda, tau, sigma;
+    
+    dim_array = mxGetDimensions(prhs[0]);
+    number_of_dims = mxGetNumberOfDimensions(prhs[0]);
+    
+    /*Handling Matlab input data*/
+    Input  = (float *) mxGetData(prhs[0]);
+    lambda =  (float) mxGetScalar(prhs[1]); /* regularization parameter */
+    sigma = (float) mxGetScalar(prhs[2]); /* Edge-preserving parameter */
+    iter_numb = 300; /* iterations number */
+    tau = 0.025; /* marching step parameter */
+    penaltytype = 1; /* Huber penalty by default */
+    
+    if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); }
+    if ((nrhs < 3) || (nrhs > 6)) mexErrMsgTxt("At least 3 parameters is required, all parameters are: Image(2D/3D), Regularisation parameter, Edge-preserving parameter, iterations number, time-marching constant, penalty type - Huber, PM or Tukey");
+    if ((nrhs == 4) || (nrhs == 5) || (nrhs == 6))  iter_numb = (int) mxGetScalar(prhs[3]); /* iterations number */
+    if ((nrhs == 5) || (nrhs == 6))  tau =  (float) mxGetScalar(prhs[4]); /* marching step parameter */
+    if (nrhs == 6)  {
+        char *penalty_type;
+        penalty_type = mxArrayToString(prhs[5]); /* Huber, PM or Tukey 'Huber' is the default */
+        if ((strcmp(penalty_type, "Huber") != 0) && (strcmp(penalty_type, "PM") != 0) && (strcmp(penalty_type, "Tukey") != 0)) mexErrMsgTxt("Choose penalty: 'Huber', 'PM' or 'Tukey',");
+        if (strcmp(penalty_type, "Huber") == 0)  penaltytype = 1;  /* enable 'Huber' penalty */
+        if (strcmp(penalty_type, "PM") == 0)  penaltytype = 2;  /* enable Perona-Malik penalty */
+        if (strcmp(penalty_type, "Tukey") == 0)  penaltytype = 3;  /* enable Tikey Biweight penalty */
+        mxFree(penalty_type);
+    }    
+    
+    /*Handling Matlab output data*/
+    dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2];
+    
+    /* output arrays*/
+    if (number_of_dims == 2) {
+        dimZ = 1; /*2D case*/
+        /* output image/volume */
+        Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
+    }
+    if (number_of_dims == 3) Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
+    
+    Diffusion_CPU_main(Input, Output, lambda, sigma, iter_numb, tau, penaltytype, dimX, dimY, dimZ);
+}
\ No newline at end of file
diff --git a/src/Matlab/mex_compile/regularisers_CPU/NonlDiff_Inp.c b/src/Matlab/mex_compile/regularisers_CPU/NonlDiff_Inp.c
new file mode 100644
index 0000000..9833392
--- /dev/null
+++ b/src/Matlab/mex_compile/regularisers_CPU/NonlDiff_Inp.c
@@ -0,0 +1,103 @@
+/*
+ * This work is part of the Core Imaging Library developed by
+ * Visual Analytics and Imaging System Group of the Science Technology
+ * Facilities Council, STFC
+ *
+ * Copyright 2017 Daniil Kazantsev
+ * Copyright 2017 Srikanth Nagella, Edoardo Pasca
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ * http://www.apache.org/licenses/LICENSE-2.0
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+#include "matrix.h"
+#include "mex.h"
+#include "Diffusion_Inpaint_core.h"
+
+/* C-OMP implementation of linear and nonlinear diffusion [1,2] for inpainting task (2D/3D case)
+ * The minimisation is performed using explicit scheme. 
+ *
+ * Input Parameters:
+ * 1. Image/volume to inpaint
+ * 2. Inpainting Mask of the same size as (1) in 'unsigned char' format  (ones mark the region to inpaint, zeros belong to the data)
+ * 3. lambda - regularization parameter
+ * 4. Edge-preserving parameter (sigma), when sigma equals to zero nonlinear diffusion -> linear diffusion
+ * 5. Number of iterations, for explicit scheme >= 150 is recommended 
+ * 6. tau - time-marching step for explicit scheme
+ * 7. Penalty type: 1 - Huber, 2 - Perona-Malik, 3 - Tukey Biweight
+ *
+ * Output:
+ * [1] Inpainted image/volume 
+ *
+ * This function is based on the paper by
+ * [1] Perona, P. and Malik, J., 1990. Scale-space and edge detection using anisotropic diffusion. IEEE Transactions on pattern analysis and machine intelligence, 12(7), pp.629-639.
+ * [2] Black, M.J., Sapiro, G., Marimont, D.H. and Heeger, D., 1998. Robust anisotropic diffusion. IEEE Transactions on image processing, 7(3), pp.421-432.
+ */
+
+void mexFunction(
+        int nlhs, mxArray *plhs[],
+        int nrhs, const mxArray *prhs[])
+        
+{
+    int number_of_dims, iter_numb, penaltytype, i, inpaint_elements;
+    mwSize dimX, dimY, dimZ;
+    const mwSize *dim_array;   
+    const mwSize *dim_array2;   
+    
+    float *Input, *Output=NULL, lambda, tau, sigma;
+    unsigned char *Mask;
+    
+    dim_array = mxGetDimensions(prhs[0]);
+    dim_array2 = mxGetDimensions(prhs[1]);
+    number_of_dims = mxGetNumberOfDimensions(prhs[0]);
+    
+    /*Handling Matlab input data*/
+    Input  = (float *) mxGetData(prhs[0]);
+    Mask  = (unsigned char *) mxGetData(prhs[1]); /* MASK */
+    lambda =  (float) mxGetScalar(prhs[2]); /* regularization parameter */
+    sigma = (float) mxGetScalar(prhs[3]); /* Edge-preserving parameter */
+    iter_numb = 300; /* iterations number */
+    tau = 0.025; /* marching step parameter */
+    penaltytype = 1; /* Huber penalty by default */    
+  
+    if ((nrhs < 4) || (nrhs > 7)) mexErrMsgTxt("At least 4 parameters is required, all parameters are: Image(2D/3D), Mask(2D/3D), Regularisation parameter, Edge-preserving parameter, iterations number, time-marching constant, penalty type - Huber, PM or Tukey");
+    if ((nrhs == 5) || (nrhs == 6) || (nrhs == 7))  iter_numb = (int) mxGetScalar(prhs[4]); /* iterations number */
+    if ((nrhs == 6) || (nrhs == 7))  tau =  (float) mxGetScalar(prhs[5]); /* marching step parameter */
+    if (nrhs == 7)  {
+        char *penalty_type;
+        penalty_type = mxArrayToString(prhs[6]); /* Huber, PM or Tukey 'Huber' is the default */
+        if ((strcmp(penalty_type, "Huber") != 0) && (strcmp(penalty_type, "PM") != 0) && (strcmp(penalty_type, "Tukey") != 0)) mexErrMsgTxt("Choose penalty: 'Huber', 'PM' or 'Tukey',");
+        if (strcmp(penalty_type, "Huber") == 0)  penaltytype = 1;  /* enable 'Huber' penalty */
+        if (strcmp(penalty_type, "PM") == 0)  penaltytype = 2;  /* enable Perona-Malik penalty */
+        if (strcmp(penalty_type, "Tukey") == 0)  penaltytype = 3;  /* enable Tikey Biweight penalty */
+        mxFree(penalty_type);
+    }    
+    
+    if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); }
+    if (mxGetClassID(prhs[1]) != mxUINT8_CLASS) {mexErrMsgTxt("The mask must be in uint8 precision");}
+    
+    dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2];
+    
+    /* output arrays*/
+    if (number_of_dims == 2) {
+        dimZ = 1; /*2D case*/
+        /* output image/volume */
+        if ((dimX != dim_array2[0]) || (dimY != dim_array2[1])) mexErrMsgTxt("Input image and the provided mask are of different dimensions!");
+        Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
+    }    
+    if (number_of_dims == 3) {
+        if ((dimX != dim_array2[0]) || (dimY != dim_array2[1]) || (dimZ != dim_array2[2])) mexErrMsgTxt("Input image and the provided mask are of different dimensions!");
+        Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
+    }    
+    
+    inpaint_elements = 0;
+    for (i=0; i<(int)(dimY*dimX*dimZ); i++) if (Mask[i] == 1) inpaint_elements++;
+    if (inpaint_elements == 0) mexErrMsgTxt("The mask is full of zeros, nothing to inpaint");        
+    Diffusion_Inpaint_CPU_main(Input, Mask, Output, lambda, sigma, iter_numb, tau, penaltytype, dimX, dimY, dimZ);
+}
\ No newline at end of file
diff --git a/src/Matlab/mex_compile/regularisers_CPU/NonlocalMarching_Inpaint.c b/src/Matlab/mex_compile/regularisers_CPU/NonlocalMarching_Inpaint.c
new file mode 100644
index 0000000..b3f2c98
--- /dev/null
+++ b/src/Matlab/mex_compile/regularisers_CPU/NonlocalMarching_Inpaint.c
@@ -0,0 +1,84 @@
+/*
+ * This work is part of the Core Imaging Library developed by
+ * Visual Analytics and Imaging System Group of the Science Technology
+ * Facilities Council, STFC
+ *
+ * Copyright 2017 Daniil Kazantsev
+ * Copyright 2017 Srikanth Nagella, Edoardo Pasca
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ * http://www.apache.org/licenses/LICENSE-2.0
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+#include "matrix.h"
+#include "mex.h"
+#include "NonlocalMarching_Inpaint_core.h"
+
+/* C-OMP implementation of Nonlocal Vertical Marching inpainting method (2D case)
+ * The method is heuristic but computationally efficent (especially for larger images).
+ * It developed specifically to smoothly inpaint horizontal or inclined missing data regions in sinograms
+ * The method WILL not work satisfactory if you have lengthy vertical stripes of missing data
+ *
+ * Input:
+ * 1. 2D image or sinogram [REQUIRED]
+ * 2. Mask of the same size as A in 'unsigned char' format  (ones mark the region to inpaint, zeros belong to the data) [REQUIRED]
+ * 3. Linear increment to increase searching window size in iterations, values from 1-3 is a good choice [OPTIONAL, default 1]
+ * 4. Number of iterations [OPTIONAL, default - calculate based on the mask]
+ *
+ * Output:
+ * 1. Inpainted sinogram  
+ * 2. updated mask
+ * Reference: TBA
+ */
+
+void mexFunction(
+        int nlhs, mxArray *plhs[],
+        int nrhs, const mxArray *prhs[])
+        
+{
+    int number_of_dims, iterations, SW_increment;
+    mwSize dimX, dimY, dimZ;
+    const mwSize *dim_array;
+    const mwSize *dim_array2;
+    
+    float *Input, *Output=NULL;
+    unsigned char *Mask, *Mask_upd=NULL;
+    
+    dim_array = mxGetDimensions(prhs[0]);
+    dim_array2 = mxGetDimensions(prhs[1]);
+    number_of_dims = mxGetNumberOfDimensions(prhs[0]);
+    
+    /*Handling Matlab input data*/
+    Input  = (float *) mxGetData(prhs[0]);
+    Mask  = (unsigned char *) mxGetData(prhs[1]); /* MASK */    
+    SW_increment = 1;
+    iterations = 0;
+            
+    if ((nrhs < 2) || (nrhs > 4)) mexErrMsgTxt("At least 4 parameters is required, all parameters are: Image(2D/3D), Mask(2D/3D), Linear increment, Iterations number");
+    if ((nrhs == 3) || (nrhs == 4))  SW_increment =  (int) mxGetScalar(prhs[2]); /* linear increment */
+    if ((nrhs == 4))  iterations =  (int) mxGetScalar(prhs[3]); /* iterations number */
+       
+    if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); }
+    if (mxGetClassID(prhs[1]) != mxUINT8_CLASS) {mexErrMsgTxt("The mask must be in uint8 precision");}    
+    
+    dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2];
+    
+    /* output arrays*/
+    if (number_of_dims == 2) {
+        dimZ = 1; /*2D case*/
+        /* output image/volume */
+        if ((dimX != dim_array2[0]) || (dimY != dim_array2[1])) mexErrMsgTxt("Input image and the provided mask are of different dimensions!");
+        Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
+        Mask_upd = (unsigned char*)mxGetPr(plhs[1] = mxCreateNumericArray(2, dim_array, mxUINT8_CLASS, mxREAL));
+    }    
+    if (number_of_dims == 3) {
+        mexErrMsgTxt("Currently 2D supported only");        
+    }           
+    NonlocalMarching_Inpaint_main(Input, Mask, Output, Mask_upd, SW_increment, iterations, 0, dimX, dimY, dimZ);
+}
\ No newline at end of file
diff --git a/src/Matlab/mex_compile/regularisers_CPU/Nonlocal_TV.c b/src/Matlab/mex_compile/regularisers_CPU/Nonlocal_TV.c
new file mode 100644
index 0000000..014c0a0
--- /dev/null
+++ b/src/Matlab/mex_compile/regularisers_CPU/Nonlocal_TV.c
@@ -0,0 +1,88 @@
+/*
+ * This work is part of the Core Imaging Library developed by
+ * Visual Analytics and Imaging System Group of the Science Technology
+ * Facilities Council, STFC and Diamond Light Source Ltd. 
+ *
+ * Copyright 2017 Daniil Kazantsev
+ * Copyright 2017 Srikanth Nagella, Edoardo Pasca
+ * Copyright 2018 Diamond Light Source Ltd. 
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ * http://www.apache.org/licenses/LICENSE-2.0
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#include "matrix.h"
+#include "mex.h"
+#include "Nonlocal_TV_core.h"
+
+#define EPS 1.0000e-9
+
+/* Matlab wrapper for C-OMP implementation of non-local regulariser
+ * Weights and associated indices must be given as an input.
+ * Gauss-Seidel fixed point iteration requires ~ 3 iterations, so the main effort
+ * goes in pre-calculation of weights and selection of patches
+ *
+ *
+ * Input Parameters:
+ * 1. 2D/3D grayscale image/volume
+ * 2. AR_i - indeces of i neighbours
+ * 3. AR_j - indeces of j neighbours
+ * 4. AR_k - indeces of k neighbours (0 - for 2D case)
+ * 5. Weights_ij(k) - associated weights 
+ * 6. regularisation parameter
+ * 7. iterations number 
+ 
+ * Output:
+ * 1. denoised image/volume 	
+ * Elmoataz, Abderrahim, Olivier Lezoray, and Sébastien Bougleux. "Nonlocal discrete regularization on weighted graphs: a framework for image and manifold processing." IEEE Trans. Image Processing 17, no. 7 (2008): 1047-1060.
+ */
+
+void mexFunction(
+        int nlhs, mxArray *plhs[],
+        int nrhs, const mxArray *prhs[])
+{
+    long number_of_dims,  dimX, dimY, dimZ;
+    int IterNumb, NumNeighb = 0;
+    unsigned short *H_i, *H_j, *H_k;
+    const int  *dim_array;
+    const int  *dim_array2;
+    float *A_orig, *Output=NULL, *Weights, lambda;
+    
+    dim_array = mxGetDimensions(prhs[0]);
+    dim_array2 = mxGetDimensions(prhs[1]);
+    number_of_dims = mxGetNumberOfDimensions(prhs[0]);
+    
+    /*Handling Matlab input data*/
+    A_orig  = (float *) mxGetData(prhs[0]); /* a 2D image or a set of 2D images (3D stack) */
+    H_i  = (unsigned short *) mxGetData(prhs[1]); /* indeces of i neighbours */
+    H_j  = (unsigned short *) mxGetData(prhs[2]); /* indeces of j neighbours */
+    H_k  = (unsigned short *) mxGetData(prhs[3]); /* indeces of k neighbours */
+    Weights = (float *) mxGetData(prhs[4]); /* weights for patches */
+    lambda = (float) mxGetScalar(prhs[5]); /* regularisation parameter */
+    IterNumb = (int) mxGetScalar(prhs[6]); /* the number of iterations */
+ 
+    dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2];   
+         
+    /*****2D INPUT *****/
+    if (number_of_dims == 2) {
+        dimZ = 0;   
+        NumNeighb = dim_array2[2];
+        Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));  
+        }
+    /*****3D INPUT *****/
+    /****************************************************/
+    if (number_of_dims == 3) {
+        NumNeighb = dim_array2[3];
+        Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
+    }
+    
+    /* run the main function here */
+    Nonlocal_TV_CPU_main(A_orig, Output, H_i, H_j, H_k, Weights, dimX, dimY, dimZ, NumNeighb, lambda, IterNumb);
+}
diff --git a/src/Matlab/mex_compile/regularisers_CPU/PatchSelect.c b/src/Matlab/mex_compile/regularisers_CPU/PatchSelect.c
new file mode 100644
index 0000000..f942539
--- /dev/null
+++ b/src/Matlab/mex_compile/regularisers_CPU/PatchSelect.c
@@ -0,0 +1,92 @@
+/*
+ * This work is part of the Core Imaging Library developed by
+ * Visual Analytics and Imaging System Group of the Science Technology
+ * Facilities Council, STFC and Diamond Light Source Ltd. 
+ *
+ * Copyright 2017 Daniil Kazantsev
+ * Copyright 2017 Srikanth Nagella, Edoardo Pasca
+ * Copyright 2018 Diamond Light Source Ltd. 
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ * http://www.apache.org/licenses/LICENSE-2.0
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#include "matrix.h"
+#include "mex.h"
+#include "PatchSelect_core.h"
+
+/* C-OMP implementation of non-local weight pre-calculation for non-local priors
+ * Weights and associated indices are stored into pre-allocated arrays and passed
+ * to the regulariser
+ *
+ *
+ * Input Parameters:
+ * 1. 2D/3D grayscale image/volume
+ * 2. Searching window (half-size of the main bigger searching window, e.g. 11)
+ * 3. Similarity window (half-size of the patch window, e.g. 2)
+ * 4. The number of neighbours to take (the most prominent after sorting neighbours will be taken)
+ * 5. noise-related parameter to calculate non-local weights
+ *
+ * Output [2D]:
+ * 1. AR_i - indeces of i neighbours
+ * 2. AR_j - indeces of j neighbours
+ * 3. Weights_ij - associated weights
+ *
+ * Output [3D]:
+ * 1. AR_i - indeces of i neighbours
+ * 2. AR_j - indeces of j neighbours
+ * 3. AR_k - indeces of j neighbours
+ * 4. Weights_ijk - associated weights
+ */
+/**************************************************/
+void mexFunction(
+        int nlhs, mxArray *plhs[],
+        int nrhs, const mxArray *prhs[])
+{
+    int number_of_dims,  SearchWindow, SimilarWin, NumNeighb;
+    mwSize dimX, dimY, dimZ;
+    unsigned short *H_i=NULL, *H_j=NULL, *H_k=NULL;
+    const int  *dim_array;
+    float *A, *Weights = NULL, h;
+    int dim_array2[3]; /* for 2D data */
+    int dim_array3[4]; /* for 3D data */
+    
+    dim_array = mxGetDimensions(prhs[0]);
+    number_of_dims = mxGetNumberOfDimensions(prhs[0]);
+    
+    /*Handling Matlab input data*/
+    A  = (float *) mxGetData(prhs[0]); /* a 2D or 3D image/volume */
+    SearchWindow = (int) mxGetScalar(prhs[1]);    /* Large Searching window */
+    SimilarWin = (int) mxGetScalar(prhs[2]);    /* Similarity window (patch-search)*/
+    NumNeighb = (int) mxGetScalar(prhs[3]); /* the total number of neighbours to take */
+    h = (float) mxGetScalar(prhs[4]); /* NLM parameter */
+
+    dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2];
+    dim_array2[0] = dimX; dim_array2[1] = dimY; dim_array2[2] = NumNeighb;  /* 2D case */
+    dim_array3[0] = dimX; dim_array3[1] = dimY; dim_array3[2] = dimZ; dim_array3[3] = NumNeighb;  /* 3D case */
+    
+    /****************2D INPUT ***************/
+    if (number_of_dims == 2) {
+        dimZ = 0;               
+        H_i = (unsigned short*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array2, mxUINT16_CLASS, mxREAL));
+        H_j = (unsigned short*)mxGetPr(plhs[1] = mxCreateNumericArray(3, dim_array2, mxUINT16_CLASS, mxREAL));
+        Weights = (float*)mxGetPr(plhs[2] = mxCreateNumericArray(3, dim_array2, mxSINGLE_CLASS, mxREAL));
+        }
+    /****************3D INPUT ***************/
+    if (number_of_dims == 3) {        
+        H_i = (unsigned short*)mxGetPr(plhs[0] = mxCreateNumericArray(4, dim_array3, mxUINT16_CLASS, mxREAL));
+        H_j = (unsigned short*)mxGetPr(plhs[1] = mxCreateNumericArray(4, dim_array3, mxUINT16_CLASS, mxREAL));
+        H_k = (unsigned short*)mxGetPr(plhs[2] = mxCreateNumericArray(4, dim_array3, mxUINT16_CLASS, mxREAL));
+        Weights = (float*)mxGetPr(plhs[3] = mxCreateNumericArray(4, dim_array3, mxSINGLE_CLASS, mxREAL));        
+    }
+    
+    PatchSelect_CPU_main(A, H_i, H_j, H_k, Weights, (long)(dimX), (long)(dimY), (long)(dimZ), SearchWindow, SimilarWin, NumNeighb, h, 0); 
+    
+ }
diff --git a/src/Matlab/mex_compile/regularisers_CPU/ROF_TV.c b/src/Matlab/mex_compile/regularisers_CPU/ROF_TV.c
new file mode 100644
index 0000000..55ef2b1
--- /dev/null
+++ b/src/Matlab/mex_compile/regularisers_CPU/ROF_TV.c
@@ -0,0 +1,77 @@
+
+/*
+ * This work is part of the Core Imaging Library developed by
+ * Visual Analytics and Imaging System Group of the Science Technology
+ * Facilities Council, STFC
+ *
+ * Copyright 2017 Daniil Kazantsev
+ * Copyright 2017 Srikanth Nagella, Edoardo Pasca
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ * http://www.apache.org/licenses/LICENSE-2.0
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+#include "matrix.h"
+#include "mex.h"
+#include "ROF_TV_core.h"
+
+/* ROF-TV denoising/regularization model [1] (2D/3D case)
+ * (MEX wrapper for MATLAB)
+ * 
+ * Input Parameters:
+ * 1. Noisy image/volume [REQUIRED]
+ * 2. lambda - regularization parameter [REQUIRED]
+ * 3. Number of iterations, for explicit scheme >= 150 is recommended  [REQUIRED]
+ * 4. tau - marching step for explicit scheme, ~1 is recommended [REQUIRED]
+ *
+ * Output:
+ * [1] Regularized image/volume 
+ *
+ * This function is based on the paper by
+ * [1] Rudin, Osher, Fatemi, "Nonlinear Total Variation based noise removal algorithms"
+ *
+ * D. Kazantsev, 2016-18
+ */
+
+void mexFunction(
+        int nlhs, mxArray *plhs[],
+        int nrhs, const mxArray *prhs[])
+        
+{
+    int number_of_dims, iter_numb;
+    mwSize dimX, dimY, dimZ;
+    const mwSize *dim_array_i;
+    float *Input, *Output=NULL, lambda, tau;    
+    
+    dim_array_i = mxGetDimensions(prhs[0]);
+    number_of_dims = mxGetNumberOfDimensions(prhs[0]);
+    
+    /*Handling Matlab input data*/
+    Input  = (float *) mxGetData(prhs[0]);
+    lambda =  (float) mxGetScalar(prhs[1]); /* regularization parameter */
+    iter_numb =  (int) mxGetScalar(prhs[2]); /* iterations number */
+    tau =  (float) mxGetScalar(prhs[3]); /* marching step parameter */  
+    
+    if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); }
+    if(nrhs != 4) mexErrMsgTxt("Four inputs reqired: Image(2D,3D), regularization parameter, iterations number,  marching step constant");
+    /*Handling Matlab output data*/
+    dimX = dim_array_i[0]; dimY = dim_array_i[1]; dimZ = dim_array_i[2];        
+    
+    /* output arrays*/
+    if (number_of_dims == 2) {
+        dimZ = 1; /*2D case*/
+        /* output image/volume */
+        Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array_i, mxSINGLE_CLASS, mxREAL));          
+    }    
+    if (number_of_dims == 3) {
+        Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array_i, mxSINGLE_CLASS, mxREAL));
+    }
+     
+    TV_ROF_CPU_main(Input, Output, lambda, iter_numb, tau, dimX, dimY, dimZ);    
+}
\ No newline at end of file
diff --git a/src/Matlab/mex_compile/regularisers_CPU/SB_TV.c b/src/Matlab/mex_compile/regularisers_CPU/SB_TV.c
new file mode 100644
index 0000000..8636322
--- /dev/null
+++ b/src/Matlab/mex_compile/regularisers_CPU/SB_TV.c
@@ -0,0 +1,91 @@
+/*
+ * This work is part of the Core Imaging Library developed by
+ * Visual Analytics and Imaging System Group of the Science Technology
+ * Facilities Council, STFC
+ *
+ * Copyright 2017 Daniil Kazantsev
+ * Copyright 2017 Srikanth Nagella, Edoardo Pasca
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ * http://www.apache.org/licenses/LICENSE-2.0
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+#include "matrix.h"
+#include "mex.h"
+#include "SB_TV_core.h"
+
+/* C-OMP implementation of Split Bregman - TV denoising-regularisation model (2D/3D) [1]
+*
+* Input Parameters:
+* 1. Noisy image/volume
+* 2. lambda - regularisation parameter
+* 3. Number of iterations [OPTIONAL parameter]
+* 4. eplsilon - tolerance constant [OPTIONAL parameter]
+* 5. TV-type: 'iso' or 'l1' [OPTIONAL parameter]
+* 6. print information: 0 (off) or 1 (on)  [OPTIONAL parameter]
+*
+* Output:
+* 1. Filtered/regularized image
+*
+* This function is based on the Matlab's code and paper by
+* [1]. Goldstein, T. and Osher, S., 2009. The split Bregman method for L1-regularized problems. SIAM journal on imaging sciences, 2(2), pp.323-343.
+*/
+
+void mexFunction(
+        int nlhs, mxArray *plhs[],
+        int nrhs, const mxArray *prhs[])
+        
+{
+    int number_of_dims, iter, methTV, printswitch;
+    mwSize dimX, dimY, dimZ;
+    const mwSize *dim_array;
+    
+    float *Input, *Output=NULL, lambda, epsil;
+    
+    number_of_dims = mxGetNumberOfDimensions(prhs[0]);
+    dim_array = mxGetDimensions(prhs[0]);
+    
+    /*Handling Matlab input data*/
+    if ((nrhs < 2) || (nrhs > 6)) mexErrMsgTxt("At least 2 parameters is required, all parameters are: Image(2D/3D), Regularization parameter, Regularization parameter, iterations number, tolerance, penalty type ('iso' or 'l1'), print switch");
+    
+    Input  = (float *) mxGetData(prhs[0]); /*noisy image (2D/3D) */
+    lambda =  (float) mxGetScalar(prhs[1]); /* regularization parameter */
+    iter = 100; /* default iterations number */
+    epsil = 0.0001; /* default tolerance constant */
+    methTV = 0;  /* default isotropic TV penalty */
+    printswitch = 0; /*default print is switched, off - 0 */
+    
+    if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); }
+    
+    if ((nrhs == 3) || (nrhs == 4) || (nrhs == 5) || (nrhs == 6))  iter = (int) mxGetScalar(prhs[2]); /* iterations number */
+    if ((nrhs == 4) || (nrhs == 5) || (nrhs == 6))  epsil =  (float) mxGetScalar(prhs[3]); /* tolerance constant */
+    if ((nrhs == 5) || (nrhs == 6))  {
+        char *penalty_type;
+        penalty_type = mxArrayToString(prhs[4]); /* choosing TV penalty: 'iso' or 'l1', 'iso' is the default */
+        if ((strcmp(penalty_type, "l1") != 0) && (strcmp(penalty_type, "iso") != 0)) mexErrMsgTxt("Choose TV type: 'iso' or 'l1',");
+        if (strcmp(penalty_type, "l1") == 0)  methTV = 1;  /* enable 'l1' penalty */
+        mxFree(penalty_type);
+    }
+    if (nrhs == 6)  {
+        printswitch = (int) mxGetScalar(prhs[5]);
+        if ((printswitch != 0) && (printswitch != 1)) mexErrMsgTxt("Print can be enabled by choosing 1 or off - 0");
+    }
+    
+    /*Handling Matlab output data*/
+    dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2];
+    
+    if (number_of_dims == 2) {
+        dimZ = 1; /*2D case*/
+        Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
+    }
+    if (number_of_dims == 3) Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
+    
+    /* running the function */
+    SB_TV_CPU_main(Input, Output, lambda, iter, epsil, methTV, printswitch, dimX, dimY, dimZ);
+}
diff --git a/src/Matlab/mex_compile/regularisers_CPU/TGV.c b/src/Matlab/mex_compile/regularisers_CPU/TGV.c
new file mode 100644
index 0000000..aa4eed4
--- /dev/null
+++ b/src/Matlab/mex_compile/regularisers_CPU/TGV.c
@@ -0,0 +1,83 @@
+/*
+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 "mex.h"
+#include "TGV_core.h"
+
+/* C-OMP implementation of Primal-Dual denoising method for 
+ * Total Generilized Variation (TGV)-L2 model [1] (2D/3D)
+ *
+ * Input Parameters:
+ * 1. Noisy image/volume (2D/3D)
+ * 2. lambda - regularisation parameter
+ * 3. parameter to control the first-order term (alpha1)
+ * 4. parameter to control the second-order term (alpha0)
+ * 5. Number of Chambolle-Pock (Primal-Dual) iterations
+ * 6. Lipshitz constant (default is 12)
+ *
+ * Output:
+ * Filtered/regulariaed image 
+ *
+ * References:
+ * [1] K. Bredies "Total Generalized Variation"
+ */
+
+void mexFunction(
+        int nlhs, mxArray *plhs[],
+        int nrhs, const mxArray *prhs[])
+        
+{
+    int number_of_dims, iter;
+    mwSize dimX, dimY, dimZ;
+    const mwSize *dim_array;
+    
+    float *Input, *Output=NULL, lambda, alpha0, alpha1, L2;
+    
+    number_of_dims = mxGetNumberOfDimensions(prhs[0]);
+    dim_array = mxGetDimensions(prhs[0]);
+    
+    /*Handling Matlab input data*/
+    if ((nrhs < 2) || (nrhs > 6)) mexErrMsgTxt("At least 2 parameters is required, all parameters are: Image(2D), Regularisation parameter, alpha0, alpha1, iterations number, Lipshitz Constant");
+    
+    Input  = (float *) mxGetData(prhs[0]); /*noisy image/volume */
+    lambda =  (float) mxGetScalar(prhs[1]); /* regularisation parameter */
+    alpha1 =  1.0f; /* parameter to control the first-order term */ 
+    alpha0 =  0.5f; /* parameter to control the second-order term */
+    iter =  300; /* Iterations number */      
+    L2 =  12.0f; /* Lipshitz constant */
+    
+    if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); }   
+    if ((nrhs == 3) || (nrhs == 4) || (nrhs == 5) || (nrhs == 6))  alpha1 =  (float) mxGetScalar(prhs[2]); /* parameter to control the first-order term */ 
+    if ((nrhs == 4) || (nrhs == 5) || (nrhs == 6))  alpha0 =  (float) mxGetScalar(prhs[3]);  /* parameter to control the second-order term */
+    if ((nrhs == 5) || (nrhs == 6))  iter =  (int) mxGetScalar(prhs[4]); /* Iterations number */      
+    if (nrhs == 6)  L2 =  (float) mxGetScalar(prhs[5]); /* Lipshitz constant */
+    
+    /*Handling Matlab output data*/
+    dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2];
+    
+    if (number_of_dims == 2) {
+        dimZ = 1; /*2D case*/
+        Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));        
+    }
+    if (number_of_dims == 3) {
+        Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
+    }       
+    /* running the function */
+    TGV_main(Input, Output, lambda, alpha1, alpha0, iter, L2, dimX, dimY, dimZ);        
+}
diff --git a/src/Matlab/mex_compile/regularisers_CPU/TNV.c b/src/Matlab/mex_compile/regularisers_CPU/TNV.c
new file mode 100644
index 0000000..acea75d
--- /dev/null
+++ b/src/Matlab/mex_compile/regularisers_CPU/TNV.c
@@ -0,0 +1,74 @@
+/*
+ * This work is part of the Core Imaging Library developed by
+ * Visual Analytics and Imaging System Group of the Science Technology
+ * Facilities Council, STFC
+ *
+ * Copyright 2017 Daniil Kazantsev
+ * Copyright 2017 Srikanth Nagella, Edoardo Pasca
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ * http://www.apache.org/licenses/LICENSE-2.0
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+#include "matrix.h"
+#include "mex.h"
+#include "TNV_core.h"
+/*
+ * C-OMP implementation of Total Nuclear Variation regularisation model (2D + channels) [1]
+ * The code is modified from the implementation by Joan Duran <joan.duran@uib.es> see
+ * "denoisingPDHG_ipol.cpp" in Joans Collaborative Total Variation package
+ *
+ * Input Parameters:
+ * 1. Noisy volume of 2D + channel dimension, i.e. 3D volume
+ * 2. lambda - regularisation parameter
+ * 3. Number of iterations [OPTIONAL parameter]
+ * 4. eplsilon - tolerance constant [OPTIONAL parameter]
+ * 5. print information: 0 (off) or 1 (on)  [OPTIONAL parameter]
+ *
+ * Output:
+ * 1. Filtered/regularized image
+ *
+ * [1]. Duran, J., Moeller, M., Sbert, C. and Cremers, D., 2016. Collaborative total variation: a general framework for vectorial TV models. SIAM Journal on Imaging Sciences, 9(1), pp.116-151.
+ */
+void mexFunction(
+        int nlhs, mxArray *plhs[],
+        int nrhs, const mxArray *prhs[])
+        
+{
+    int number_of_dims, iter;
+    mwSize dimX, dimY, dimZ;
+    const mwSize *dim_array;
+    float *Input, *Output=NULL, lambda, epsil;
+    
+    number_of_dims = mxGetNumberOfDimensions(prhs[0]);
+    dim_array = mxGetDimensions(prhs[0]);
+    
+    /*Handling Matlab input data*/
+    if ((nrhs < 2) || (nrhs > 4)) mexErrMsgTxt("At least 2 parameters is required, all parameters are: Image(2D + channels), Regularisation parameter, Regularization parameter, iterations number, tolerance");
+    
+    Input  = (float *) mxGetData(prhs[0]); /* noisy sequence of channels (2D + channels) */
+    lambda =  (float) mxGetScalar(prhs[1]); /* regularization parameter */
+    iter = 1000; /* default iterations number */
+    epsil = 1.00e-05; /* default tolerance constant */
+    
+    if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); }
+    
+    if ((nrhs == 3) || (nrhs == 4))  iter = (int) mxGetScalar(prhs[2]); /* iterations number */
+    if (nrhs == 4)  epsil =  (float) mxGetScalar(prhs[3]); /* tolerance constant */
+    
+    /*Handling Matlab output data*/
+    dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2];
+    
+    if (number_of_dims == 2) mexErrMsgTxt("The input must be 3D: [X,Y,Channels]");
+    if (number_of_dims == 3) {
+        Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
+        /* running the function */
+        TNV_CPU_main(Input, Output, lambda, iter, epsil, dimX, dimY, dimZ);
+    }
+}
\ No newline at end of file
diff --git a/src/Matlab/mex_compile/regularisers_CPU/TV_energy.c b/src/Matlab/mex_compile/regularisers_CPU/TV_energy.c
new file mode 100644
index 0000000..d457f46
--- /dev/null
+++ b/src/Matlab/mex_compile/regularisers_CPU/TV_energy.c
@@ -0,0 +1,72 @@
+/*
+ * This work is part of the Core Imaging Library developed by
+ * Visual Analytics and Imaging System Group of the Science Technology
+ * Facilities Council, STFC
+ *
+ * Copyright 2017 Daniil Kazantsev
+ * Copyright 2017 Srikanth Nagella, Edoardo Pasca
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ * http://www.apache.org/licenses/LICENSE-2.0
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+#include "matrix.h"
+#include "mex.h"
+#include "utils.h"
+/*
+ * Function to calculate TV energy value with respect to the denoising variational problem
+ * 
+ * Input:
+ * 1. Denoised Image/volume
+ * 2. Original (noisy) Image/volume
+ * 3. lambda - regularisation parameter 
+ * 
+ * Output:
+ * 1. Energy function value
+ * 
+ */
+void mexFunction(
+        int nlhs, mxArray *plhs[],
+        int nrhs, const mxArray *prhs[])
+        
+{
+    int number_of_dims, type;
+    
+    mwSize dimX, dimY, dimZ;
+    const mwSize *dim_array;
+    float *Input, *Input0, lambda;
+    
+    number_of_dims = mxGetNumberOfDimensions(prhs[0]);
+    dim_array = mxGetDimensions(prhs[0]);
+    
+    /*Handling Matlab input data*/
+    if ((nrhs != 4)) mexErrMsgTxt("4 inputs: Two images or volumes of the same size required, estimated and the original (noisy), regularisation parameter, type");
+    
+    Input  = (float *) mxGetData(prhs[0]); /* Denoised Image/volume */
+    Input0  = (float *) mxGetData(prhs[1]); /* Original (noisy) Image/volume */
+    lambda =  (float) mxGetScalar(prhs[2]); /* regularisation parameter */
+    type =  (int) mxGetScalar(prhs[3]); /* type of energy */
+    
+    if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); }
+    if (mxGetClassID(prhs[1]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); }
+    
+    /*output energy function value */
+    plhs[0] = mxCreateNumericMatrix(1, 1, mxSINGLE_CLASS, mxREAL);
+    float *funcvalA = (float *) mxGetData(plhs[0]);
+    
+    /*Handling Matlab output data*/
+    dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2];
+    
+    if (number_of_dims == 2) {
+		TV_energy2D(Input, Input0, funcvalA, lambda, type, dimX, dimY);
+		}
+    if (number_of_dims == 3) {
+        TV_energy3D(Input, Input0, funcvalA, lambda, type, dimX, dimY, dimZ);
+    }
+}
diff --git a/src/Matlab/mex_compile/regularisers_GPU/Diffusion_4thO_GPU.cpp b/src/Matlab/mex_compile/regularisers_GPU/Diffusion_4thO_GPU.cpp
new file mode 100644
index 0000000..0cc042b
--- /dev/null
+++ b/src/Matlab/mex_compile/regularisers_GPU/Diffusion_4thO_GPU.cpp
@@ -0,0 +1,77 @@
+/*
+ * This work is part of the Core Imaging Library developed by
+ * Visual Analytics and Imaging System Group of the Science Technology
+ * Facilities Council, STFC
+ *
+ * Copyright 2017 Daniil Kazantsev
+ * Copyright 2017 Srikanth Nagella, Edoardo Pasca
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ * http://www.apache.org/licenses/LICENSE-2.0
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+#include "matrix.h"
+#include "mex.h"
+#include "Diffus_4thO_GPU_core.h"
+
+/* CUDA implementation of fourth-order diffusion scheme [1] for piecewise-smooth recovery (2D/3D case)
+ * The minimisation is performed using explicit scheme. 
+ *
+ * Input Parameters:
+ * 1. Noisy image/volume [REQUIRED]
+ * 2. lambda - regularization parameter [REQUIRED]
+ * 3. Edge-preserving parameter (sigma) [REQUIRED]
+ * 4. Number of iterations, for explicit scheme >= 150 is recommended [OPTIONAL, default 300]
+ * 5. tau - time-marching step for the explicit scheme [OPTIONAL, default 0.015]
+ *
+ * Output:
+ * [1] Regularized image/volume 
+ *
+ * This function is based on the paper by
+ * [1] Hajiaboli, M.R., 2011. An anisotropic fourth-order diffusion filter for image noise removal. International Journal of Computer Vision, 92(2), pp.177-191.
+ */
+
+void mexFunction(
+        int nlhs, mxArray *plhs[],
+        int nrhs, const mxArray *prhs[])
+        
+{
+    int number_of_dims, iter_numb;
+    mwSize dimX, dimY, dimZ;
+    const mwSize *dim_array;
+    float *Input, *Output=NULL, lambda, tau, sigma;
+    
+    dim_array = mxGetDimensions(prhs[0]);
+    number_of_dims = mxGetNumberOfDimensions(prhs[0]);
+    
+    /*Handling Matlab input data*/
+    Input  = (float *) mxGetData(prhs[0]);
+    lambda =  (float) mxGetScalar(prhs[1]); /* regularization parameter */
+    sigma = (float) mxGetScalar(prhs[2]); /* Edge-preserving parameter */
+    iter_numb = 300; /* iterations number */
+    tau = 0.01; /* marching step parameter */
+    
+    if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); }
+    if ((nrhs < 3) || (nrhs > 5)) mexErrMsgTxt("At least 3 parameters is required, all parameters are: Image(2D/3D), Regularisation parameter, Edge-preserving parameter, iterations number, time-marching constant");
+    if ((nrhs == 4) || (nrhs == 5))  iter_numb = (int) mxGetScalar(prhs[3]); /* iterations number */
+    if (nrhs == 5)  tau =  (float) mxGetScalar(prhs[4]); /* marching step parameter */
+    
+    /*Handling Matlab output data*/
+    dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2];
+    
+    /* output arrays*/
+    if (number_of_dims == 2) {
+        dimZ = 1; /*2D case*/
+        /* output image/volume */
+        Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
+    }
+    if (number_of_dims == 3) Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
+    
+    Diffus4th_GPU_main(Input, Output, lambda, sigma, iter_numb, tau, dimX, dimY, dimZ);
+}
\ No newline at end of file
diff --git a/src/Matlab/mex_compile/regularisers_GPU/FGP_TV_GPU.cpp b/src/Matlab/mex_compile/regularisers_GPU/FGP_TV_GPU.cpp
new file mode 100644
index 0000000..c174e75
--- /dev/null
+++ b/src/Matlab/mex_compile/regularisers_GPU/FGP_TV_GPU.cpp
@@ -0,0 +1,97 @@
+/*
+ * This work is part of the Core Imaging Library developed by
+ * Visual Analytics and Imaging System Group of the Science Technology
+ * Facilities Council, STFC
+ *
+ * Copyright 2017 Daniil Kazantsev
+ * Copyright 2017 Srikanth Nagella, Edoardo Pasca
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ * http://www.apache.org/licenses/LICENSE-2.0
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+#include "matrix.h"
+#include "mex.h"
+#include "TV_FGP_GPU_core.h"
+
+/* GPU (CUDA) implementation of FGP-TV [1] denoising/regularization model (2D/3D case)
+ *
+ * Input Parameters:
+ * 1. Noisy image/volume
+ * 2. lambdaPar - regularization parameter
+ * 3. Number of iterations
+ * 4. eplsilon: tolerance constant
+ * 5. TV-type: methodTV - 'iso' (0) or 'l1' (1)
+ * 6. nonneg: 'nonnegativity (0 is OFF by default)
+ * 7. print information: 0 (off) or 1 (on)
+ *
+ * Output:
+ * [1] Filtered/regularized image
+ *
+ * This function is based on the Matlab's code and paper by
+ * [1] Amir Beck and Marc Teboulle, "Fast Gradient-Based Algorithms for Constrained Total Variation Image Denoising and Deblurring Problems"
+ */
+
+void mexFunction(
+        int nlhs, mxArray *plhs[],
+        int nrhs, const mxArray *prhs[])
+        
+{
+    int number_of_dims, iter, methTV, printswitch, nonneg;
+    mwSize dimX, dimY, dimZ;
+    const mwSize *dim_array;
+    
+    float *Input, *Output=NULL, lambda, epsil;
+    
+    number_of_dims = mxGetNumberOfDimensions(prhs[0]);
+    dim_array = mxGetDimensions(prhs[0]);
+    
+    /*Handling Matlab input data*/
+    if ((nrhs < 2) || (nrhs > 7)) mexErrMsgTxt("At least 2 parameters is required, all parameters are: Image(2D/3D), Regularization parameter. The full list of parameters: Image(2D/3D), Regularization parameter, iterations number, tolerance, penalty type ('iso' or 'l1'), nonnegativity switch, print switch");
+    
+    Input  = (float *) mxGetData(prhs[0]); /*noisy image (2D/3D) */
+    lambda =  (float) mxGetScalar(prhs[1]); /* regularization parameter */
+    iter = 300; /* default iterations number */
+    epsil = 0.0001; /* default tolerance constant */
+    methTV = 0;  /* default isotropic TV penalty */
+    nonneg = 0; /* default nonnegativity switch, off - 0 */
+    printswitch = 0; /*default print is switched, off - 0 */
+    
+    if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); }
+    
+    if ((nrhs == 3) || (nrhs == 4) || (nrhs == 5) || (nrhs == 6) || (nrhs == 7))  iter = (int) mxGetScalar(prhs[2]); /* iterations number */
+    if ((nrhs == 4) || (nrhs == 5) || (nrhs == 6) || (nrhs == 7))  epsil =  (float) mxGetScalar(prhs[3]); /* tolerance constant */
+    if ((nrhs == 5) || (nrhs == 6) || (nrhs == 7))  {
+        char *penalty_type;
+        penalty_type = mxArrayToString(prhs[4]); /* choosing TV penalty: 'iso' or 'l1', 'iso' is the default */
+        if ((strcmp(penalty_type, "l1") != 0) && (strcmp(penalty_type, "iso") != 0)) mexErrMsgTxt("Choose TV type: 'iso' or 'l1',");
+        if (strcmp(penalty_type, "l1") == 0)  methTV = 1;  /* enable 'l1' penalty */
+        mxFree(penalty_type);
+    }
+    if ((nrhs == 6) || (nrhs == 7))  {
+        nonneg = (int) mxGetScalar(prhs[5]);
+        if ((nonneg != 0) && (nonneg != 1)) mexErrMsgTxt("Nonnegativity constraint can be enabled by choosing 1 or off - 0");
+    }
+    if (nrhs == 7)  {
+        printswitch = (int) mxGetScalar(prhs[6]);
+        if ((printswitch != 0) && (printswitch != 1)) mexErrMsgTxt("Print can be enabled by choosing 1 or off - 0");
+    }
+    
+    /*Handling Matlab output data*/
+    dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2];
+    
+    if (number_of_dims == 2) {
+        dimZ = 1; /*2D case*/
+        Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
+    }
+    if (number_of_dims == 3) Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
+    
+    /* running the function */
+    TV_FGP_GPU_main(Input, Output, lambda, iter, epsil, methTV, nonneg, printswitch, dimX, dimY, dimZ);    
+}
\ No newline at end of file
diff --git a/src/Matlab/mex_compile/regularisers_GPU/FGP_dTV_GPU.cpp b/src/Matlab/mex_compile/regularisers_GPU/FGP_dTV_GPU.cpp
new file mode 100644
index 0000000..3f5a4b3
--- /dev/null
+++ b/src/Matlab/mex_compile/regularisers_GPU/FGP_dTV_GPU.cpp
@@ -0,0 +1,113 @@
+/*
+ * This work is part of the Core Imaging Library developed by
+ * Visual Analytics and Imaging System Group of the Science Technology
+ * Facilities Council, STFC
+ *
+ * Copyright 2017 Daniil Kazantsev
+ * Copyright 2017 Srikanth Nagella, Edoardo Pasca
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ * http://www.apache.org/licenses/LICENSE-2.0
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+#include "matrix.h"
+#include "mex.h"
+#include "dTV_FGP_GPU_core.h"
+
+/* CUDA implementation of FGP-dTV [1,2] denoising/regularization model (2D/3D case)
+ * which employs structural similarity of the level sets of two images/volumes, see [1,2]
+ * The current implementation updates image 1 while image 2 is being fixed.
+ *
+ * Input Parameters:
+ * 1. Noisy image/volume [REQUIRED]
+ * 2. Additional reference image/volume of the same dimensions as (1) [REQUIRED]
+ * 3. lambdaPar - regularization parameter [REQUIRED]
+ * 4. Number of iterations [OPTIONAL]
+ * 5. eplsilon: tolerance constant [OPTIONAL]
+ * 6. eta: smoothing constant to calculate gradient of the reference [OPTIONAL] * 
+ * 7. TV-type: methodTV - 'iso' (0) or 'l1' (1) [OPTIONAL]
+ * 8. nonneg: 'nonnegativity (0 is OFF by default) [OPTIONAL]
+ * 9. print information: 0 (off) or 1 (on) [OPTIONAL]
+ *
+ * Output:
+ * [1] Filtered/regularized image/volume
+ *
+ * This function is based on the Matlab's codes and papers by
+ * [1] Amir Beck and Marc Teboulle, "Fast Gradient-Based Algorithms for Constrained Total Variation Image Denoising and Deblurring Problems"
+ * [2] M. J. Ehrhardt and M. M. Betcke, Multi-Contrast MRI Reconstruction with Structure-Guided Total Variation, SIAM Journal on Imaging Sciences 9(3), pp. 1084–1106
+ */
+void mexFunction(
+        int nlhs, mxArray *plhs[],
+        int nrhs, const mxArray *prhs[])
+        
+{
+    int number_of_dims, iter, methTV, printswitch, nonneg;
+    mwSize dimX, dimY, dimZ;
+    const mwSize *dim_array;
+    const mwSize *dim_array2;
+    
+    float *Input, *InputRef, *Output=NULL, lambda, epsil, eta;
+    
+    number_of_dims = mxGetNumberOfDimensions(prhs[0]);
+    dim_array = mxGetDimensions(prhs[0]);
+    dim_array2 = mxGetDimensions(prhs[1]);
+    
+    /*Handling Matlab input data*/
+    if ((nrhs < 3) || (nrhs > 9)) mexErrMsgTxt("At least 3 parameters is required, all parameters are: Image(2D/3D), Reference(2D/3D), Regularization parameter, iterations number, tolerance, smoothing constant, penalty type ('iso' or 'l1'), nonnegativity switch, print switch");
+    
+    Input  = (float *) mxGetData(prhs[0]); /*noisy image (2D/3D) */
+    InputRef  = (float *) mxGetData(prhs[1]); /* reference image (2D/3D) */
+    lambda =  (float) mxGetScalar(prhs[2]); /* regularization parameter */
+    iter = 300; /* default iterations number */
+    epsil = 0.0001; /* default tolerance constant */
+    eta = 0.01; /* default smoothing constant */
+    methTV = 0;  /* default isotropic TV penalty */
+    nonneg = 0; /* default nonnegativity switch, off - 0 */
+    printswitch = 0; /*default print is switched, off - 0 */
+    
+        
+    if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); }
+    if (mxGetClassID(prhs[1]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); }
+    
+    /*Handling Matlab output data*/
+    dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2];
+    if (number_of_dims == 2) { if ((dimX != dim_array2[0]) || (dimY != dim_array2[1])) mexErrMsgTxt("The input images have different dimensionalities");}
+    if (number_of_dims == 3) { if ((dimX != dim_array2[0]) || (dimY != dim_array2[1]) || (dimZ != dim_array2[2])) mexErrMsgTxt("The input volumes have different dimensionalities");}   
+    
+    
+    if ((nrhs == 4) || (nrhs == 5) || (nrhs == 6) || (nrhs == 7) || (nrhs == 8) || (nrhs == 9))  iter = (int) mxGetScalar(prhs[3]); /* iterations number */
+    if ((nrhs == 5) || (nrhs == 6) || (nrhs == 7) || (nrhs == 8) || (nrhs == 9))  epsil =  (float) mxGetScalar(prhs[4]); /* tolerance constant */
+    if ((nrhs == 6) || (nrhs == 7) || (nrhs == 8) || (nrhs == 9))  {
+    eta =  (float) mxGetScalar(prhs[5]); /* smoothing constant for the gradient of InputRef */
+    }
+    if ((nrhs == 7) || (nrhs == 8) || (nrhs == 9))  {        
+        char *penalty_type;
+        penalty_type = mxArrayToString(prhs[6]); /* choosing TV penalty: 'iso' or 'l1', 'iso' is the default */
+        if ((strcmp(penalty_type, "l1") != 0) && (strcmp(penalty_type, "iso") != 0)) mexErrMsgTxt("Choose TV type: 'iso' or 'l1',");
+        if (strcmp(penalty_type, "l1") == 0)  methTV = 1;  /* enable 'l1' penalty */
+        mxFree(penalty_type);
+    }    
+    if ((nrhs == 8) || (nrhs == 9))  {
+        nonneg = (int) mxGetScalar(prhs[7]);
+        if ((nonneg != 0) && (nonneg != 1)) mexErrMsgTxt("Nonnegativity constraint can be enabled by choosing 1 or off - 0");
+    }
+    if (nrhs == 9)  {
+        printswitch = (int) mxGetScalar(prhs[8]);
+        if ((printswitch != 0) && (printswitch != 1)) mexErrMsgTxt("Print can be enabled by choosing 1 or off - 0");
+    }    
+   
+    if (number_of_dims == 2) {
+        dimZ = 1; /*2D case*/
+        Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
+    }
+    if (number_of_dims == 3) Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
+    
+    /* running the function */
+    dTV_FGP_GPU_main(Input, InputRef, Output, lambda, iter, epsil, eta, methTV, nonneg, printswitch, dimX, dimY, dimZ);
+}
\ No newline at end of file
diff --git a/src/Matlab/mex_compile/regularisers_GPU/LLT_ROF_GPU.cpp b/src/Matlab/mex_compile/regularisers_GPU/LLT_ROF_GPU.cpp
new file mode 100644
index 0000000..e8da4ce
--- /dev/null
+++ b/src/Matlab/mex_compile/regularisers_GPU/LLT_ROF_GPU.cpp
@@ -0,0 +1,83 @@
+/*
+ * This work is part of the Core Imaging Library developed by
+ * Visual Analytics and Imaging System Group of the Science Technology
+ * Facilities Council, STFC
+ *
+ * Copyright 2017 Daniil Kazantsev
+ * Copyright 2017 Srikanth Nagella, Edoardo Pasca
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ * http://www.apache.org/licenses/LICENSE-2.0
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+#include "matrix.h"
+#include "mex.h"
+#include "LLT_ROF_GPU_core.h"
+
+/* CUDA implementation of Lysaker, Lundervold and Tai (LLT) model [1] combined with Rudin-Osher-Fatemi [2] TV regularisation penalty.
+* 
+* This penalty can deliver visually pleasant piecewise-smooth recovery if regularisation parameters are selected well. 
+* The rule of thumb for selection is to start with lambdaLLT = 0 (just the ROF-TV model) and then proceed to increase 
+* lambdaLLT starting with smaller values. 
+*
+* Input Parameters:
+* 1. U0 - original noise image/volume
+* 2. lambdaROF - ROF-related regularisation parameter
+* 3. lambdaLLT - LLT-related regularisation parameter
+* 4. tau - time-marching step 
+* 5. iter - iterations number (for both models)
+*
+* Output:
+* Filtered/regularised image
+*
+* References: 
+* [1] Lysaker, M., Lundervold, A. and Tai, X.C., 2003. Noise removal using fourth-order partial differential equation with applications to medical magnetic resonance images in space and time. IEEE Transactions on image processing, 12(12), pp.1579-1590.
+* [2] Rudin, Osher, Fatemi, "Nonlinear Total Variation based noise removal algorithms"
+*/
+
+void mexFunction(
+        int nlhs, mxArray *plhs[],
+        int nrhs, const mxArray *prhs[])
+        
+{
+    int number_of_dims, iterationsNumb;
+    mwSize dimX, dimY, dimZ;
+    const mwSize *dim_array;
+    
+    float *Input, *Output=NULL, lambdaROF, lambdaLLT, tau;
+    
+    dim_array = mxGetDimensions(prhs[0]);
+    number_of_dims = mxGetNumberOfDimensions(prhs[0]);
+    
+    if ((nrhs < 3) || (nrhs > 5)) mexErrMsgTxt("At least 3 parameters is required, all parameters are: Image(2D/3D), Regularisation parameter (ROF), Regularisation parameter (LTT), iterations number, time-marching parameter");
+    
+    /*Handling Matlab input data*/
+    Input  = (float *) mxGetData(prhs[0]);
+    lambdaROF =  (float) mxGetScalar(prhs[1]); /* ROF regularization parameter */
+    lambdaLLT =  (float) mxGetScalar(prhs[2]); /* ROF regularization parameter */    
+    iterationsNumb = 250;
+    tau =  0.0025;
+    
+    if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); }   
+    if ((nrhs == 4) || (nrhs == 5)) iterationsNumb =  (int) mxGetScalar(prhs[3]); /* iterations number */    
+    if (nrhs == 5) tau =  (float) mxGetScalar(prhs[4]); /* marching step parameter */  
+        
+    /*Handling Matlab output data*/
+    dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2];
+    
+    /* output arrays*/
+    if (number_of_dims == 2) {
+        dimZ = 1; /*2D case*/
+        /* output image/volume */
+        Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));                        
+    }    
+    if (number_of_dims == 3) Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));   
+  
+    LLT_ROF_GPU_main(Input, Output, lambdaROF, lambdaLLT, iterationsNumb, tau, dimX, dimY, dimZ);    
+}
\ No newline at end of file
diff --git a/src/Matlab/mex_compile/regularisers_GPU/NonlDiff_GPU.cpp b/src/Matlab/mex_compile/regularisers_GPU/NonlDiff_GPU.cpp
new file mode 100644
index 0000000..1cd0cdc
--- /dev/null
+++ b/src/Matlab/mex_compile/regularisers_GPU/NonlDiff_GPU.cpp
@@ -0,0 +1,92 @@
+/*
+ * This work is part of the Core Imaging Library developed by
+ * Visual Analytics and Imaging System Group of the Science Technology
+ * Facilities Council, STFC
+ *
+ * Copyright 2017 Daniil Kazantsev
+ * Copyright 2017 Srikanth Nagella, Edoardo Pasca
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ * http://www.apache.org/licenses/LICENSE-2.0
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+#include "matrix.h"
+#include "mex.h"
+#include <stdio.h>
+#include <string.h>
+#include "NonlDiff_GPU_core.h"
+
+/* CUDA implementation of linear and nonlinear diffusion with the regularisation model [1,2] (2D/3D case)
+ * The minimisation is performed using explicit scheme. 
+ *
+ * Input Parameters:
+ * 1. Noisy image/volume 
+ * 2. lambda - regularization parameter
+ * 3. Edge-preserving parameter (sigma), when sigma equals to zero nonlinear diffusion -> linear diffusion
+ * 4. Number of iterations, for explicit scheme >= 150 is recommended 
+ * 5. tau - time-marching step for explicit scheme
+ * 6. Penalty type: 1 - Huber, 2 - Perona-Malik, 3 - Tukey Biweight
+ *
+ * Output:
+ * [1] Regularized image/volume 
+ *
+ * This function is based on the paper by
+ * [1] Perona, P. and Malik, J., 1990. Scale-space and edge detection using anisotropic diffusion. IEEE Transactions on pattern analysis and machine intelligence, 12(7), pp.629-639.
+ * [2] Black, M.J., Sapiro, G., Marimont, D.H. and Heeger, D., 1998. Robust anisotropic diffusion. IEEE Transactions on image processing, 7(3), pp.421-432.
+ */
+
+void mexFunction(
+        int nlhs, mxArray *plhs[],
+        int nrhs, const mxArray *prhs[])
+        
+{
+    int number_of_dims, iter_numb, penaltytype;
+    mwSize dimX, dimY, dimZ;
+    const mwSize *dim_array;
+    
+    float *Input, *Output=NULL, lambda, tau, sigma;
+    
+    dim_array = mxGetDimensions(prhs[0]);
+    number_of_dims = mxGetNumberOfDimensions(prhs[0]);
+    
+    /*Handling Matlab input data*/
+    Input  = (float *) mxGetData(prhs[0]);
+    lambda =  (float) mxGetScalar(prhs[1]); /* regularization parameter */
+    sigma = (float) mxGetScalar(prhs[2]); /* Edge-preserving parameter */
+    iter_numb = 300; /* iterations number */
+    tau = 0.025; /* marching step parameter */
+    penaltytype = 1; /* Huber penalty by default */
+    
+    if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); }
+    if ((nrhs < 3) || (nrhs > 6)) mexErrMsgTxt("At least 3 parameters is required, all parameters are: Image(2D/3D), Regularisation parameter, Edge-preserving parameter, iterations number, time-marching constant, penalty type - Huber, PM or Tukey");
+    if ((nrhs == 4) || (nrhs == 5) || (nrhs == 6))  iter_numb = (int) mxGetScalar(prhs[3]); /* iterations number */
+    if ((nrhs == 5) || (nrhs == 6))  tau =  (float) mxGetScalar(prhs[4]); /* marching step parameter */
+    if (nrhs == 6)  {
+        char *penalty_type;
+        penalty_type = mxArrayToString(prhs[5]); /* Huber, PM or Tukey 'Huber' is the default */
+        if ((strcmp(penalty_type, "Huber") != 0) && (strcmp(penalty_type, "PM") != 0) && (strcmp(penalty_type, "Tukey") != 0)) mexErrMsgTxt("Choose penalty: 'Huber', 'PM' or 'Tukey',");
+        if (strcmp(penalty_type, "Huber") == 0)  penaltytype = 1;  /* enable 'Huber' penalty */
+        if (strcmp(penalty_type, "PM") == 0)  penaltytype = 2;  /* enable Perona-Malik penalty */
+        if (strcmp(penalty_type, "Tukey") == 0)  penaltytype = 3;  /* enable Tikey Biweight penalty */
+        mxFree(penalty_type);
+    }    
+    
+    /*Handling Matlab output data*/
+    dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2];
+    
+    /* output arrays*/
+    if (number_of_dims == 2) {
+        dimZ = 1; /*2D case*/
+        /* output image/volume */
+        Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
+    }
+    if (number_of_dims == 3) Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
+    
+    NonlDiff_GPU_main(Input, Output, lambda, sigma, iter_numb, tau, penaltytype, dimX, dimY, dimZ);
+}
\ No newline at end of file
diff --git a/src/Matlab/mex_compile/regularisers_GPU/ROF_TV_GPU.cpp b/src/Matlab/mex_compile/regularisers_GPU/ROF_TV_GPU.cpp
new file mode 100644
index 0000000..bd01d55
--- /dev/null
+++ b/src/Matlab/mex_compile/regularisers_GPU/ROF_TV_GPU.cpp
@@ -0,0 +1,74 @@
+/*
+ * This work is part of the Core Imaging Library developed by
+ * Visual Analytics and Imaging System Group of the Science Technology
+ * Facilities Council, STFC
+ *
+ * Copyright 2017 Daniil Kazantsev
+ * Copyright 2017 Srikanth Nagella, Edoardo Pasca
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ * http://www.apache.org/licenses/LICENSE-2.0
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+#include "matrix.h"
+#include "mex.h"
+#include "TV_ROF_GPU_core.h"
+
+/* ROF-TV denoising/regularization model [1] (2D/3D case)
+ * (MEX wrapper for MATLAB)
+ * 
+ * Input Parameters:
+ * 1. Noisy image/volume [REQUIRED]
+ * 2. lambda - regularization parameter [REQUIRED]
+ * 3. Number of iterations, for explicit scheme >= 150 is recommended  [REQUIRED]
+ * 4. tau - marching step for explicit scheme, ~1 is recommended [REQUIRED]
+ *
+ * Output:
+ * [1] Regularized image/volume 
+ *
+ * This function is based on the paper by
+ * [1] Rudin, Osher, Fatemi, "Nonlinear Total Variation based noise removal algorithms"
+ *
+ * D. Kazantsev, 2016-18
+ */
+void mexFunction(
+        int nlhs, mxArray *plhs[],
+        int nrhs, const mxArray *prhs[])
+        
+{
+    int number_of_dims, iter_numb;
+    mwSize dimX, dimY, dimZ;
+    const mwSize *dim_array;
+    
+    float *Input, *Output=NULL, lambda, tau;
+    
+    dim_array = mxGetDimensions(prhs[0]);
+    number_of_dims = mxGetNumberOfDimensions(prhs[0]);
+    
+    /*Handling Matlab input data*/
+    Input  = (float *) mxGetData(prhs[0]);
+    lambda =  (float) mxGetScalar(prhs[1]); /* regularization parameter */
+    iter_numb =  (int) mxGetScalar(prhs[2]); /* iterations number */
+    tau =  (float) mxGetScalar(prhs[3]); /* marching step parameter */  
+    
+    if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); }
+    if(nrhs != 4) mexErrMsgTxt("Four inputs reqired: Image(2D,3D), regularization parameter, iterations number,  marching step constant");
+    /*Handling Matlab output data*/
+    dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2];
+    
+    /* output arrays*/
+    if (number_of_dims == 2) {
+        dimZ = 1; /*2D case*/
+        /* output image/volume */
+        Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));                        
+    }    
+    if (number_of_dims == 3) Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
+    
+    TV_ROF_GPU_main(Input, Output, lambda, iter_numb, tau, dimX, dimY, dimZ);    
+}
\ No newline at end of file
diff --git a/src/Matlab/mex_compile/regularisers_GPU/SB_TV_GPU.cpp b/src/Matlab/mex_compile/regularisers_GPU/SB_TV_GPU.cpp
new file mode 100644
index 0000000..9d1328f
--- /dev/null
+++ b/src/Matlab/mex_compile/regularisers_GPU/SB_TV_GPU.cpp
@@ -0,0 +1,91 @@
+/*
+ * This work is part of the Core Imaging Library developed by
+ * Visual Analytics and Imaging System Group of the Science Technology
+ * Facilities Council, STFC
+ *
+ * Copyright 2017 Daniil Kazantsev
+ * Copyright 2017 Srikanth Nagella, Edoardo Pasca
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ * http://www.apache.org/licenses/LICENSE-2.0
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+#include "matrix.h"
+#include "mex.h"
+#include "TV_SB_GPU_core.h"
+
+/* CUDA mex-file for implementation of Split Bregman - TV denoising-regularisation model (2D/3D) [1]
+*
+* Input Parameters:
+* 1. Noisy image/volume
+* 2. lambda - regularisation parameter
+* 3. Number of iterations [OPTIONAL parameter]
+* 4. eplsilon - tolerance constant [OPTIONAL parameter]
+* 5. TV-type: 'iso' or 'l1' [OPTIONAL parameter]
+* 6. print information: 0 (off) or 1 (on)  [OPTIONAL parameter]
+*
+* Output:
+* 1. Filtered/regularized image
+*
+* This function is based on the Matlab's code and paper by
+* [1]. Goldstein, T. and Osher, S., 2009. The split Bregman method for L1-regularized problems. SIAM journal on imaging sciences, 2(2), pp.323-343.
+*/
+
+void mexFunction(
+        int nlhs, mxArray *plhs[],
+        int nrhs, const mxArray *prhs[])
+        
+{
+    int number_of_dims, iter, methTV, printswitch;
+    mwSize dimX, dimY, dimZ;
+    const mwSize *dim_array;
+    
+    float *Input, *Output=NULL, lambda, epsil;
+    
+    number_of_dims = mxGetNumberOfDimensions(prhs[0]);
+    dim_array = mxGetDimensions(prhs[0]);
+    
+    /*Handling Matlab input data*/
+    if ((nrhs < 2) || (nrhs > 6)) mexErrMsgTxt("At least 2 parameters is required, all parameters are: Image(2D/3D), Regularization parameter, Regularization parameter, iterations number, tolerance, penalty type ('iso' or 'l1'), print switch");
+    
+    Input  = (float *) mxGetData(prhs[0]); /*noisy image (2D/3D) */
+    lambda =  (float) mxGetScalar(prhs[1]); /* regularization parameter */
+    iter = 100; /* default iterations number */
+    epsil = 0.0001; /* default tolerance constant */
+    methTV = 0;  /* default isotropic TV penalty */
+    printswitch = 0; /*default print is switched, off - 0 */
+    
+    if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); }
+    
+    if ((nrhs == 3) || (nrhs == 4) || (nrhs == 5) || (nrhs == 6))  iter = (int) mxGetScalar(prhs[2]); /* iterations number */
+    if ((nrhs == 4) || (nrhs == 5) || (nrhs == 6))  epsil =  (float) mxGetScalar(prhs[3]); /* tolerance constant */
+    if ((nrhs == 5) || (nrhs == 6))  {
+        char *penalty_type;
+        penalty_type = mxArrayToString(prhs[4]); /* choosing TV penalty: 'iso' or 'l1', 'iso' is the default */
+        if ((strcmp(penalty_type, "l1") != 0) && (strcmp(penalty_type, "iso") != 0)) mexErrMsgTxt("Choose TV type: 'iso' or 'l1',");
+        if (strcmp(penalty_type, "l1") == 0)  methTV = 1;  /* enable 'l1' penalty */
+        mxFree(penalty_type);
+    }
+    if (nrhs == 6)  {
+        printswitch = (int) mxGetScalar(prhs[5]);
+        if ((printswitch != 0) && (printswitch != 1)) mexErrMsgTxt("Print can be enabled by choosing 1 or off - 0");
+    }
+    
+    /*Handling Matlab output data*/
+    dimX = dim_array[0]; dimY = dim_array[1]; dimZ = dim_array[2];
+    
+    if (number_of_dims == 2) {
+        dimZ = 1; /*2D case*/
+        Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
+    }
+    if (number_of_dims == 3) Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));
+    
+    /* running the function */
+    TV_SB_GPU_main(Input, Output, lambda, iter, epsil, methTV, printswitch, dimX, dimY, dimZ);
+}
diff --git a/src/Matlab/mex_compile/regularisers_GPU/TGV_GPU.cpp b/src/Matlab/mex_compile/regularisers_GPU/TGV_GPU.cpp
new file mode 100644
index 0000000..edb551d
--- /dev/null
+++ b/src/Matlab/mex_compile/regularisers_GPU/TGV_GPU.cpp
@@ -0,0 +1,79 @@
+/*
+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 "mex.h"
+#include "TGV_GPU_core.h"
+
+/* CUDA implementation of Primal-Dual denoising method for 
+ * Total Generilized Variation (TGV)-L2 model [1] (2D case only)
+ *
+ * Input Parameters:
+ * 1. Noisy image (2D) (required)
+ * 2. lambda - regularisation parameter (required)
+ * 3. parameter to control the first-order term (alpha1) (default - 1)
+ * 4. parameter to control the second-order term (alpha0) (default - 0.5)
+ * 5. Number of Chambolle-Pock (Primal-Dual) iterations (default is 300)
+ * 6. Lipshitz constant (default is 12)
+ *
+ * Output:
+ * Filtered/regulariaed image 
+ *
+ * References:
+ * [1] K. Bredies "Total Generalized Variation"
+ */
+
+void mexFunction(
+        int nlhs, mxArray *plhs[],
+        int nrhs, const mxArray *prhs[])
+        
+{
+    int number_of_dims, iter;
+    mwSize dimX, dimY;
+    const mwSize *dim_array;
+    float *Input, *Output=NULL, lambda, alpha0, alpha1, L2;
+    
+    number_of_dims = mxGetNumberOfDimensions(prhs[0]);
+    dim_array = mxGetDimensions(prhs[0]);
+    
+    /*Handling Matlab input data*/
+    if ((nrhs < 2) || (nrhs > 6)) mexErrMsgTxt("At least 2 parameters is required, all parameters are: Image(2D), Regularisation parameter, alpha0, alpha1, iterations number, Lipshitz Constant");
+    
+    Input  = (float *) mxGetData(prhs[0]); /*noisy image (2D) */
+    lambda =  (float) mxGetScalar(prhs[1]); /* regularisation parameter */
+    alpha1 =  1.0f; /* parameter to control the first-order term */ 
+    alpha0 =  0.5f; /* parameter to control the second-order term */
+    iter =  300; /* Iterations number */      
+    L2 =  12.0f; /* Lipshitz constant */
+    
+    if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) {mexErrMsgTxt("The input image must be in a single precision"); }   
+    if ((nrhs == 3) || (nrhs == 4) || (nrhs == 5) || (nrhs == 6))  alpha1 =  (float) mxGetScalar(prhs[2]); /* parameter to control the first-order term */ 
+    if ((nrhs == 4) || (nrhs == 5) || (nrhs == 6))  alpha0 =  (float) mxGetScalar(prhs[3]);  /* parameter to control the second-order term */
+    if ((nrhs == 5) || (nrhs == 6))  iter =  (int) mxGetScalar(prhs[4]); /* Iterations number */      
+    if (nrhs == 6)  L2 =  (float) mxGetScalar(prhs[5]); /* Lipshitz constant */
+    
+    /*Handling Matlab output data*/
+    dimX = dim_array[0]; dimY = dim_array[1];
+    
+    if (number_of_dims == 2) {
+        Output = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));
+        /* running the function */
+        TGV_GPU_main(Input, Output, lambda, alpha1, alpha0, iter, L2, dimX, dimY);        
+    }
+    if (number_of_dims == 3) {mexErrMsgTxt("Only 2D images accepted");}       
+}
diff --git a/src/Matlab/supp/RMSE.m b/src/Matlab/supp/RMSE.m
new file mode 100644
index 0000000..002f776
--- /dev/null
+++ b/src/Matlab/supp/RMSE.m
@@ -0,0 +1,7 @@
+function err = RMSE(signal1, signal2)

+%RMSE Root Mean Squared Error

+

+err = sum((signal1 - signal2).^2)/length(signal1);  % MSE

+err = sqrt(err);                                    % RMSE

+

+end
\ No newline at end of file
diff --git a/src/Matlab/supp/my_red_yellowMAP.mat b/src/Matlab/supp/my_red_yellowMAP.mat
new file mode 100644
index 0000000..c2a5b87
--- /dev/null
+++ b/src/Matlab/supp/my_red_yellowMAP.mat