Merge pull request #104 from vais-ral/newdirstructure

New directory structure, Merged other changes. The build script checks old and new structure.

29 files changed, 0 insertions, 6559 deletions

diff --git a/Wrappers/Python/CMakeLists.txt b/Wrappers/Python/CMakeLists.txt
deleted file mode 100644
index c2ef855..0000000
--- a/Wrappers/Python/CMakeLists.txt
+++ /dev/null
@@ -1,141 +0,0 @@
-#   Copyright 2018 Edoardo Pasca
-cmake_minimum_required (VERSION 3.0)
-
-project(regulariserPython)
-#https://stackoverflow.com/questions/13298504/using-cmake-with-setup-py
-
-# The version number.
-
-#set (CIL_VERSION $ENV{CIL_VERSION} CACHE INTERNAL "Core Imaging Library version" FORCE)
-
-# conda orchestrated build
-message("CIL_VERSION: ${CIL_VERSION}")
-#include (GenerateExportHeader)
-
-find_package(PythonInterp REQUIRED)
-if (PYTHONINTERP_FOUND)
-  message ("Current Python " ${PYTHON_VERSION_STRING} " found " ${PYTHON_EXECUTABLE})
-endif()
-
-	
-## Build the regularisers package as a library
-message("Creating Regularisers as shared library")
-
-message("CMAKE_CXX_FLAGS ${CMAKE_CXX_FLAGS}")
-
-set(CMAKE_BUILD_TYPE "Release")
-
-if(WIN32)
-  set (FLAGS "/DWIN32 /EHsc /openmp /DCCPiCore_EXPORTS")
-  set (CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} /NODEFAULTLIB:MSVCRT.lib")
-  
-  set (EXTRA_LIBRARIES)
-		
-  message("library lib: ${LIBRARY_LIB}")
-  
-elseif(UNIX)
-   set (FLAGS "-fopenmp -O2 -funsigned-char -Wall  -Wl,--no-undefined  -DCCPiReconstructionIterative_EXPORTS -std=c++0x")  
-   set (EXTRA_LIBRARIES 
-		"gomp"
-		)
-endif()
-
-# GPU regularisers
-if (BUILD_CUDA)
-    find_package(CUDA)
-    if (CUDA_FOUND)
-      message("CUDA FOUND")
-      set (SETUP_GPU_WRAPPERS "extra_libraries += ['cilregcuda']\n\
-setup( \n\
-        name='ccpi', \n\
-        description='CCPi Core Imaging Library - Image regularisers GPU',\n\
-        version=cil_version,\n\
-        cmdclass = {'build_ext': build_ext},\n\
-        ext_modules = [Extension('ccpi.filters.gpu_regularisers',\n\
-                                  sources=[ \n\
-                                          os.path.join('.' , 'src', 'gpu_regularisers.pyx' ),\n\
-                                            ],\n\
-                                 include_dirs=extra_include_dirs, \n\
-                                 library_dirs=extra_library_dirs, \n\
-                                 extra_compile_args=extra_compile_args, \n\
-                                 libraries=extra_libraries ), \n\
-        ],\n\
-        zip_safe = False,	\n\
-        packages = {'ccpi','ccpi.filters'},\n\
-    )")
-    else()
-      message("CUDA NOT FOUND")
-      set(SETUP_GPU_WRAPPERS "#CUDA NOT FOUND")
-    endif()
-endif()
-configure_file("${CMAKE_CURRENT_SOURCE_DIR}/setup-regularisers.py.in" "${CMAKE_CURRENT_BINARY_DIR}/setup-regularisers.py")
-
-
-find_package(PythonInterp)
-find_package(PythonLibs)
-if (PYTHONINTERP_FOUND)
-  message(STATUS "Found PYTHON_EXECUTABLE=${PYTHON_EXECUTABLE}")
-  message(STATUS "Python version ${PYTHON_VERSION_STRING}")
-endif()
-if (PYTHONLIBS_FOUND)
-  message(STATUS "Found PYTHON_INCLUDE_DIRS=${PYTHON_INCLUDE_DIRS}")
-  message(STATUS "Found PYTHON_LIBRARIES=${PYTHON_LIBRARIES}")
-endif()
-
-if (PYTHONINTERP_FOUND)
-    message("Python found " ${PYTHON_EXECUTABLE})
-    set(SETUP_PY_IN "${CMAKE_CURRENT_SOURCE_DIR}/setup-regularisers.py.in")
-    set(SETUP_PY    "${CMAKE_CURRENT_BINARY_DIR}/setup-regularisers.py")
-    #set(DEPS        "${CMAKE_CURRENT_SOURCE_DIR}/module/__init__.py")
-    set (DEPS       "${CMAKE_BINARY_DIR}/Core/")
-    set(OUTPUT      "${CMAKE_CURRENT_BINARY_DIR}/build/timestamp")
-
-    configure_file(${SETUP_PY_IN} ${SETUP_PY})
-
-    message("Core binary dir " ${CMAKE_BINARY_DIR}/Core/${CMAKE_BUILD_TYPE})
-    
-    if (CONDA_BUILD)
-      add_custom_command(OUTPUT ${OUTPUT}
-                       COMMAND ${CMAKE_COMMAND} -E copy_directory ${CMAKE_CURRENT_SOURCE_DIR}/src ${CMAKE_CURRENT_BINARY_DIR}/src
-                       COMMAND ${CMAKE_COMMAND} -E copy_directory ${CMAKE_CURRENT_SOURCE_DIR}/ccpi ${CMAKE_CURRENT_BINARY_DIR}/ccpi
-                       COMMAND ${CMAKE_COMMAND} -E env CIL_VERSION=${CIL_VERSION}
-                                                       PREFIX=${CMAKE_SOURCE_DIR}/Core 
-                                                       LIBRARY_INC=${CMAKE_SOURCE_DIR}/Core 
-                                                       LIBRARY_LIB=${CMAKE_BINARY_DIR}/Core
-                                                       ${PYTHON_EXECUTABLE} ${SETUP_PY} install
-                       COMMAND ${CMAKE_COMMAND} -E touch ${OUTPUT}
-                       DEPENDS cilreg)
-
-    else()
-      if (WIN32)
-        add_custom_command(OUTPUT ${OUTPUT}
-                       COMMAND ${CMAKE_COMMAND} -E copy_directory ${CMAKE_CURRENT_SOURCE_DIR}/src ${CMAKE_CURRENT_BINARY_DIR}/src
-                       COMMAND ${CMAKE_COMMAND} -E copy_directory ${CMAKE_CURRENT_SOURCE_DIR}/ccpi ${CMAKE_CURRENT_BINARY_DIR}/ccpi
-                       COMMAND ${CMAKE_COMMAND} -E env CIL_VERSION=${CIL_VERSION}
-                                                       PREFIX=${CMAKE_SOURCE_DIR}/Core 
-                                                       LIBRARY_INC=${CMAKE_SOURCE_DIR}/Core 
-                                                       LIBRARY_LIB=${CMAKE_BINARY_DIR}/Core/${CMAKE_BUILD_TYPE}
-                                                       ${PYTHON_EXECUTABLE} ${SETUP_PY} build_ext --inplace
-                       COMMAND ${CMAKE_COMMAND} -E touch ${OUTPUT}
-                       DEPENDS cilreg)
-      else()
-        add_custom_command(OUTPUT ${OUTPUT}
-                       COMMAND ${CMAKE_COMMAND} -E copy_directory ${CMAKE_CURRENT_SOURCE_DIR}/src ${CMAKE_CURRENT_BINARY_DIR}/src
-                       COMMAND ${CMAKE_COMMAND} -E copy_directory ${CMAKE_CURRENT_SOURCE_DIR}/ccpi ${CMAKE_CURRENT_BINARY_DIR}/ccpi
-                       COMMAND ${CMAKE_COMMAND} -E env CIL_VERSION=${CIL_VERSION}
-                                                       PREFIX=${CMAKE_SOURCE_DIR}/Core 
-                                                       LIBRARY_INC=${CMAKE_SOURCE_DIR}/Core 
-                                                       LIBRARY_LIB=${CMAKE_BINARY_DIR}/Core
-                                                       ${PYTHON_EXECUTABLE} ${SETUP_PY} build_ext --inplace
-                       COMMAND ${CMAKE_COMMAND} -E touch ${OUTPUT}
-                       DEPENDS cilreg)
-      endif()
-      install(DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/ccpi 
-              DESTINATION ${PYTHON_DEST})
-    endif()
-    
-    
-    add_custom_target(PythonWrapper ALL DEPENDS ${OUTPUT})
-
-    #install(CODE "execute_process(COMMAND ${PYTHON} ${SETUP_PY} install)")
-endif()
diff --git a/Wrappers/Python/ccpi/__init__.py b/Wrappers/Python/ccpi/__init__.py
deleted file mode 100644
index e69de29..0000000
--- a/Wrappers/Python/ccpi/__init__.py
+++ /dev/null
diff --git a/Wrappers/Python/ccpi/filters/__init__.py b/Wrappers/Python/ccpi/filters/__init__.py
deleted file mode 100644
index e69de29..0000000
--- a/Wrappers/Python/ccpi/filters/__init__.py
+++ /dev/null
diff --git a/Wrappers/Python/ccpi/filters/regularisers.py b/Wrappers/Python/ccpi/filters/regularisers.py
deleted file mode 100644
index 588ea32..0000000
--- a/Wrappers/Python/ccpi/filters/regularisers.py
+++ /dev/null
@@ -1,214 +0,0 @@
-"""
-script which assigns a proper device core function based on a flag ('cpu' or 'gpu')
-"""
-
-from ccpi.filters.cpu_regularisers import TV_ROF_CPU, TV_FGP_CPU, TV_SB_CPU, dTV_FGP_CPU, TNV_CPU, NDF_CPU, Diff4th_CPU, TGV_CPU, LLT_ROF_CPU, PATCHSEL_CPU, NLTV_CPU
-try:
-    from ccpi.filters.gpu_regularisers import TV_ROF_GPU, TV_FGP_GPU, TV_SB_GPU, dTV_FGP_GPU, NDF_GPU, Diff4th_GPU, TGV_GPU, LLT_ROF_GPU, PATCHSEL_GPU
-    gpu_enabled = True
-except ImportError:
-    gpu_enabled = False    
-from ccpi.filters.cpu_regularisers import NDF_INPAINT_CPU, NVM_INPAINT_CPU
-
-def ROF_TV(inputData, regularisation_parameter, iterations,
-                     time_marching_parameter,device='cpu'):
-    if device == 'cpu':
-        return TV_ROF_CPU(inputData,
-                     regularisation_parameter,
-                     iterations, 
-                     time_marching_parameter)
-    elif device == 'gpu' and gpu_enabled:
-        return TV_ROF_GPU(inputData,
-                     regularisation_parameter,
-                     iterations, 
-                     time_marching_parameter)
-    else:
-        if not gpu_enabled and device == 'gpu':
-            raise ValueError ('GPU is not available')
-        raise ValueError('Unknown device {0}. Expecting gpu or cpu'\
-                         .format(device))
-
-def FGP_TV(inputData, regularisation_parameter,iterations,
-                     tolerance_param, methodTV, nonneg, printM, device='cpu'):
-    if device == 'cpu':
-        return TV_FGP_CPU(inputData,
-                     regularisation_parameter,
-                     iterations, 
-                     tolerance_param,
-                     methodTV,
-                     nonneg,
-                     printM)
-    elif device == 'gpu' and gpu_enabled:
-        return TV_FGP_GPU(inputData,
-                     regularisation_parameter,
-                     iterations, 
-                     tolerance_param,
-                     methodTV,
-                     nonneg,
-                     printM)
-    else:
-        if not gpu_enabled and device == 'gpu':
-            raise ValueError ('GPU is not available')
-        raise ValueError('Unknown device {0}. Expecting gpu or cpu'\
-                         .format(device))
-def SB_TV(inputData, regularisation_parameter, iterations,
-                     tolerance_param, methodTV, printM, device='cpu'):
-    if device == 'cpu':
-        return TV_SB_CPU(inputData,
-                     regularisation_parameter,
-                     iterations, 
-                     tolerance_param,
-                     methodTV,
-                     printM)
-    elif device == 'gpu' and gpu_enabled:
-        return TV_SB_GPU(inputData,
-                     regularisation_parameter,
-                     iterations, 
-                     tolerance_param,
-                     methodTV,
-                     printM)
-    else:
-        if not gpu_enabled and device == 'gpu':
-            raise ValueError ('GPU is not available')
-        raise ValueError('Unknown device {0}. Expecting gpu or cpu'\
-                         .format(device))
-def FGP_dTV(inputData, refdata, regularisation_parameter, iterations,
-                     tolerance_param, eta_const, methodTV, nonneg, printM, device='cpu'):
-    if device == 'cpu':
-        return dTV_FGP_CPU(inputData,
-                     refdata,
-                     regularisation_parameter,
-                     iterations, 
-                     tolerance_param,
-                     eta_const,
-                     methodTV,
-                     nonneg,
-                     printM)
-    elif device == 'gpu' and gpu_enabled:
-        return dTV_FGP_GPU(inputData,
-                     refdata,
-                     regularisation_parameter,
-                     iterations, 
-                     tolerance_param,
-                     eta_const,
-                     methodTV,
-                     nonneg,
-                     printM)
-    else:
-        if not gpu_enabled and device == 'gpu':
-            raise ValueError ('GPU is not available')
-        raise ValueError('Unknown device {0}. Expecting gpu or cpu'\
-                         .format(device))
-def TNV(inputData, regularisation_parameter, iterations, tolerance_param):
-        return TNV_CPU(inputData,
-                     regularisation_parameter,
-                     iterations, 
-                     tolerance_param)
-def NDF(inputData, regularisation_parameter, edge_parameter, iterations,
-                     time_marching_parameter, penalty_type, device='cpu'):
-    if device == 'cpu':
-        return NDF_CPU(inputData,
-                     regularisation_parameter,
-                     edge_parameter,
-                     iterations, 
-                     time_marching_parameter,
-                     penalty_type)
-    elif device == 'gpu' and gpu_enabled:
-        return NDF_GPU(inputData,
-                     regularisation_parameter,
-                     edge_parameter,
-                     iterations, 
-                     time_marching_parameter,
-                     penalty_type)
-    else:
-        if not gpu_enabled and device == 'gpu':
-    	    raise ValueError ('GPU is not available')
-        raise ValueError('Unknown device {0}. Expecting gpu or cpu'\
-                         .format(device))
-def Diff4th(inputData, regularisation_parameter, edge_parameter, iterations,
-                     time_marching_parameter, device='cpu'):
-    if device == 'cpu':
-        return Diff4th_CPU(inputData,
-                     regularisation_parameter,
-                     edge_parameter,
-                     iterations, 
-                     time_marching_parameter)
-    elif device == 'gpu' and gpu_enabled:
-        return Diff4th_GPU(inputData,
-                     regularisation_parameter,
-                     edge_parameter,
-                     iterations, 
-                     time_marching_parameter)
-    else:
-        if not gpu_enabled and device == 'gpu':
-            raise ValueError ('GPU is not available')
-        raise ValueError('Unknown device {0}. Expecting gpu or cpu'\
-                         .format(device))
-        
-def PatchSelect(inputData, searchwindow, patchwindow, neighbours, edge_parameter, device='cpu'):
-    if device == 'cpu':
-        return PATCHSEL_CPU(inputData,
-                     searchwindow,
-                     patchwindow,
-                     neighbours, 
-                     edge_parameter)
-    elif device == 'gpu' and gpu_enabled:
-        return PATCHSEL_GPU(inputData,
-                     searchwindow,
-                     patchwindow,
-                     neighbours, 
-                     edge_parameter)
-    else:
-        if not gpu_enabled and device == 'gpu':
-            raise ValueError ('GPU is not available')
-        raise ValueError('Unknown device {0}. Expecting gpu or cpu'\
-                         .format(device))
-
-def NLTV(inputData, H_i, H_j, H_k, Weights, regularisation_parameter, iterations):
-    return NLTV_CPU(inputData,
-                     H_i,
-                     H_j,
-                     H_k, 
-                     Weights,
-                     regularisation_parameter,
-                     iterations)
-
-def TGV(inputData, regularisation_parameter, alpha1, alpha0, iterations,
-                     LipshitzConst, device='cpu'):
-    if device == 'cpu':
-        return TGV_CPU(inputData, 
-					regularisation_parameter, 
-					alpha1, 
-					alpha0, 
-					iterations,
-                    LipshitzConst)
-    elif device == 'gpu' and gpu_enabled:
-        return TGV_GPU(inputData, 
-					regularisation_parameter, 
-					alpha1, 
-					alpha0, 
-					iterations,
-                    LipshitzConst)
-    else:
-        if not gpu_enabled and device == 'gpu':
-            raise ValueError ('GPU is not available')
-        raise ValueError('Unknown device {0}. Expecting gpu or cpu'\
-                         .format(device))
-def LLT_ROF(inputData, regularisation_parameterROF, regularisation_parameterLLT, iterations,
-                     time_marching_parameter, device='cpu'):
-    if device == 'cpu':
-        return LLT_ROF_CPU(inputData, regularisation_parameterROF, regularisation_parameterLLT, iterations, time_marching_parameter)
-    elif device == 'gpu' and gpu_enabled:
-        return LLT_ROF_GPU(inputData, regularisation_parameterROF, regularisation_parameterLLT, iterations, time_marching_parameter)
-    else:
-        if not gpu_enabled and device == 'gpu':
-            raise ValueError ('GPU is not available')
-        raise ValueError('Unknown device {0}. Expecting gpu or cpu'\
-                         .format(device))
-def NDF_INP(inputData, maskData, regularisation_parameter, edge_parameter, iterations,
-                     time_marching_parameter, penalty_type):
-        return NDF_INPAINT_CPU(inputData, maskData, regularisation_parameter, 
-        edge_parameter, iterations, time_marching_parameter, penalty_type)
-        
-def NVM_INP(inputData, maskData, SW_increment, iterations):
-        return NVM_INPAINT_CPU(inputData, maskData, SW_increment, iterations)
diff --git a/Wrappers/Python/ccpi/supp/__init__.py b/Wrappers/Python/ccpi/supp/__init__.py
deleted file mode 100644
index e69de29..0000000
--- a/Wrappers/Python/ccpi/supp/__init__.py
+++ /dev/null
diff --git a/Wrappers/Python/ccpi/supp/qualitymetrics.py b/Wrappers/Python/ccpi/supp/qualitymetrics.py
deleted file mode 100644
index f44d832..0000000
--- a/Wrappers/Python/ccpi/supp/qualitymetrics.py
+++ /dev/null
@@ -1,65 +0,0 @@
-#!/usr/bin/env python2
-# -*- coding: utf-8 -*-
-"""
-A class for some standard image quality metrics
-"""
-import numpy as np
-
-class QualityTools:
-    def __init__(self, im1, im2):
-        if im1.size != im2.size:
-            print ('Error: Sizes of images/volumes are different')
-            raise SystemExit
-        self.im1 = im1 # image or volume - 1
-        self.im2 = im2 # image or volume - 2
-    def nrmse(self):
-        """ Normalised Root Mean Square Error """
-        rmse = np.sqrt(np.sum((self.im2 - self.im1) ** 2) / float(self.im1.size))
-        max_val = max(np.max(self.im1), np.max(self.im2))
-        min_val = min(np.min(self.im1), np.min(self.im2))
-        return 1 - (rmse / (max_val - min_val))
-    def rmse(self):
-        """ Root Mean Square Error """
-        rmse = np.sqrt(np.sum((self.im1 - self.im2) ** 2) / float(self.im1.size))
-        return rmse
-    def ssim(self, window, k=(0.01, 0.03), l=255):
-        from scipy.signal import fftconvolve
-        """See https://ece.uwaterloo.ca/~z70wang/research/ssim/"""
-        # Check if the window is smaller than the images.
-        for a, b in zip(window.shape, self.im1.shape):
-            if a > b:
-                return None, None
-        # Values in k must be positive according to the base implementation.
-        for ki in k:
-            if ki < 0:
-                return None, None
-    
-        c1 = (k[0] * l) ** 2
-        c2 = (k[1] * l) ** 2
-        window = window/np.sum(window)
-    
-        mu1 = fftconvolve(self.im1, window, mode='valid')
-        mu2 = fftconvolve(self.im2, window, mode='valid')
-        mu1_sq = mu1 * mu1
-        mu2_sq = mu2 * mu2
-        mu1_mu2 = mu1 * mu2
-        sigma1_sq = fftconvolve(self.im1 * self.im1, window, mode='valid') - mu1_sq
-        sigma2_sq = fftconvolve(self.im2 * self.im2, window, mode='valid') - mu2_sq
-        sigma12 = fftconvolve(self.im1 * self.im2, window, mode='valid') - mu1_mu2
-    
-        if c1 > 0 and c2 > 0:
-            num = (2 * mu1_mu2 + c1) * (2 * sigma12 + c2)
-            den = (mu1_sq + mu2_sq + c1) * (sigma1_sq + sigma2_sq + c2)
-            ssim_map = num / den
-        else:
-            num1 = 2 * mu1_mu2 + c1
-            num2 = 2 * sigma12 + c2
-            den1 = mu1_sq + mu2_sq + c1
-            den2 = sigma1_sq + sigma2_sq + c2
-            ssim_map = np.ones(np.shape(mu1))
-            index = (den1 * den2) > 0
-            ssim_map[index] = (num1[index] * num2[index]) / (den1[index] * den2[index])
-            index = (den1 != 0) & (den2 == 0)
-            ssim_map[index] = num1[index] / den1[index]    
-        mssim = ssim_map.mean()
-        return mssim, ssim_map
diff --git a/Wrappers/Python/conda-recipe/bld.bat b/Wrappers/Python/conda-recipe/bld.bat
deleted file mode 100644
index 6c84355..0000000
--- a/Wrappers/Python/conda-recipe/bld.bat
+++ /dev/null
@@ -1,20 +0,0 @@
-IF NOT DEFINED CIL_VERSION (
-ECHO CIL_VERSION Not Defined.
-exit 1
-)
-
-mkdir "%SRC_DIR%\ccpi"
-ROBOCOPY /E "%RECIPE_DIR%\..\.." "%SRC_DIR%\ccpi"
-ROBOCOPY /E "%RECIPE_DIR%\..\..\..\Core" "%SRC_DIR%\Core"
-::cd %SRC_DIR%\ccpi\Python
-cd %SRC_DIR%
-
-:: issue cmake to create setup.py
-cmake -G "NMake Makefiles" %RECIPE_DIR%\..\..\..\ -DBUILD_PYTHON_WRAPPERS=ON -DCONDA_BUILD=ON -DBUILD_CUDA=OFF -DCMAKE_BUILD_TYPE="Release" -DLIBRARY_LIB="%CONDA_PREFIX%\lib" -DLIBRARY_INC="%CONDA_PREFIX%" -DCMAKE_INSTALL_PREFIX="%PREFIX%\Library" 
-
-::%PYTHON% setup-regularisers.py build_ext
-::if errorlevel 1 exit 1
-::%PYTHON% setup-regularisers.py install
-::if errorlevel 1 exit 1
-nmake install
-if errorlevel 1 exit 1
\ No newline at end of file
diff --git a/Wrappers/Python/conda-recipe/build.sh b/Wrappers/Python/conda-recipe/build.sh
deleted file mode 100644
index 39c0f2c..0000000
--- a/Wrappers/Python/conda-recipe/build.sh
+++ /dev/null
@@ -1,17 +0,0 @@
-
-mkdir "$SRC_DIR/ccpi"
-cp -rv "$RECIPE_DIR/../.." "$SRC_DIR/ccpi"
-cp -rv "$RECIPE_DIR/../../../Core" "$SRC_DIR/Core"
-
-cd $SRC_DIR
-##cuda=off
-
-cmake -G "Unix Makefiles" $RECIPE_DIR/../../../ -DBUILD_PYTHON_WRAPPER=ON -DCONDA_BUILD=ON -DBUILD_CUDA=ON -DCMAKE_BUILD_TYPE="Release" -DLIBRARY_LIB=$CONDA_PREFIX/lib -DLIBRARY_INC=$CONDA_PREFIX -DCMAKE_INSTALL_PREFIX=$PREFIX
-
-
-make install
-
-#$PYTHON setup-regularisers.py build_ext
-#$PYTHON setup-regularisers.py install
-
-
diff --git a/Wrappers/Python/conda-recipe/conda_build_config.yaml b/Wrappers/Python/conda-recipe/conda_build_config.yaml
deleted file mode 100644
index fbe82dc..0000000
--- a/Wrappers/Python/conda-recipe/conda_build_config.yaml
+++ /dev/null
@@ -1,9 +0,0 @@
-python:
-  - 2.7 # [not win]
-  - 3.5
-  - 3.6
-#  - 3.7
-numpy:
-  - 1.12
-  - 1.14
-  - 1.15
diff --git a/Wrappers/Python/conda-recipe/lena_gray_512.tif b/Wrappers/Python/conda-recipe/lena_gray_512.tif
deleted file mode 100644
index f80cafc..0000000
--- a/Wrappers/Python/conda-recipe/lena_gray_512.tif
+++ /dev/null
diff --git a/Wrappers/Python/conda-recipe/meta.yaml b/Wrappers/Python/conda-recipe/meta.yaml
deleted file mode 100644
index 7435b2b..0000000
--- a/Wrappers/Python/conda-recipe/meta.yaml
+++ /dev/null
@@ -1,40 +0,0 @@
-package:
-  name: ccpi-regulariser
-  version: {{CIL_VERSION}}
-  
-build:
-  preserve_egg_dir: False
-  number: 0
-  script_env:
-    - CIL_VERSION
-  
-test:
-  files:
-    - lena_gray_512.tif
-  requires:
-    - pillow=4.1.1
-
-requirements:
-  build:
-    - python
-    - numpy {{ numpy }}
-    - setuptools
-    - cython
-    - vc 14 # [win and py36] 
-    - vc 14 # [win and py35] 
-    - vc 9  # [win and py27]
-    - cmake 
-
-  run:
-    - {{ pin_compatible('numpy', max_pin='x.x') }}
-    - python
-    - numpy
-    - vc 14 # [win and py36] 
-    - vc 14 # [win and py35] 
-    - vc 9  # [win and py27]
-    - libgcc-ng
-
-about:
-  home: http://www.ccpi.ac.uk
-  license:  BSD license
-  summary: 'CCPi Core Imaging Library Quantification Toolbox'
diff --git a/Wrappers/Python/conda-recipe/run_test.py b/Wrappers/Python/conda-recipe/run_test.py
deleted file mode 100755
index 21f3216..0000000
--- a/Wrappers/Python/conda-recipe/run_test.py
+++ /dev/null
@@ -1,819 +0,0 @@
-import unittest

-import numpy as np

-import os

-import timeit

-from ccpi.filters.regularisers import ROF_TV, FGP_TV, SB_TV, TGV, LLT_ROF, FGP_dTV, NDF, Diff4th

-from PIL import Image

-

-class TiffReader(object):

-    def imread(self, filename):

-        return np.asarray(Image.open(filename))

-###############################################################################

-def printParametersToString(pars):

-        txt = r''

-        for key, value in pars.items():

-            if key== 'algorithm' :

-                txt += "{0} = {1}".format(key, value.__name__)

-            elif key == 'input':

-                txt += "{0} = {1}".format(key, np.shape(value))

-            elif key == 'refdata':

-                txt += "{0} = {1}".format(key, np.shape(value))

-            else:

-                txt += "{0} = {1}".format(key, value)

-            txt += '\n'

-        return txt

-def nrmse(im1, im2):

-    rmse = np.sqrt(np.sum((im2 - im1) ** 2) / float(im1.size))

-    max_val = max(np.max(im1), np.max(im2))

-    min_val = min(np.min(im1), np.min(im2))

-    return 1 - (rmse / (max_val - min_val))

-    

-def rmse(im1, im2):

-    rmse = np.sqrt(np.sum((im1 - im2) ** 2) / float(im1.size))

-    return rmse

-###############################################################################

-

-class TestRegularisers(unittest.TestCase):

-    

-

-    def test_ROF_TV_CPU_vs_GPU(self):

-        #print ("tomas debug test function")

-        print(__name__)

-        filename = os.path.join("lena_gray_512.tif")

-        plt = TiffReader()

-        # read image

-        Im = plt.imread(filename)                     

-        Im = np.asarray(Im, dtype='float32')

-        

-        Im = Im/255

-        perc = 0.05

-        u0 = Im + np.random.normal(loc = 0 ,

-                                          scale = perc * Im , 

-                                          size = np.shape(Im))

-        u_ref = Im + np.random.normal(loc = 0 ,

-                                          scale = 0.01 * Im , 

-                                          size = np.shape(Im))

-        

-        # map the u0 u0->u0>0

-        # f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1)

-        u0 = u0.astype('float32')

-        u_ref = u_ref.astype('float32')

-        

-        print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")

-        print ("____________ROF-TV bench___________________")

-        print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")

-        

-        # set parameters

-        pars = {'algorithm': ROF_TV, \

-        'input' : u0,\

-        'regularisation_parameter':0.04,\

-        'number_of_iterations': 2500,\

-        'time_marching_parameter': 0.00002

-        }

-        print ("#############ROF TV CPU####################")

-        start_time = timeit.default_timer()

-        rof_cpu = ROF_TV(pars['input'],

-                     pars['regularisation_parameter'],

-                     pars['number_of_iterations'],

-                     pars['time_marching_parameter'],'cpu')

-        rms = rmse(Im, rof_cpu)

-        pars['rmse'] = rms

-        

-        txtstr = printParametersToString(pars)

-        txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)

-        print (txtstr)

-        print ("##############ROF TV GPU##################")

-        start_time = timeit.default_timer()

-        try:

-            rof_gpu = ROF_TV(pars['input'], 

-                             pars['regularisation_parameter'],

-                             pars['number_of_iterations'], 

-                             pars['time_marching_parameter'],'gpu')

-        except ValueError as ve:

-            self.skipTest("Results not comparable. GPU computing error.")

-

-        rms = rmse(Im, rof_gpu)

-        pars['rmse'] = rms

-        pars['algorithm'] = ROF_TV

-        txtstr = printParametersToString(pars)

-        txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)

-        print (txtstr)

-        print ("--------Compare the results--------")

-        tolerance = 1e-04

-        diff_im = np.zeros(np.shape(rof_cpu))

-        diff_im = abs(rof_cpu - rof_gpu)

-        diff_im[diff_im > tolerance] = 1

-        self.assertLessEqual(diff_im.sum() , 1)

-        

-    def test_FGP_TV_CPU_vs_GPU(self):

-        print(__name__)

-        filename = os.path.join("lena_gray_512.tif")

-        plt = TiffReader()

-        # read image

-        Im = plt.imread(filename)                     

-        Im = np.asarray(Im, dtype='float32')

-        

-        Im = Im/255

-        perc = 0.05

-        u0 = Im + np.random.normal(loc = 0 ,

-                                          scale = perc * Im , 

-                                          size = np.shape(Im))

-        u_ref = Im + np.random.normal(loc = 0 ,

-                                          scale = 0.01 * Im , 

-                                          size = np.shape(Im))

-        

-        # map the u0 u0->u0>0

-        # f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1)

-        u0 = u0.astype('float32')

-        u_ref = u_ref.astype('float32')

-        

-        print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")

-        print ("____________FGP-TV bench___________________")

-        print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")

-        

-        

-        # set parameters

-        pars = {'algorithm' : FGP_TV, \

-                'input' : u0,\

-                'regularisation_parameter':0.04, \

-                'number_of_iterations' :1200 ,\

-                'tolerance_constant':0.00001,\

-                'methodTV': 0 ,\

-                'nonneg': 0 ,\

-                'printingOut': 0 

-                }

-                

-        print ("#############FGP TV CPU####################")

-        start_time = timeit.default_timer()

-        fgp_cpu = FGP_TV(pars['input'], 

-                      pars['regularisation_parameter'],

-                      pars['number_of_iterations'],

-                      pars['tolerance_constant'], 

-                      pars['methodTV'],

-                      pars['nonneg'],

-                      pars['printingOut'],'cpu')  

-                     

-                     

-        rms = rmse(Im, fgp_cpu)

-        pars['rmse'] = rms

-        

-        txtstr = printParametersToString(pars)

-        txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)

-        print (txtstr)

-        

-        print ("##############FGP TV GPU##################")

-        start_time = timeit.default_timer()

-        try:

-            fgp_gpu = FGP_TV(pars['input'], 

-                      pars['regularisation_parameter'],

-                      pars['number_of_iterations'],

-                      pars['tolerance_constant'], 

-                      pars['methodTV'],

-                      pars['nonneg'],

-                      pars['printingOut'],'gpu')

-

-        except ValueError as ve:

-            self.skipTest("Results not comparable. GPU computing error.")

-

-        rms = rmse(Im, fgp_gpu)

-        pars['rmse'] = rms

-        pars['algorithm'] = FGP_TV

-        txtstr = printParametersToString(pars)

-        txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)

-        print (txtstr)

-        

-        print ("--------Compare the results--------")

-        tolerance = 1e-05

-        diff_im = np.zeros(np.shape(fgp_cpu))

-        diff_im = abs(fgp_cpu - fgp_gpu)

-        diff_im[diff_im > tolerance] = 1

-

-        self.assertLessEqual(diff_im.sum() , 1)

-

-    def test_SB_TV_CPU_vs_GPU(self):

-        print(__name__)

-        filename = os.path.join("lena_gray_512.tif")

-        plt = TiffReader()

-        # read image

-        Im = plt.imread(filename)                     

-        Im = np.asarray(Im, dtype='float32')

-        

-        Im = Im/255

-        perc = 0.05

-        u0 = Im + np.random.normal(loc = 0 ,

-                                          scale = perc * Im , 

-                                          size = np.shape(Im))

-        u_ref = Im + np.random.normal(loc = 0 ,

-                                          scale = 0.01 * Im , 

-                                          size = np.shape(Im))

-        

-        # map the u0 u0->u0>0

-        # f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1)

-        u0 = u0.astype('float32')

-        u_ref = u_ref.astype('float32')

-        

-        print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")

-        print ("____________SB-TV bench___________________")

-        print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")

-        

-        

-        # set parameters

-        pars = {'algorithm' : SB_TV, \

-                'input' : u0,\

-                'regularisation_parameter':0.04, \

-                'number_of_iterations' :150 ,\

-                'tolerance_constant':1e-05,\

-                'methodTV': 0 ,\

-                'printingOut': 0 

-                }

-                

-        print ("#############SB-TV CPU####################")

-        start_time = timeit.default_timer()

-        sb_cpu = SB_TV(pars['input'], 

-                      pars['regularisation_parameter'],

-                      pars['number_of_iterations'],

-                      pars['tolerance_constant'], 

-                      pars['methodTV'],

-                      pars['printingOut'],'cpu')  

-                     

-                     

-        rms = rmse(Im, sb_cpu)

-        pars['rmse'] = rms

-        

-        txtstr = printParametersToString(pars)

-        txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)

-        print (txtstr)

-        

-        print ("##############SB TV GPU##################")

-        start_time = timeit.default_timer()

-        try:

-            

-            sb_gpu = SB_TV(pars['input'], 

-                      pars['regularisation_parameter'],

-                      pars['number_of_iterations'],

-                      pars['tolerance_constant'], 

-                      pars['methodTV'],

-                      pars['printingOut'],'gpu')

-

-        except ValueError as ve:

-            self.skipTest("Results not comparable. GPU computing error.")

-

-        rms = rmse(Im, sb_gpu)

-        pars['rmse'] = rms

-        pars['algorithm'] = SB_TV

-        txtstr = printParametersToString(pars)

-        txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)

-        print (txtstr)

-        print ("--------Compare the results--------")

-        tolerance = 1e-05

-        diff_im = np.zeros(np.shape(sb_cpu))

-        diff_im = abs(sb_cpu - sb_gpu)

-        diff_im[diff_im > tolerance] = 1

-        self.assertLessEqual(diff_im.sum(), 1)

-

-    def test_TGV_CPU_vs_GPU(self):

-        print(__name__)

-        filename = os.path.join("lena_gray_512.tif")

-        plt = TiffReader()

-        # read image

-        Im = plt.imread(filename)                     

-        Im = np.asarray(Im, dtype='float32')

-        

-        Im = Im/255

-        perc = 0.05

-        u0 = Im + np.random.normal(loc = 0 ,

-                                          scale = perc * Im , 

-                                          size = np.shape(Im))

-        u_ref = Im + np.random.normal(loc = 0 ,

-                                          scale = 0.01 * Im , 

-                                          size = np.shape(Im))

-        

-        # map the u0 u0->u0>0

-        # f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1)

-        u0 = u0.astype('float32')

-        u_ref = u_ref.astype('float32')

-        

-        print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")

-        print ("____________TGV bench___________________")

-        print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")

-        

-        

-        # set parameters

-        pars = {'algorithm' : TGV, \

-                'input' : u0,\

-                'regularisation_parameter':0.04, \

-                'alpha1':1.0,\

-                'alpha0':2.0,\

-                'number_of_iterations' :250 ,\

-                'LipshitzConstant' :12 ,\

-                }

-                

-        print ("#############TGV CPU####################")

-        start_time = timeit.default_timer()

-        tgv_cpu = TGV(pars['input'], 

-                      pars['regularisation_parameter'],

-                      pars['alpha1'],

-                      pars['alpha0'],

-                      pars['number_of_iterations'],

-                      pars['LipshitzConstant'],'cpu')

-                     

-        rms = rmse(Im, tgv_cpu)

-        pars['rmse'] = rms

-        

-        txtstr = printParametersToString(pars)

-        txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)

-        print (txtstr)

-        

-        print ("##############TGV GPU##################")

-        start_time = timeit.default_timer()

-        try:

-            tgv_gpu = TGV(pars['input'], 

-                      pars['regularisation_parameter'],

-                      pars['alpha1'],

-                      pars['alpha0'],

-                      pars['number_of_iterations'],

-                      pars['LipshitzConstant'],'gpu')

-

-        except ValueError as ve:

-            self.skipTest("Results not comparable. GPU computing error.")

-

-        rms = rmse(Im, tgv_gpu)

-        pars['rmse'] = rms

-        pars['algorithm'] = TGV

-        txtstr = printParametersToString(pars)

-        txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)

-        print (txtstr)

-        print ("--------Compare the results--------")

-        tolerance = 1e-05

-        diff_im = np.zeros(np.shape(tgv_gpu))

-        diff_im = abs(tgv_cpu - tgv_gpu)

-        diff_im[diff_im > tolerance] = 1

-        self.assertLessEqual(diff_im.sum() , 1)

-

-    def test_LLT_ROF_CPU_vs_GPU(self):

-        print(__name__)

-        filename = os.path.join("lena_gray_512.tif")

-        plt = TiffReader()

-        # read image

-        Im = plt.imread(filename)                     

-        Im = np.asarray(Im, dtype='float32')

-        

-        Im = Im/255

-        perc = 0.05

-        u0 = Im + np.random.normal(loc = 0 ,

-                                          scale = perc * Im , 

-                                          size = np.shape(Im))

-        u_ref = Im + np.random.normal(loc = 0 ,

-                                          scale = 0.01 * Im , 

-                                          size = np.shape(Im))

-        

-        # map the u0 u0->u0>0

-        # f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1)

-        u0 = u0.astype('float32')

-        u_ref = u_ref.astype('float32')

-        

-        print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")

-        print ("____________LLT-ROF bench___________________")

-        print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")

-        

-        

-        # set parameters

-        pars = {'algorithm' : LLT_ROF, \

-                'input' : u0,\

-                'regularisation_parameterROF':0.04, \

-                'regularisation_parameterLLT':0.01, \

-                'number_of_iterations' :1000 ,\

-                'time_marching_parameter' :0.0001 ,\

-                }

-                

-        print ("#############LLT- ROF CPU####################")

-        start_time = timeit.default_timer()

-        lltrof_cpu = LLT_ROF(pars['input'], 

-                      pars['regularisation_parameterROF'],

-                      pars['regularisation_parameterLLT'],

-                      pars['number_of_iterations'],

-                      pars['time_marching_parameter'],'cpu')

-        

-        rms = rmse(Im, lltrof_cpu)

-        pars['rmse'] = rms

-        

-        txtstr = printParametersToString(pars)

-        txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)

-        print (txtstr)

-        print ("#############LLT- ROF GPU####################")

-        start_time = timeit.default_timer()

-        try:

-            lltrof_gpu = LLT_ROF(pars['input'], 

-                      pars['regularisation_parameterROF'],

-                      pars['regularisation_parameterLLT'],

-                      pars['number_of_iterations'],

-                      pars['time_marching_parameter'],'gpu')

-        

-        except ValueError as ve:

-            self.skipTest("Results not comparable. GPU computing error.")

-

-        rms = rmse(Im, lltrof_gpu)

-        pars['rmse'] = rms

-        pars['algorithm'] = LLT_ROF

-        txtstr = printParametersToString(pars)

-        txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)

-        print (txtstr)

-        print ("--------Compare the results--------")

-        tolerance = 1e-04

-        diff_im = np.zeros(np.shape(lltrof_gpu))

-        diff_im = abs(lltrof_cpu - lltrof_gpu)

-        diff_im[diff_im > tolerance] = 1

-        self.assertLessEqual(diff_im.sum(), 1)

-

-    def test_NDF_CPU_vs_GPU(self):

-        print(__name__)

-        filename = os.path.join("lena_gray_512.tif")

-        plt = TiffReader()

-        # read image

-        Im = plt.imread(filename)                     

-        Im = np.asarray(Im, dtype='float32')

-        

-        Im = Im/255

-        perc = 0.05

-        u0 = Im + np.random.normal(loc = 0 ,

-                                          scale = perc * Im , 

-                                          size = np.shape(Im))

-        u_ref = Im + np.random.normal(loc = 0 ,

-                                          scale = 0.01 * Im , 

-                                          size = np.shape(Im))

-        

-        # map the u0 u0->u0>0

-        # f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1)

-        u0 = u0.astype('float32')

-        u_ref = u_ref.astype('float32')

-        

-        print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")

-        print ("_______________NDF bench___________________")

-        print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")

-        

-        

-        # set parameters

-        pars = {'algorithm' : NDF, \

-                'input' : u0,\

-                'regularisation_parameter':0.06, \

-                'edge_parameter':0.04,\

-                'number_of_iterations' :1000 ,\

-                'time_marching_parameter':0.025,\

-                'penalty_type':  1

-                }

-                

-        print ("#############NDF CPU####################")

-        start_time = timeit.default_timer()

-        ndf_cpu = NDF(pars['input'], 

-                      pars['regularisation_parameter'],

-                      pars['edge_parameter'], 

-                      pars['number_of_iterations'],

-                      pars['time_marching_parameter'], 

-                      pars['penalty_type'],'cpu')

-                     

-        rms = rmse(Im, ndf_cpu)

-        pars['rmse'] = rms

-        

-        txtstr = printParametersToString(pars)

-        txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)

-        print (txtstr)

-        

-        print ("##############NDF GPU##################")

-        start_time = timeit.default_timer()

-        try:

-            ndf_gpu = NDF(pars['input'], 

-                      pars['regularisation_parameter'],

-                      pars['edge_parameter'], 

-                      pars['number_of_iterations'],

-                      pars['time_marching_parameter'], 

-                      pars['penalty_type'],'gpu')

-                     

-        except ValueError as ve:

-            self.skipTest("Results not comparable. GPU computing error.")

-        rms = rmse(Im, ndf_gpu)

-        pars['rmse'] = rms

-        pars['algorithm'] = NDF

-        txtstr = printParametersToString(pars)

-        txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)

-        print (txtstr)

-        print ("--------Compare the results--------")

-        tolerance = 1e-05

-        diff_im = np.zeros(np.shape(ndf_cpu))

-        diff_im = abs(ndf_cpu - ndf_gpu)

-        diff_im[diff_im > tolerance] = 1

-        self.assertLessEqual(diff_im.sum(), 1)

-

-        

-    def test_Diff4th_CPU_vs_GPU(self):

-        filename = os.path.join("lena_gray_512.tif")

-        plt = TiffReader()

-        # read image

-        Im = plt.imread(filename)                     

-        Im = np.asarray(Im, dtype='float32')

-        

-        Im = Im/255

-        perc = 0.05

-        u0 = Im + np.random.normal(loc = 0 ,

-                                          scale = perc * Im , 

-                                          size = np.shape(Im))

-        u_ref = Im + np.random.normal(loc = 0 ,

-                                          scale = 0.01 * Im , 

-                                          size = np.shape(Im))

-        

-        # map the u0 u0->u0>0

-        # f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1)

-        u0 = u0.astype('float32')

-        u_ref = u_ref.astype('float32')

-        

-        print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")

-        print ("___Anisotropic Diffusion 4th Order (2D)____")

-        print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")

-        

-        # set parameters

-        pars = {'algorithm' : Diff4th, \

-        'input' : u0,\

-        'regularisation_parameter':3.5, \

-        'edge_parameter':0.02,\

-        'number_of_iterations' :500 ,\

-        'time_marching_parameter':0.001

-        }

-        

-        print ("#############Diff4th CPU####################")

-        start_time = timeit.default_timer()

-        diff4th_cpu = Diff4th(pars['input'], 

-                      pars['regularisation_parameter'],

-                      pars['edge_parameter'], 

-                      pars['number_of_iterations'],

-                      pars['time_marching_parameter'],'cpu')

-                     

-        rms = rmse(Im, diff4th_cpu)

-        pars['rmse'] = rms

-

-        txtstr = printParametersToString(pars)

-        txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)

-        print (txtstr)

-        print ("##############Diff4th GPU##################")

-        start_time = timeit.default_timer()

-        try:

-            diff4th_gpu = Diff4th(pars['input'], 

-                      pars['regularisation_parameter'],

-                      pars['edge_parameter'], 

-                      pars['number_of_iterations'],

-                      pars['time_marching_parameter'], 'gpu')

-                     

-        except ValueError as ve:

-            self.skipTest("Results not comparable. GPU computing error.")

-        rms = rmse(Im, diff4th_gpu)

-        pars['rmse'] = rms

-        pars['algorithm'] = Diff4th

-        txtstr = printParametersToString(pars)

-        txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)

-        print (txtstr)

-        print ("--------Compare the results--------")

-        tolerance = 1e-05

-        diff_im = np.zeros(np.shape(diff4th_cpu))

-        diff_im = abs(diff4th_cpu - diff4th_gpu)

-        diff_im[diff_im > tolerance] = 1

-        self.assertLessEqual(diff_im.sum() , 1)

-

-    def test_FDGdTV_CPU_vs_GPU(self):

-        filename = os.path.join("lena_gray_512.tif")

-        plt = TiffReader()

-        # read image

-        Im = plt.imread(filename)                     

-        Im = np.asarray(Im, dtype='float32')

-        

-        Im = Im/255

-        perc = 0.05

-        u0 = Im + np.random.normal(loc = 0 ,

-                                          scale = perc * Im , 

-                                          size = np.shape(Im))

-        u_ref = Im + np.random.normal(loc = 0 ,

-                                          scale = 0.01 * Im , 

-                                          size = np.shape(Im))

-        

-        # map the u0 u0->u0>0

-        # f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1)

-        u0 = u0.astype('float32')

-        u_ref = u_ref.astype('float32')

-        

-

-        print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")

-        print ("____________FGP-dTV bench___________________")

-        print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")

-        

-        # set parameters

-        pars = {'algorithm' : FGP_dTV, \

-                'input' : u0,\

-                'refdata' : u_ref,\

-                'regularisation_parameter':0.04, \

-                'number_of_iterations' :1000 ,\

-                'tolerance_constant':1e-07,\

-                'eta_const':0.2,\

-                'methodTV': 0 ,\

-                'nonneg': 0 ,\

-                'printingOut': 0 

-                }

-                

-        print ("#############FGP dTV CPU####################")

-        start_time = timeit.default_timer()

-        fgp_dtv_cpu = FGP_dTV(pars['input'], 

-                      pars['refdata'], 

-                      pars['regularisation_parameter'],

-                      pars['number_of_iterations'],

-                      pars['tolerance_constant'], 

-                      pars['eta_const'], 

-                      pars['methodTV'],

-                      pars['nonneg'],

-                      pars['printingOut'],'cpu')

-                     

-                     

-        rms = rmse(Im, fgp_dtv_cpu)

-        pars['rmse'] = rms

-        

-        txtstr = printParametersToString(pars)

-        txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)

-        print (txtstr)

-        print ("##############FGP dTV GPU##################")

-        start_time = timeit.default_timer()

-        try:

-            fgp_dtv_gpu = FGP_dTV(pars['input'], 

-                      pars['refdata'], 

-                      pars['regularisation_parameter'],

-                      pars['number_of_iterations'],

-                      pars['tolerance_constant'], 

-                      pars['eta_const'], 

-                      pars['methodTV'],

-                      pars['nonneg'],

-                      pars['printingOut'],'gpu')

-        except ValueError as ve:

-            self.skipTest("Results not comparable. GPU computing error.")

-        rms = rmse(Im, fgp_dtv_gpu)

-        pars['rmse'] = rms

-        pars['algorithm'] = FGP_dTV

-        txtstr = printParametersToString(pars)

-        txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)

-        print (txtstr)

-        print ("--------Compare the results--------")

-        tolerance = 1e-05

-        diff_im = np.zeros(np.shape(fgp_dtv_cpu))

-        diff_im = abs(fgp_dtv_cpu - fgp_dtv_gpu)

-        diff_im[diff_im > tolerance] = 1

-        self.assertLessEqual(diff_im.sum(), 1)

-

-    def test_cpu_ROF_TV(self):

-        #filename = os.path.join(".." , ".." , ".." , "data" ,"testLena.npy")

-        

-        filename = os.path.join("lena_gray_512.tif")

-

-        plt = TiffReader()

-        # read image

-        Im = plt.imread(filename)                     

-        Im = np.asarray(Im, dtype='float32')

-        Im = Im/255

-        

-        """

-        # read noiseless image

-        Im = plt.imread(filename)

-        Im = np.asarray(Im, dtype='float32')

-        """

-        tolerance = 1e-05

-        rms_rof_exp = 8.313131464999238e-05 #expected value for ROF model

-

-        # set parameters for ROF-TV

-        pars_rof_tv = {'algorithm': ROF_TV, \

-                            'input' : Im,\

-                            'regularisation_parameter':0.04,\

-                            'number_of_iterations': 50,\

-                            'time_marching_parameter': 0.00001

-                            }

-        print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")

-        print ("_________testing ROF-TV (2D, CPU)__________")

-        print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")

-        rof_cpu = ROF_TV(pars_rof_tv['input'],

-             pars_rof_tv['regularisation_parameter'],

-             pars_rof_tv['number_of_iterations'],

-             pars_rof_tv['time_marching_parameter'],'cpu')

-        rms_rof = rmse(Im, rof_cpu)

-        

-        # now compare obtained rms with the expected value

-        self.assertLess(abs(rms_rof-rms_rof_exp) , tolerance)

-    def test_cpu_FGP_TV(self):

-        #filename = os.path.join(".." , ".." , ".." , "data" ,"testLena.npy")

-        

-        filename = os.path.join("lena_gray_512.tif")

-

-        plt = TiffReader()

-        # read image

-        Im = plt.imread(filename)                     

-        Im = np.asarray(Im, dtype='float32')

-        Im = Im/255

-        """

-        # read noiseless image

-        Im = plt.imread(filename)

-        Im = np.asarray(Im, dtype='float32')

-        """

-        tolerance = 1e-05

-        rms_fgp_exp = 0.019152347 #expected value for FGP model

-        

-        pars_fgp_tv = {'algorithm' : FGP_TV, \

-                            'input' : Im,\

-                            'regularisation_parameter':0.04, \

-                            'number_of_iterations' :50 ,\

-                            'tolerance_constant':1e-06,\

-                            'methodTV': 0 ,\

-                            'nonneg': 0 ,\

-                            'printingOut': 0 

-                            }

-        print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")

-        print ("_________testing FGP-TV (2D, CPU)__________")

-        print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")

-        fgp_cpu = FGP_TV(pars_fgp_tv['input'], 

-              pars_fgp_tv['regularisation_parameter'],

-              pars_fgp_tv['number_of_iterations'],

-              pars_fgp_tv['tolerance_constant'], 

-              pars_fgp_tv['methodTV'],

-              pars_fgp_tv['nonneg'],

-              pars_fgp_tv['printingOut'],'cpu')  

-        rms_fgp = rmse(Im, fgp_cpu)

-        # now compare obtained rms with the expected value

-        self.assertLess(abs(rms_fgp-rms_fgp_exp) , tolerance)

-

-    def test_gpu_ROF(self):

-        #filename = os.path.join(".." , ".." , ".." , "data" ,"testLena.npy")

-        filename = os.path.join("lena_gray_512.tif")

-

-        plt = TiffReader()

-        # read image

-        Im = plt.imread(filename)

-        Im = np.asarray(Im, dtype='float32')

-        Im = Im/255

-        

-        tolerance = 1e-05

-        rms_rof_exp = 8.313131464999238e-05 #expected value for ROF model

-        

-        # set parameters for ROF-TV

-        pars_rof_tv = {'algorithm': ROF_TV, \

-                            'input' : Im,\

-                            'regularisation_parameter':0.04,\

-                            'number_of_iterations': 50,\

-                            'time_marching_parameter': 0.00001

-                            }

-        print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")

-        print ("_________testing ROF-TV (2D, GPU)__________")

-        print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")

-        try:

-            rof_gpu = ROF_TV(pars_rof_tv['input'],

-             pars_rof_tv['regularisation_parameter'],

-             pars_rof_tv['number_of_iterations'],

-             pars_rof_tv['time_marching_parameter'],'gpu')

-        except ValueError as ve:

-            self.skipTest("Results not comparable. GPU computing error.")

-

-        rms_rof = rmse(Im, rof_gpu)

-        # now compare obtained rms with the expected value

-        self.assertLess(abs(rms_rof-rms_rof_exp) , tolerance)

-    

-    def test_gpu_FGP(self):

-        #filename = os.path.join(".." , ".." , ".." , "data" ,"testLena.npy")

-        filename = os.path.join("lena_gray_512.tif")

-

-        plt = TiffReader()

-        # read image

-        Im = plt.imread(filename)                     

-        Im = np.asarray(Im, dtype='float32')

-        Im = Im/255

-        tolerance = 1e-05

-        

-        rms_fgp_exp = 0.019152347 #expected value for FGP model

-        

-        # set parameters for FGP-TV

-        pars_fgp_tv = {'algorithm' : FGP_TV, \

-                            'input' : Im,\

-                            'regularisation_parameter':0.04, \

-                            'number_of_iterations' :50 ,\

-                            'tolerance_constant':1e-06,\

-                            'methodTV': 0 ,\

-                            'nonneg': 0 ,\

-                            'printingOut': 0 

-                            }

-        print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")

-        print ("_________testing FGP-TV (2D, GPU)__________")

-        print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")

-        try:

-            fgp_gpu = FGP_TV(pars_fgp_tv['input'], 

-              pars_fgp_tv['regularisation_parameter'],

-              pars_fgp_tv['number_of_iterations'],

-              pars_fgp_tv['tolerance_constant'], 

-              pars_fgp_tv['methodTV'],

-              pars_fgp_tv['nonneg'],

-              pars_fgp_tv['printingOut'],'gpu')  

-        except ValueError as ve:

-            self.skipTest("Results not comparable. GPU computing error.")

-        rms_fgp = rmse(Im, fgp_gpu)

-        # now compare obtained rms with the expected value

-

-        self.assertLess(abs(rms_fgp-rms_fgp_exp) , tolerance)

-

-if __name__ == '__main__':

-    unittest.main()

diff --git a/Wrappers/Python/demos/SoftwareX_supp/Demo_RealData_Recon_SX.py b/Wrappers/Python/demos/SoftwareX_supp/Demo_RealData_Recon_SX.py
deleted file mode 100644
index 01491d9..0000000
--- a/Wrappers/Python/demos/SoftwareX_supp/Demo_RealData_Recon_SX.py
+++ /dev/null
@@ -1,231 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-"""
-This demo scripts support the following publication: 
-"CCPi-Regularisation Toolkit for computed tomographic image reconstruction with 
-proximal splitting algorithms" by Daniil Kazantsev, Edoardo Pasca, Martin J. Turner,
- Philip J. Withers; Software X, 2019
-____________________________________________________________________________
-* Reads real tomographic data (stored at Zenodo)
---- https://doi.org/10.5281/zenodo.2578893
-* Reconstructs using TomoRec software
-* Saves reconstructed images 
-____________________________________________________________________________
->>>>> Dependencies: <<<<<
-1. ASTRA toolbox: conda install -c astra-toolbox astra-toolbox
-2. TomoRec: conda install -c dkazanc tomorec
-or install from https://github.com/dkazanc/TomoRec
-3. libtiff if one needs to save tiff images:
-    install pip install libtiff
-
-@author: Daniil Kazantsev, e:mail daniil.kazantsev@diamond.ac.uk
-GPLv3 license (ASTRA toolbox)
-"""
-import numpy as np
-import matplotlib.pyplot as plt
-import h5py
-from tomorec.supp.suppTools import normaliser
-import time
-
-# load dendritic projection data
-h5f = h5py.File('data/DendrData_3D.h5','r')
-dataRaw = h5f['dataRaw'][:]
-flats = h5f['flats'][:]
-darks = h5f['darks'][:]
-angles_rad = h5f['angles_rad'][:]
-h5f.close()
-#%%
-# normalise the data [detectorsVert, Projections, detectorsHoriz]
-data_norm = normaliser(dataRaw, flats, darks, log='log')
-del dataRaw, darks, flats
-
-intens_max = 2.3
-plt.figure() 
-plt.subplot(131)
-plt.imshow(data_norm[:,150,:],vmin=0, vmax=intens_max)
-plt.title('2D Projection (analytical)')
-plt.subplot(132)
-plt.imshow(data_norm[300,:,:],vmin=0, vmax=intens_max)
-plt.title('Sinogram view')
-plt.subplot(133)
-plt.imshow(data_norm[:,:,600],vmin=0, vmax=intens_max)
-plt.title('Tangentogram view')
-plt.show()
-
-detectorHoriz = np.size(data_norm,2)
-det_y_crop = [i for i in range(0,detectorHoriz-22)]
-N_size = 950 # reconstruction domain
-time_label = int(time.time())
-#%%
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("%%%%%%%%%%%%Reconstructing with FBP method %%%%%%%%%%%%%%%%%")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-from tomorec.methodsDIR import RecToolsDIR
-
-RectoolsDIR = RecToolsDIR(DetectorsDimH = np.size(det_y_crop),  # DetectorsDimH # detector dimension (horizontal)
-                    DetectorsDimV = 100,  # DetectorsDimV # detector dimension (vertical) for 3D case only
-                    AnglesVec = angles_rad, # array of angles in radians
-                    ObjSize = N_size, # a scalar to define reconstructed object dimensions
-                    device='gpu')
-
-FBPrec = RectoolsDIR.FBP(data_norm[0:100,:,det_y_crop])
-
-sliceSel = 50
-max_val = 0.003
-plt.figure() 
-plt.subplot(131)
-plt.imshow(FBPrec[sliceSel,:,:],vmin=0, vmax=max_val, cmap="gray")
-plt.title('FBP Reconstruction, axial view')
-
-plt.subplot(132)
-plt.imshow(FBPrec[:,sliceSel,:],vmin=0, vmax=max_val, cmap="gray")
-plt.title('FBP Reconstruction, coronal view')
-
-plt.subplot(133)
-plt.imshow(FBPrec[:,:,sliceSel],vmin=0, vmax=max_val, cmap="gray")
-plt.title('FBP Reconstruction, sagittal view')
-plt.show()
-
-# saving to tiffs (16bit)
-"""
-from libtiff import TIFF
-FBPrec += np.abs(np.min(FBPrec))
-multiplier = (int)(65535/(np.max(FBPrec)))
-
-# saving to tiffs (16bit)
-for i in range(0,np.size(FBPrec,0)):
-    tiff = TIFF.open('Dendr_FBP'+'_'+str(i)+'.tiff', mode='w')
-    tiff.write_image(np.uint16(FBPrec[i,:,:]*multiplier))
-    tiff.close()
-"""
-#%%
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("Reconstructing with ADMM method using TomoRec software")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-# initialise TomoRec ITERATIVE reconstruction class ONCE
-from tomorec.methodsIR import RecToolsIR
-RectoolsIR = RecToolsIR(DetectorsDimH =  np.size(det_y_crop),  # DetectorsDimH # detector dimension (horizontal)
-                    DetectorsDimV = 100,  # DetectorsDimV # detector dimension (vertical) for 3D case only
-                    AnglesVec = angles_rad, # array of angles in radians
-                    ObjSize = N_size, # a scalar to define reconstructed object dimensions
-                    datafidelity='LS',# data fidelity, choose LS, PWLS (wip), GH (wip), Student (wip)
-                    nonnegativity='ENABLE', # enable nonnegativity constraint (set to 'ENABLE')
-                    OS_number = None, # the number of subsets, NONE/(or > 1) ~ classical / ordered subsets
-                    tolerance = 1e-08, # tolerance to stop outer iterations earlier
-                    device='gpu')
-#%%
-print ("Reconstructing with ADMM method using SB-TV penalty")
-RecADMM_reg_sbtv = RectoolsIR.ADMM(data_norm[0:100,:,det_y_crop],
-                              rho_const = 2000.0, \
-                              iterationsADMM = 15, \
-                              regularisation = 'SB_TV', \
-                              regularisation_parameter = 0.00085,\
-                              regularisation_iterations = 50)
-
-sliceSel = 50
-max_val = 0.003
-plt.figure() 
-plt.subplot(131)
-plt.imshow(RecADMM_reg_sbtv[sliceSel,:,:],vmin=0, vmax=max_val, cmap="gray")
-plt.title('3D ADMM-SB-TV Reconstruction, axial view')
-
-plt.subplot(132)
-plt.imshow(RecADMM_reg_sbtv[:,sliceSel,:],vmin=0, vmax=max_val, cmap="gray")
-plt.title('3D ADMM-SB-TV Reconstruction, coronal view')
-
-plt.subplot(133)
-plt.imshow(RecADMM_reg_sbtv[:,:,sliceSel],vmin=0, vmax=max_val, cmap="gray")
-plt.title('3D ADMM-SB-TV Reconstruction, sagittal view')
-plt.show()
-
-
-# saving to tiffs (16bit)
-"""
-from libtiff import TIFF
-multiplier = (int)(65535/(np.max(RecADMM_reg_sbtv)))
-for i in range(0,np.size(RecADMM_reg_sbtv,0)):
-    tiff = TIFF.open('Dendr_ADMM_SBTV'+'_'+str(i)+'.tiff', mode='w')
-    tiff.write_image(np.uint16(RecADMM_reg_sbtv[i,:,:]*multiplier))
-    tiff.close()
-"""
-# Saving recpnstructed data with a unique time label
-np.save('Dendr_ADMM_SBTV'+str(time_label)+'.npy', RecADMM_reg_sbtv)
-del RecADMM_reg_sbtv
-#%%
-print ("Reconstructing with ADMM method using ROF-LLT penalty")
-RecADMM_reg_rofllt = RectoolsIR.ADMM(data_norm[0:100,:,det_y_crop],
-                              rho_const = 2000.0, \
-                              iterationsADMM = 15, \
-                              regularisation = 'LLT_ROF', \
-                              regularisation_parameter = 0.0009,\
-                              regularisation_parameter2 = 0.0007,\
-                              time_marching_parameter = 0.001,\
-                              regularisation_iterations = 550)
-
-sliceSel = 50
-max_val = 0.003
-plt.figure() 
-plt.subplot(131)
-plt.imshow(RecADMM_reg_rofllt[sliceSel,:,:],vmin=0, vmax=max_val)
-plt.title('3D ADMM-ROFLLT Reconstruction, axial view')
-
-plt.subplot(132)
-plt.imshow(RecADMM_reg_rofllt[:,sliceSel,:],vmin=0, vmax=max_val)
-plt.title('3D ADMM-ROFLLT Reconstruction, coronal view')
-
-plt.subplot(133)
-plt.imshow(RecADMM_reg_rofllt[:,:,sliceSel],vmin=0, vmax=max_val)
-plt.title('3D ADMM-ROFLLT Reconstruction, sagittal view')
-plt.show()
-
-# saving to tiffs (16bit)
-"""
-from libtiff import TIFF
-multiplier = (int)(65535/(np.max(RecADMM_reg_rofllt)))
-for i in range(0,np.size(RecADMM_reg_rofllt,0)):
-    tiff = TIFF.open('Dendr_ADMM_ROFLLT'+'_'+str(i)+'.tiff', mode='w')
-    tiff.write_image(np.uint16(RecADMM_reg_rofllt[i,:,:]*multiplier))
-    tiff.close()
-"""
-
-# Saving recpnstructed data with a unique time label
-np.save('Dendr_ADMM_ROFLLT'+str(time_label)+'.npy', RecADMM_reg_rofllt)
-del RecADMM_reg_rofllt
-#%%
-print ("Reconstructing with ADMM method using TGV penalty")
-RecADMM_reg_tgv = RectoolsIR.ADMM(data_norm[0:100,:,det_y_crop],
-                              rho_const = 2000.0, \
-                              iterationsADMM = 15, \
-                              regularisation = 'TGV', \
-                              regularisation_parameter = 0.01,\
-                              regularisation_iterations = 500)
-
-sliceSel = 50
-max_val = 0.003
-plt.figure() 
-plt.subplot(131)
-plt.imshow(RecADMM_reg_tgv[sliceSel,:,:],vmin=0, vmax=max_val)
-plt.title('3D ADMM-TGV Reconstruction, axial view')
-
-plt.subplot(132)
-plt.imshow(RecADMM_reg_tgv[:,sliceSel,:],vmin=0, vmax=max_val)
-plt.title('3D ADMM-TGV Reconstruction, coronal view')
-
-plt.subplot(133)
-plt.imshow(RecADMM_reg_tgv[:,:,sliceSel],vmin=0, vmax=max_val)
-plt.title('3D ADMM-TGV Reconstruction, sagittal view')
-plt.show()
-
-# saving to tiffs (16bit)
-"""
-from libtiff import TIFF
-multiplier = (int)(65535/(np.max(RecADMM_reg_tgv)))
-for i in range(0,np.size(RecADMM_reg_tgv,0)):
-    tiff = TIFF.open('Dendr_ADMM_TGV'+'_'+str(i)+'.tiff', mode='w')
-    tiff.write_image(np.uint16(RecADMM_reg_tgv[i,:,:]*multiplier))
-    tiff.close()
-"""
-# Saving recpnstructed data with a unique time label
-np.save('Dendr_ADMM_TGV'+str(time_label)+'.npy', RecADMM_reg_tgv)
-del RecADMM_reg_tgv
-#%%
\ No newline at end of file
diff --git a/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_ParOptimis_SX.py b/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_ParOptimis_SX.py
deleted file mode 100644
index 59ffc0e..0000000
--- a/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_ParOptimis_SX.py
+++ /dev/null
@@ -1,161 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-"""
-This demo scripts support the following publication: 
-"CCPi-Regularisation Toolkit for computed tomographic image reconstruction with 
-proximal splitting algorithms" by Daniil Kazantsev, Edoardo Pasca, Martin J. Turner,
- Philip J. Withers; Software X, 2019
-____________________________________________________________________________
-* Reads data which is previosly generated by TomoPhantom software (Zenodo link)
---- https://doi.org/10.5281/zenodo.2578893
-* Optimises for the regularisation parameters which later used in the script:
-Demo_SimulData_Recon_SX.py
-____________________________________________________________________________
->>>>> Dependencies: <<<<<
->>>>> Dependencies: <<<<<
-1. ASTRA toolbox: conda install -c astra-toolbox astra-toolbox
-2. TomoRec: conda install -c dkazanc tomorec
-or install from https://github.com/dkazanc/TomoRec
-
-@author: Daniil Kazantsev, e:mail daniil.kazantsev@diamond.ac.uk
-GPLv3 license (ASTRA toolbox)
-"""
-#import timeit
-import matplotlib.pyplot as plt
-import numpy as np
-import h5py
-from ccpi.supp.qualitymetrics import QualityTools
-
-# loading the data 
-h5f = h5py.File('data/TomoSim_data1550671417.h5','r')
-phantom = h5f['phantom'][:]
-projdata_norm = h5f['projdata_norm'][:]
-proj_angles = h5f['proj_angles'][:]
-h5f.close()
-
-[Vert_det, AnglesNum, Horiz_det] = np.shape(projdata_norm)
-N_size = Vert_det
-
-sliceSel = 128
-#plt.gray()
-plt.figure() 
-plt.subplot(131)
-plt.imshow(phantom[sliceSel,:,:],vmin=0, vmax=1)
-plt.title('3D Phantom, axial view')
-
-plt.subplot(132)
-plt.imshow(phantom[:,sliceSel,:],vmin=0, vmax=1)
-plt.title('3D Phantom, coronal view')
-
-plt.subplot(133)
-plt.imshow(phantom[:,:,sliceSel],vmin=0, vmax=1)
-plt.title('3D Phantom, sagittal view')
-plt.show()
-
-intens_max = 240
-plt.figure() 
-plt.subplot(131)
-plt.imshow(projdata_norm[:,sliceSel,:],vmin=0, vmax=intens_max)
-plt.title('2D Projection (erroneous)')
-plt.subplot(132)
-plt.imshow(projdata_norm[sliceSel,:,:],vmin=0, vmax=intens_max)
-plt.title('Sinogram view')
-plt.subplot(133)
-plt.imshow(projdata_norm[:,:,sliceSel],vmin=0, vmax=intens_max)
-plt.title('Tangentogram view')
-plt.show()
-#%%
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("Reconstructing with ADMM method using TomoRec software")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-# initialise TomoRec ITERATIVE reconstruction class ONCE
-from tomorec.methodsIR import RecToolsIR
-RectoolsIR = RecToolsIR(DetectorsDimH = Horiz_det,  # DetectorsDimH # detector dimension (horizontal)
-                    DetectorsDimV = Vert_det,  # DetectorsDimV # detector dimension (vertical) for 3D case only
-                    AnglesVec = proj_angles, # array of angles in radians
-                    ObjSize = N_size, # a scalar to define reconstructed object dimensions
-                    datafidelity='LS',# data fidelity, choose LS, PWLS (wip), GH (wip), Student (wip)
-                    nonnegativity='ENABLE', # enable nonnegativity constraint (set to 'ENABLE')
-                    OS_number = None, # the number of subsets, NONE/(or > 1) ~ classical / ordered subsets
-                    tolerance = 1e-08, # tolerance to stop outer iterations earlier
-                    device='gpu')
-#%%
-param_space = 30
-reg_param_sb_vec = np.linspace(0.03,0.15,param_space,dtype='float32') # a vector of parameters
-erros_vec_sbtv = np.zeros((param_space)) # a vector of errors
-
-print ("Reconstructing with ADMM method using SB-TV penalty")
-for i in range(0,param_space):
-    RecADMM_reg_sbtv = RectoolsIR.ADMM(projdata_norm,
-                                  rho_const = 2000.0, \
-                                  iterationsADMM = 15, \
-                                  regularisation = 'SB_TV', \
-                                  regularisation_parameter = reg_param_sb_vec[i],\
-                                  regularisation_iterations = 50)
-    # calculate errors 
-    Qtools = QualityTools(phantom, RecADMM_reg_sbtv)
-    erros_vec_sbtv[i] = Qtools.rmse()
-    print("RMSE for regularisation parameter {} for ADMM-SB-TV is {}".format(reg_param_sb_vec[i],erros_vec_sbtv[i]))
-
-plt.figure() 
-plt.plot(erros_vec_sbtv)
-
-# Saving generated data with a unique time label
-h5f = h5py.File('Optim_admm_sbtv.h5', 'w')
-h5f.create_dataset('reg_param_sb_vec', data=reg_param_sb_vec)
-h5f.create_dataset('erros_vec_sbtv', data=erros_vec_sbtv)
-h5f.close()
-#%%
-param_space = 30
-reg_param_rofllt_vec = np.linspace(0.03,0.15,param_space,dtype='float32') # a vector of parameters
-erros_vec_rofllt = np.zeros((param_space)) # a vector of errors
-
-print ("Reconstructing with ADMM method using ROF-LLT penalty")
-for i in range(0,param_space):
-    RecADMM_reg_rofllt = RectoolsIR.ADMM(projdata_norm,
-                                  rho_const = 2000.0, \
-                                  iterationsADMM = 15, \
-                                  regularisation = 'LLT_ROF', \
-                                  regularisation_parameter = reg_param_rofllt_vec[i],\
-                                  regularisation_parameter2 = 0.005,\
-                                  regularisation_iterations = 600)
-    # calculate errors 
-    Qtools = QualityTools(phantom, RecADMM_reg_rofllt)
-    erros_vec_rofllt[i] = Qtools.rmse()
-    print("RMSE for regularisation parameter {} for ADMM-ROF-LLT is {}".format(reg_param_rofllt_vec[i],erros_vec_rofllt[i]))
-
-plt.figure() 
-plt.plot(erros_vec_rofllt)
-
-# Saving generated data with a unique time label
-h5f = h5py.File('Optim_admm_rofllt.h5', 'w')
-h5f.create_dataset('reg_param_rofllt_vec', data=reg_param_rofllt_vec)
-h5f.create_dataset('erros_vec_rofllt', data=erros_vec_rofllt)
-h5f.close()
-#%%
-param_space = 30
-reg_param_tgv_vec = np.linspace(0.03,0.15,param_space,dtype='float32') # a vector of parameters
-erros_vec_tgv = np.zeros((param_space)) # a vector of errors
-
-print ("Reconstructing with ADMM method using TGV penalty")
-for i in range(0,param_space):
-    RecADMM_reg_tgv = RectoolsIR.ADMM(projdata_norm,
-                                  rho_const = 2000.0, \
-                                  iterationsADMM = 15, \
-                                  regularisation = 'TGV', \
-                                  regularisation_parameter = reg_param_tgv_vec[i],\
-                                  regularisation_iterations = 600)
-    # calculate errors 
-    Qtools = QualityTools(phantom, RecADMM_reg_tgv)
-    erros_vec_tgv[i] = Qtools.rmse()
-    print("RMSE for regularisation parameter {} for ADMM-TGV is {}".format(reg_param_tgv_vec[i],erros_vec_tgv[i]))
-
-plt.figure() 
-plt.plot(erros_vec_tgv)
-
-# Saving generated data with a unique time label
-h5f = h5py.File('Optim_admm_tgv.h5', 'w')
-h5f.create_dataset('reg_param_tgv_vec', data=reg_param_tgv_vec)
-h5f.create_dataset('erros_vec_tgv', data=erros_vec_tgv)
-h5f.close()
-#%%
\ No newline at end of file
diff --git a/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_Recon_SX.py b/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_Recon_SX.py
deleted file mode 100644
index 93b0cef..0000000
--- a/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_Recon_SX.py
+++ /dev/null
@@ -1,309 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-"""
-This demo scripts support the following publication: 
-"CCPi-Regularisation Toolkit for computed tomographic image reconstruction with 
-proximal splitting algorithms" by Daniil Kazantsev, Edoardo Pasca, Martin J. Turner,
- Philip J. Withers; Software X, 2019
-____________________________________________________________________________
-* Reads data which is previously generated by TomoPhantom software (Zenodo link)
---- https://doi.org/10.5281/zenodo.2578893
-* Reconstruct using optimised regularisation parameters (see Demo_SimulData_ParOptimis_SX.py)
-____________________________________________________________________________
->>>>> Dependencies: <<<<<
-1. ASTRA toolbox: conda install -c astra-toolbox astra-toolbox
-2. TomoRec: conda install -c dkazanc tomorec
-or install from https://github.com/dkazanc/TomoRec
-
-@author: Daniil Kazantsev, e:mail daniil.kazantsev@diamond.ac.uk
-GPLv3 license (ASTRA toolbox)
-"""
-#import timeit
-import matplotlib.pyplot as plt
-import matplotlib.gridspec as gridspec
-import numpy as np
-import h5py
-from ccpi.supp.qualitymetrics import QualityTools
-from scipy.signal import gaussian
-
-# loading the data 
-h5f = h5py.File('data/TomoSim_data1550671417.h5','r')
-phantom = h5f['phantom'][:]
-projdata_norm = h5f['projdata_norm'][:]
-proj_angles = h5f['proj_angles'][:]
-h5f.close()
-
-[Vert_det, AnglesNum, Horiz_det] = np.shape(projdata_norm)
-N_size = Vert_det
-
-# loading optmisation parameters (the result of running Demo_SimulData_ParOptimis_SX)
-h5f = h5py.File('optim_param/Optim_admm_sbtv.h5','r')
-reg_param_sb_vec = h5f['reg_param_sb_vec'][:]
-erros_vec_sbtv = h5f['erros_vec_sbtv'][:]
-h5f.close()
-
-h5f = h5py.File('optim_param/Optim_admm_rofllt.h5','r')
-reg_param_rofllt_vec = h5f['reg_param_rofllt_vec'][:]
-erros_vec_rofllt = h5f['erros_vec_rofllt'][:]
-h5f.close()
-
-h5f = h5py.File('optim_param/Optim_admm_tgv.h5','r')
-reg_param_tgv_vec = h5f['reg_param_tgv_vec'][:]
-erros_vec_tgv = h5f['erros_vec_tgv'][:]
-h5f.close()
-
-index_minSBTV = min(xrange(len(erros_vec_sbtv)), key=erros_vec_sbtv.__getitem__)
-index_minROFLLT = min(xrange(len(erros_vec_rofllt)), key=erros_vec_rofllt.__getitem__)
-index_minTGV = min(xrange(len(erros_vec_tgv)), key=erros_vec_tgv.__getitem__)
-# assign optimal regularisation parameters:
-optimReg_sbtv = reg_param_sb_vec[index_minSBTV]
-optimReg_rofllt = reg_param_rofllt_vec[index_minROFLLT]
-optimReg_tgv = reg_param_tgv_vec[index_minTGV]
-#%%
-# plot loaded data
-sliceSel = 128
-#plt.figure() 
-fig, (ax1, ax2) = plt.subplots(figsize=(15, 5), ncols=2)
-plt.rcParams.update({'xtick.labelsize': 'x-small'})
-plt.rcParams.update({'ytick.labelsize':'x-small'})
-plt.subplot(121)
-one = plt.imshow(phantom[sliceSel,:,:],vmin=0, vmax=1, interpolation='none', cmap="PuOr")
-fig.colorbar(one, ax=ax1)
-plt.title('3D Phantom, axial (X-Y) view')
-plt.subplot(122)
-two = plt.imshow(phantom[:,sliceSel,:],vmin=0, vmax=1,interpolation='none', cmap="PuOr")
-fig.colorbar(two, ax=ax2)
-plt.title('3D Phantom, coronal (Y-Z) view')
-"""
-plt.subplot(133)
-plt.imshow(phantom[:,:,sliceSel],vmin=0, vmax=1, cmap="PuOr")
-plt.title('3D Phantom, sagittal view')
-
-"""
-plt.show()
-#%%
-intens_max = 220
-plt.figure() 
-plt.rcParams.update({'xtick.labelsize': 'x-small'})
-plt.rcParams.update({'ytick.labelsize':'x-small'})
-plt.subplot(131)
-plt.imshow(projdata_norm[:,sliceSel,:],vmin=0, vmax=intens_max, cmap="PuOr")
-plt.xlabel('X-detector', fontsize=16)
-plt.ylabel('Z-detector', fontsize=16)
-plt.title('2D Projection (X-Z) view', fontsize=19)
-plt.subplot(132)
-plt.imshow(projdata_norm[sliceSel,:,:],vmin=0, vmax=intens_max, cmap="PuOr")
-plt.xlabel('X-detector', fontsize=16)
-plt.ylabel('Projection angle', fontsize=16)
-plt.title('Sinogram (X-Y) view', fontsize=19)
-plt.subplot(133)
-plt.imshow(projdata_norm[:,:,sliceSel],vmin=0, vmax=intens_max, cmap="PuOr")
-plt.xlabel('Projection angle', fontsize=16)
-plt.ylabel('Z-detector', fontsize=16)
-plt.title('Vertical (Y-Z) view', fontsize=19)
-plt.show()
-#plt.savefig('projdata.pdf', format='pdf', dpi=1200)
-#%%
-# initialise TomoRec DIRECT reconstruction class ONCE
-from tomorec.methodsDIR import RecToolsDIR
-RectoolsDIR = RecToolsDIR(DetectorsDimH = Horiz_det,  # DetectorsDimH # detector dimension (horizontal)
-                    DetectorsDimV = Vert_det,  # DetectorsDimV # detector dimension (vertical) for 3D case only
-                    AnglesVec = proj_angles, # array of angles in radians
-                    ObjSize = N_size, # a scalar to define reconstructed object dimensions
-                    device = 'gpu')
-#%%
-print ("Reconstruction using FBP from TomoRec")
-recFBP= RectoolsDIR.FBP(projdata_norm) # FBP reconstruction
-#%%
-x0, y0 = 0, 127 # These are in _pixel_ coordinates!!
-x1, y1 = 255, 127
-
-sliceSel = int(0.5*N_size)
-max_val = 1
-plt.figure(figsize = (20,5))
-gs1 = gridspec.GridSpec(1, 3)
-gs1.update(wspace=0.1, hspace=0.05) # set the spacing between axes. 
-ax1 = plt.subplot(gs1[0])
-plt.imshow(recFBP[sliceSel,:,:],vmin=0, vmax=max_val, cmap="PuOr")
-ax1.plot([x0, x1], [y0, y1], 'ko-', linestyle='--')
-plt.colorbar(ax=ax1)
-plt.title('FBP Reconstruction, axial (X-Y) view', fontsize=19)
-ax1.set_aspect('equal')
-ax3 = plt.subplot(gs1[1])
-plt.plot(phantom[sliceSel,sliceSel,0:N_size],color='k',linewidth=2)
-plt.plot(recFBP[sliceSel,sliceSel,0:N_size],linestyle='--',color='g')
-plt.title('Profile', fontsize=19)
-ax2 = plt.subplot(gs1[2])
-plt.imshow(recFBP[:,sliceSel,:],vmin=0, vmax=max_val, cmap="PuOr")
-plt.title('FBP Reconstruction, coronal (Y-Z) view', fontsize=19)
-ax2.set_aspect('equal')
-plt.show()
-#plt.savefig('FBP_phantom.pdf', format='pdf', dpi=1600)
-
-# calculate errors 
-Qtools = QualityTools(phantom, recFBP)
-RMSE_fbp = Qtools.rmse()
-print("Root Mean Square Error for FBP is {}".format(RMSE_fbp))
-
-# SSIM measure
-Qtools = QualityTools(phantom[128,:,:]*255, recFBP[128,:,:]*235)
-win = np.array([gaussian(11, 1.5)])
-win2d = win * (win.T)
-ssim_fbp = Qtools.ssim(win2d)
-print("Mean SSIM for FBP is {}".format(ssim_fbp[0]))
-#%%
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("Reconstructing with ADMM method using TomoRec software")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-# initialise TomoRec ITERATIVE reconstruction class ONCE
-from tomorec.methodsIR import RecToolsIR
-RectoolsIR = RecToolsIR(DetectorsDimH = Horiz_det,  # DetectorsDimH # detector dimension (horizontal)
-                    DetectorsDimV = Vert_det,  # DetectorsDimV # detector dimension (vertical) for 3D case only
-                    AnglesVec = proj_angles, # array of angles in radians
-                    ObjSize = N_size, # a scalar to define reconstructed object dimensions
-                    datafidelity='LS',# data fidelity, choose LS, PWLS (wip), GH (wip), Student (wip)
-                    nonnegativity='ENABLE', # enable nonnegativity constraint (set to 'ENABLE')
-                    OS_number = None, # the number of subsets, NONE/(or > 1) ~ classical / ordered subsets
-                    tolerance = 1e-08, # tolerance to stop outer iterations earlier
-                    device='gpu')
-#%%
-print ("Reconstructing with ADMM method using SB-TV penalty")
-RecADMM_reg_sbtv = RectoolsIR.ADMM(projdata_norm,
-                                  rho_const = 2000.0, \
-                                  iterationsADMM = 25, \
-                                  regularisation = 'SB_TV', \
-                                  regularisation_parameter = optimReg_sbtv,\
-                                  regularisation_iterations = 50)
-
-sliceSel = int(0.5*N_size)
-max_val = 1
-plt.figure(figsize = (20,3))
-gs1 = gridspec.GridSpec(1, 4)
-gs1.update(wspace=0.02, hspace=0.01) # set the spacing between axes. 
-ax1 = plt.subplot(gs1[0])
-plt.plot(reg_param_sb_vec, erros_vec_sbtv, color='k',linewidth=2)
-plt.xlabel('Regularisation parameter', fontsize=16)
-plt.ylabel('RMSE value', fontsize=16)
-plt.title('Regularisation selection', fontsize=19)
-ax2 = plt.subplot(gs1[1])
-plt.imshow(RecADMM_reg_sbtv[sliceSel,:,:],vmin=0, vmax=max_val, cmap="PuOr")
-ax2.plot([x0, x1], [y0, y1], 'ko-', linestyle='--')
-plt.title('ADMM-SBTV (X-Y) view', fontsize=19)
-#ax2.set_aspect('equal')
-ax3 = plt.subplot(gs1[2])
-plt.plot(phantom[sliceSel,sliceSel,0:N_size],color='k',linewidth=2)
-plt.plot(RecADMM_reg_sbtv[sliceSel,sliceSel,0:N_size],linestyle='--',color='g')
-plt.title('Profile', fontsize=19)
-ax4 = plt.subplot(gs1[3])
-plt.imshow(RecADMM_reg_sbtv[:,sliceSel,:],vmin=0, vmax=max_val, cmap="PuOr")
-plt.title('ADMM-SBTV (Y-Z) view', fontsize=19)
-plt.colorbar(ax=ax4)
-plt.show()
-plt.savefig('SBTV_phantom.pdf', format='pdf', dpi=1600)
-
-# calculate errors 
-Qtools = QualityTools(phantom, RecADMM_reg_sbtv)
-RMSE_admm_sbtv = Qtools.rmse()
-print("Root Mean Square Error for ADMM-SB-TV is {}".format(RMSE_admm_sbtv))
-
-# SSIM measure
-Qtools = QualityTools(phantom[128,:,:]*255, RecADMM_reg_sbtv[128,:,:]*235)
-win = np.array([gaussian(11, 1.5)])
-win2d = win * (win.T)
-ssim_admm_sbtv = Qtools.ssim(win2d)
-print("Mean SSIM ADMM-SBTV is {}".format(ssim_admm_sbtv[0]))
-#%%
-print ("Reconstructing with ADMM method using ROFLLT penalty")
-RecADMM_reg_rofllt = RectoolsIR.ADMM(projdata_norm,
-                                  rho_const = 2000.0, \
-                                  iterationsADMM = 25, \
-                                  regularisation = 'LLT_ROF', \
-                                  regularisation_parameter = optimReg_rofllt,\
-                                  regularisation_parameter2 = 0.0085,\
-                                  regularisation_iterations = 600)
-
-sliceSel = int(0.5*N_size)
-max_val = 1
-plt.figure(figsize = (20,3))
-gs1 = gridspec.GridSpec(1, 4)
-gs1.update(wspace=0.02, hspace=0.01) # set the spacing between axes. 
-ax1 = plt.subplot(gs1[0])
-plt.plot(reg_param_rofllt_vec, erros_vec_rofllt, color='k',linewidth=2)
-plt.xlabel('Regularisation parameter', fontsize=16)
-plt.ylabel('RMSE value', fontsize=16)
-plt.title('Regularisation selection', fontsize=19)
-ax2 = plt.subplot(gs1[1])
-plt.imshow(RecADMM_reg_rofllt[sliceSel,:,:],vmin=0, vmax=max_val, cmap="PuOr")
-ax2.plot([x0, x1], [y0, y1], 'ko-', linestyle='--')
-plt.title('ADMM-ROFLLT (X-Y) view', fontsize=19)
-#ax2.set_aspect('equal')
-ax3 = plt.subplot(gs1[2])
-plt.plot(phantom[sliceSel,sliceSel,0:N_size],color='k',linewidth=2)
-plt.plot(RecADMM_reg_rofllt[sliceSel,sliceSel,0:N_size],linestyle='--',color='g')
-plt.title('Profile', fontsize=19)
-ax4 = plt.subplot(gs1[3])
-plt.imshow(RecADMM_reg_rofllt[:,sliceSel,:],vmin=0, vmax=max_val, cmap="PuOr")
-plt.title('ADMM-ROFLLT (Y-Z) view', fontsize=19)
-plt.colorbar(ax=ax4)
-plt.show()
-#plt.savefig('ROFLLT_phantom.pdf', format='pdf', dpi=1600)
-
-# calculate errors 
-Qtools = QualityTools(phantom, RecADMM_reg_rofllt)
-RMSE_admm_rofllt = Qtools.rmse()
-print("Root Mean Square Error for ADMM-ROF-LLT is {}".format(RMSE_admm_rofllt))
-
-# SSIM measure
-Qtools = QualityTools(phantom[128,:,:]*255, RecADMM_reg_rofllt[128,:,:]*235)
-win = np.array([gaussian(11, 1.5)])
-win2d = win * (win.T)
-ssim_admm_rifllt = Qtools.ssim(win2d)
-print("Mean SSIM ADMM-ROFLLT is {}".format(ssim_admm_rifllt[0]))
-#%%
-print ("Reconstructing with ADMM method using TGV penalty")
-RecADMM_reg_tgv = RectoolsIR.ADMM(projdata_norm,
-                                  rho_const = 2000.0, \
-                                  iterationsADMM = 25, \
-                                  regularisation = 'TGV', \
-                                  regularisation_parameter = optimReg_tgv,\
-                                  regularisation_iterations = 600)
-#%%
-sliceSel = int(0.5*N_size)
-max_val = 1
-plt.figure(figsize = (20,3))
-gs1 = gridspec.GridSpec(1, 4)
-gs1.update(wspace=0.02, hspace=0.01) # set the spacing between axes. 
-ax1 = plt.subplot(gs1[0])
-plt.plot(reg_param_tgv_vec, erros_vec_tgv, color='k',linewidth=2)
-plt.xlabel('Regularisation parameter', fontsize=16)
-plt.ylabel('RMSE value', fontsize=16)
-plt.title('Regularisation selection', fontsize=19)
-ax2 = plt.subplot(gs1[1])
-plt.imshow(RecADMM_reg_tgv[sliceSel,:,:],vmin=0, vmax=max_val, cmap="PuOr")
-ax2.plot([x0, x1], [y0, y1], 'ko-', linestyle='--')
-plt.title('ADMM-TGV (X-Y) view', fontsize=19)
-#ax2.set_aspect('equal')
-ax3 = plt.subplot(gs1[2])
-plt.plot(phantom[sliceSel,sliceSel,0:N_size],color='k',linewidth=2)
-plt.plot(RecADMM_reg_tgv[sliceSel,sliceSel,0:N_size],linestyle='--',color='g')
-plt.title('Profile', fontsize=19)
-ax4 = plt.subplot(gs1[3])
-plt.imshow(RecADMM_reg_tgv[:,sliceSel,:],vmin=0, vmax=max_val, cmap="PuOr")
-plt.title('ADMM-TGV (Y-Z) view', fontsize=19)
-plt.colorbar(ax=ax4)
-plt.show()
-#plt.savefig('TGV_phantom.pdf', format='pdf', dpi=1600)
-
-# calculate errors 
-Qtools = QualityTools(phantom, RecADMM_reg_tgv)
-RMSE_admm_tgv = Qtools.rmse()
-print("Root Mean Square Error for ADMM-TGV is {}".format(RMSE_admm_tgv))
-
-# SSIM measure
-#Create a 2d gaussian for the window parameter
-Qtools = QualityTools(phantom[128,:,:]*255, RecADMM_reg_tgv[128,:,:]*235)
-win = np.array([gaussian(11, 1.5)])
-win2d = win * (win.T)
-ssim_admm_tgv = Qtools.ssim(win2d)
-print("Mean SSIM ADMM-TGV is {}".format(ssim_admm_tgv[0]))
-#%%
\ No newline at end of file
diff --git a/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_SX.py b/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_SX.py
deleted file mode 100644
index cdf4325..0000000
--- a/Wrappers/Python/demos/SoftwareX_supp/Demo_SimulData_SX.py
+++ /dev/null
@@ -1,117 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-"""
-This demo scripts support the following publication: 
-"CCPi-Regularisation Toolkit for computed tomographic image reconstruction with 
-proximal splitting algorithms" by Daniil Kazantsev, Edoardo Pasca, Martin J. Turner,
- Philip J. Withers; Software X, 2019
-____________________________________________________________________________
-* Runs TomoPhantom software to simulate tomographic projection data with
-some imaging errors and noise
-* Saves the data into hdf file to be uploaded in reconstruction scripts
-__________________________________________________________________________
-
->>>>> Dependencies: <<<<<
-1. TomoPhantom software for phantom and data generation
-
-@author: Daniil Kazantsev, e:mail daniil.kazantsev@diamond.ac.uk
-Apache 2.0 license
-"""
-import timeit
-import os
-import matplotlib.pyplot as plt
-import numpy as np
-import tomophantom
-from tomophantom import TomoP3D
-from tomophantom.supp.flatsgen import flats
-from tomophantom.supp.normraw import normaliser_sim
-
-print ("Building 3D phantom using TomoPhantom software")
-tic=timeit.default_timer()
-model = 16 # select a model number from the library
-N_size = 256 # Define phantom dimensions using a scalar value (cubic phantom)
-path = os.path.dirname(tomophantom.__file__)
-path_library3D = os.path.join(path, "Phantom3DLibrary.dat")
-#This will generate a N_size x N_size x N_size phantom (3D)
-phantom_tm = TomoP3D.Model(model, N_size, path_library3D)
-toc=timeit.default_timer()
-Run_time = toc - tic
-print("Phantom has been built in {} seconds".format(Run_time))
-
-sliceSel = int(0.5*N_size)
-#plt.gray()
-plt.figure() 
-plt.subplot(131)
-plt.imshow(phantom_tm[sliceSel,:,:],vmin=0, vmax=1)
-plt.title('3D Phantom, axial view')
-
-plt.subplot(132)
-plt.imshow(phantom_tm[:,sliceSel,:],vmin=0, vmax=1)
-plt.title('3D Phantom, coronal view')
-
-plt.subplot(133)
-plt.imshow(phantom_tm[:,:,sliceSel],vmin=0, vmax=1)
-plt.title('3D Phantom, sagittal view')
-plt.show()
-
-# Projection geometry related parameters:
-Horiz_det = int(np.sqrt(2)*N_size) # detector column count (horizontal)
-Vert_det = N_size # detector row count (vertical) (no reason for it to be > N)
-angles_num = int(0.35*np.pi*N_size); # angles number
-angles = np.linspace(0.0,179.9,angles_num,dtype='float32') # in degrees
-angles_rad = angles*(np.pi/180.0)
-#%%
-print ("Building 3D analytical projection data with TomoPhantom")
-projData3D_analyt= TomoP3D.ModelSino(model, N_size, Horiz_det, Vert_det, angles, path_library3D)
-
-intens_max = N_size
-sliceSel = int(0.5*N_size)
-plt.figure() 
-plt.subplot(131)
-plt.imshow(projData3D_analyt[:,sliceSel,:],vmin=0, vmax=intens_max)
-plt.title('2D Projection (analytical)')
-plt.subplot(132)
-plt.imshow(projData3D_analyt[sliceSel,:,:],vmin=0, vmax=intens_max)
-plt.title('Sinogram view')
-plt.subplot(133)
-plt.imshow(projData3D_analyt[:,:,sliceSel],vmin=0, vmax=intens_max)
-plt.title('Tangentogram view')
-plt.show()
-#%%
-print ("Simulate flat fields, add noise and normalise projections...")
-flatsnum = 20 # generate 20 flat fields
-flatsSIM = flats(Vert_det, Horiz_det, maxheight = 0.1, maxthickness = 3, sigma_noise = 0.2, sigmasmooth = 3, flatsnum=flatsnum)
-
-plt.figure() 
-plt.imshow(flatsSIM[0,:,:],vmin=0, vmax=1)
-plt.title('A selected simulated flat-field')
-#%%
-# Apply normalisation of data and add noise
-flux_intensity = 60000 # controls the level of noise 
-sigma_flats = 0.01 # contro the level of noise in flats (higher creates more ring artifacts)
-projData3D_norm = normaliser_sim(projData3D_analyt, flatsSIM, sigma_flats, flux_intensity)
-
-intens_max = N_size
-sliceSel = int(0.5*N_size)
-plt.figure() 
-plt.subplot(131)
-plt.imshow(projData3D_norm[:,sliceSel,:],vmin=0, vmax=intens_max)
-plt.title('2D Projection (erroneous)')
-plt.subplot(132)
-plt.imshow(projData3D_norm[sliceSel,:,:],vmin=0, vmax=intens_max)
-plt.title('Sinogram view')
-plt.subplot(133)
-plt.imshow(projData3D_norm[:,:,sliceSel],vmin=0, vmax=intens_max)
-plt.title('Tangentogram view')
-plt.show()
-#%%
-import h5py
-import time
-time_label = int(time.time())
-# Saving generated data with a unique time label
-h5f = h5py.File('TomoSim_data'+str(time_label)+'.h5', 'w')
-h5f.create_dataset('phantom', data=phantom_tm)
-h5f.create_dataset('projdata_norm', data=projData3D_norm)
-h5f.create_dataset('proj_angles', data=angles_rad)
-h5f.close()
-#%%
\ No newline at end of file
diff --git a/Wrappers/Python/demos/SoftwareX_supp/Readme.md b/Wrappers/Python/demos/SoftwareX_supp/Readme.md
deleted file mode 100644
index 54e83f1..0000000
--- a/Wrappers/Python/demos/SoftwareX_supp/Readme.md
+++ /dev/null
@@ -1,26 +0,0 @@
-
-# SoftwareX publication [1] supporting files
-
-## Decription:
-The scripts here support publication in SoftwareX journal [1] to ensure reproducibility of the research. The scripts linked with data shared at Zenodo. 
-
-## Data:
-Data is shared at Zenodo [here](https://doi.org/10.5281/zenodo.2578893)
-
-## Dependencies:
-1. [ASTRA toolbox](https://github.com/astra-toolbox/astra-toolbox): `conda install -c astra-toolbox astra-toolbox`
-2. [TomoRec](https://github.com/dkazanc/TomoRec): `conda install -c dkazanc tomorec`
-3. [Tomophantom](https://github.com/dkazanc/TomoPhantom): `conda install tomophantom -c ccpi`
-
-## Files description: 
-- `Demo_SimulData_SX.py` - simulates 3D projection data using [Tomophantom](https://github.com/dkazanc/TomoPhantom) software. One can skip this module if the data is taken from [Zenodo](https://doi.org/10.5281/zenodo.2578893)
-- `Demo_SimulData_ParOptimis_SX.py` - runs computationally extensive calculations for optimal regularisation parameters, the result are saved into directory `optim_param`. This script can be also skipped. 
-- `Demo_SimulData_Recon_SX.py` - using established regularisation parameters, one runs iterative reconstruction
-- `Demo_RealData_Recon_SX.py` - runs real data reconstructions. Can be quite intense on memory so reduce the size of the reconstructed volume if needed. 
-
-### References:
-[1] "CCPi-Regularisation Toolkit for computed tomographic image reconstruction with proximal splitting algorithms" by Daniil Kazantsev, Edoardo Pasca, Martin J. Turner and Philip J. Withers; SoftwareX, 2019. 
-
-### Acknowledgments:
-CCPi-RGL software is a product of the [CCPi](https://www.ccpi.ac.uk/) group, STFC SCD software developers and Diamond Light Source (DLS). Any relevant questions/comments can be e-mailed to Daniil Kazantsev at dkazanc@hotmail.com
-
diff --git a/Wrappers/Python/demos/SoftwareX_supp/optim_param/Optim_admm_rofllt.h5 b/Wrappers/Python/demos/SoftwareX_supp/optim_param/Optim_admm_rofllt.h5
deleted file mode 100644
index 63bc4fd..0000000
--- a/Wrappers/Python/demos/SoftwareX_supp/optim_param/Optim_admm_rofllt.h5
+++ /dev/null
diff --git a/Wrappers/Python/demos/SoftwareX_supp/optim_param/Optim_admm_sbtv.h5 b/Wrappers/Python/demos/SoftwareX_supp/optim_param/Optim_admm_sbtv.h5
deleted file mode 100644
index 03c0c14..0000000
--- a/Wrappers/Python/demos/SoftwareX_supp/optim_param/Optim_admm_sbtv.h5
+++ /dev/null
diff --git a/Wrappers/Python/demos/SoftwareX_supp/optim_param/Optim_admm_tgv.h5 b/Wrappers/Python/demos/SoftwareX_supp/optim_param/Optim_admm_tgv.h5
deleted file mode 100644
index 056d915..0000000
--- a/Wrappers/Python/demos/SoftwareX_supp/optim_param/Optim_admm_tgv.h5
+++ /dev/null
diff --git a/Wrappers/Python/demos/demo_cpu_inpainters.py b/Wrappers/Python/demos/demo_cpu_inpainters.py
deleted file mode 100644
index c61ea50..0000000
--- a/Wrappers/Python/demos/demo_cpu_inpainters.py
+++ /dev/null
@@ -1,194 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-"""
-Demonstration of CPU inpainters
-@authors: Daniil Kazantsev, Edoardo Pasca
-"""
-
-import matplotlib.pyplot as plt
-import numpy as np
-import os
-import timeit
-from scipy import io
-from ccpi.filters.regularisers import NDF_INP, NVM_INP
-from ccpi.supp.qualitymetrics import QualityTools
-###############################################################################
-def printParametersToString(pars):
-        txt = r''
-        for key, value in pars.items():
-            if key== 'algorithm' :
-                txt += "{0} = {1}".format(key, value.__name__)
-            elif key == 'input':
-                txt += "{0} = {1}".format(key, np.shape(value))
-            elif key == 'maskData':
-                txt += "{0} = {1}".format(key, np.shape(value))
-            else:
-                txt += "{0} = {1}".format(key, value)
-            txt += '\n'
-        return txt
-###############################################################################
-
-# read sinogram and the mask
-filename = os.path.join(".." , ".." , ".." , "data" ,"SinoInpaint.mat")
-sino = io.loadmat(filename)
-sino_full = sino.get('Sinogram')
-Mask = sino.get('Mask')
-[angles_dim,detectors_dim] = sino_full.shape
-sino_full = sino_full/np.max(sino_full)
-#apply mask to sinogram
-sino_cut = sino_full*(1-Mask)
-#sino_cut_new = np.zeros((angles_dim,detectors_dim),'float32')
-#sino_cut_new = sino_cut.copy(order='c')
-#sino_cut_new[:] = sino_cut[:]
-sino_cut_new = np.ascontiguousarray(sino_cut, dtype=np.float32);
-#mask = np.zeros((angles_dim,detectors_dim),'uint8')
-#mask =Mask.copy(order='c')
-#mask[:] = Mask[:]
-mask = np.ascontiguousarray(Mask, dtype=np.uint8);
-
-plt.figure(1)
-plt.subplot(121)
-plt.imshow(sino_cut_new,vmin=0.0, vmax=1)
-plt.title('Missing Data sinogram')
-plt.subplot(122)
-plt.imshow(mask)
-plt.title('Mask')
-plt.show()
-#%%
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("___Inpainting using linear diffusion (2D)__")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-
-## plot 
-fig = plt.figure(2)
-plt.suptitle('Performance of linear inpainting using the CPU')
-a=fig.add_subplot(1,2,1)
-a.set_title('Missing data sinogram')
-imgplot = plt.imshow(sino_cut_new,cmap="gray")
-
-# set parameters
-pars = {'algorithm' : NDF_INP, \
-        'input' : sino_cut_new,\
-        'maskData' : mask,\
-        'regularisation_parameter':5000,\
-        'edge_parameter':0,\
-        'number_of_iterations' :5000 ,\
-        'time_marching_parameter':0.000075,\
-        'penalty_type':0
-        }
-        
-start_time = timeit.default_timer()
-ndf_inp_linear = NDF_INP(pars['input'],
-              pars['maskData'],
-              pars['regularisation_parameter'],
-              pars['edge_parameter'], 
-              pars['number_of_iterations'],
-              pars['time_marching_parameter'], 
-              pars['penalty_type'])
-
-Qtools = QualityTools(sino_full, ndf_inp_linear)
-pars['rmse'] = Qtools.rmse()
-
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(1,2,2)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
-# place a text box in upper left in axes coords
-a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
-         verticalalignment='top', bbox=props)
-imgplot = plt.imshow(ndf_inp_linear, cmap="gray")
-plt.title('{}'.format('Linear diffusion inpainting results'))
-#%%
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("_Inpainting using nonlinear diffusion (2D)_")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-
-## plot 
-fig = plt.figure(3)
-plt.suptitle('Performance of nonlinear diffusion inpainting using the CPU')
-a=fig.add_subplot(1,2,1)
-a.set_title('Missing data sinogram')
-imgplot = plt.imshow(sino_cut_new,cmap="gray")
-
-# set parameters
-pars = {'algorithm' : NDF_INP, \
-        'input' : sino_cut_new,\
-        'maskData' : mask,\
-        'regularisation_parameter':80,\
-        'edge_parameter':0.00009,\
-        'number_of_iterations' :1500 ,\
-        'time_marching_parameter':0.000008,\
-        'penalty_type':1
-        }
-        
-start_time = timeit.default_timer()
-ndf_inp_nonlinear = NDF_INP(pars['input'],
-              pars['maskData'],
-              pars['regularisation_parameter'],
-              pars['edge_parameter'], 
-              pars['number_of_iterations'],
-              pars['time_marching_parameter'], 
-              pars['penalty_type'])
-             
-
-Qtools = QualityTools(sino_full, ndf_inp_nonlinear)
-pars['rmse'] = Qtools.rmse()
-
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(1,2,2)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
-# place a text box in upper left in axes coords
-a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
-         verticalalignment='top', bbox=props)
-imgplot = plt.imshow(ndf_inp_nonlinear, cmap="gray")
-plt.title('{}'.format('Nonlinear diffusion inpainting results'))
-#%%
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("Inpainting using nonlocal vertical marching")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-
-## plot 
-fig = plt.figure(4)
-plt.suptitle('Performance of NVM inpainting using the CPU')
-a=fig.add_subplot(1,2,1)
-a.set_title('Missing data sinogram')
-imgplot = plt.imshow(sino_cut,cmap="gray")
-
-# set parameters
-pars = {'algorithm' : NVM_INP, \
-        'input' : sino_cut_new,\
-        'maskData' : mask,\
-        'SW_increment': 1,\
-        'number_of_iterations' : 150
-        }
-        
-start_time = timeit.default_timer()
-(nvm_inp, mask_upd) = NVM_INP(pars['input'],
-              pars['maskData'],
-              pars['SW_increment'],
-              pars['number_of_iterations'])
-             
-
-Qtools = QualityTools(sino_full, nvm_inp)
-pars['rmse'] = Qtools.rmse()
-
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(1,2,2)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
-# place a text box in upper left in axes coords
-a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
-         verticalalignment='top', bbox=props)
-imgplot = plt.imshow(nvm_inp, cmap="gray")
-plt.title('{}'.format('Nonlocal Vertical Marching inpainting results'))
-#%%
diff --git a/Wrappers/Python/demos/demo_cpu_regularisers.py b/Wrappers/Python/demos/demo_cpu_regularisers.py
deleted file mode 100644
index b8dadf5..0000000
--- a/Wrappers/Python/demos/demo_cpu_regularisers.py
+++ /dev/null
@@ -1,572 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-"""
-Created on Thu Feb 22 11:39:43 2018
-
-Demonstration of CPU regularisers 
-
-@authors: Daniil Kazantsev, Edoardo Pasca
-"""
-
-import matplotlib.pyplot as plt
-import numpy as np
-import os
-import timeit
-from ccpi.filters.regularisers import ROF_TV, FGP_TV, SB_TV, TGV, LLT_ROF, FGP_dTV, TNV, NDF, Diff4th
-from ccpi.filters.regularisers import PatchSelect, NLTV
-from ccpi.supp.qualitymetrics import QualityTools
-###############################################################################
-def printParametersToString(pars):
-        txt = r''
-        for key, value in pars.items():
-            if key== 'algorithm' :
-                txt += "{0} = {1}".format(key, value.__name__)
-            elif key == 'input':
-                txt += "{0} = {1}".format(key, np.shape(value))
-            elif key == 'refdata':
-                txt += "{0} = {1}".format(key, np.shape(value))
-            else:
-                txt += "{0} = {1}".format(key, value)
-            txt += '\n'
-        return txt
-###############################################################################
-#%%
-filename = os.path.join(".." , ".." , ".." , "data" ,"lena_gray_512.tif")
-
-# read image
-Im = plt.imread(filename)
-Im = np.asarray(Im, dtype='float32')
-
-Im = Im/255.0
-perc = 0.05
-u0 = Im + np.random.normal(loc = 0 ,
-                                  scale = perc * Im , 
-                                  size = np.shape(Im))
-u_ref = Im + np.random.normal(loc = 0 ,
-                                  scale = 0.01 * Im , 
-                                  size = np.shape(Im))
-(N,M) = np.shape(u0)
-# map the u0 u0->u0>0
-# f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1)
-u0 = u0.astype('float32')
-u_ref = u_ref.astype('float32')
-
-# change dims to check that modules work with non-squared images
-"""
-M = M-100
-u_ref2 = np.zeros([N,M],dtype='float32')
-u_ref2[:,0:M] = u_ref[:,0:M]
-u_ref = u_ref2
-del u_ref2
-
-u02 = np.zeros([N,M],dtype='float32')
-u02[:,0:M] = u0[:,0:M]
-u0 = u02
-del u02
-
-Im2 = np.zeros([N,M],dtype='float32')
-Im2[:,0:M] = Im[:,0:M]
-Im = Im2
-del Im2
-"""
-#%%
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("_______________ROF-TV (2D)_________________")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-
-## plot 
-fig = plt.figure()
-plt.suptitle('Performance of ROF-TV regulariser using the CPU')
-a=fig.add_subplot(1,2,1)
-a.set_title('Noisy Image')
-imgplot = plt.imshow(u0,cmap="gray")
-
-# set parameters
-pars = {'algorithm': ROF_TV, \
-        'input' : u0,\
-        'regularisation_parameter':0.02,\
-        'number_of_iterations': 2000,\
-        'time_marching_parameter': 0.0025
-        }
-print ("#############ROF TV CPU####################")
-start_time = timeit.default_timer()
-rof_cpu = ROF_TV(pars['input'],
-             pars['regularisation_parameter'],
-             pars['number_of_iterations'],
-             pars['time_marching_parameter'],'cpu')
-
-Qtools = QualityTools(Im, rof_cpu)
-pars['rmse'] = Qtools.rmse()
-
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(1,2,2)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
-# place a text box in upper left in axes coords
-a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
-         verticalalignment='top', bbox=props)
-imgplot = plt.imshow(rof_cpu, cmap="gray")
-plt.title('{}'.format('CPU results'))
-
-#%%
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("_______________FGP-TV (2D)__________________")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-
-## plot 
-fig = plt.figure()
-plt.suptitle('Performance of FGP-TV regulariser using the CPU')
-a=fig.add_subplot(1,2,1)
-a.set_title('Noisy Image')
-imgplot = plt.imshow(u0,cmap="gray")
-
-# set parameters
-pars = {'algorithm' : FGP_TV, \
-        'input' : u0,\
-        'regularisation_parameter':0.04, \
-        'number_of_iterations' :2000 ,\
-        'tolerance_constant':1e-06,\
-        'methodTV': 0 ,\
-        'nonneg': 0 ,\
-        'printingOut': 0 
-        }
-        
-print ("#############FGP TV CPU####################")
-start_time = timeit.default_timer()
-fgp_cpu = FGP_TV(pars['input'], 
-              pars['regularisation_parameter'],
-              pars['number_of_iterations'],
-              pars['tolerance_constant'], 
-              pars['methodTV'],
-              pars['nonneg'],
-              pars['printingOut'],'cpu')  
-             
-
-Qtools = QualityTools(Im, fgp_cpu)
-pars['rmse'] = Qtools.rmse()
-
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(1,2,2)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
-# place a text box in upper left in axes coords
-a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
-         verticalalignment='top', bbox=props)
-imgplot = plt.imshow(fgp_cpu, cmap="gray")
-plt.title('{}'.format('CPU results'))
-
-#%%
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("_______________SB-TV (2D)__________________")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-
-## plot 
-fig = plt.figure()
-plt.suptitle('Performance of SB-TV regulariser using the CPU')
-a=fig.add_subplot(1,2,1)
-a.set_title('Noisy Image')
-imgplot = plt.imshow(u0,cmap="gray")
-
-# set parameters
-pars = {'algorithm' : SB_TV, \
-        'input' : u0,\
-        'regularisation_parameter':0.04, \
-        'number_of_iterations' :150 ,\
-        'tolerance_constant':1e-06,\
-        'methodTV': 0 ,\
-        'printingOut': 0 
-        }
-        
-print ("#############SB TV CPU####################")
-start_time = timeit.default_timer()
-sb_cpu = SB_TV(pars['input'], 
-              pars['regularisation_parameter'],
-              pars['number_of_iterations'],
-              pars['tolerance_constant'], 
-              pars['methodTV'],
-              pars['printingOut'],'cpu')  
-             
-Qtools = QualityTools(Im, sb_cpu)
-pars['rmse'] = Qtools.rmse()
-
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(1,2,2)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
-# place a text box in upper left in axes coords
-a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
-         verticalalignment='top', bbox=props)
-imgplot = plt.imshow(sb_cpu, cmap="gray")
-plt.title('{}'.format('CPU results'))
-#%%
-
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("_____Total Generalised Variation (2D)______")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-
-## plot 
-fig = plt.figure()
-plt.suptitle('Performance of TGV regulariser using the CPU')
-a=fig.add_subplot(1,2,1)
-a.set_title('Noisy Image')
-imgplot = plt.imshow(u0,cmap="gray")
-
-# set parameters
-pars = {'algorithm' : TGV, \
-        'input' : u0,\
-        'regularisation_parameter':0.04, \
-        'alpha1':1.0,\
-        'alpha0':2.0,\
-        'number_of_iterations' :1350 ,\
-        'LipshitzConstant' :12 ,\
-        }
-        
-print ("#############TGV CPU####################")
-start_time = timeit.default_timer()
-tgv_cpu = TGV(pars['input'], 
-              pars['regularisation_parameter'],
-              pars['alpha1'],
-              pars['alpha0'],
-              pars['number_of_iterations'],
-              pars['LipshitzConstant'],'cpu')
-             
-             
-Qtools = QualityTools(Im, tgv_cpu)
-pars['rmse'] = Qtools.rmse()
-
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(1,2,2)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
-# place a text box in upper left in axes coords
-a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
-         verticalalignment='top', bbox=props)
-imgplot = plt.imshow(tgv_cpu, cmap="gray")
-plt.title('{}'.format('CPU results'))
-
-#%%
-
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("______________LLT- ROF (2D)________________")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-
-## plot 
-fig = plt.figure()
-plt.suptitle('Performance of LLT-ROF regulariser using the CPU')
-a=fig.add_subplot(1,2,1)
-a.set_title('Noisy Image')
-imgplot = plt.imshow(u0,cmap="gray")
-
-# set parameters
-pars = {'algorithm' : LLT_ROF, \
-        'input' : u0,\
-        'regularisation_parameterROF':0.04, \
-        'regularisation_parameterLLT':0.01, \
-        'number_of_iterations' :500 ,\
-        'time_marching_parameter' :0.0025 ,\
-        }
-        
-print ("#############LLT- ROF CPU####################")
-start_time = timeit.default_timer()
-lltrof_cpu = LLT_ROF(pars['input'], 
-              pars['regularisation_parameterROF'],
-              pars['regularisation_parameterLLT'],
-              pars['number_of_iterations'],
-              pars['time_marching_parameter'],'cpu')
-
-Qtools = QualityTools(Im, lltrof_cpu)
-pars['rmse'] = Qtools.rmse()
-
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(1,2,2)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
-# place a text box in upper left in axes coords
-a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
-         verticalalignment='top', bbox=props)
-imgplot = plt.imshow(lltrof_cpu, cmap="gray")
-plt.title('{}'.format('CPU results'))
-
-#%%
-
-
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("________________NDF (2D)___________________")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-
-## plot 
-fig = plt.figure()
-plt.suptitle('Performance of NDF regulariser using the CPU')
-a=fig.add_subplot(1,2,1)
-a.set_title('Noisy Image')
-imgplot = plt.imshow(u0,cmap="gray")
-
-# set parameters
-pars = {'algorithm' : NDF, \
-        'input' : u0,\
-        'regularisation_parameter':0.025, \
-        'edge_parameter':0.015,\
-        'number_of_iterations' :500 ,\
-        'time_marching_parameter':0.025,\
-        'penalty_type':1
-        }
-        
-print ("#############NDF CPU################")
-start_time = timeit.default_timer()
-ndf_cpu = NDF(pars['input'], 
-              pars['regularisation_parameter'],
-              pars['edge_parameter'], 
-              pars['number_of_iterations'],
-              pars['time_marching_parameter'], 
-              pars['penalty_type'],'cpu')  
-             
-Qtools = QualityTools(Im, ndf_cpu)
-pars['rmse'] = Qtools.rmse()
-
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(1,2,2)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
-# place a text box in upper left in axes coords
-a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
-         verticalalignment='top', bbox=props)
-imgplot = plt.imshow(ndf_cpu, cmap="gray")
-plt.title('{}'.format('CPU results'))
-
-#%%
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("___Anisotropic Diffusion 4th Order (2D)____")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-
-## plot 
-fig = plt.figure()
-plt.suptitle('Performance of Diff4th regulariser using the CPU')
-a=fig.add_subplot(1,2,1)
-a.set_title('Noisy Image')
-imgplot = plt.imshow(u0,cmap="gray")
-
-# set parameters
-pars = {'algorithm' : Diff4th, \
-        'input' : u0,\
-        'regularisation_parameter':3.5, \
-        'edge_parameter':0.02,\
-        'number_of_iterations' :500 ,\
-        'time_marching_parameter':0.0015
-        }
-        
-print ("#############Diff4th CPU################")
-start_time = timeit.default_timer()
-diff4_cpu = Diff4th(pars['input'], 
-              pars['regularisation_parameter'],
-              pars['edge_parameter'], 
-              pars['number_of_iterations'],
-              pars['time_marching_parameter'],'cpu')
-             
-Qtools = QualityTools(Im, diff4_cpu)
-pars['rmse'] = Qtools.rmse()
-
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(1,2,2)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
-# place a text box in upper left in axes coords
-a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
-         verticalalignment='top', bbox=props)
-imgplot = plt.imshow(diff4_cpu, cmap="gray")
-plt.title('{}'.format('CPU results'))
-
-#%%
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("___Nonlocal patches pre-calculation____")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-start_time = timeit.default_timer()
-# set parameters
-pars = {'algorithm' : PatchSelect, \
-        'input' : u0,\
-        'searchwindow': 7, \
-        'patchwindow': 2,\
-        'neighbours' : 15 ,\
-        'edge_parameter':0.18}
-
-H_i, H_j, Weights = PatchSelect(pars['input'], 
-              pars['searchwindow'],
-              pars['patchwindow'], 
-              pars['neighbours'],
-              pars['edge_parameter'],'cpu')
-              
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-"""
-plt.figure()
-plt.imshow(Weights[0,:,:],cmap="gray",interpolation="nearest",vmin=0, vmax=1)
-plt.show()
-"""
-#%%
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("___Nonlocal Total Variation penalty____")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-## plot 
-fig = plt.figure()
-plt.suptitle('Performance of NLTV regulariser using the CPU')
-a=fig.add_subplot(1,2,1)
-a.set_title('Noisy Image')
-imgplot = plt.imshow(u0,cmap="gray")
-
-pars2 = {'algorithm' : NLTV, \
-        'input' : u0,\
-        'H_i': H_i, \
-        'H_j': H_j,\
-        'H_k' : 0,\
-        'Weights' : Weights,\
-        'regularisation_parameter': 0.04,\
-        'iterations': 3
-        }
-start_time = timeit.default_timer()
-nltv_cpu = NLTV(pars2['input'], 
-              pars2['H_i'],
-              pars2['H_j'], 
-              pars2['H_k'],
-              pars2['Weights'],
-              pars2['regularisation_parameter'],
-              pars2['iterations'])
-
-Qtools = QualityTools(Im, nltv_cpu)
-pars['rmse'] = Qtools.rmse()
-
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(1,2,2)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
-# place a text box in upper left in axes coords
-a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
-         verticalalignment='top', bbox=props)
-imgplot = plt.imshow(nltv_cpu, cmap="gray")
-plt.title('{}'.format('CPU results'))
-
-#%%
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("_____________FGP-dTV (2D)__________________")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-
-## plot 
-fig = plt.figure()
-plt.suptitle('Performance of FGP-dTV regulariser using the CPU')
-a=fig.add_subplot(1,2,1)
-a.set_title('Noisy Image')
-imgplot = plt.imshow(u0,cmap="gray")
-
-# set parameters
-pars = {'algorithm' : FGP_dTV, \
-        'input' : u0,\
-        'refdata' : u_ref,\
-        'regularisation_parameter':0.04, \
-        'number_of_iterations' :2000 ,\
-        'tolerance_constant':1e-06,\
-        'eta_const':0.2,\
-        'methodTV': 0 ,\
-        'nonneg': 0 ,\
-        'printingOut': 0 
-        }
-        
-print ("#############FGP dTV CPU####################")
-start_time = timeit.default_timer()
-fgp_dtv_cpu = FGP_dTV(pars['input'], 
-              pars['refdata'], 
-              pars['regularisation_parameter'],
-              pars['number_of_iterations'],
-              pars['tolerance_constant'], 
-              pars['eta_const'], 
-              pars['methodTV'],
-              pars['nonneg'],
-              pars['printingOut'],'cpu')
-             
-Qtools = QualityTools(Im, fgp_dtv_cpu)
-pars['rmse'] = Qtools.rmse()
-
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(1,2,2)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
-# place a text box in upper left in axes coords
-a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
-         verticalalignment='top', bbox=props)
-imgplot = plt.imshow(fgp_dtv_cpu, cmap="gray")
-plt.title('{}'.format('CPU results'))
-#%%
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("__________Total nuclear Variation__________")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-
-## plot 
-fig = plt.figure()
-plt.suptitle('Performance of TNV regulariser using the CPU')
-a=fig.add_subplot(1,2,1)
-a.set_title('Noisy Image')
-imgplot = plt.imshow(u0,cmap="gray")
-
-channelsNo = 5
-noisyVol = np.zeros((channelsNo,N,M),dtype='float32')
-idealVol = np.zeros((channelsNo,N,M),dtype='float32')
-
-for i in range (channelsNo):
-    noisyVol[i,:,:] = Im + np.random.normal(loc = 0 , scale = perc * Im , size = np.shape(Im))
-    idealVol[i,:,:] = Im
-
-# set parameters
-pars = {'algorithm' : TNV, \
-        'input' : noisyVol,\
-        'regularisation_parameter': 0.04, \
-        'number_of_iterations' : 200 ,\
-        'tolerance_constant':1e-05
-        }
-        
-print ("#############TNV CPU#################")
-start_time = timeit.default_timer()
-tnv_cpu = TNV(pars['input'],           
-              pars['regularisation_parameter'],
-              pars['number_of_iterations'],
-              pars['tolerance_constant'])
-
-Qtools = QualityTools(idealVol, tnv_cpu)
-pars['rmse'] = Qtools.rmse()
-
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(1,2,2)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
-# place a text box in upper left in axes coords
-a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
-         verticalalignment='top', bbox=props)
-imgplot = plt.imshow(tnv_cpu[3,:,:], cmap="gray")
-plt.title('{}'.format('CPU results'))
diff --git a/Wrappers/Python/demos/demo_cpu_regularisers3D.py b/Wrappers/Python/demos/demo_cpu_regularisers3D.py
deleted file mode 100644
index df8af27..0000000
--- a/Wrappers/Python/demos/demo_cpu_regularisers3D.py
+++ /dev/null
@@ -1,463 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-"""
-Created on Thu Feb 22 11:39:43 2018
-
-Demonstration of 3D CPU regularisers 
-
-@authors: Daniil Kazantsev, Edoardo Pasca
-"""
-
-import matplotlib.pyplot as plt
-import numpy as np
-import os
-import timeit
-from ccpi.filters.regularisers import ROF_TV, FGP_TV, SB_TV, TGV, LLT_ROF, FGP_dTV, NDF, Diff4th
-from ccpi.supp.qualitymetrics import QualityTools
-###############################################################################
-def printParametersToString(pars):
-        txt = r''
-        for key, value in pars.items():
-            if key== 'algorithm' :
-                txt += "{0} = {1}".format(key, value.__name__)
-            elif key == 'input':
-                txt += "{0} = {1}".format(key, np.shape(value))
-            elif key == 'refdata':
-                txt += "{0} = {1}".format(key, np.shape(value))
-            else:
-                txt += "{0} = {1}".format(key, value)
-            txt += '\n'
-        return txt
-###############################################################################
-#%%
-filename = os.path.join(".." , ".." , ".." , "data" ,"lena_gray_512.tif")
-
-# read image
-Im = plt.imread(filename)
-Im = np.asarray(Im, dtype='float32')
-
-Im = Im/255
-perc = 0.05
-u0 = Im + np.random.normal(loc = 0 ,
-                                  scale = perc * Im , 
-                                  size = np.shape(Im))
-u_ref = Im + np.random.normal(loc = 0 ,
-                                  scale = 0.01 * Im , 
-                                  size = np.shape(Im))
-(N,M) = np.shape(u0)
-# map the u0 u0->u0>0
-# f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1)
-u0 = u0.astype('float32')
-u_ref = u_ref.astype('float32')
-
-# change dims to check that modules work with non-squared images
-"""
-M = M-100
-u_ref2 = np.zeros([N,M],dtype='float32')
-u_ref2[:,0:M] = u_ref[:,0:M]
-u_ref = u_ref2
-del u_ref2
-
-u02 = np.zeros([N,M],dtype='float32')
-u02[:,0:M] = u0[:,0:M]
-u0 = u02
-del u02
-
-Im2 = np.zeros([N,M],dtype='float32')
-Im2[:,0:M] = Im[:,0:M]
-Im = Im2
-del Im2
-"""
-slices = 15
-
-noisyVol = np.zeros((slices,N,M),dtype='float32')
-noisyRef = np.zeros((slices,N,M),dtype='float32')
-idealVol = np.zeros((slices,N,M),dtype='float32')
-
-for i in range (slices):
-    noisyVol[i,:,:] = Im + np.random.normal(loc = 0 , scale = perc * Im , size = np.shape(Im))
-    noisyRef[i,:,:] = Im + np.random.normal(loc = 0 , scale = 0.01 * Im , size = np.shape(Im))
-    idealVol[i,:,:] = Im
-
-#%%
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("_______________ROF-TV (3D)_________________")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-
-## plot 
-fig = plt.figure()
-plt.suptitle('Performance of ROF-TV regulariser using the CPU')
-a=fig.add_subplot(1,2,1)
-a.set_title('Noisy 15th slice of a volume')
-imgplot = plt.imshow(noisyVol[10,:,:],cmap="gray")
-
-# set parameters
-pars = {'algorithm': ROF_TV, \
-        'input' : noisyVol,\
-        'regularisation_parameter':0.04,\
-        'number_of_iterations': 500,\
-        'time_marching_parameter': 0.0025
-        }
-print ("#############ROF TV CPU####################")
-start_time = timeit.default_timer()
-rof_cpu3D = ROF_TV(pars['input'],
-             pars['regularisation_parameter'],
-             pars['number_of_iterations'],
-             pars['time_marching_parameter'],'cpu')
-
-Qtools = QualityTools(idealVol, rof_cpu3D)
-pars['rmse'] = Qtools.rmse()
-
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(1,2,2)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
-# place a text box in upper left in axes coords
-a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
-         verticalalignment='top', bbox=props)
-imgplot = plt.imshow(rof_cpu3D[10,:,:], cmap="gray")
-plt.title('{}'.format('Recovered volume on the CPU using ROF-TV'))
-
-#%%
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("_______________FGP-TV (3D)__________________")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-
-## plot 
-fig = plt.figure()
-plt.suptitle('Performance of FGP-TV regulariser using the CPU')
-a=fig.add_subplot(1,2,1)
-a.set_title('Noisy Image')
-imgplot = plt.imshow(noisyVol[10,:,:],cmap="gray")
-
-# set parameters
-pars = {'algorithm' : FGP_TV, \
-        'input' : noisyVol,\
-        'regularisation_parameter':0.04, \
-        'number_of_iterations' :300 ,\
-        'tolerance_constant':0.00001,\
-        'methodTV': 0 ,\
-        'nonneg': 0 ,\
-        'printingOut': 0 
-        }
-        
-print ("#############FGP TV CPU####################")
-start_time = timeit.default_timer()
-fgp_cpu3D = FGP_TV(pars['input'], 
-              pars['regularisation_parameter'],
-              pars['number_of_iterations'],
-              pars['tolerance_constant'], 
-              pars['methodTV'],
-              pars['nonneg'],
-              pars['printingOut'],'cpu')  
-             
-Qtools = QualityTools(idealVol, fgp_cpu3D)
-pars['rmse'] = Qtools.rmse()
-
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(1,2,2)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
-# place a text box in upper left in axes coords
-a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
-         verticalalignment='top', bbox=props)
-imgplot = plt.imshow(fgp_cpu3D[10,:,:], cmap="gray")
-plt.title('{}'.format('Recovered volume on the CPU using FGP-TV'))
-
-#%%
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("_______________SB-TV (3D)_________________")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-
-## plot 
-fig = plt.figure()
-plt.suptitle('Performance of SB-TV regulariser using the CPU')
-a=fig.add_subplot(1,2,1)
-a.set_title('Noisy Image')
-imgplot = plt.imshow(noisyVol[10,:,:],cmap="gray")
-
-# set parameters
-pars = {'algorithm' : SB_TV, \
-        'input' : noisyVol,\
-        'regularisation_parameter':0.04, \
-        'number_of_iterations' :150 ,\
-        'tolerance_constant':0.00001,\
-        'methodTV': 0 ,\
-        'printingOut': 0 
-        }
-        
-print ("#############SB TV CPU####################")
-start_time = timeit.default_timer()
-sb_cpu3D = SB_TV(pars['input'], 
-              pars['regularisation_parameter'],
-              pars['number_of_iterations'],
-              pars['tolerance_constant'], 
-              pars['methodTV'],
-              pars['printingOut'],'cpu')
-             
-
-Qtools = QualityTools(idealVol, sb_cpu3D)
-pars['rmse'] = Qtools.rmse()
-
-
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(1,2,2)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
-# place a text box in upper left in axes coords
-a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
-         verticalalignment='top', bbox=props)
-imgplot = plt.imshow(sb_cpu3D[10,:,:], cmap="gray")
-plt.title('{}'.format('Recovered volume on the CPU using SB-TV'))
-
-#%%
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("_______________LLT-ROF (3D)_________________")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-
-## plot 
-fig = plt.figure()
-plt.suptitle('Performance of LLT-ROF regulariser using the CPU')
-a=fig.add_subplot(1,2,1)
-a.set_title('Noisy Image')
-imgplot = plt.imshow(noisyVol[10,:,:],cmap="gray")
-
-# set parameters
-pars = {'algorithm' : LLT_ROF, \
-        'input' : noisyVol,\
-        'regularisation_parameterROF':0.04, \
-        'regularisation_parameterLLT':0.015, \
-        'number_of_iterations' :300 ,\
-        'time_marching_parameter' :0.0025 ,\
-        }
-
-print ("#############LLT ROF CPU####################")
-start_time = timeit.default_timer()
-lltrof_cpu3D = LLT_ROF(pars['input'], 
-              pars['regularisation_parameterROF'],
-              pars['regularisation_parameterLLT'],
-              pars['number_of_iterations'],
-              pars['time_marching_parameter'],'cpu')
-
-
-Qtools = QualityTools(idealVol, lltrof_cpu3D)
-pars['rmse'] = Qtools.rmse()
-
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(1,2,2)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
-# place a text box in upper left in axes coords
-a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
-         verticalalignment='top', bbox=props)
-imgplot = plt.imshow(lltrof_cpu3D[10,:,:], cmap="gray")
-plt.title('{}'.format('Recovered volume on the CPU using LLT-ROF'))
-
-#%%
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("_______________TGV (3D)_________________")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-
-## plot 
-fig = plt.figure()
-plt.suptitle('Performance of TGV regulariser using the CPU')
-a=fig.add_subplot(1,2,1)
-a.set_title('Noisy Image')
-imgplot = plt.imshow(noisyVol[10,:,:],cmap="gray")
-
-# set parameters
-pars = {'algorithm' : TGV, \
-        'input' : noisyVol,\
-        'regularisation_parameter':0.04, \
-        'alpha1':1.0,\
-        'alpha0':2.0,\
-        'number_of_iterations' :250 ,\
-        'LipshitzConstant' :12 ,\
-        }
-
-print ("#############TGV CPU####################")
-start_time = timeit.default_timer()
-tgv_cpu3D = TGV(pars['input'], 
-              pars['regularisation_parameter'],
-              pars['alpha1'],
-              pars['alpha0'],
-              pars['number_of_iterations'],
-              pars['LipshitzConstant'],'cpu')
-             
-
-Qtools = QualityTools(idealVol, tgv_cpu3D)
-pars['rmse'] = Qtools.rmse()
-
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(1,2,2)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
-# place a text box in upper left in axes coords
-a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
-         verticalalignment='top', bbox=props)
-imgplot = plt.imshow(tgv_cpu3D[10,:,:], cmap="gray")
-plt.title('{}'.format('Recovered volume on the CPU using TGV'))
-
-#%%
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("________________NDF (3D)___________________")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-
-## plot 
-fig = plt.figure()
-plt.suptitle('Performance of NDF regulariser using the CPU')
-a=fig.add_subplot(1,2,1)
-a.set_title('Noisy volume')
-imgplot = plt.imshow(noisyVol[10,:,:],cmap="gray")
-
-# set parameters
-pars = {'algorithm' : NDF, \
-        'input' : noisyVol,\
-        'regularisation_parameter':0.025, \
-        'edge_parameter':0.015,\
-        'number_of_iterations' :500 ,\
-        'time_marching_parameter':0.025,\
-        'penalty_type':  1
-        }
-        
-print ("#############NDF CPU################")
-start_time = timeit.default_timer()
-ndf_cpu3D = NDF(pars['input'], 
-              pars['regularisation_parameter'],
-              pars['edge_parameter'], 
-              pars['number_of_iterations'],
-              pars['time_marching_parameter'], 
-              pars['penalty_type'])  
-             
-
-Qtools = QualityTools(idealVol, ndf_cpu3D)
-pars['rmse'] = Qtools.rmse()
-
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(1,2,2)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
-# place a text box in upper left in axes coords
-a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
-         verticalalignment='top', bbox=props)
-imgplot = plt.imshow(ndf_cpu3D[10,:,:], cmap="gray")
-plt.title('{}'.format('Recovered volume on the CPU using NDF iterations'))
-
-#%%
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("___Anisotropic Diffusion 4th Order (2D)____")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-
-## plot 
-fig = plt.figure()
-plt.suptitle('Performance of Diff4th regulariser using the CPU')
-a=fig.add_subplot(1,2,1)
-a.set_title('Noisy volume')
-imgplot = plt.imshow(noisyVol[10,:,:],cmap="gray")
-
-# set parameters
-pars = {'algorithm' : Diff4th, \
-        'input' : noisyVol,\
-        'regularisation_parameter':3.5, \
-        'edge_parameter':0.02,\
-        'number_of_iterations' :300 ,\
-        'time_marching_parameter':0.0015
-        }
-        
-print ("#############Diff4th CPU################")
-start_time = timeit.default_timer()
-diff4th_cpu3D = Diff4th(pars['input'], 
-              pars['regularisation_parameter'],
-              pars['edge_parameter'], 
-              pars['number_of_iterations'],
-              pars['time_marching_parameter'])  
-             
-
-Qtools = QualityTools(idealVol, diff4th_cpu3D)
-pars['rmse'] = Qtools.rmse()
-
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(1,2,2)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
-# place a text box in upper left in axes coords
-a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
-         verticalalignment='top', bbox=props)
-imgplot = plt.imshow(diff4th_cpu3D[10,:,:], cmap="gray")
-plt.title('{}'.format('Recovered volume on the CPU using DIFF4th iterations'))
-
-#%%
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("_______________FGP-dTV (3D)__________________")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-
-## plot 
-fig = plt.figure()
-plt.suptitle('Performance of FGP-dTV regulariser using the CPU')
-a=fig.add_subplot(1,2,1)
-a.set_title('Noisy Image')
-imgplot = plt.imshow(noisyVol[10,:,:],cmap="gray")
-
-# set parameters
-pars = {'algorithm' : FGP_dTV,\
-        'input' : noisyVol,\
-        'refdata' : noisyRef,\
-        'regularisation_parameter':0.04, \
-        'number_of_iterations' :300 ,\
-        'tolerance_constant':0.00001,\
-        'eta_const':0.2,\
-        'methodTV': 0 ,\
-        'nonneg': 0 ,\
-        'printingOut': 0 
-        }
-        
-print ("#############FGP dTV CPU####################")
-start_time = timeit.default_timer()
-fgp_dTV_cpu3D = FGP_dTV(pars['input'],
-              pars['refdata'], 
-              pars['regularisation_parameter'],
-              pars['number_of_iterations'],
-              pars['tolerance_constant'], 
-              pars['eta_const'],
-              pars['methodTV'],
-              pars['nonneg'],
-              pars['printingOut'],'cpu')
-             
-
-Qtools = QualityTools(idealVol, fgp_dTV_cpu3D)
-pars['rmse'] = Qtools.rmse()
-
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(1,2,2)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
-# place a text box in upper left in axes coords
-a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
-         verticalalignment='top', bbox=props)
-imgplot = plt.imshow(fgp_dTV_cpu3D[10,:,:], cmap="gray")
-plt.title('{}'.format('Recovered volume on the CPU using FGP-dTV'))
-#%%
diff --git a/Wrappers/Python/demos/demo_cpu_vs_gpu_regularisers.py b/Wrappers/Python/demos/demo_cpu_vs_gpu_regularisers.py
deleted file mode 100644
index 6c4ab5e..0000000
--- a/Wrappers/Python/demos/demo_cpu_vs_gpu_regularisers.py
+++ /dev/null
@@ -1,794 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-"""
-Created on Thu Feb 22 11:39:43 2018
-
-Demonstration of CPU implementation against the GPU one
-
-@authors: Daniil Kazantsev, Edoardo Pasca
-"""
-
-import matplotlib.pyplot as plt
-import numpy as np
-import os
-import timeit
-from ccpi.filters.regularisers import ROF_TV, FGP_TV, SB_TV, TGV, LLT_ROF, FGP_dTV, NDF, Diff4th
-from ccpi.filters.regularisers import PatchSelect
-from ccpi.supp.qualitymetrics import QualityTools
-###############################################################################
-def printParametersToString(pars):
-        txt = r''
-        for key, value in pars.items():
-            if key== 'algorithm' :
-                txt += "{0} = {1}".format(key, value.__name__)
-            elif key == 'input':
-                txt += "{0} = {1}".format(key, np.shape(value))
-            elif key == 'refdata':
-                txt += "{0} = {1}".format(key, np.shape(value))
-            else:
-                txt += "{0} = {1}".format(key, value)
-            txt += '\n'
-        return txt
-###############################################################################
-
-filename = os.path.join(".." , ".." , ".." , "data" ,"lena_gray_512.tif")
-
-# read image
-Im = plt.imread(filename)                     
-Im = np.asarray(Im, dtype='float32')
-
-Im = Im/255
-perc = 0.05
-u0 = Im + np.random.normal(loc = 0 ,
-                                  scale = perc * Im , 
-                                  size = np.shape(Im))
-u_ref = Im + np.random.normal(loc = 0 ,
-                                  scale = 0.01 * Im , 
-                                  size = np.shape(Im))
-
-# map the u0 u0->u0>0
-# f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1)
-u0 = u0.astype('float32')
-u_ref = u_ref.astype('float32')
-
-#%%
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("____________ROF-TV bench___________________")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-
-## plot 
-fig = plt.figure()
-plt.suptitle('Comparison of ROF-TV regulariser using CPU and GPU implementations')
-a=fig.add_subplot(1,4,1)
-a.set_title('Noisy Image')
-imgplot = plt.imshow(u0,cmap="gray")
-
-# set parameters
-pars = {'algorithm': ROF_TV, \
-        'input' : u0,\
-        'regularisation_parameter':0.04,\
-        'number_of_iterations': 4500,\
-        'time_marching_parameter': 0.00002
-        }
-print ("#############ROF TV CPU####################")
-start_time = timeit.default_timer()
-rof_cpu = ROF_TV(pars['input'],
-             pars['regularisation_parameter'],
-             pars['number_of_iterations'],
-             pars['time_marching_parameter'],'cpu')
-
-Qtools = QualityTools(Im, rof_cpu)
-pars['rmse'] = Qtools.rmse()
-
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(1,4,2)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
-# place a text box in upper left in axes coords
-a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
-         verticalalignment='top', bbox=props)
-imgplot = plt.imshow(rof_cpu, cmap="gray")
-plt.title('{}'.format('CPU results'))
-
-print ("##############ROF TV GPU##################")
-start_time = timeit.default_timer()
-rof_gpu = ROF_TV(pars['input'], 
-                     pars['regularisation_parameter'],
-                     pars['number_of_iterations'], 
-                     pars['time_marching_parameter'],'gpu')
-
-Qtools = QualityTools(Im, rof_gpu)
-pars['rmse'] = Qtools.rmse()
-
-pars['algorithm'] = ROF_TV
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(1,4,3)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
-# place a text box in upper left in axes coords
-a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
-         verticalalignment='top', bbox=props)
-imgplot = plt.imshow(rof_gpu, cmap="gray")
-plt.title('{}'.format('GPU results'))
-
-
-print ("--------Compare the results--------")
-tolerance = 1e-05
-diff_im = np.zeros(np.shape(rof_cpu))
-diff_im = abs(rof_cpu - rof_gpu)
-diff_im[diff_im > tolerance] = 1
-a=fig.add_subplot(1,4,4)
-imgplot = plt.imshow(diff_im, vmin=0, vmax=1, cmap="gray")
-plt.title('{}'.format('Pixels larger threshold difference'))
-if (diff_im.sum() > 1):
-    print ("Arrays do not match!")
-else:
-    print ("Arrays match")
-
-#%%
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("____________FGP-TV bench___________________")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-
-## plot 
-fig = plt.figure()
-plt.suptitle('Comparison of FGP-TV regulariser using CPU and GPU implementations')
-a=fig.add_subplot(1,4,1)
-a.set_title('Noisy Image')
-imgplot = plt.imshow(u0,cmap="gray")
-
-# set parameters
-pars = {'algorithm' : FGP_TV, \
-        'input' : u0,\
-        'regularisation_parameter':0.04, \
-        'number_of_iterations' :1200 ,\
-        'tolerance_constant':0.00001,\
-        'methodTV': 0 ,\
-        'nonneg': 0 ,\
-        'printingOut': 0 
-        }
-        
-print ("#############FGP TV CPU####################")
-start_time = timeit.default_timer()
-fgp_cpu = FGP_TV(pars['input'], 
-              pars['regularisation_parameter'],
-              pars['number_of_iterations'],
-              pars['tolerance_constant'], 
-              pars['methodTV'],
-              pars['nonneg'],
-              pars['printingOut'],'cpu')  
-             
-
-Qtools = QualityTools(Im, fgp_cpu)
-pars['rmse'] = Qtools.rmse()
-
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(1,4,2)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
-# place a text box in upper left in axes coords
-a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
-         verticalalignment='top', bbox=props)
-imgplot = plt.imshow(fgp_cpu, cmap="gray")
-plt.title('{}'.format('CPU results'))
-
-
-print ("##############FGP TV GPU##################")
-start_time = timeit.default_timer()
-fgp_gpu = FGP_TV(pars['input'], 
-              pars['regularisation_parameter'],
-              pars['number_of_iterations'],
-              pars['tolerance_constant'], 
-              pars['methodTV'],
-              pars['nonneg'],
-              pars['printingOut'],'gpu')
-
-Qtools = QualityTools(Im, fgp_gpu)
-pars['rmse'] = Qtools.rmse()
-
-pars['algorithm'] = FGP_TV
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(1,4,3)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
-# place a text box in upper left in axes coords
-a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
-         verticalalignment='top', bbox=props)
-imgplot = plt.imshow(fgp_gpu, cmap="gray")
-plt.title('{}'.format('GPU results'))
-
-
-print ("--------Compare the results--------")
-tolerance = 1e-05
-diff_im = np.zeros(np.shape(fgp_cpu))
-diff_im = abs(fgp_cpu - fgp_gpu)
-diff_im[diff_im > tolerance] = 1
-a=fig.add_subplot(1,4,4)
-imgplot = plt.imshow(diff_im, vmin=0, vmax=1, cmap="gray")
-plt.title('{}'.format('Pixels larger threshold difference'))
-if (diff_im.sum() > 1):
-    print ("Arrays do not match!")
-else:
-    print ("Arrays match")
-
-#%%
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("____________SB-TV bench___________________")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-
-## plot 
-fig = plt.figure()
-plt.suptitle('Comparison of SB-TV regulariser using CPU and GPU implementations')
-a=fig.add_subplot(1,4,1)
-a.set_title('Noisy Image')
-imgplot = plt.imshow(u0,cmap="gray")
-
-# set parameters
-pars = {'algorithm' : SB_TV, \
-        'input' : u0,\
-        'regularisation_parameter':0.04, \
-        'number_of_iterations' :150 ,\
-        'tolerance_constant':1e-05,\
-        'methodTV': 0 ,\
-        'printingOut': 0 
-        }
-        
-print ("#############SB-TV CPU####################")
-start_time = timeit.default_timer()
-sb_cpu = SB_TV(pars['input'], 
-              pars['regularisation_parameter'],
-              pars['number_of_iterations'],
-              pars['tolerance_constant'], 
-              pars['methodTV'],
-              pars['printingOut'],'cpu')  
-             
-
-Qtools = QualityTools(Im, sb_cpu)
-pars['rmse'] = Qtools.rmse()
-
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(1,4,2)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
-# place a text box in upper left in axes coords
-a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
-         verticalalignment='top', bbox=props)
-imgplot = plt.imshow(sb_cpu, cmap="gray")
-plt.title('{}'.format('CPU results'))
-
-
-print ("##############SB TV GPU##################")
-start_time = timeit.default_timer()
-sb_gpu = SB_TV(pars['input'], 
-              pars['regularisation_parameter'],
-              pars['number_of_iterations'],
-              pars['tolerance_constant'], 
-              pars['methodTV'],
-              pars['printingOut'],'gpu')
-
-Qtools = QualityTools(Im, sb_gpu)
-pars['rmse'] = Qtools.rmse()
-pars['algorithm'] = SB_TV
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(1,4,3)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
-# place a text box in upper left in axes coords
-a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
-         verticalalignment='top', bbox=props)
-imgplot = plt.imshow(sb_gpu, cmap="gray")
-plt.title('{}'.format('GPU results'))
-
-print ("--------Compare the results--------")
-tolerance = 1e-05
-diff_im = np.zeros(np.shape(sb_cpu))
-diff_im = abs(sb_cpu - sb_gpu)
-diff_im[diff_im > tolerance] = 1
-a=fig.add_subplot(1,4,4)
-imgplot = plt.imshow(diff_im, vmin=0, vmax=1, cmap="gray")
-plt.title('{}'.format('Pixels larger threshold difference'))
-if (diff_im.sum() > 1):
-    print ("Arrays do not match!")
-else:
-    print ("Arrays match")
-#%%
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("____________TGV bench___________________")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-
-## plot 
-fig = plt.figure()
-plt.suptitle('Comparison of TGV regulariser using CPU and GPU implementations')
-a=fig.add_subplot(1,4,1)
-a.set_title('Noisy Image')
-imgplot = plt.imshow(u0,cmap="gray")
-
-# set parameters
-pars = {'algorithm' : TGV, \
-        'input' : u0,\
-        'regularisation_parameter':0.04, \
-        'alpha1':1.0,\
-        'alpha0':2.0,\
-        'number_of_iterations' :400 ,\
-        'LipshitzConstant' :12 ,\
-        }
-        
-print ("#############TGV CPU####################")
-start_time = timeit.default_timer()
-tgv_cpu = TGV(pars['input'], 
-              pars['regularisation_parameter'],
-              pars['alpha1'],
-              pars['alpha0'],
-              pars['number_of_iterations'],
-              pars['LipshitzConstant'],'cpu')
-             
-Qtools = QualityTools(Im, tgv_cpu)
-pars['rmse'] = Qtools.rmse()
-
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(1,4,2)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
-# place a text box in upper left in axes coords
-a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
-         verticalalignment='top', bbox=props)
-imgplot = plt.imshow(tgv_cpu, cmap="gray")
-plt.title('{}'.format('CPU results'))
-
-print ("##############TGV GPU##################")
-start_time = timeit.default_timer()
-tgv_gpu = TGV(pars['input'], 
-              pars['regularisation_parameter'],
-              pars['alpha1'],
-              pars['alpha0'],
-              pars['number_of_iterations'],
-              pars['LipshitzConstant'],'gpu')
-                                   
-Qtools = QualityTools(Im, tgv_gpu)
-pars['rmse'] = Qtools.rmse()
-pars['algorithm'] = TGV
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(1,4,3)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
-# place a text box in upper left in axes coords
-a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
-         verticalalignment='top', bbox=props)
-imgplot = plt.imshow(tgv_gpu, cmap="gray")
-plt.title('{}'.format('GPU results'))
-
-print ("--------Compare the results--------")
-tolerance = 1e-05
-diff_im = np.zeros(np.shape(tgv_gpu))
-diff_im = abs(tgv_cpu - tgv_gpu)
-diff_im[diff_im > tolerance] = 1
-a=fig.add_subplot(1,4,4)
-imgplot = plt.imshow(diff_im, vmin=0, vmax=1, cmap="gray")
-plt.title('{}'.format('Pixels larger threshold difference'))
-if (diff_im.sum() > 1):
-    print ("Arrays do not match!")
-else:
-    print ("Arrays match")
-#%%
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("____________LLT-ROF bench___________________")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-
-## plot 
-fig = plt.figure()
-plt.suptitle('Comparison of LLT-ROF regulariser using CPU and GPU implementations')
-a=fig.add_subplot(1,4,1)
-a.set_title('Noisy Image')
-imgplot = plt.imshow(u0,cmap="gray")
-
-# set parameters
-pars = {'algorithm' : LLT_ROF, \
-        'input' : u0,\
-        'regularisation_parameterROF':0.04, \
-        'regularisation_parameterLLT':0.01, \
-        'number_of_iterations' :4500 ,\
-        'time_marching_parameter' :0.00002 ,\
-        }
-        
-print ("#############LLT- ROF CPU####################")
-start_time = timeit.default_timer()
-lltrof_cpu = LLT_ROF(pars['input'], 
-              pars['regularisation_parameterROF'],
-              pars['regularisation_parameterLLT'],
-              pars['number_of_iterations'],
-              pars['time_marching_parameter'],'cpu')
-
-Qtools = QualityTools(Im, lltrof_cpu)
-pars['rmse'] = Qtools.rmse()
-
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(1,4,2)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
-# place a text box in upper left in axes coords
-a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
-         verticalalignment='top', bbox=props)
-imgplot = plt.imshow(lltrof_cpu, cmap="gray")
-plt.title('{}'.format('CPU results'))
-
-print ("#############LLT- ROF GPU####################")
-start_time = timeit.default_timer()
-lltrof_gpu = LLT_ROF(pars['input'], 
-              pars['regularisation_parameterROF'],
-              pars['regularisation_parameterLLT'],
-              pars['number_of_iterations'],
-              pars['time_marching_parameter'],'gpu')
-
-Qtools = QualityTools(Im, lltrof_gpu)
-pars['rmse'] = Qtools.rmse()
-
-pars['algorithm'] = LLT_ROF
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(1,4,3)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
-# place a text box in upper left in axes coords
-a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
-         verticalalignment='top', bbox=props)
-imgplot = plt.imshow(lltrof_gpu, cmap="gray")
-plt.title('{}'.format('GPU results'))
-
-print ("--------Compare the results--------")
-tolerance = 1e-05
-diff_im = np.zeros(np.shape(lltrof_gpu))
-diff_im = abs(lltrof_cpu - lltrof_gpu)
-diff_im[diff_im > tolerance] = 1
-a=fig.add_subplot(1,4,4)
-imgplot = plt.imshow(diff_im, vmin=0, vmax=1, cmap="gray")
-plt.title('{}'.format('Pixels larger threshold difference'))
-if (diff_im.sum() > 1):
-    print ("Arrays do not match!")
-else:
-    print ("Arrays match")
-#%%
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("_______________NDF bench___________________")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-
-## plot 
-fig = plt.figure()
-plt.suptitle('Comparison of NDF regulariser using CPU and GPU implementations')
-a=fig.add_subplot(1,4,1)
-a.set_title('Noisy Image')
-imgplot = plt.imshow(u0,cmap="gray")
-
-# set parameters
-pars = {'algorithm' : NDF, \
-        'input' : u0,\
-        'regularisation_parameter':0.06, \
-        'edge_parameter':0.04,\
-        'number_of_iterations' :1000 ,\
-        'time_marching_parameter':0.025,\
-        'penalty_type':  1
-        }
-        
-print ("#############NDF CPU####################")
-start_time = timeit.default_timer()
-ndf_cpu = NDF(pars['input'], 
-              pars['regularisation_parameter'],
-              pars['edge_parameter'], 
-              pars['number_of_iterations'],
-              pars['time_marching_parameter'], 
-              pars['penalty_type'],'cpu')
-
-Qtools = QualityTools(Im, ndf_cpu)
-pars['rmse'] = Qtools.rmse()
-
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(1,4,2)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
-# place a text box in upper left in axes coords
-a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
-         verticalalignment='top', bbox=props)
-imgplot = plt.imshow(ndf_cpu, cmap="gray")
-plt.title('{}'.format('CPU results'))
-
-
-print ("##############NDF GPU##################")
-start_time = timeit.default_timer()
-ndf_gpu = NDF(pars['input'], 
-              pars['regularisation_parameter'],
-              pars['edge_parameter'], 
-              pars['number_of_iterations'],
-              pars['time_marching_parameter'], 
-              pars['penalty_type'],'gpu')
-
-Qtools = QualityTools(Im, ndf_gpu)
-pars['rmse'] = Qtools.rmse()
-pars['algorithm'] = NDF
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(1,4,3)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
-# place a text box in upper left in axes coords
-a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
-         verticalalignment='top', bbox=props)
-imgplot = plt.imshow(ndf_gpu, cmap="gray")
-plt.title('{}'.format('GPU results'))
-
-print ("--------Compare the results--------")
-tolerance = 1e-05
-diff_im = np.zeros(np.shape(ndf_cpu))
-diff_im = abs(ndf_cpu - ndf_gpu)
-diff_im[diff_im > tolerance] = 1
-a=fig.add_subplot(1,4,4)
-imgplot = plt.imshow(diff_im, vmin=0, vmax=1, cmap="gray")
-plt.title('{}'.format('Pixels larger threshold difference'))
-if (diff_im.sum() > 1):
-    print ("Arrays do not match!")
-else:
-    print ("Arrays match")
-
-#%%
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("___Anisotropic Diffusion 4th Order (2D)____")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-
-## plot 
-fig = plt.figure()
-plt.suptitle('Comparison of Diff4th regulariser using CPU and GPU implementations')
-a=fig.add_subplot(1,4,1)
-a.set_title('Noisy Image')
-imgplot = plt.imshow(u0,cmap="gray")
-
-# set parameters
-pars = {'algorithm' : Diff4th, \
-        'input' : u0,\
-        'regularisation_parameter':3.5, \
-        'edge_parameter':0.02,\
-        'number_of_iterations' :500 ,\
-        'time_marching_parameter':0.001
-        }
-
-print ("#############Diff4th CPU####################")
-start_time = timeit.default_timer()
-diff4th_cpu = Diff4th(pars['input'], 
-              pars['regularisation_parameter'],
-              pars['edge_parameter'], 
-              pars['number_of_iterations'],
-              pars['time_marching_parameter'],'cpu')
-
-Qtools = QualityTools(Im, diff4th_cpu)
-pars['rmse'] = Qtools.rmse()
-
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(1,4,2)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
-# place a text box in upper left in axes coords
-a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
-         verticalalignment='top', bbox=props)
-imgplot = plt.imshow(diff4th_cpu, cmap="gray")
-plt.title('{}'.format('CPU results'))
-
-print ("##############Diff4th GPU##################")
-start_time = timeit.default_timer()
-diff4th_gpu = Diff4th(pars['input'], 
-              pars['regularisation_parameter'],
-              pars['edge_parameter'], 
-              pars['number_of_iterations'],
-              pars['time_marching_parameter'], 'gpu')
-
-Qtools = QualityTools(Im, diff4th_gpu)
-pars['rmse'] = Qtools.rmse()
-pars['algorithm'] = Diff4th
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(1,4,3)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
-# place a text box in upper left in axes coords
-a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
-         verticalalignment='top', bbox=props)
-imgplot = plt.imshow(diff4th_gpu, cmap="gray")
-plt.title('{}'.format('GPU results'))
-
-print ("--------Compare the results--------")
-tolerance = 1e-05
-diff_im = np.zeros(np.shape(diff4th_cpu))
-diff_im = abs(diff4th_cpu - diff4th_gpu)
-diff_im[diff_im > tolerance] = 1
-a=fig.add_subplot(1,4,4)
-imgplot = plt.imshow(diff_im, vmin=0, vmax=1, cmap="gray")
-plt.title('{}'.format('Pixels larger threshold difference'))
-if (diff_im.sum() > 1):
-    print ("Arrays do not match!")
-else:
-    print ("Arrays match")
-
-#%%
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("____________FGP-dTV bench___________________")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-
-## plot 
-fig = plt.figure()
-plt.suptitle('Comparison of FGP-dTV regulariser using CPU and GPU implementations')
-a=fig.add_subplot(1,4,1)
-a.set_title('Noisy Image')
-imgplot = plt.imshow(u0,cmap="gray")
-
-# set parameters
-pars = {'algorithm' : FGP_dTV, \
-        'input' : u0,\
-        'refdata' : u_ref,\
-        'regularisation_parameter':0.04, \
-        'number_of_iterations' :1000 ,\
-        'tolerance_constant':1e-07,\
-        'eta_const':0.2,\
-        'methodTV': 0 ,\
-        'nonneg': 0 ,\
-        'printingOut': 0 
-        }
-        
-print ("#############FGP dTV CPU####################")
-start_time = timeit.default_timer()
-fgp_dtv_cpu = FGP_dTV(pars['input'], 
-              pars['refdata'], 
-              pars['regularisation_parameter'],
-              pars['number_of_iterations'],
-              pars['tolerance_constant'], 
-              pars['eta_const'], 
-              pars['methodTV'],
-              pars['nonneg'],
-              pars['printingOut'],'cpu')
-
-Qtools = QualityTools(Im, fgp_dtv_cpu)
-pars['rmse'] = Qtools.rmse()
-
-
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(1,4,2)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
-# place a text box in upper left in axes coords
-a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
-         verticalalignment='top', bbox=props)
-imgplot = plt.imshow(fgp_dtv_cpu, cmap="gray")
-plt.title('{}'.format('CPU results'))
-
-print ("##############FGP dTV GPU##################")
-start_time = timeit.default_timer()
-fgp_dtv_gpu = FGP_dTV(pars['input'], 
-              pars['refdata'], 
-              pars['regularisation_parameter'],
-              pars['number_of_iterations'],
-              pars['tolerance_constant'], 
-              pars['eta_const'], 
-              pars['methodTV'],
-              pars['nonneg'],
-              pars['printingOut'],'gpu')
-Qtools = QualityTools(Im, fgp_dtv_gpu)
-pars['rmse'] = Qtools.rmse()
-pars['algorithm'] = FGP_dTV
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(1,4,3)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
-# place a text box in upper left in axes coords
-a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
-         verticalalignment='top', bbox=props)
-imgplot = plt.imshow(fgp_dtv_gpu, cmap="gray")
-plt.title('{}'.format('GPU results'))
-
-
-print ("--------Compare the results--------")
-tolerance = 1e-05
-diff_im = np.zeros(np.shape(fgp_dtv_cpu))
-diff_im = abs(fgp_dtv_cpu - fgp_dtv_gpu)
-diff_im[diff_im > tolerance] = 1
-a=fig.add_subplot(1,4,4)
-imgplot = plt.imshow(diff_im, vmin=0, vmax=1, cmap="gray")
-plt.title('{}'.format('Pixels larger threshold difference'))
-if (diff_im.sum() > 1):
-    print ("Arrays do not match!")
-else:
-    print ("Arrays match")
-#%%
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("____Non-local regularisation bench_________")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-
-## plot 
-fig = plt.figure()
-plt.suptitle('Comparison of Nonlocal TV regulariser using CPU and GPU implementations')
-a=fig.add_subplot(1,2,1)
-a.set_title('Noisy Image')
-imgplot = plt.imshow(u0,cmap="gray")
-
-pars = {'algorithm' : PatchSelect, \
-        'input' : u0,\
-        'searchwindow': 7, \
-        'patchwindow': 2,\
-        'neighbours' : 15 ,\
-        'edge_parameter':0.18}
-
-print ("############## Nonlocal Patches on CPU##################")
-start_time = timeit.default_timer()
-H_i, H_j, WeightsCPU = PatchSelect(pars['input'], 
-              pars['searchwindow'],
-              pars['patchwindow'], 
-              pars['neighbours'],
-              pars['edge_parameter'],'cpu')
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-
-print ("############## Nonlocal Patches on GPU##################")
-start_time = timeit.default_timer()
-start_time = timeit.default_timer()
-H_i, H_j, WeightsGPU = PatchSelect(pars['input'], 
-              pars['searchwindow'],
-              pars['patchwindow'], 
-              pars['neighbours'],
-              pars['edge_parameter'],'gpu')
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-
-print ("--------Compare the results--------")
-tolerance = 1e-05
-diff_im = np.zeros(np.shape(u0))
-diff_im = abs(WeightsCPU[0,:,:] - WeightsGPU[0,:,:])
-diff_im[diff_im > tolerance] = 1
-a=fig.add_subplot(1,2,2)
-imgplot = plt.imshow(diff_im, vmin=0, vmax=1, cmap="gray")
-plt.title('{}'.format('Pixels larger threshold difference'))
-if (diff_im.sum() > 1):
-    print ("Arrays do not match!")
-else:
-    print ("Arrays match")
-#%%
\ No newline at end of file
diff --git a/Wrappers/Python/demos/demo_gpu_regularisers.py b/Wrappers/Python/demos/demo_gpu_regularisers.py
deleted file mode 100644
index 54a1c14..0000000
--- a/Wrappers/Python/demos/demo_gpu_regularisers.py
+++ /dev/null
@@ -1,512 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-"""
-Created on Thu Feb 22 11:39:43 2018
-
-Demonstration of GPU regularisers
-
-@authors: Daniil Kazantsev, Edoardo Pasca
-"""
-
-import matplotlib.pyplot as plt
-import numpy as np
-import os
-import timeit
-from ccpi.filters.regularisers import ROF_TV, FGP_TV, SB_TV, TGV, LLT_ROF, FGP_dTV, NDF, Diff4th
-from ccpi.filters.regularisers import PatchSelect, NLTV
-from ccpi.supp.qualitymetrics import QualityTools
-###############################################################################
-def printParametersToString(pars):
-        txt = r''
-        for key, value in pars.items():
-            if key== 'algorithm' :
-                txt += "{0} = {1}".format(key, value.__name__)
-            elif key == 'input':
-                txt += "{0} = {1}".format(key, np.shape(value))
-            elif key == 'refdata':
-                txt += "{0} = {1}".format(key, np.shape(value))
-            else:
-                txt += "{0} = {1}".format(key, value)
-            txt += '\n'
-        return txt
-###############################################################################
-#%%
-filename = os.path.join(".." , ".." , ".." , "data" ,"lena_gray_512.tif")
-
-# read image
-Im = plt.imread(filename)                     
-Im = np.asarray(Im, dtype='float32')
-
-Im = Im/255
-perc = 0.05
-u0 = Im + np.random.normal(loc = 0 ,
-                                  scale = perc * Im , 
-                                  size = np.shape(Im))
-u_ref = Im + np.random.normal(loc = 0 ,
-                                  scale = 0.01 * Im , 
-                                  size = np.shape(Im))
-(N,M) = np.shape(u0)
-# map the u0 u0->u0>0
-# f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1)
-u0 = u0.astype('float32')
-u_ref = u_ref.astype('float32')
-"""
-M = M-100
-u_ref2 = np.zeros([N,M],dtype='float32')
-u_ref2[:,0:M] = u_ref[:,0:M]
-u_ref = u_ref2
-del u_ref2
-
-u02 = np.zeros([N,M],dtype='float32')
-u02[:,0:M] = u0[:,0:M]
-u0 = u02
-del u02
-
-Im2 = np.zeros([N,M],dtype='float32')
-Im2[:,0:M] = Im[:,0:M]
-Im = Im2
-del Im2
-"""
-#%%
-
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("____________ROF-TV regulariser_____________")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-
-## plot 
-fig = plt.figure()
-plt.suptitle('Performance of the ROF-TV regulariser using the GPU')
-a=fig.add_subplot(1,2,1)
-a.set_title('Noisy Image')
-imgplot = plt.imshow(u0,cmap="gray")
-
-# set parameters
-pars = {'algorithm': ROF_TV, \
-        'input' : u0,\
-        'regularisation_parameter':0.04,\
-        'number_of_iterations': 1200,\
-        'time_marching_parameter': 0.0025
-        }
-print ("##############ROF TV GPU##################")
-start_time = timeit.default_timer()
-rof_gpu = ROF_TV(pars['input'], 
-                     pars['regularisation_parameter'],
-                     pars['number_of_iterations'], 
-                     pars['time_marching_parameter'],'gpu')
-
-Qtools = QualityTools(Im, rof_gpu)
-pars['rmse'] = Qtools.rmse()
-pars['algorithm'] = ROF_TV
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(1,2,2)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
-# place a text box in upper left in axes coords
-a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
-         verticalalignment='top', bbox=props)
-imgplot = plt.imshow(rof_gpu, cmap="gray")
-plt.title('{}'.format('GPU results'))
-
-#%%
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("____________FGP-TV regulariser_____________")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-
-## plot 
-fig = plt.figure()
-plt.suptitle('Performance of the FGP-TV regulariser using the GPU')
-a=fig.add_subplot(1,2,1)
-a.set_title('Noisy Image')
-imgplot = plt.imshow(u0,cmap="gray")
-
-# set parameters
-pars = {'algorithm' : FGP_TV, \
-        'input' : u0,\
-        'regularisation_parameter':0.04, \
-        'number_of_iterations' :1200 ,\
-        'tolerance_constant':1e-06,\
-        'methodTV': 0 ,\
-        'nonneg': 0 ,\
-        'printingOut': 0 
-        }
-
-print ("##############FGP TV GPU##################")
-start_time = timeit.default_timer()
-fgp_gpu = FGP_TV(pars['input'], 
-              pars['regularisation_parameter'],
-              pars['number_of_iterations'],
-              pars['tolerance_constant'], 
-              pars['methodTV'],
-              pars['nonneg'],
-              pars['printingOut'],'gpu')
-Qtools = QualityTools(Im, fgp_gpu)
-pars['rmse'] = Qtools.rmse()
-pars['algorithm'] = FGP_TV
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(1,2,2)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
-# place a text box in upper left in axes coords
-a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
-         verticalalignment='top', bbox=props)
-imgplot = plt.imshow(fgp_gpu, cmap="gray")
-plt.title('{}'.format('GPU results'))
-
-#%%
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("____________SB-TV regulariser______________")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-
-## plot 
-fig = plt.figure()
-plt.suptitle('Performance of the SB-TV regulariser using the GPU')
-a=fig.add_subplot(1,2,1)
-a.set_title('Noisy Image')
-imgplot = plt.imshow(u0,cmap="gray")
-
-# set parameters
-pars = {'algorithm' : SB_TV, \
-        'input' : u0,\
-        'regularisation_parameter':0.04, \
-        'number_of_iterations' :150 ,\
-        'tolerance_constant':1e-06,\
-        'methodTV': 0 ,\
-        'printingOut': 0 
-        }
-
-print ("##############SB TV GPU##################")
-start_time = timeit.default_timer()
-sb_gpu = SB_TV(pars['input'], 
-              pars['regularisation_parameter'],
-              pars['number_of_iterations'],
-              pars['tolerance_constant'], 
-              pars['methodTV'],
-              pars['printingOut'],'gpu')
-
-Qtools = QualityTools(Im, sb_gpu)
-pars['rmse'] = Qtools.rmse()
-pars['algorithm'] = SB_TV
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(1,2,2)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
-# place a text box in upper left in axes coords
-a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
-         verticalalignment='top', bbox=props)
-imgplot = plt.imshow(sb_gpu, cmap="gray")
-plt.title('{}'.format('GPU results'))
-#%%
-
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("_____Total Generalised Variation (2D)______")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-
-## plot 
-fig = plt.figure()
-plt.suptitle('Performance of TGV regulariser using the GPU')
-a=fig.add_subplot(1,2,1)
-a.set_title('Noisy Image')
-imgplot = plt.imshow(u0,cmap="gray")
-
-# set parameters
-pars = {'algorithm' : TGV, \
-        'input' : u0,\
-        'regularisation_parameter':0.04, \
-        'alpha1':1.0,\
-        'alpha0':2.0,\
-        'number_of_iterations' :1250 ,\
-        'LipshitzConstant' :12 ,\
-        }
-        
-print ("#############TGV CPU####################")
-start_time = timeit.default_timer()
-tgv_gpu = TGV(pars['input'], 
-              pars['regularisation_parameter'],
-              pars['alpha1'],
-              pars['alpha0'],
-              pars['number_of_iterations'],
-              pars['LipshitzConstant'],'gpu')  
-
-Qtools = QualityTools(Im, tgv_gpu)
-pars['rmse'] = Qtools.rmse()
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(1,2,2)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
-# place a text box in upper left in axes coords
-a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
-         verticalalignment='top', bbox=props)
-imgplot = plt.imshow(tgv_gpu, cmap="gray")
-plt.title('{}'.format('GPU results'))
-
-#%%
-
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("______________LLT- ROF (2D)________________")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-
-## plot 
-fig = plt.figure()
-plt.suptitle('Performance of LLT-ROF regulariser using the GPU')
-a=fig.add_subplot(1,2,1)
-a.set_title('Noisy Image')
-imgplot = plt.imshow(u0,cmap="gray")
-
-# set parameters
-pars = {'algorithm' : LLT_ROF, \
-        'input' : u0,\
-        'regularisation_parameterROF':0.04, \
-        'regularisation_parameterLLT':0.01, \
-        'number_of_iterations' :500 ,\
-        'time_marching_parameter' :0.0025 ,\
-        }
-        
-print ("#############LLT- ROF GPU####################")
-start_time = timeit.default_timer()
-lltrof_gpu = LLT_ROF(pars['input'], 
-              pars['regularisation_parameterROF'],
-              pars['regularisation_parameterLLT'],
-              pars['number_of_iterations'],
-              pars['time_marching_parameter'],'gpu')
-             
-Qtools = QualityTools(Im, lltrof_gpu)
-pars['rmse'] = Qtools.rmse()
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(1,2,2)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
-# place a text box in upper left in axes coords
-a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
-         verticalalignment='top', bbox=props)
-imgplot = plt.imshow(lltrof_gpu, cmap="gray")
-plt.title('{}'.format('GPU results'))
-
-#%%
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("_______________NDF regulariser_____________")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-
-## plot 
-fig = plt.figure()
-plt.suptitle('Performance of the NDF regulariser using the GPU')
-a=fig.add_subplot(1,2,1)
-a.set_title('Noisy Image')
-imgplot = plt.imshow(u0,cmap="gray")
-
-# set parameters
-pars = {'algorithm' : NDF, \
-        'input' : u0,\
-        'regularisation_parameter':0.025, \
-        'edge_parameter':0.015,\
-        'number_of_iterations' :500 ,\
-        'time_marching_parameter':0.025,\
-        'penalty_type':  1
-        }
-
-print ("##############NDF GPU##################")
-start_time = timeit.default_timer()
-ndf_gpu = NDF(pars['input'], 
-              pars['regularisation_parameter'],
-              pars['edge_parameter'], 
-              pars['number_of_iterations'],
-              pars['time_marching_parameter'], 
-              pars['penalty_type'],'gpu')  
-
-Qtools = QualityTools(Im, ndf_gpu)
-pars['rmse'] = Qtools.rmse()
-pars['algorithm'] = NDF
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(1,2,2)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
-# place a text box in upper left in axes coords
-a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
-         verticalalignment='top', bbox=props)
-imgplot = plt.imshow(ndf_gpu, cmap="gray")
-plt.title('{}'.format('GPU results'))
-
-#%%
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("___Anisotropic Diffusion 4th Order (2D)____")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-
-## plot 
-fig = plt.figure()
-plt.suptitle('Performance of Diff4th regulariser using the GPU')
-a=fig.add_subplot(1,2,1)
-a.set_title('Noisy Image')
-imgplot = plt.imshow(u0,cmap="gray")
-
-# set parameters
-pars = {'algorithm' : Diff4th, \
-        'input' : u0,\
-        'regularisation_parameter':3.5, \
-        'edge_parameter':0.02,\
-        'number_of_iterations' :500 ,\
-        'time_marching_parameter':0.0015
-        }
-        
-print ("#############DIFF4th CPU################")
-start_time = timeit.default_timer()
-diff4_gpu = Diff4th(pars['input'], 
-              pars['regularisation_parameter'],
-              pars['edge_parameter'], 
-              pars['number_of_iterations'],
-              pars['time_marching_parameter'],'gpu')
-
-Qtools = QualityTools(Im, diff4_gpu)
-pars['algorithm'] = Diff4th
-pars['rmse'] = Qtools.rmse()
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(1,2,2)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
-# place a text box in upper left in axes coords
-a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
-         verticalalignment='top', bbox=props)
-imgplot = plt.imshow(diff4_gpu, cmap="gray")
-plt.title('{}'.format('GPU results'))
-
-#%%
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("___Nonlocal patches pre-calculation____")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-start_time = timeit.default_timer()
-# set parameters
-pars = {'algorithm' : PatchSelect, \
-        'input' : u0,\
-        'searchwindow': 7, \
-        'patchwindow': 2,\
-        'neighbours' : 15 ,\
-        'edge_parameter':0.18}
-
-H_i, H_j, Weights = PatchSelect(pars['input'], 
-              pars['searchwindow'],
-              pars['patchwindow'], 
-              pars['neighbours'],
-              pars['edge_parameter'],'gpu')
-              
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-"""
-plt.figure()
-plt.imshow(Weights[0,:,:],cmap="gray",interpolation="nearest",vmin=0, vmax=1)
-plt.show()
-"""
-#%%
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("___Nonlocal Total Variation penalty____")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-## plot 
-fig = plt.figure()
-plt.suptitle('Performance of NLTV regulariser using the CPU')
-a=fig.add_subplot(1,2,1)
-a.set_title('Noisy Image')
-imgplot = plt.imshow(u0,cmap="gray")
-
-pars2 = {'algorithm' : NLTV, \
-        'input' : u0,\
-        'H_i': H_i, \
-        'H_j': H_j,\
-        'H_k' : 0,\
-        'Weights' : Weights,\
-        'regularisation_parameter': 0.02,\
-        'iterations': 3
-        }
-start_time = timeit.default_timer()
-nltv_cpu = NLTV(pars2['input'], 
-              pars2['H_i'],
-              pars2['H_j'], 
-              pars2['H_k'],
-              pars2['Weights'],
-              pars2['regularisation_parameter'],
-              pars2['iterations'])
-
-Qtools = QualityTools(Im, nltv_cpu)
-pars['rmse'] = Qtools.rmse()
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(1,2,2)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
-# place a text box in upper left in axes coords
-a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
-         verticalalignment='top', bbox=props)
-imgplot = plt.imshow(nltv_cpu, cmap="gray")
-plt.title('{}'.format('CPU results'))
-#%%
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("____________FGP-dTV bench___________________")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-
-## plot 
-fig = plt.figure()
-plt.suptitle('Performance of the FGP-dTV regulariser using the GPU')
-a=fig.add_subplot(1,2,1)
-a.set_title('Noisy Image')
-imgplot = plt.imshow(u0,cmap="gray")
-
-# set parameters
-pars = {'algorithm' : FGP_dTV, \
-        'input' : u0,\
-        'refdata' : u_ref,\
-        'regularisation_parameter':0.04, \
-        'number_of_iterations' :2000 ,\
-        'tolerance_constant':1e-06,\
-        'eta_const':0.2,\
-        'methodTV': 0 ,\
-        'nonneg': 0 ,\
-        'printingOut': 0 
-        }
-
-print ("##############FGP dTV GPU##################")
-start_time = timeit.default_timer()
-fgp_dtv_gpu = FGP_dTV(pars['input'], 
-              pars['refdata'], 
-              pars['regularisation_parameter'],
-              pars['number_of_iterations'],
-              pars['tolerance_constant'], 
-              pars['eta_const'], 
-              pars['methodTV'],
-              pars['nonneg'],
-              pars['printingOut'],'gpu')
-
-Qtools = QualityTools(Im, fgp_dtv_gpu)
-pars['rmse'] = Qtools.rmse()
-pars['algorithm'] = FGP_dTV
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(1,2,2)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
-# place a text box in upper left in axes coords
-a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
-         verticalalignment='top', bbox=props)
-imgplot = plt.imshow(fgp_dtv_gpu, cmap="gray")
-plt.title('{}'.format('GPU results'))
diff --git a/Wrappers/Python/demos/demo_gpu_regularisers3D.py b/Wrappers/Python/demos/demo_gpu_regularisers3D.py
deleted file mode 100644
index d50c08e..0000000
--- a/Wrappers/Python/demos/demo_gpu_regularisers3D.py
+++ /dev/null
@@ -1,455 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-"""
-Created on Thu Feb 22 11:39:43 2018
-
-Demonstration of GPU regularisers
-
-@authors: Daniil Kazantsev, Edoardo Pasca
-"""
-
-import matplotlib.pyplot as plt
-import numpy as np
-import os
-import timeit
-from ccpi.filters.regularisers import ROF_TV, FGP_TV, SB_TV, TGV, LLT_ROF, FGP_dTV, NDF, Diff4th
-from ccpi.supp.qualitymetrics import QualityTools
-###############################################################################
-def printParametersToString(pars):
-        txt = r''
-        for key, value in pars.items():
-            if key== 'algorithm' :
-                txt += "{0} = {1}".format(key, value.__name__)
-            elif key == 'input':
-                txt += "{0} = {1}".format(key, np.shape(value))
-            elif key == 'refdata':
-                txt += "{0} = {1}".format(key, np.shape(value))
-            else:
-                txt += "{0} = {1}".format(key, value)
-            txt += '\n'
-        return txt
-###############################################################################
-#%%
-filename = os.path.join(".." , ".." , ".." , "data" ,"lena_gray_512.tif")
-
-# read image
-Im = plt.imread(filename)                     
-Im = np.asarray(Im, dtype='float32')
-
-Im = Im/255
-perc = 0.05
-u0 = Im + np.random.normal(loc = 0 ,
-                                  scale = perc * Im , 
-                                  size = np.shape(Im))
-u_ref = Im + np.random.normal(loc = 0 ,
-                                  scale = 0.01 * Im , 
-                                  size = np.shape(Im))
-(N,M) = np.shape(u0)
-# map the u0 u0->u0>0
-# f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1)
-u0 = u0.astype('float32')
-u_ref = u_ref.astype('float32')
-"""
-M = M-100
-u_ref2 = np.zeros([N,M],dtype='float32')
-u_ref2[:,0:M] = u_ref[:,0:M]
-u_ref = u_ref2
-del u_ref2
-
-u02 = np.zeros([N,M],dtype='float32')
-u02[:,0:M] = u0[:,0:M]
-u0 = u02
-del u02
-
-Im2 = np.zeros([N,M],dtype='float32')
-Im2[:,0:M] = Im[:,0:M]
-Im = Im2
-del Im2
-"""
-
-
-slices = 20
-
-filename = os.path.join(".." , ".." , ".." , "data" ,"lena_gray_512.tif")
-Im = plt.imread(filename)
-Im = np.asarray(Im, dtype='float32')
-
-Im = Im/255
-perc = 0.05
-
-noisyVol = np.zeros((slices,N,N),dtype='float32')
-noisyRef = np.zeros((slices,N,N),dtype='float32')
-idealVol = np.zeros((slices,N,N),dtype='float32')
-
-for i in range (slices):
-    noisyVol[i,:,:] = Im + np.random.normal(loc = 0 , scale = perc * Im , size = np.shape(Im))
-    noisyRef[i,:,:] = Im + np.random.normal(loc = 0 , scale = 0.01 * Im , size = np.shape(Im))
-    idealVol[i,:,:] = Im
-
-#%%
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("_______________ROF-TV (3D)_________________")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-
-## plot 
-fig = plt.figure()
-plt.suptitle('Performance of ROF-TV regulariser using the GPU')
-a=fig.add_subplot(1,2,1)
-a.set_title('Noisy 15th slice of a volume')
-imgplot = plt.imshow(noisyVol[10,:,:],cmap="gray")
-
-# set parameters
-pars = {'algorithm': ROF_TV, \
-        'input' : noisyVol,\
-        'regularisation_parameter':0.04,\
-        'number_of_iterations': 500,\
-        'time_marching_parameter': 0.0025        
-        }
-print ("#############ROF TV GPU####################")
-start_time = timeit.default_timer()
-rof_gpu3D = ROF_TV(pars['input'],
-             pars['regularisation_parameter'],
-             pars['number_of_iterations'],
-             pars['time_marching_parameter'],'gpu')
-
-Qtools = QualityTools(idealVol, rof_gpu3D)
-pars['rmse'] = Qtools.rmse()
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(1,2,2)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
-# place a text box in upper left in axes coords
-a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
-         verticalalignment='top', bbox=props)
-imgplot = plt.imshow(rof_gpu3D[10,:,:], cmap="gray")
-plt.title('{}'.format('Recovered volume on the GPU using ROF-TV'))
-#%%
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("_______________FGP-TV (3D)__________________")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-
-## plot 
-fig = plt.figure()
-plt.suptitle('Performance of FGP-TV regulariser using the GPU')
-a=fig.add_subplot(1,2,1)
-a.set_title('Noisy Image')
-imgplot = plt.imshow(noisyVol[10,:,:],cmap="gray")
-
-# set parameters
-pars = {'algorithm' : FGP_TV, \
-        'input' : noisyVol,\
-        'regularisation_parameter':0.04, \
-        'number_of_iterations' :300 ,\
-        'tolerance_constant':0.00001,\
-        'methodTV': 0 ,\
-        'nonneg': 0 ,\
-        'printingOut': 0 
-        }
-
-print ("#############FGP TV GPU####################")
-start_time = timeit.default_timer()
-fgp_gpu3D = FGP_TV(pars['input'], 
-              pars['regularisation_parameter'],
-              pars['number_of_iterations'],
-              pars['tolerance_constant'], 
-              pars['methodTV'],
-              pars['nonneg'],
-              pars['printingOut'],'gpu')
-
-Qtools = QualityTools(idealVol, fgp_gpu3D)
-pars['rmse'] = Qtools.rmse()
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(1,2,2)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
-# place a text box in upper left in axes coords
-a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
-         verticalalignment='top', bbox=props)
-imgplot = plt.imshow(fgp_gpu3D[10,:,:], cmap="gray")
-plt.title('{}'.format('Recovered volume on the GPU using FGP-TV'))
-
-#%%
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("_______________SB-TV (3D)__________________")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-
-## plot 
-fig = plt.figure()
-plt.suptitle('Performance of SB-TV regulariser using the GPU')
-a=fig.add_subplot(1,2,1)
-a.set_title('Noisy Image')
-imgplot = plt.imshow(noisyVol[10,:,:],cmap="gray")
-
-# set parameters
-pars = {'algorithm' : SB_TV, \
-        'input' : noisyVol,\
-        'regularisation_parameter':0.04, \
-        'number_of_iterations' :100 ,\
-        'tolerance_constant':1e-05,\
-        'methodTV': 0 ,\
-        'printingOut': 0 
-        }
-
-print ("#############SB TV GPU####################")
-start_time = timeit.default_timer()
-sb_gpu3D = SB_TV(pars['input'], 
-              pars['regularisation_parameter'],
-              pars['number_of_iterations'],
-              pars['tolerance_constant'], 
-              pars['methodTV'],
-              pars['printingOut'],'gpu')
-
-Qtools = QualityTools(idealVol, sb_gpu3D)
-pars['rmse'] = Qtools.rmse()
-
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(1,2,2)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
-# place a text box in upper left in axes coords
-a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
-         verticalalignment='top', bbox=props)
-imgplot = plt.imshow(sb_gpu3D[10,:,:], cmap="gray")
-plt.title('{}'.format('Recovered volume on the GPU using SB-TV'))
-#%%
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("_______________LLT-ROF (3D)_________________")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-
-## plot 
-fig = plt.figure()
-plt.suptitle('Performance of LLT-ROF regulariser using the GPU')
-a=fig.add_subplot(1,2,1)
-a.set_title('Noisy Image')
-imgplot = plt.imshow(noisyVol[10,:,:],cmap="gray")
-
-# set parameters
-pars = {'algorithm' : LLT_ROF, \
-        'input' : noisyVol,\
-        'regularisation_parameterROF':0.04, \
-        'regularisation_parameterLLT':0.015, \
-        'number_of_iterations' :300 ,\
-        'time_marching_parameter' :0.0025 ,\
-        }
-
-print ("#############LLT ROF CPU####################")
-start_time = timeit.default_timer()
-lltrof_gpu3D = LLT_ROF(pars['input'], 
-              pars['regularisation_parameterROF'],
-              pars['regularisation_parameterLLT'],
-              pars['number_of_iterations'],
-              pars['time_marching_parameter'],'gpu')
-
-Qtools = QualityTools(idealVol, lltrof_gpu3D)
-pars['rmse'] = Qtools.rmse()
-
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(1,2,2)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
-# place a text box in upper left in axes coords
-a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
-         verticalalignment='top', bbox=props)
-imgplot = plt.imshow(lltrof_gpu3D[10,:,:], cmap="gray")
-plt.title('{}'.format('Recovered volume on the GPU using LLT-ROF'))
-
-#%%
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("_______________TGV (3D)_________________")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-
-## plot 
-fig = plt.figure()
-plt.suptitle('Performance of TGV regulariser using the GPU')
-a=fig.add_subplot(1,2,1)
-a.set_title('Noisy Image')
-imgplot = plt.imshow(noisyVol[10,:,:],cmap="gray")
-
-# set parameters
-pars = {'algorithm' : TGV, \
-        'input' : noisyVol,\
-        'regularisation_parameter':0.04, \
-        'alpha1':1.0,\
-        'alpha0':2.0,\
-        'number_of_iterations' :600 ,\
-        'LipshitzConstant' :12 ,\
-        }
-
-print ("#############TGV GPU####################")
-start_time = timeit.default_timer()
-tgv_gpu3D = TGV(pars['input'], 
-              pars['regularisation_parameter'],
-              pars['alpha1'],
-              pars['alpha0'],
-              pars['number_of_iterations'],
-              pars['LipshitzConstant'],'gpu')
-
-Qtools = QualityTools(idealVol, tgv_gpu3D)
-pars['rmse'] = Qtools.rmse()
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(1,2,2)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
-# place a text box in upper left in axes coords
-a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
-         verticalalignment='top', bbox=props)
-imgplot = plt.imshow(tgv_gpu3D[10,:,:], cmap="gray")
-plt.title('{}'.format('Recovered volume on the GPU using TGV'))
-#%%
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("_______________NDF-TV (3D)_________________")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-
-## plot 
-fig = plt.figure()
-plt.suptitle('Performance of NDF regulariser using the GPU')
-a=fig.add_subplot(1,2,1)
-a.set_title('Noisy Image')
-imgplot = plt.imshow(noisyVol[10,:,:],cmap="gray")
-
-# set parameters
-pars = {'algorithm' : NDF, \
-        'input' : noisyVol,\
-        'regularisation_parameter':0.025, \
-        'edge_parameter':0.015,\
-        'number_of_iterations' :500 ,\
-        'time_marching_parameter':0.025,\
-        'penalty_type':  1
-        }
-
-print ("#############NDF GPU####################")
-start_time = timeit.default_timer()
-ndf_gpu3D = NDF(pars['input'], 
-              pars['regularisation_parameter'],
-              pars['edge_parameter'], 
-              pars['number_of_iterations'],
-              pars['time_marching_parameter'], 
-              pars['penalty_type'],'gpu')
-
-Qtools = QualityTools(idealVol, ndf_gpu3D)
-pars['rmse'] = Qtools.rmse()
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(1,2,2)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
-# place a text box in upper left in axes coords
-a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
-         verticalalignment='top', bbox=props)
-imgplot = plt.imshow(ndf_gpu3D[10,:,:], cmap="gray")
-plt.title('{}'.format('Recovered volume on the GPU using NDF'))
-
-#%%
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("___Anisotropic Diffusion 4th Order (3D)____")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-
-## plot 
-fig = plt.figure()
-plt.suptitle('Performance of DIFF4th regulariser using the GPU')
-a=fig.add_subplot(1,2,1)
-a.set_title('Noisy Image')
-imgplot = plt.imshow(noisyVol[10,:,:],cmap="gray")
-
-# set parameters
-pars = {'algorithm' : Diff4th, \
-        'input' : noisyVol,\
-        'regularisation_parameter':3.5, \
-        'edge_parameter':0.02,\
-        'number_of_iterations' :300 ,\
-        'time_marching_parameter':0.0015
-        }
-        
-print ("#############DIFF4th CPU################")
-start_time = timeit.default_timer()
-diff4_gpu3D = Diff4th(pars['input'], 
-              pars['regularisation_parameter'],
-              pars['edge_parameter'], 
-              pars['number_of_iterations'],
-              pars['time_marching_parameter'],'gpu')
-
-Qtools = QualityTools(idealVol, diff4_gpu3D)
-pars['rmse'] = Qtools.rmse()
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(1,2,2)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
-# place a text box in upper left in axes coords
-a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
-         verticalalignment='top', bbox=props)
-imgplot = plt.imshow(diff4_gpu3D[10,:,:], cmap="gray")
-plt.title('{}'.format('GPU results'))
-
-#%%
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-print ("_______________FGP-dTV (3D)________________")
-print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
-
-## plot 
-fig = plt.figure()
-plt.suptitle('Performance of FGP-dTV regulariser using the GPU')
-a=fig.add_subplot(1,2,1)
-a.set_title('Noisy Image')
-imgplot = plt.imshow(noisyVol[10,:,:],cmap="gray")
-
-# set parameters
-pars = {'algorithm' : FGP_dTV, \
-        'input' : noisyVol,\
-        'refdata' : noisyRef,\
-        'regularisation_parameter':0.04, \
-        'number_of_iterations' :300 ,\
-        'tolerance_constant':0.00001,\
-        'eta_const':0.2,\
-        'methodTV': 0 ,\
-        'nonneg': 0 ,\
-        'printingOut': 0 
-        }
-
-print ("#############FGP TV GPU####################")
-start_time = timeit.default_timer()
-fgp_dTV_gpu3D = FGP_dTV(pars['input'],
-              pars['refdata'], 
-              pars['regularisation_parameter'],
-              pars['number_of_iterations'],
-              pars['tolerance_constant'], 
-              pars['eta_const'],
-              pars['methodTV'],
-              pars['nonneg'],
-              pars['printingOut'],'gpu')
-
-Qtools = QualityTools(idealVol, fgp_dTV_gpu3D)
-pars['rmse'] = Qtools.rmse()
-
-txtstr = printParametersToString(pars)
-txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
-print (txtstr)
-a=fig.add_subplot(1,2,2)
-
-# these are matplotlib.patch.Patch properties
-props = dict(boxstyle='round', facecolor='wheat', alpha=0.75)
-# place a text box in upper left in axes coords
-a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
-         verticalalignment='top', bbox=props)
-imgplot = plt.imshow(fgp_dTV_gpu3D[10,:,:], cmap="gray")
-plt.title('{}'.format('Recovered volume on the GPU using FGP-dTV'))
-#%%
diff --git a/Wrappers/Python/setup-regularisers.py.in b/Wrappers/Python/setup-regularisers.py.in
deleted file mode 100644
index 16fb380..0000000
--- a/Wrappers/Python/setup-regularisers.py.in
+++ /dev/null
@@ -1,75 +0,0 @@
-#!/usr/bin/env python
-
-import setuptools
-from distutils.core import setup
-from distutils.extension import Extension
-from Cython.Distutils import build_ext
-
-import os
-import sys
-import numpy
-import platform	
-
-cil_version=os.environ['CIL_VERSION']
-if  cil_version == '':
-    print("Please set the environmental variable CIL_VERSION")
-    sys.exit(1)
-	
-library_include_path = ""
-library_lib_path = ""
-try:
-    library_include_path = os.environ['LIBRARY_INC']
-    library_lib_path = os.environ['LIBRARY_LIB']
-except:
-    library_include_path = os.environ['PREFIX']+'/include'
-    pass
-    
-extra_include_dirs = [numpy.get_include(), library_include_path]
-#extra_library_dirs = [os.path.join(library_include_path, "..", "lib")]
-extra_compile_args = []
-extra_library_dirs = [library_lib_path]
-extra_compile_args = []
-extra_link_args = []
-extra_libraries = ['cilreg']
-
-print ("extra_library_dirs " , extra_library_dirs)
-
-extra_include_dirs += [os.path.join(".." , ".." , "Core"),
-                       os.path.join(".." , ".." , "Core",  "regularisers_CPU"),
-                       os.path.join(".." , ".." , "Core",  "inpainters_CPU"),
-                       os.path.join(".." , ".." , "Core",  "regularisers_GPU" , "TV_FGP" ) , 
-                       os.path.join(".." , ".." , "Core",  "regularisers_GPU" , "TV_ROF" ) , 
-                       os.path.join(".." , ".." , "Core",  "regularisers_GPU" , "TV_SB" ) ,
-                       os.path.join(".." , ".." , "Core",  "regularisers_GPU" , "TGV" ) ,
-                       os.path.join(".." , ".." , "Core",  "regularisers_GPU" , "LLTROF" ) ,
-                       os.path.join(".." , ".." , "Core",  "regularisers_GPU" , "NDF" ) ,
-                       os.path.join(".." , ".." , "Core",  "regularisers_GPU" , "dTV_FGP" ) , 
-                       os.path.join(".." , ".." , "Core",  "regularisers_GPU" , "DIFF4th" ) , 
-                       os.path.join(".." , ".." , "Core",  "regularisers_GPU" , "PatchSelect" ) ,
-						   "."]
-
-if platform.system() == 'Windows':				   
-    extra_compile_args[0:] = ['/DWIN32','/EHsc','/DBOOST_ALL_NO_LIB' , '/openmp' ]   
-else:
-    extra_compile_args = ['-fopenmp','-O2', '-funsigned-char', '-Wall', '-std=c++0x']
-    extra_libraries += [@EXTRA_OMP_LIB@]
-    
-setup(
-    name='ccpi',
-	description='CCPi Core Imaging Library - Image regularisers',
-	version=cil_version,
-    cmdclass = {'build_ext': build_ext},
-    ext_modules = [Extension("ccpi.filters.cpu_regularisers",
-                             sources=[os.path.join("." , "src", "cpu_regularisers.pyx" ) ],
-                             include_dirs=extra_include_dirs, 
-							 library_dirs=extra_library_dirs, 
-							 extra_compile_args=extra_compile_args, 
-							 libraries=extra_libraries ), 
-    
-    ],
-	zip_safe = False,	
-	packages = {'ccpi','ccpi.filters', 'ccpi.supp'},
-)
-
-
-@SETUP_GPU_WRAPPERS@
diff --git a/Wrappers/Python/src/cpu_regularisers.pyx b/Wrappers/Python/src/cpu_regularisers.pyx
deleted file mode 100644
index 11a0617..0000000
--- a/Wrappers/Python/src/cpu_regularisers.pyx
+++ /dev/null
@@ -1,685 +0,0 @@
-# distutils: language=c++
-"""
-Copyright 2018 CCPi
-Licensed under the Apache License, Version 2.0 (the "License");
-you may not use this file except in compliance with the License.
-You may obtain a copy of the License at
-    http://www.apache.org/licenses/LICENSE-2.0
-Unless required by applicable law or agreed to in writing, software
-distributed under the License is distributed on an "AS IS" BASIS,
-WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-See the License for the specific language governing permissions and
-limitations under the License.
-
-Author: Edoardo Pasca, Daniil Kazantsev
-"""
-
-import cython
-import numpy as np
-cimport numpy as np
-
-cdef extern float TV_ROF_CPU_main(float *Input, float *Output, float lambdaPar, int iterationsNumb, float tau, int dimX, int dimY, int dimZ);
-cdef extern float TV_FGP_CPU_main(float *Input, float *Output, float lambdaPar, int iterationsNumb, float epsil, int methodTV, int nonneg, int printM, int dimX, int dimY, int dimZ);
-cdef extern float SB_TV_CPU_main(float *Input, float *Output, float lambdaPar, int iterationsNumb, float epsil, int methodTV, int printM, int dimX, int dimY, int dimZ);
-cdef extern float LLT_ROF_CPU_main(float *Input, float *Output, float lambdaROF, float lambdaLLT, int iterationsNumb, float tau, int dimX, int dimY, int dimZ);
-cdef extern float TGV_main(float *Input, float *Output, float lambdaPar, float alpha1, float alpha0, int iterationsNumb, float L2, int dimX, int dimY, int dimZ);
-cdef extern float Diffusion_CPU_main(float *Input, float *Output, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int penaltytype, int dimX, int dimY, int dimZ);
-cdef extern float Diffus4th_CPU_main(float *Input, float *Output, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int dimX, int dimY, int dimZ);
-cdef extern float TNV_CPU_main(float *Input, float *u, float lambdaPar, int maxIter, float tol, int dimX, int dimY, int dimZ);
-cdef extern float dTV_FGP_CPU_main(float *Input, float *InputRef, float *Output, float lambdaPar, int iterationsNumb, float epsil, float eta, int methodTV, int nonneg, int printM, int dimX, int dimY, int dimZ);
-cdef extern float PatchSelect_CPU_main(float *Input, unsigned short *H_i, unsigned short *H_j, unsigned short *H_k, float *Weights, int dimX, int dimY, int dimZ, int SearchWindow, int SimilarWin, int NumNeighb, float h, int switchM);
-cdef extern float Nonlocal_TV_CPU_main(float *A_orig, float *Output, unsigned short *H_i, unsigned short *H_j, unsigned short *H_k, float *Weights, int dimX, int dimY, int dimZ, int NumNeighb, float lambdaReg, int IterNumb);
-
-cdef extern float Diffusion_Inpaint_CPU_main(float *Input, unsigned char *Mask, float *Output, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int penaltytype, int dimX, int dimY, int dimZ);
-cdef extern float NonlocalMarching_Inpaint_main(float *Input, unsigned char *M, float *Output, unsigned char *M_upd, int SW_increment, int iterationsNumb, int trigger, int dimX, int dimY, int dimZ);
-cdef extern float TV_energy2D(float *U, float *U0, float *E_val, float lambdaPar, int type, int dimX, int dimY);
-cdef extern float TV_energy3D(float *U, float *U0, float *E_val, float lambdaPar, int type, int dimX, int dimY, int dimZ);
-#****************************************************************#
-#********************** Total-variation ROF *********************#
-#****************************************************************#
-def TV_ROF_CPU(inputData, regularisation_parameter, iterationsNumb, marching_step_parameter):
-    if inputData.ndim == 2:
-        return TV_ROF_2D(inputData, regularisation_parameter, iterationsNumb, marching_step_parameter)
-    elif inputData.ndim == 3:
-        return TV_ROF_3D(inputData, regularisation_parameter, iterationsNumb, marching_step_parameter)
-
-def TV_ROF_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
-                     float regularisation_parameter,
-                     int iterationsNumb,                     
-                     float marching_step_parameter):
-    cdef long dims[2]
-    dims[0] = inputData.shape[0]
-    dims[1] = inputData.shape[1]
-    
-    cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
-            np.zeros([dims[0],dims[1]], dtype='float32')
-                   
-    # Run ROF iterations for 2D data 
-    TV_ROF_CPU_main(&inputData[0,0], &outputData[0,0], regularisation_parameter, iterationsNumb, marching_step_parameter, dims[1], dims[0], 1)
-    
-    return outputData
-            
-def TV_ROF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
-                     float regularisation_parameter,
-                     int iterationsNumb,
-                     float marching_step_parameter):
-    cdef long dims[3]
-    dims[0] = inputData.shape[0]
-    dims[1] = inputData.shape[1]
-    dims[2] = inputData.shape[2]
-    
-    cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
-            np.zeros([dims[0],dims[1],dims[2]], dtype='float32')
-           
-    # Run ROF iterations for 3D data 
-    TV_ROF_CPU_main(&inputData[0,0,0], &outputData[0,0,0], regularisation_parameter, iterationsNumb, marching_step_parameter, dims[2], dims[1], dims[0])
-
-    return outputData
-
-#****************************************************************#
-#********************** Total-variation FGP *********************#
-#****************************************************************#
-#******** Total-variation Fast-Gradient-Projection (FGP)*********#
-def TV_FGP_CPU(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, printM):
-    if inputData.ndim == 2:
-        return TV_FGP_2D(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, printM)
-    elif inputData.ndim == 3:
-        return TV_FGP_3D(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, printM)
-
-def TV_FGP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
-                     float regularisation_parameter,
-                     int iterationsNumb, 
-                     float tolerance_param,
-                     int methodTV,
-                     int nonneg,
-                     int printM):
-                         
-    cdef long dims[2]
-    dims[0] = inputData.shape[0]
-    dims[1] = inputData.shape[1]
-    
-    cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
-            np.zeros([dims[0],dims[1]], dtype='float32')
-                   
-    #/* Run FGP-TV iterations for 2D data */
-    TV_FGP_CPU_main(&inputData[0,0], &outputData[0,0], regularisation_parameter, 
-                       iterationsNumb, 
-                       tolerance_param,
-                       methodTV,
-                       nonneg,
-                       printM,
-                       dims[1],dims[0],1)
-    
-    return outputData        
-            
-def TV_FGP_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
-                     float regularisation_parameter,
-                     int iterationsNumb, 
-                     float tolerance_param,
-                     int methodTV,
-                     int nonneg,
-                     int printM):
-    cdef long dims[3]
-    dims[0] = inputData.shape[0]
-    dims[1] = inputData.shape[1]
-    dims[2] = inputData.shape[2]
-    
-    cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
-            np.zeros([dims[0], dims[1], dims[2]], dtype='float32')
-           
-    #/* Run FGP-TV iterations for 3D data */
-    TV_FGP_CPU_main(&inputData[0,0,0], &outputData[0,0,0], regularisation_parameter,
-                       iterationsNumb, 
-                       tolerance_param,
-                       methodTV,
-                       nonneg,
-                       printM,
-                       dims[2], dims[1], dims[0])
-    return outputData 
-
-#***************************************************************#
-#********************** Total-variation SB *********************#
-#***************************************************************#
-#*************** Total-variation Split Bregman (SB)*************#
-def TV_SB_CPU(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, printM):
-    if inputData.ndim == 2:
-        return TV_SB_2D(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, printM)
-    elif inputData.ndim == 3:
-        return TV_SB_3D(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, printM)
-
-def TV_SB_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
-                     float regularisation_parameter,
-                     int iterationsNumb, 
-                     float tolerance_param,
-                     int methodTV,
-                     int printM):
-                         
-    cdef long dims[2]
-    dims[0] = inputData.shape[0]
-    dims[1] = inputData.shape[1]
-    
-    cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
-            np.zeros([dims[0],dims[1]], dtype='float32')
-                   
-    #/* Run SB-TV iterations for 2D data */
-    SB_TV_CPU_main(&inputData[0,0], &outputData[0,0], regularisation_parameter, 
-                       iterationsNumb, 
-                       tolerance_param,
-                       methodTV,
-                       printM,
-                       dims[1],dims[0],1)
-    
-    return outputData        
-            
-def TV_SB_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
-                     float regularisation_parameter,
-                     int iterationsNumb, 
-                     float tolerance_param,
-                     int methodTV,
-                     int printM):
-    cdef long dims[3]
-    dims[0] = inputData.shape[0]
-    dims[1] = inputData.shape[1]
-    dims[2] = inputData.shape[2]
-    
-    cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
-            np.zeros([dims[0], dims[1], dims[2]], dtype='float32')
-           
-    #/* Run SB-TV iterations for 3D data */
-    SB_TV_CPU_main(&inputData[0,0,0], &outputData[0,0,0], regularisation_parameter,
-                       iterationsNumb, 
-                       tolerance_param,
-                       methodTV,
-                       printM,
-                       dims[2], dims[1], dims[0])
-    return outputData 
-
-#***************************************************************#
-#***************** Total Generalised Variation *****************#
-#***************************************************************#
-def TGV_CPU(inputData, regularisation_parameter, alpha1, alpha0, iterations, LipshitzConst):
-    if inputData.ndim == 2:
-        return TGV_2D(inputData, regularisation_parameter, alpha1, alpha0, 
-                      iterations, LipshitzConst)
-    elif inputData.ndim == 3:
-        return TGV_3D(inputData, regularisation_parameter, alpha1, alpha0, 
-                      iterations, LipshitzConst)
-
-def TGV_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
-                     float regularisation_parameter,
-                     float alpha1,
-                     float alpha0,
-                     int iterationsNumb, 
-                     float LipshitzConst):
-                         
-    cdef long dims[2]
-    dims[0] = inputData.shape[0]
-    dims[1] = inputData.shape[1]
-    
-    cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
-            np.zeros([dims[0],dims[1]], dtype='float32')
-                   
-    #/* Run TGV iterations for 2D data */
-    TGV_main(&inputData[0,0], &outputData[0,0], regularisation_parameter, 
-                       alpha1,
-                       alpha0,
-                       iterationsNumb, 
-                       LipshitzConst,
-                       dims[1],dims[0],1)
-    return outputData
-def TGV_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
-                     float regularisation_parameter,
-                     float alpha1,
-                     float alpha0,
-                     int iterationsNumb, 
-                     float LipshitzConst):
-                         
-    cdef long dims[3]
-    dims[0] = inputData.shape[0]
-    dims[1] = inputData.shape[1]
-    dims[2] = inputData.shape[2]
-    
-    cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
-            np.zeros([dims[0], dims[1], dims[2]], dtype='float32')
-                   
-    #/* Run TGV iterations for 3D data */
-    TGV_main(&inputData[0,0,0], &outputData[0,0,0], regularisation_parameter, 
-                       alpha1,
-                       alpha0,
-                       iterationsNumb, 
-                       LipshitzConst,
-                       dims[2], dims[1], dims[0])
-    return outputData
-
-#***************************************************************#
-#******************* ROF - LLT regularisation ******************#
-#***************************************************************#
-def LLT_ROF_CPU(inputData, regularisation_parameterROF, regularisation_parameterLLT, iterations, time_marching_parameter):
-    if inputData.ndim == 2:
-        return LLT_ROF_2D(inputData, regularisation_parameterROF, regularisation_parameterLLT, iterations, time_marching_parameter)
-    elif inputData.ndim == 3:
-        return LLT_ROF_3D(inputData, regularisation_parameterROF, regularisation_parameterLLT, iterations, time_marching_parameter)
-
-def LLT_ROF_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
-                     float regularisation_parameterROF,
-                     float regularisation_parameterLLT,
-                     int iterations, 
-                     float time_marching_parameter):
-                         
-    cdef long dims[2]
-    dims[0] = inputData.shape[0]
-    dims[1] = inputData.shape[1]
-    
-    cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
-            np.zeros([dims[0],dims[1]], dtype='float32')
-                   
-    #/* Run ROF-LLT iterations for 2D data */
-    LLT_ROF_CPU_main(&inputData[0,0], &outputData[0,0], regularisation_parameterROF, regularisation_parameterLLT, iterations, time_marching_parameter, dims[1],dims[0],1)
-    return outputData
-
-def LLT_ROF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
-                     float regularisation_parameterROF,
-                     float regularisation_parameterLLT,
-                     int iterations, 
-                     float time_marching_parameter):
-						 
-    cdef long dims[3]
-    dims[0] = inputData.shape[0]
-    dims[1] = inputData.shape[1]
-    dims[2] = inputData.shape[2]
-    
-    cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
-            np.zeros([dims[0], dims[1], dims[2]], dtype='float32')
-           
-    #/* Run ROF-LLT iterations for 3D data */
-    LLT_ROF_CPU_main(&inputData[0,0,0], &outputData[0,0,0], regularisation_parameterROF, regularisation_parameterLLT, iterations, time_marching_parameter, dims[2], dims[1], dims[0])
-    return outputData 
-
-#****************************************************************#
-#**************Directional Total-variation FGP ******************#
-#****************************************************************#
-#******** Directional TV Fast-Gradient-Projection (FGP)*********#
-def dTV_FGP_CPU(inputData, refdata, regularisation_parameter, iterationsNumb, tolerance_param, eta_const, methodTV, nonneg, printM):
-    if inputData.ndim == 2:
-        return dTV_FGP_2D(inputData, refdata, regularisation_parameter, iterationsNumb, tolerance_param, eta_const, methodTV, nonneg, printM)
-    elif inputData.ndim == 3:
-        return dTV_FGP_3D(inputData, refdata, regularisation_parameter, iterationsNumb, tolerance_param, eta_const, methodTV, nonneg, printM)
-
-def dTV_FGP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
-               np.ndarray[np.float32_t, ndim=2, mode="c"] refdata,
-                     float regularisation_parameter,
-                     int iterationsNumb, 
-                     float tolerance_param,
-                     float eta_const,
-                     int methodTV,
-                     int nonneg,
-                     int printM):
-                         
-    cdef long dims[2]
-    dims[0] = inputData.shape[0]
-    dims[1] = inputData.shape[1]
-    
-    cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
-            np.zeros([dims[0],dims[1]], dtype='float32')
-                   
-    #/* Run FGP-dTV iterations for 2D data */
-    dTV_FGP_CPU_main(&inputData[0,0], &refdata[0,0], &outputData[0,0], regularisation_parameter, 
-                       iterationsNumb, 
-                       tolerance_param,
-                       eta_const,
-                       methodTV,                       
-                       nonneg,
-                       printM,
-                       dims[1], dims[0], 1)
-    
-    return outputData        
-            
-def dTV_FGP_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
-               np.ndarray[np.float32_t, ndim=3, mode="c"] refdata,
-                     float regularisation_parameter,
-                     int iterationsNumb, 
-                     float tolerance_param,
-                     float eta_const,
-                     int methodTV,
-                     int nonneg,
-                     int printM):
-    cdef long dims[3]
-    dims[0] = inputData.shape[0]
-    dims[1] = inputData.shape[1]
-    dims[2] = inputData.shape[2]
-    
-    cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
-            np.zeros([dims[0], dims[1], dims[2]], dtype='float32')
-           
-    #/* Run FGP-dTV iterations for 3D data */
-    dTV_FGP_CPU_main(&inputData[0,0,0], &refdata[0,0,0], &outputData[0,0,0], regularisation_parameter,
-                       iterationsNumb, 
-                       tolerance_param,
-                       eta_const,
-                       methodTV,
-                       nonneg,
-                       printM,
-                       dims[2], dims[1], dims[0])
-    return outputData
-    
-#****************************************************************#
-#*********************Total Nuclear Variation********************#
-#****************************************************************#
-def TNV_CPU(inputData, regularisation_parameter, iterationsNumb, tolerance_param):
-    if inputData.ndim == 2:
-        return 
-    elif inputData.ndim == 3:
-        return TNV_3D(inputData, regularisation_parameter, iterationsNumb, tolerance_param)
-
-def TNV_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
-                     float regularisation_parameter,
-                     int iterationsNumb,
-                     float tolerance_param):
-    cdef long dims[3]
-    dims[0] = inputData.shape[0]
-    dims[1] = inputData.shape[1]
-    dims[2] = inputData.shape[2]
-    
-    cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
-            np.zeros([dims[0],dims[1],dims[2]], dtype='float32')
-           
-    # Run TNV iterations for 3D (X,Y,Channels) data 
-    TNV_CPU_main(&inputData[0,0,0], &outputData[0,0,0], regularisation_parameter, iterationsNumb, tolerance_param, dims[2], dims[1], dims[0])
-    return outputData
-#****************************************************************#
-#***************Nonlinear (Isotropic) Diffusion******************#
-#****************************************************************#
-def NDF_CPU(inputData, regularisation_parameter, edge_parameter, iterationsNumb,time_marching_parameter, penalty_type):
-    if inputData.ndim == 2:
-        return NDF_2D(inputData, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type)
-    elif inputData.ndim == 3:
-        return NDF_3D(inputData, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type)
-
-def NDF_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
-                     float regularisation_parameter,
-                     float edge_parameter,
-                     int iterationsNumb,                     
-                     float time_marching_parameter,
-                     int penalty_type):
-    cdef long dims[2]
-    dims[0] = inputData.shape[0]
-    dims[1] = inputData.shape[1]
-    
-    cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
-            np.zeros([dims[0],dims[1]], dtype='float32')   
-    
-    # Run Nonlinear Diffusion iterations for 2D data 
-    Diffusion_CPU_main(&inputData[0,0], &outputData[0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type, dims[1], dims[0], 1)
-    return outputData
-            
-def NDF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
-                     float regularisation_parameter,
-                     float edge_parameter,
-                     int iterationsNumb,                     
-                     float time_marching_parameter,
-                     int penalty_type):
-    cdef long dims[3]
-    dims[0] = inputData.shape[0]
-    dims[1] = inputData.shape[1]
-    dims[2] = inputData.shape[2]
-    
-    cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
-            np.zeros([dims[0],dims[1],dims[2]], dtype='float32')
-    
-    # Run Nonlinear Diffusion iterations for  3D data 
-    Diffusion_CPU_main(&inputData[0,0,0], &outputData[0,0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type, dims[2], dims[1], dims[0])
-
-    return outputData
-
-#****************************************************************#
-#*************Anisotropic Fourth-Order diffusion*****************#
-#****************************************************************#
-def Diff4th_CPU(inputData, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter):
-    if inputData.ndim == 2:
-        return Diff4th_2D(inputData, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter)
-    elif inputData.ndim == 3:
-        return Diff4th_3D(inputData, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter)
-
-def Diff4th_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
-                     float regularisation_parameter,
-                     float edge_parameter,
-                     int iterationsNumb,                     
-                     float time_marching_parameter):
-    cdef long dims[2]
-    dims[0] = inputData.shape[0]
-    dims[1] = inputData.shape[1]
-    
-    cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
-            np.zeros([dims[0],dims[1]], dtype='float32')   
-    
-    # Run Anisotropic Fourth-Order diffusion for 2D data 
-    Diffus4th_CPU_main(&inputData[0,0], &outputData[0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, dims[1], dims[0], 1)
-    return outputData
-          
-def Diff4th_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
-                     float regularisation_parameter,
-                     float edge_parameter,
-                     int iterationsNumb,
-                     float time_marching_parameter):
-    cdef long dims[3]
-    dims[0] = inputData.shape[0]
-    dims[1] = inputData.shape[1]
-    dims[2] = inputData.shape[2]
-    
-    cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
-            np.zeros([dims[0],dims[1],dims[2]], dtype='float32')
-    
-    # Run Anisotropic Fourth-Order diffusion for  3D data 
-    Diffus4th_CPU_main(&inputData[0,0,0], &outputData[0,0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, dims[2], dims[1], dims[0])
-
-    return outputData
-
-#****************************************************************#
-#***************Patch-based weights calculation******************#
-#****************************************************************#
-def PATCHSEL_CPU(inputData, searchwindow, patchwindow, neighbours, edge_parameter):
-    if inputData.ndim == 2:
-        return PatchSel_2D(inputData, searchwindow, patchwindow, neighbours, edge_parameter)
-    elif inputData.ndim == 3:
-        return 1
-def PatchSel_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
-                     int searchwindow,
-                     int patchwindow,
-                     int neighbours,
-                     float edge_parameter):
-    cdef long dims[3]
-    dims[0] = neighbours
-    dims[1] = inputData.shape[0]
-    dims[2] = inputData.shape[1]
-    
-    
-    cdef np.ndarray[np.float32_t, ndim=3, mode="c"] Weights = \
-            np.zeros([dims[0], dims[1],dims[2]], dtype='float32')
-    
-    cdef np.ndarray[np.uint16_t, ndim=3, mode="c"] H_i = \
-            np.zeros([dims[0], dims[1],dims[2]], dtype='uint16')
-            
-    cdef np.ndarray[np.uint16_t, ndim=3, mode="c"] H_j = \
-            np.zeros([dims[0], dims[1],dims[2]], dtype='uint16')
-
-    # Run patch-based weight selection function
-    PatchSelect_CPU_main(&inputData[0,0], &H_j[0,0,0], &H_i[0,0,0], &H_i[0,0,0], &Weights[0,0,0], dims[2], dims[1], 0, searchwindow, patchwindow,  neighbours,  edge_parameter, 1)
-    return H_i, H_j, Weights
-"""
-def PatchSel_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
-                     int searchwindow,
-                     int patchwindow,
-                     int neighbours,
-                     float edge_parameter):
-    cdef long dims[4]
-    dims[0] = inputData.shape[0]
-    dims[1] = inputData.shape[1]
-    dims[2] = inputData.shape[2]
-    dims[3] = neighbours
-    
-    cdef np.ndarray[np.float32_t, ndim=4, mode="c"] Weights = \
-            np.zeros([dims[3],dims[0],dims[1],dims[2]], dtype='float32')
-    
-    cdef np.ndarray[np.uint16_t, ndim=4, mode="c"] H_i = \
-            np.zeros([dims[3],dims[0],dims[1],dims[2]], dtype='uint16')
-            
-    cdef np.ndarray[np.uint16_t, ndim=4, mode="c"] H_j = \
-            np.zeros([dims[3],dims[0],dims[1],dims[2]], dtype='uint16')
-            
-    cdef np.ndarray[np.uint16_t, ndim=4, mode="c"] H_k = \
-            np.zeros([dims[3],dims[0],dims[1],dims[2]], dtype='uint16')
-
-    # Run patch-based weight selection function
-    PatchSelect_CPU_main(&inputData[0,0,0], &H_i[0,0,0,0], &H_j[0,0,0,0], &H_k[0,0,0,0], &Weights[0,0,0,0], dims[2], dims[1], dims[0], searchwindow, patchwindow,  neighbours, edge_parameter, 1)
-    return H_i, H_j, H_k, Weights
-"""
-
-#****************************************************************#
-#***************Non-local Total Variation******************#
-#****************************************************************#
-def NLTV_CPU(inputData, H_i, H_j, H_k, Weights, regularisation_parameter, iterations):
-    if inputData.ndim == 2:
-        return NLTV_2D(inputData, H_i, H_j, Weights, regularisation_parameter, iterations)
-    elif inputData.ndim == 3:
-        return 1
-def NLTV_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
-                     np.ndarray[np.uint16_t, ndim=3, mode="c"] H_i,
-                     np.ndarray[np.uint16_t, ndim=3, mode="c"] H_j,
-                     np.ndarray[np.float32_t, ndim=3, mode="c"] Weights,
-                     float regularisation_parameter,
-                     int iterations):
-
-    cdef long dims[2]
-    dims[0] = inputData.shape[0]
-    dims[1] = inputData.shape[1]
-    neighbours = H_i.shape[0]
-    
-    cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
-            np.zeros([dims[0],dims[1]], dtype='float32')
-    
-    # Run nonlocal TV regularisation
-    Nonlocal_TV_CPU_main(&inputData[0,0], &outputData[0,0], &H_i[0,0,0], &H_j[0,0,0], &H_i[0,0,0], &Weights[0,0,0], dims[1], dims[0], 0, neighbours, regularisation_parameter, iterations)
-    return outputData
-
-#*********************Inpainting WITH****************************#
-#***************Nonlinear (Isotropic) Diffusion******************#
-#****************************************************************#
-def NDF_INPAINT_CPU(inputData, maskData, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type):
-    if inputData.ndim == 2:
-        return NDF_INP_2D(inputData, maskData, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type)
-    elif inputData.ndim == 3:
-        return NDF_INP_3D(inputData, maskData, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type)
-
-def NDF_INP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
-                     np.ndarray[np.uint8_t, ndim=2, mode="c"] maskData,
-                     float regularisation_parameter,
-                     float edge_parameter,
-                     int iterationsNumb,
-                     float time_marching_parameter,
-                     int penalty_type):
-
-    cdef long dims[2]
-    dims[0] = inputData.shape[0]
-    dims[1] = inputData.shape[1]
-
-
-    cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
-            np.zeros([dims[0],dims[1]], dtype='float32')
-    
-    # Run Inpaiting by Diffusion iterations for 2D data 
-    Diffusion_Inpaint_CPU_main(&inputData[0,0], &maskData[0,0], &outputData[0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type, dims[1], dims[0], 1)
-    return outputData
-            
-def NDF_INP_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
-                     np.ndarray[np.uint8_t, ndim=3, mode="c"] maskData,
-                     float regularisation_parameter,
-                     float edge_parameter,
-                     int iterationsNumb,
-                     float time_marching_parameter,
-                     int penalty_type):
-    cdef long dims[3]
-    dims[0] = inputData.shape[0]
-    dims[1] = inputData.shape[1]
-    dims[2] = inputData.shape[2]
-    
-    cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
-            np.zeros([dims[0],dims[1],dims[2]], dtype='float32')
-    
-    # Run Inpaiting by Diffusion iterations for 3D data 
-    Diffusion_Inpaint_CPU_main(&inputData[0,0,0], &maskData[0,0,0], &outputData[0,0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type, dims[2], dims[1], dims[0])
-
-    return outputData
-#*********************Inpainting WITH****************************#
-#***************Nonlocal Vertical Marching method****************#
-#****************************************************************#
-def NVM_INPAINT_CPU(inputData, maskData, SW_increment, iterationsNumb):
-    if inputData.ndim == 2:
-        return NVM_INP_2D(inputData, maskData, SW_increment, iterationsNumb)
-    elif inputData.ndim == 3:
-        return 
-
-def NVM_INP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
-               np.ndarray[np.uint8_t, ndim=2, mode="c"] maskData,
-                     int SW_increment,
-                     int iterationsNumb):
-    cdef long dims[2]
-    dims[0] = inputData.shape[0]
-    dims[1] = inputData.shape[1]
-    
-    cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
-            np.zeros([dims[0],dims[1]], dtype='float32')   
-    
-    cdef np.ndarray[np.uint8_t, ndim=2, mode="c"] maskData_upd = \
-            np.zeros([dims[0],dims[1]], dtype='uint8')
-    
-    # Run Inpaiting by Nonlocal vertical marching method for 2D data 
-    NonlocalMarching_Inpaint_main(&inputData[0,0], &maskData[0,0], &outputData[0,0], 
-                                  &maskData_upd[0,0],
-                                  SW_increment, iterationsNumb, 1, dims[1], dims[0], 1)
-    
-    return (outputData, maskData_upd)
-
-
-#****************************************************************#
-#***************Calculation of TV-energy functional**************#
-#****************************************************************#
-def TV_ENERGY(inputData, inputData0, regularisation_parameter, typeFunctional):
-    if inputData.ndim == 2:
-        return TV_ENERGY_2D(inputData, inputData0, regularisation_parameter, typeFunctional)
-    elif inputData.ndim == 3:
-        return TV_ENERGY_3D(inputData, inputData0, regularisation_parameter, typeFunctional)
-
-def TV_ENERGY_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
-                 np.ndarray[np.float32_t, ndim=2, mode="c"] inputData0, 
-                     float regularisation_parameter,
-                     int typeFunctional):
-    
-    cdef long dims[2]
-    dims[0] = inputData.shape[0]
-    dims[1] = inputData.shape[1]
-    
-    cdef np.ndarray[np.float32_t, ndim=1, mode="c"] outputData = \
-            np.zeros([1], dtype='float32')
-                   
-    # run function    
-    TV_energy2D(&inputData[0,0], &inputData0[0,0], &outputData[0], regularisation_parameter, typeFunctional, dims[1], dims[0])
-    
-    return outputData
-            
-def TV_ENERGY_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
-                 np.ndarray[np.float32_t, ndim=3, mode="c"] inputData0, 
-                     float regularisation_parameter,
-                     int typeFunctional):
-						 
-    cdef long dims[3]
-    dims[0] = inputData.shape[0]
-    dims[1] = inputData.shape[1]
-    dims[2] = inputData.shape[2]
-    
-    cdef np.ndarray[np.float32_t, ndim=1, mode="c"] outputData = \
-            np.zeros([1], dtype='float32')
-           
-    # Run function
-    TV_energy3D(&inputData[0,0,0], &inputData0[0,0,0], &outputData[0], regularisation_parameter, typeFunctional, dims[2], dims[1], dims[0])
-
-    return outputData
diff --git a/Wrappers/Python/src/gpu_regularisers.pyx b/Wrappers/Python/src/gpu_regularisers.pyx
deleted file mode 100644
index b52f669..0000000
--- a/Wrappers/Python/src/gpu_regularisers.pyx
+++ /dev/null
@@ -1,640 +0,0 @@
-# distutils: language=c++
-"""
-Copyright 2018 CCPi
-Licensed under the Apache License, Version 2.0 (the "License");
-you may not use this file except in compliance with the License.
-You may obtain a copy of the License at
-    http://www.apache.org/licenses/LICENSE-2.0
-Unless required by applicable law or agreed to in writing, software
-distributed under the License is distributed on an "AS IS" BASIS,
-WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-See the License for the specific language governing permissions and
-limitations under the License.
-
-Author: Edoardo Pasca, Daniil Kazantsev
-"""
-
-import cython
-import numpy as np
-cimport numpy as np
-
-CUDAErrorMessage = 'CUDA error'
-
-cdef extern int TV_ROF_GPU_main(float* Input, float* Output, float lambdaPar, int iter, float tau, int N, int M, int Z);
-cdef extern int TV_FGP_GPU_main(float *Input, float *Output, float lambdaPar, int iter, float epsil, int methodTV, int nonneg, int printM, int N, int M, int Z);
-cdef extern int TV_SB_GPU_main(float *Input, float *Output, float lambdaPar, int iter, float epsil, int methodTV, int printM, int N, int M, int Z);
-cdef extern int TGV_GPU_main(float *Input, float *Output, float lambdaPar, float alpha1, float alpha0, int iterationsNumb, float L2, int dimX, int dimY, int dimZ);
-cdef extern int LLT_ROF_GPU_main(float *Input, float *Output, float lambdaROF, float lambdaLLT, int iterationsNumb, float tau, int N, int M, int Z);
-cdef extern int NonlDiff_GPU_main(float *Input, float *Output, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int penaltytype, int N, int M, int Z);
-cdef extern int dTV_FGP_GPU_main(float *Input, float *InputRef, float *Output, float lambdaPar, int iterationsNumb, float epsil, float eta, int methodTV, int nonneg, int printM, int N, int M, int Z);
-cdef extern int Diffus4th_GPU_main(float *Input, float *Output, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int N, int M, int Z);
-cdef extern int PatchSelect_GPU_main(float *Input, unsigned short *H_i, unsigned short *H_j, float *Weights, int N, int M, int SearchWindow, int SimilarWin, int NumNeighb, float h);
-
-# Total-variation Rudin-Osher-Fatemi (ROF)
-def TV_ROF_GPU(inputData,
-                     regularisation_parameter,
-                     iterations, 
-                     time_marching_parameter):
-    if inputData.ndim == 2:
-        return ROFTV2D(inputData, 
-                     regularisation_parameter,
-                     iterations,
-                     time_marching_parameter)
-    elif inputData.ndim == 3:
-        return ROFTV3D(inputData, 
-                     regularisation_parameter,
-                     iterations, 
-                     time_marching_parameter)
-                     
-# Total-variation Fast-Gradient-Projection (FGP)
-def TV_FGP_GPU(inputData,
-                     regularisation_parameter,
-                     iterations, 
-                     tolerance_param,
-                     methodTV,
-                     nonneg,
-                     printM):
-    if inputData.ndim == 2:
-        return FGPTV2D(inputData,
-                     regularisation_parameter,
-                     iterations, 
-                     tolerance_param,
-                     methodTV,
-                     nonneg,
-                     printM)
-    elif inputData.ndim == 3:
-        return FGPTV3D(inputData,
-                     regularisation_parameter,
-                     iterations, 
-                     tolerance_param,
-                     methodTV,
-                     nonneg,
-                     printM)
-# Total-variation Split Bregman (SB)
-def TV_SB_GPU(inputData,
-                     regularisation_parameter,
-                     iterations, 
-                     tolerance_param,
-                     methodTV,
-                     printM):
-    if inputData.ndim == 2:
-        return SBTV2D(inputData,
-                     regularisation_parameter,
-                     iterations, 
-                     tolerance_param,
-                     methodTV,
-                     printM)
-    elif inputData.ndim == 3:
-        return SBTV3D(inputData,
-                     regularisation_parameter,
-                     iterations, 
-                     tolerance_param,
-                     methodTV,
-                     printM)
-# LLT-ROF model
-def LLT_ROF_GPU(inputData, regularisation_parameterROF, regularisation_parameterLLT, iterations, time_marching_parameter):
-    if inputData.ndim == 2:
-        return LLT_ROF_GPU2D(inputData, regularisation_parameterROF, regularisation_parameterLLT, iterations, time_marching_parameter)
-    elif inputData.ndim == 3:
-        return LLT_ROF_GPU3D(inputData, regularisation_parameterROF, regularisation_parameterLLT, iterations, time_marching_parameter)
-# Total Generilised Variation (TGV)
-def TGV_GPU(inputData, regularisation_parameter, alpha1, alpha0, iterations, LipshitzConst):
-    if inputData.ndim == 2:
-        return TGV2D(inputData, regularisation_parameter, alpha1, alpha0, iterations, LipshitzConst)
-    elif inputData.ndim == 3:
-        return TGV3D(inputData, regularisation_parameter, alpha1, alpha0, iterations, LipshitzConst)
-# Directional Total-variation Fast-Gradient-Projection (FGP)
-def dTV_FGP_GPU(inputData,
-                     refdata,
-                     regularisation_parameter,
-                     iterations, 
-                     tolerance_param,
-                     eta_const,
-                     methodTV,
-                     nonneg,
-                     printM):
-    if inputData.ndim == 2:
-        return FGPdTV2D(inputData,
-                     refdata,
-                     regularisation_parameter,
-                     iterations, 
-                     tolerance_param,
-                     eta_const,
-                     methodTV,
-                     nonneg,
-                     printM)
-    elif inputData.ndim == 3:
-        return FGPdTV3D(inputData,
-                     refdata,
-                     regularisation_parameter,
-                     iterations, 
-                     tolerance_param,
-                     eta_const,
-                     methodTV,
-                     nonneg,
-                     printM)
-# Nonlocal Isotropic Diffusion (NDF)
-def NDF_GPU(inputData,
-                     regularisation_parameter,
-                     edge_parameter,
-                     iterations, 
-                     time_marching_parameter,
-                     penalty_type):
-    if inputData.ndim == 2:
-        return NDF_GPU_2D(inputData,
-                     regularisation_parameter,
-                     edge_parameter,
-                     iterations, 
-                     time_marching_parameter,
-                     penalty_type)
-    elif inputData.ndim == 3:
-        return NDF_GPU_3D(inputData,
-                     regularisation_parameter,
-                     edge_parameter,
-                     iterations, 
-                     time_marching_parameter,
-                     penalty_type)
-# Anisotropic Fourth-Order diffusion
-def Diff4th_GPU(inputData,
-                     regularisation_parameter,
-                     edge_parameter,
-                     iterations, 
-                     time_marching_parameter):
-    if inputData.ndim == 2:
-        return Diff4th_2D(inputData,
-                     regularisation_parameter,
-                     edge_parameter,
-                     iterations, 
-                     time_marching_parameter)
-    elif inputData.ndim == 3:
-        return Diff4th_3D(inputData,
-                     regularisation_parameter,
-                     edge_parameter,
-                     iterations, 
-                     time_marching_parameter)
-                     
-#****************************************************************#
-#********************** Total-variation ROF *********************#
-#****************************************************************#
-def ROFTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
-                     float regularisation_parameter,
-                     int iterations, 
-                     float time_marching_parameter):
-    
-    cdef long dims[2]
-    dims[0] = inputData.shape[0]
-    dims[1] = inputData.shape[1]
-
-    cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
-		    np.zeros([dims[0],dims[1]], dtype='float32')
-          
-    # Running CUDA code here
-    if (TV_ROF_GPU_main(
-            &inputData[0,0], &outputData[0,0], 
-                       regularisation_parameter,
-                       iterations , 
-                       time_marching_parameter, 
-                       dims[1], dims[0], 1)==0):
-        return outputData;
-    else:
-        raise ValueError(CUDAErrorMessage);
-    
-def ROFTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
-                     float regularisation_parameter,
-                     int iterations, 
-                     float time_marching_parameter):
-    
-    cdef long dims[3]
-    dims[0] = inputData.shape[0]
-    dims[1] = inputData.shape[1]
-    dims[2] = inputData.shape[2]
-
-    cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
-		    np.zeros([dims[0],dims[1],dims[2]], dtype='float32')
-          
-    # Running CUDA code here    
-    if (TV_ROF_GPU_main(
-            &inputData[0,0,0], &outputData[0,0,0], 
-                       regularisation_parameter,
-                       iterations , 
-                       time_marching_parameter, 
-                       dims[2], dims[1], dims[0])==0):
-        return outputData;
-    else:
-        raise ValueError(CUDAErrorMessage);
-#****************************************************************#
-#********************** Total-variation FGP *********************#
-#****************************************************************#
-#******** Total-variation Fast-Gradient-Projection (FGP)*********#
-def FGPTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
-                     float regularisation_parameter,
-                     int iterations, 
-                     float tolerance_param,
-                     int methodTV,
-                     int nonneg,
-                     int printM):
-    
-    cdef long dims[2]
-    dims[0] = inputData.shape[0]
-    dims[1] = inputData.shape[1]
-
-    cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
-		    np.zeros([dims[0],dims[1]], dtype='float32')
-          
-    # Running CUDA code here    
-    if (TV_FGP_GPU_main(&inputData[0,0], &outputData[0,0],
-                       regularisation_parameter, 
-                       iterations, 
-                       tolerance_param,
-                       methodTV,
-                       nonneg,
-                       printM,
-                       dims[1], dims[0], 1)==0):
-        return outputData;
-    else:
-        raise ValueError(CUDAErrorMessage);
-
-    
-def FGPTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
-                     float regularisation_parameter,
-                     int iterations, 
-                     float tolerance_param,
-                     int methodTV,
-                     int nonneg,
-                     int printM):
-    
-    cdef long dims[3]
-    dims[0] = inputData.shape[0]
-    dims[1] = inputData.shape[1]
-    dims[2] = inputData.shape[2]
-
-    cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
-		    np.zeros([dims[0],dims[1],dims[2]], dtype='float32')
-          
-    # Running CUDA code here    
-    if (TV_FGP_GPU_main(&inputData[0,0,0], &outputData[0,0,0],
-                       regularisation_parameter , 
-                       iterations, 
-                       tolerance_param,
-                       methodTV,
-                       nonneg,
-                       printM,
-                       dims[2], dims[1], dims[0])==0):
-        return outputData;
-    else:
-        raise ValueError(CUDAErrorMessage);
-
-#***************************************************************#
-#********************** Total-variation SB *********************#
-#***************************************************************#
-#*************** Total-variation Split Bregman (SB)*************#
-def SBTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
-                     float regularisation_parameter,
-                     int iterations, 
-                     float tolerance_param,
-                     int methodTV,
-                     int printM):
-    
-    cdef long dims[2]
-    dims[0] = inputData.shape[0]
-    dims[1] = inputData.shape[1]
-
-    cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
-		    np.zeros([dims[0],dims[1]], dtype='float32')
-          
-    # Running CUDA code here    
-    if (TV_SB_GPU_main(&inputData[0,0], &outputData[0,0],
-                       regularisation_parameter, 
-                       iterations, 
-                       tolerance_param,
-                       methodTV,
-                       printM,
-                       dims[1], dims[0], 1)==0):
-        return outputData;
-    else:
-        raise ValueError(CUDAErrorMessage);
-
-    
-def SBTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
-                     float regularisation_parameter,
-                     int iterations, 
-                     float tolerance_param,
-                     int methodTV,
-                     int printM):
-    
-    cdef long dims[3]
-    dims[0] = inputData.shape[0]
-    dims[1] = inputData.shape[1]
-    dims[2] = inputData.shape[2]
-
-    cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
-		    np.zeros([dims[0],dims[1],dims[2]], dtype='float32')
-          
-    # Running CUDA code here    
-    if (TV_SB_GPU_main(&inputData[0,0,0], &outputData[0,0,0],
-                       regularisation_parameter , 
-                       iterations, 
-                       tolerance_param,
-                       methodTV,
-                       printM,
-                       dims[2], dims[1], dims[0])==0):
-        return outputData;
-    else:
-        raise ValueError(CUDAErrorMessage);
-
-
-#***************************************************************#
-#************************ LLT-ROF model ************************#
-#***************************************************************#
-#************Joint LLT-ROF model for higher order **************#
-def LLT_ROF_GPU2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
-                     float regularisation_parameterROF,
-                     float regularisation_parameterLLT,
-                     int iterations, 
-                     float time_marching_parameter):
-    
-    cdef long dims[2]
-    dims[0] = inputData.shape[0]
-    dims[1] = inputData.shape[1]
-
-    cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
-		    np.zeros([dims[0],dims[1]], dtype='float32')
-          
-    # Running CUDA code here    
-    if (LLT_ROF_GPU_main(&inputData[0,0], &outputData[0,0],regularisation_parameterROF, regularisation_parameterLLT, iterations, time_marching_parameter, dims[1],dims[0],1)==0):
-        return outputData;
-    else:
-        raise ValueError(CUDAErrorMessage);
-
-    
-def LLT_ROF_GPU3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
-                     float regularisation_parameterROF,
-                     float regularisation_parameterLLT,
-                     int iterations, 
-                     float time_marching_parameter):
-    
-    cdef long dims[3]
-    dims[0] = inputData.shape[0]
-    dims[1] = inputData.shape[1]
-    dims[2] = inputData.shape[2]
-
-    cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
-		    np.zeros([dims[0],dims[1],dims[2]], dtype='float32')
-          
-    # Running CUDA code here    
-    if (LLT_ROF_GPU_main(&inputData[0,0,0], &outputData[0,0,0], regularisation_parameterROF, regularisation_parameterLLT, iterations, time_marching_parameter, dims[2], dims[1], dims[0])==0):
-        return outputData;
-    else:
-        raise ValueError(CUDAErrorMessage);
-
-
-#***************************************************************#
-#***************** Total Generalised Variation *****************#
-#***************************************************************#
-def TGV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
-                     float regularisation_parameter,
-                     float alpha1,
-                     float alpha0,
-                     int iterationsNumb, 
-                     float LipshitzConst):
-                         
-    cdef long dims[2]
-    dims[0] = inputData.shape[0]
-    dims[1] = inputData.shape[1]
-    
-    cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
-            np.zeros([dims[0],dims[1]], dtype='float32')
-                   
-    #/* Run TGV iterations for 2D data */
-    if (TGV_GPU_main(&inputData[0,0], &outputData[0,0], regularisation_parameter,
-                       alpha1,
-                       alpha0,
-                       iterationsNumb, 
-                       LipshitzConst,
-                       dims[1],dims[0], 1)==0):
-        return outputData
-    else:
-        raise ValueError(CUDAErrorMessage);
-
-def TGV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
-                     float regularisation_parameter,
-                     float alpha1,
-                     float alpha0,
-                     int iterationsNumb, 
-                     float LipshitzConst):
-    
-    cdef long dims[3]
-    dims[0] = inputData.shape[0]
-    dims[1] = inputData.shape[1]
-    dims[2] = inputData.shape[2]
-
-    cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
-		    np.zeros([dims[0],dims[1],dims[2]], dtype='float32')
-          
-    # Running CUDA code here    
-    if (TGV_GPU_main(
-            &inputData[0,0,0], &outputData[0,0,0], regularisation_parameter,
-                       alpha1,
-                       alpha0,
-                       iterationsNumb, 
-                       LipshitzConst,
-                       dims[2], dims[1], dims[0])==0):
-        return outputData;
-    else:
-        raise ValueError(CUDAErrorMessage);
-
-
-#****************************************************************#
-#**************Directional Total-variation FGP ******************#
-#****************************************************************#
-#******** Directional TV Fast-Gradient-Projection (FGP)*********#
-def FGPdTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
-             np.ndarray[np.float32_t, ndim=2, mode="c"] refdata,
-                     float regularisation_parameter,
-                     int iterations, 
-                     float tolerance_param,
-                     float eta_const,
-                     int methodTV,
-                     int nonneg,
-                     int printM):
-    
-    cdef long dims[2]
-    dims[0] = inputData.shape[0]
-    dims[1] = inputData.shape[1]
-
-    cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
-		    np.zeros([dims[0],dims[1]], dtype='float32')
-          
-    # Running CUDA code here    
-    if (dTV_FGP_GPU_main(&inputData[0,0], &refdata[0,0], &outputData[0,0],
-                       regularisation_parameter, 
-                       iterations, 
-                       tolerance_param,
-                       eta_const,
-                       methodTV,
-                       nonneg,
-                       printM,
-                       dims[1], dims[0], 1)==0):
-        return outputData
-    else:
-        raise ValueError(CUDAErrorMessage);
-
-    
-def FGPdTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
-             np.ndarray[np.float32_t, ndim=3, mode="c"] refdata, 
-                     float regularisation_parameter,
-                     int iterations, 
-                     float tolerance_param,
-                     float eta_const,
-                     int methodTV,
-                     int nonneg,
-                     int printM):
-    
-    cdef long dims[3]
-    dims[0] = inputData.shape[0]
-    dims[1] = inputData.shape[1]
-    dims[2] = inputData.shape[2]
-
-    cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
-		    np.zeros([dims[0],dims[1],dims[2]], dtype='float32')
-          
-    # Running CUDA code here    
-    if (dTV_FGP_GPU_main(&inputData[0,0,0], &refdata[0,0,0], &outputData[0,0,0],
-                       regularisation_parameter , 
-                       iterations, 
-                       tolerance_param,
-                       eta_const,
-                       methodTV,
-                       nonneg,
-                       printM,
-                       dims[2], dims[1], dims[0])==0):
-        return outputData;
-    else:
-        raise ValueError(CUDAErrorMessage);
-
-
-#****************************************************************#
-#***************Nonlinear (Isotropic) Diffusion******************#
-#****************************************************************#
-def NDF_GPU_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
-                     float regularisation_parameter,
-                     float edge_parameter,
-                     int iterationsNumb,                     
-                     float time_marching_parameter,
-                     int penalty_type):
-    cdef long dims[2]
-    dims[0] = inputData.shape[0]
-    dims[1] = inputData.shape[1]
-    
-    cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
-            np.zeros([dims[0],dims[1]], dtype='float32')
-    
-    #rangecheck = penalty_type < 1 and penalty_type > 3
-    #if not rangecheck:
-#        raise ValueError('Choose penalty type as 1 for Huber, 2 - Perona-Malik, 3 - Tukey Biweight')
-    
-    # Run Nonlinear Diffusion iterations for 2D data 
-    # Running CUDA code here  
-    if (NonlDiff_GPU_main(&inputData[0,0], &outputData[0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type, dims[1], dims[0], 1)==0):
-        return outputData;
-    else:
-        raise ValueError(CUDAErrorMessage);
-
-            
-def NDF_GPU_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
-                     float regularisation_parameter,
-                     float edge_parameter,
-                     int iterationsNumb,                     
-                     float time_marching_parameter,
-                     int penalty_type):
-    cdef long dims[3]
-    dims[0] = inputData.shape[0]
-    dims[1] = inputData.shape[1]
-    dims[2] = inputData.shape[2]
-    
-    cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
-            np.zeros([dims[0],dims[1],dims[2]], dtype='float32')    
-       
-    # Run Nonlinear Diffusion iterations for  3D data 
-    # Running CUDA code here  
-    if (NonlDiff_GPU_main(&inputData[0,0,0], &outputData[0,0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type, dims[2], dims[1], dims[0])==0):
-        return outputData;
-    else:
-        raise ValueError(CUDAErrorMessage);
-
-#****************************************************************#
-#************Anisotropic Fourth-Order diffusion******************#
-#****************************************************************#
-def Diff4th_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 
-                     float regularisation_parameter,
-                     float edge_parameter,
-                     int iterationsNumb,
-                     float time_marching_parameter):
-    cdef long dims[2]
-    dims[0] = inputData.shape[0]
-    dims[1] = inputData.shape[1]
-    
-    cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
-            np.zeros([dims[0],dims[1]], dtype='float32')
-    
-    # Run Anisotropic Fourth-Order diffusion for 2D data 
-    # Running CUDA code here  
-    if (Diffus4th_GPU_main(&inputData[0,0], &outputData[0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, dims[1], dims[0], 1)==0):
-        return outputData
-    else:
-        raise ValueError(CUDAErrorMessage);
-
-            
-def Diff4th_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
-                     float regularisation_parameter,
-                     float edge_parameter,
-                     int iterationsNumb,
-                     float time_marching_parameter):
-    cdef long dims[3]
-    dims[0] = inputData.shape[0]
-    dims[1] = inputData.shape[1]
-    dims[2] = inputData.shape[2]
-    
-    cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
-            np.zeros([dims[0],dims[1],dims[2]], dtype='float32')    
-       
-    # Run Anisotropic Fourth-Order diffusion for  3D data 
-    # Running CUDA code here  
-    if (Diffus4th_GPU_main(&inputData[0,0,0], &outputData[0,0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, dims[2], dims[1], dims[0])==0):
-        return outputData;
-    else:
-        raise ValueError(CUDAErrorMessage);
-
-#****************************************************************#
-#************Patch-based weights pre-selection******************#
-#****************************************************************#
-def PATCHSEL_GPU(inputData, searchwindow, patchwindow, neighbours, edge_parameter):
-    if inputData.ndim == 2:
-        return PatchSel_2D(inputData, searchwindow, patchwindow, neighbours, edge_parameter)
-    elif inputData.ndim == 3:
-        return 1
-def PatchSel_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
-                     int searchwindow,
-                     int patchwindow,
-                     int neighbours,
-                     float edge_parameter):
-    cdef long dims[3]
-    dims[0] = neighbours
-    dims[1] = inputData.shape[0]
-    dims[2] = inputData.shape[1]    
-    
-    cdef np.ndarray[np.float32_t, ndim=3, mode="c"] Weights = \
-            np.zeros([dims[0], dims[1],dims[2]], dtype='float32')
-    
-    cdef np.ndarray[np.uint16_t, ndim=3, mode="c"] H_i = \
-            np.zeros([dims[0], dims[1],dims[2]], dtype='uint16')
-            
-    cdef np.ndarray[np.uint16_t, ndim=3, mode="c"] H_j = \
-            np.zeros([dims[0], dims[1],dims[2]], dtype='uint16')
-
-    # Run patch-based weight selection function
-    if (PatchSelect_GPU_main(&inputData[0,0], &H_j[0,0,0], &H_i[0,0,0], &Weights[0,0,0], dims[2], dims[1], searchwindow, patchwindow,  neighbours,  edge_parameter)==0):
-        return H_i, H_j, Weights;
-    else:
-        raise ValueError(CUDAErrorMessage);
-