diff options
| author | Tomas Kulhanek <tomas.kulhanek@stfc.ac.uk> | 2019-02-21 02:10:14 -0500 |
|---|---|---|
| committer | Tomas Kulhanek <tomas.kulhanek@stfc.ac.uk> | 2019-02-21 02:10:14 -0500 |
| commit | 3caa686662f7d937cf7eb852dde437cd66e79a6e (patch) | |
| tree | 76088f5924ff9278e0a37140fce888cd89b84a7e /Wrappers/Python | |
| parent | 8f2e86726669b9dadb3c788e0ea681d397a2eeb7 (diff) | |
| download | regularization-3caa686662f7d937cf7eb852dde437cd66e79a6e.tar.gz regularization-3caa686662f7d937cf7eb852dde437cd66e79a6e.tar.bz2 regularization-3caa686662f7d937cf7eb852dde437cd66e79a6e.tar.xz regularization-3caa686662f7d937cf7eb852dde437cd66e79a6e.zip | |
restructured sources
Diffstat (limited to 'Wrappers/Python')
19 files changed, 0 insertions, 5668 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/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/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/demo_cpu_inpainters.py b/Wrappers/Python/demos/demo_cpu_inpainters.py deleted file mode 100644 index 3b4191b..0000000 --- a/Wrappers/Python/demos/demo_cpu_inpainters.py +++ /dev/null @@ -1,192 +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 qualitymetrics import rmse -############################################################################### -def printParametersToString(pars): - txt = r'' - for key, value in pars.items(): - if key== 'algorithm' : - txt += "{0} = {1}".format(key, value.__name__) - elif key == 'input': - txt += "{0} = {1}".format(key, np.shape(value)) - elif key == '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']) - -rms = rmse(sino_full, ndf_inp_linear) -pars['rmse'] = rms - -txtstr = printParametersToString(pars) -txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) -print (txtstr) -a=fig.add_subplot(1,2,2) - -# these are matplotlib.patch.Patch properties -props = dict(boxstyle='round', facecolor='wheat', alpha=0.75) -# place a text box in upper left in axes coords -a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(ndf_inp_linear, cmap="gray") -plt.title('{}'.format('Linear diffusion inpainting results')) -#%% -print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") -print ("_Inpainting using nonlinear diffusion (2D)_") -print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") - -## plot -fig = plt.figure(3) -plt.suptitle('Performance of nonlinear diffusion inpainting using the CPU') -a=fig.add_subplot(1,2,1) -a.set_title('Missing data sinogram') -imgplot = plt.imshow(sino_cut_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']) - -rms = rmse(sino_full, ndf_inp_nonlinear) -pars['rmse'] = rms - -txtstr = printParametersToString(pars) -txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) -print (txtstr) -a=fig.add_subplot(1,2,2) - -# these are matplotlib.patch.Patch properties -props = dict(boxstyle='round', facecolor='wheat', alpha=0.75) -# place a text box in upper left in axes coords -a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(ndf_inp_nonlinear, cmap="gray") -plt.title('{}'.format('Nonlinear diffusion inpainting results')) -#%% -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']) - -rms = rmse(sino_full, nvm_inp) -pars['rmse'] = rms - -txtstr = printParametersToString(pars) -txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) -print (txtstr) -a=fig.add_subplot(1,2,2) - -# these are matplotlib.patch.Patch properties -props = dict(boxstyle='round', facecolor='wheat', alpha=0.75) -# place a text box in upper left in axes coords -a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(nvm_inp, cmap="gray") -plt.title('{}'.format('Nonlocal Vertical Marching inpainting results')) -#%% diff --git a/Wrappers/Python/demos/demo_cpu_regularisers.py b/Wrappers/Python/demos/demo_cpu_regularisers.py deleted file mode 100644 index e6befa9..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 qualitymetrics import rmse -############################################################################### -def printParametersToString(pars): - txt = r'' - for key, value in pars.items(): - if key== 'algorithm' : - txt += "{0} = {1}".format(key, value.__name__) - elif key == 'input': - txt += "{0} = {1}".format(key, np.shape(value)) - elif key == 'refdata': - txt += "{0} = {1}".format(key, np.shape(value)) - else: - txt += "{0} = {1}".format(key, value) - txt += '\n' - return txt -############################################################################### -#%% -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.04,\ - 'number_of_iterations': 1200,\ - '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') -rms = rmse(Im, rof_cpu) -pars['rmse'] = rms - -txtstr = printParametersToString(pars) -txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) -print (txtstr) -a=fig.add_subplot(1,2,2) - -# these are matplotlib.patch.Patch properties -props = dict(boxstyle='round', facecolor='wheat', alpha=0.75) -# place a text box in upper left in axes coords -a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(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') - - -rms = rmse(Im, fgp_cpu) -pars['rmse'] = rms - -txtstr = printParametersToString(pars) -txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) -print (txtstr) -a=fig.add_subplot(1,2,2) - -# these are matplotlib.patch.Patch properties -props = dict(boxstyle='round', facecolor='wheat', alpha=0.75) -# place a text box in upper left in axes coords -a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(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') - - -rms = rmse(Im, sb_cpu) -pars['rmse'] = rms - -txtstr = printParametersToString(pars) -txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) -print (txtstr) -a=fig.add_subplot(1,2,2) - -# these are matplotlib.patch.Patch properties -props = dict(boxstyle='round', facecolor='wheat', alpha=0.75) -# place a text box in upper left in axes coords -a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(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') - - -rms = rmse(Im, tgv_cpu) -pars['rmse'] = rms - -txtstr = printParametersToString(pars) -txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) -print (txtstr) -a=fig.add_subplot(1,2,2) - -# these are matplotlib.patch.Patch properties -props = dict(boxstyle='round', facecolor='wheat', alpha=0.75) -# place a text box in upper left in axes coords -a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(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') - -rms = rmse(Im, lltrof_cpu) -pars['rmse'] = rms - -txtstr = printParametersToString(pars) -txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) -print (txtstr) -a=fig.add_subplot(1,2,2) - -# these are matplotlib.patch.Patch properties -props = dict(boxstyle='round', facecolor='wheat', alpha=0.75) -# place a text box in upper left in axes coords -a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(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') - -rms = rmse(Im, ndf_cpu) -pars['rmse'] = rms - -txtstr = printParametersToString(pars) -txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) -print (txtstr) -a=fig.add_subplot(1,2,2) - -# these are matplotlib.patch.Patch properties -props = dict(boxstyle='round', facecolor='wheat', alpha=0.75) -# place a text box in upper left in axes coords -a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(ndf_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') - -rms = rmse(Im, diff4_cpu) -pars['rmse'] = rms - -txtstr = printParametersToString(pars) -txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) -print (txtstr) -a=fig.add_subplot(1,2,2) - -# these are matplotlib.patch.Patch properties -props = dict(boxstyle='round', facecolor='wheat', alpha=0.75) -# place a text box in upper left in axes coords -a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(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']) - -rms = rmse(Im, nltv_cpu) -pars['rmse'] = rms - -txtstr = printParametersToString(pars) -txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) -print (txtstr) -a=fig.add_subplot(1,2,2) - -# these are matplotlib.patch.Patch properties -props = dict(boxstyle='round', facecolor='wheat', alpha=0.75) -# place a text box in upper left in axes coords -a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(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') - -rms = rmse(Im, fgp_dtv_cpu) -pars['rmse'] = rms - -txtstr = printParametersToString(pars) -txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) -print (txtstr) -a=fig.add_subplot(1,2,2) - -# these are matplotlib.patch.Patch properties -props = dict(boxstyle='round', facecolor='wheat', alpha=0.75) -# place a text box in upper left in axes coords -a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(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']) - -rms = rmse(idealVol, tnv_cpu) -pars['rmse'] = rms - -txtstr = printParametersToString(pars) -txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) -print (txtstr) -a=fig.add_subplot(1,2,2) - -# these are matplotlib.patch.Patch properties -props = dict(boxstyle='round', facecolor='wheat', alpha=0.75) -# place a text box in upper left in axes coords -a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(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 2d2fc22..0000000 --- a/Wrappers/Python/demos/demo_cpu_regularisers3D.py +++ /dev/null @@ -1,458 +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 qualitymetrics import rmse -############################################################################### -def printParametersToString(pars): - txt = r'' - for key, value in pars.items(): - if key== 'algorithm' : - txt += "{0} = {1}".format(key, value.__name__) - elif key == 'input': - txt += "{0} = {1}".format(key, np.shape(value)) - elif key == 'refdata': - txt += "{0} = {1}".format(key, np.shape(value)) - else: - txt += "{0} = {1}".format(key, value) - txt += '\n' - return txt -############################################################################### -#%% -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') -rms = rmse(idealVol, rof_cpu3D) -pars['rmse'] = rms - -txtstr = printParametersToString(pars) -txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) -print (txtstr) -a=fig.add_subplot(1,2,2) - -# these are matplotlib.patch.Patch properties -props = dict(boxstyle='round', facecolor='wheat', alpha=0.75) -# place a text box in upper left in axes coords -a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(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') - - -rms = rmse(idealVol, fgp_cpu3D) -pars['rmse'] = rms - -txtstr = printParametersToString(pars) -txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) -print (txtstr) -a=fig.add_subplot(1,2,2) - -# these are matplotlib.patch.Patch properties -props = dict(boxstyle='round', facecolor='wheat', alpha=0.75) -# place a text box in upper left in axes coords -a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(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') - -rms = rmse(idealVol, sb_cpu3D) -pars['rmse'] = rms - -txtstr = printParametersToString(pars) -txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) -print (txtstr) -a=fig.add_subplot(1,2,2) - -# these are matplotlib.patch.Patch properties -props = dict(boxstyle='round', facecolor='wheat', alpha=0.75) -# place a text box in upper left in axes coords -a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(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') - -rms = rmse(idealVol, lltrof_cpu3D) -pars['rmse'] = rms - -txtstr = printParametersToString(pars) -txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) -print (txtstr) -a=fig.add_subplot(1,2,2) - -# these are matplotlib.patch.Patch properties -props = dict(boxstyle='round', facecolor='wheat', alpha=0.75) -# place a text box in upper left in axes coords -a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(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') - - -rms = rmse(idealVol, tgv_cpu3D) -pars['rmse'] = rms - -txtstr = printParametersToString(pars) -txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) -print (txtstr) -a=fig.add_subplot(1,2,2) - -# these are matplotlib.patch.Patch properties -props = dict(boxstyle='round', facecolor='wheat', alpha=0.75) -# place a text box in upper left in axes coords -a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(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']) - -rms = rmse(idealVol, ndf_cpu3D) -pars['rmse'] = rms - -txtstr = printParametersToString(pars) -txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) -print (txtstr) -a=fig.add_subplot(1,2,2) - -# these are matplotlib.patch.Patch properties -props = dict(boxstyle='round', facecolor='wheat', alpha=0.75) -# place a text box in upper left in axes coords -a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(ndf_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']) - -rms = rmse(idealVol, diff4th_cpu3D) -pars['rmse'] = rms - -txtstr = printParametersToString(pars) -txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) -print (txtstr) -a=fig.add_subplot(1,2,2) - -# these are matplotlib.patch.Patch properties -props = dict(boxstyle='round', facecolor='wheat', alpha=0.75) -# place a text box in upper left in axes coords -a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(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') - - -rms = rmse(idealVol, fgp_dTV_cpu3D) -pars['rmse'] = rms - -txtstr = printParametersToString(pars) -txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) -print (txtstr) -a=fig.add_subplot(1,2,2) - -# these are matplotlib.patch.Patch properties -props = dict(boxstyle='round', facecolor='wheat', alpha=0.75) -# place a text box in upper left in axes coords -a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(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 230a761..0000000 --- a/Wrappers/Python/demos/demo_cpu_vs_gpu_regularisers.py +++ /dev/null @@ -1,790 +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 qualitymetrics import rmse -############################################################################### -def printParametersToString(pars): - txt = r'' - for key, value in pars.items(): - if key== 'algorithm' : - txt += "{0} = {1}".format(key, value.__name__) - elif key == 'input': - txt += "{0} = {1}".format(key, np.shape(value)) - elif key == 'refdata': - txt += "{0} = {1}".format(key, np.shape(value)) - else: - txt += "{0} = {1}".format(key, value) - txt += '\n' - return txt -############################################################################### - -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') -rms = rmse(Im, rof_cpu) -pars['rmse'] = rms - -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') - -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) -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') - - -rms = rmse(Im, fgp_cpu) -pars['rmse'] = rms - -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') - -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) -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') - - -rms = rmse(Im, sb_cpu) -pars['rmse'] = rms - -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') - -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) -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') - -rms = rmse(Im, tgv_cpu) -pars['rmse'] = rms - -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') - -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) -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') - -rms = rmse(Im, lltrof_cpu) -pars['rmse'] = rms - -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') - -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) -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') - -rms = rmse(Im, ndf_cpu) -pars['rmse'] = rms - -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') - -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) -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') - -rms = rmse(Im, diff4th_cpu) -pars['rmse'] = rms - -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') - -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) -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') - - -rms = rmse(Im, fgp_dtv_cpu) -pars['rmse'] = rms - -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') -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) -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 e1c6575..0000000 --- a/Wrappers/Python/demos/demo_gpu_regularisers.py +++ /dev/null @@ -1,518 +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 qualitymetrics import rmse -############################################################################### -def printParametersToString(pars): - txt = r'' - for key, value in pars.items(): - if key== 'algorithm' : - txt += "{0} = {1}".format(key, value.__name__) - elif key == 'input': - txt += "{0} = {1}".format(key, np.shape(value)) - elif key == 'refdata': - txt += "{0} = {1}".format(key, np.shape(value)) - else: - txt += "{0} = {1}".format(key, value) - txt += '\n' - return txt -############################################################################### -#%% -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') - -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) -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') - -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) -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') - -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) -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') - - -rms = rmse(Im, tgv_gpu) -pars['rmse'] = rms - -txtstr = printParametersToString(pars) -txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) -print (txtstr) -a=fig.add_subplot(1,2,2) - -# these are matplotlib.patch.Patch properties -props = dict(boxstyle='round', facecolor='wheat', alpha=0.75) -# place a text box in upper left in axes coords -a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(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') - - -rms = rmse(Im, lltrof_gpu) -pars['rmse'] = rms - -txtstr = printParametersToString(pars) -txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) -print (txtstr) -a=fig.add_subplot(1,2,2) - -# these are matplotlib.patch.Patch properties -props = dict(boxstyle='round', facecolor='wheat', alpha=0.75) -# place a text box in upper left in axes coords -a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(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') - -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) -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') - -rms = rmse(Im, diff4_gpu) -pars['rmse'] = rms - -txtstr = printParametersToString(pars) -txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) -print (txtstr) -a=fig.add_subplot(1,2,2) - -# these are matplotlib.patch.Patch properties -props = dict(boxstyle='round', facecolor='wheat', alpha=0.75) -# place a text box in upper left in axes coords -a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(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']) - -rms = rmse(Im, nltv_cpu) -pars['rmse'] = rms - -txtstr = printParametersToString(pars) -txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) -print (txtstr) -a=fig.add_subplot(1,2,2) - -# these are matplotlib.patch.Patch properties -props = dict(boxstyle='round', facecolor='wheat', alpha=0.75) -# place a text box in upper left in axes coords -a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(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') - -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) -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 b6058d2..0000000 --- a/Wrappers/Python/demos/demo_gpu_regularisers3D.py +++ /dev/null @@ -1,460 +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 qualitymetrics import rmse -############################################################################### -def printParametersToString(pars): - txt = r'' - for key, value in pars.items(): - if key== 'algorithm' : - txt += "{0} = {1}".format(key, value.__name__) - elif key == 'input': - txt += "{0} = {1}".format(key, np.shape(value)) - elif key == 'refdata': - txt += "{0} = {1}".format(key, np.shape(value)) - else: - txt += "{0} = {1}".format(key, value) - txt += '\n' - return txt -############################################################################### -#%% -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') -rms = rmse(idealVol, rof_gpu3D) -pars['rmse'] = rms - -txtstr = printParametersToString(pars) -txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) -print (txtstr) -a=fig.add_subplot(1,2,2) - -# these are matplotlib.patch.Patch properties -props = dict(boxstyle='round', facecolor='wheat', alpha=0.75) -# place a text box in upper left in axes coords -a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(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') - -rms = rmse(idealVol, fgp_gpu3D) -pars['rmse'] = rms - -txtstr = printParametersToString(pars) -txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) -print (txtstr) -a=fig.add_subplot(1,2,2) - -# these are matplotlib.patch.Patch properties -props = dict(boxstyle='round', facecolor='wheat', alpha=0.75) -# place a text box in upper left in axes coords -a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(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') - -rms = rmse(idealVol, sb_gpu3D) -pars['rmse'] = rms - -txtstr = printParametersToString(pars) -txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) -print (txtstr) -a=fig.add_subplot(1,2,2) - -# these are matplotlib.patch.Patch properties -props = dict(boxstyle='round', facecolor='wheat', alpha=0.75) -# place a text box in upper left in axes coords -a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(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') - -rms = rmse(idealVol, lltrof_gpu3D) -pars['rmse'] = rms - -txtstr = printParametersToString(pars) -txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) -print (txtstr) -a=fig.add_subplot(1,2,2) - -# these are matplotlib.patch.Patch properties -props = dict(boxstyle='round', facecolor='wheat', alpha=0.75) -# place a text box in upper left in axes coords -a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(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') - - -rms = rmse(idealVol, tgv_gpu3D) -pars['rmse'] = rms - -txtstr = printParametersToString(pars) -txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) -print (txtstr) -a=fig.add_subplot(1,2,2) - -# these are matplotlib.patch.Patch properties -props = dict(boxstyle='round', facecolor='wheat', alpha=0.75) -# place a text box in upper left in axes coords -a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(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') - -rms = rmse(idealVol, ndf_gpu3D) -pars['rmse'] = rms - -txtstr = printParametersToString(pars) -txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) -print (txtstr) -a=fig.add_subplot(1,2,2) - -# these are matplotlib.patch.Patch properties -props = dict(boxstyle='round', facecolor='wheat', alpha=0.75) -# place a text box in upper left in axes coords -a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(ndf_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') - -rms = rmse(idealVol, diff4_gpu3D) -pars['rmse'] = rms - -txtstr = printParametersToString(pars) -txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) -print (txtstr) -a=fig.add_subplot(1,2,2) - -# these are matplotlib.patch.Patch properties -props = dict(boxstyle='round', facecolor='wheat', alpha=0.75) -# place a text box in upper left in axes coords -a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(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') - -rms = rmse(idealVol, fgp_dTV_gpu3D) -pars['rmse'] = rms - -txtstr = printParametersToString(pars) -txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) -print (txtstr) -a=fig.add_subplot(1,2,2) - -# these are matplotlib.patch.Patch properties -props = dict(boxstyle='round', facecolor='wheat', alpha=0.75) -# place a text box in upper left in axes coords -a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(fgp_dTV_gpu3D[10,:,:], cmap="gray") -plt.title('{}'.format('Recovered volume on the GPU using FGP-dTV')) -#%% diff --git a/Wrappers/Python/demos/qualitymetrics.py b/Wrappers/Python/demos/qualitymetrics.py deleted file mode 100644 index 850829e..0000000 --- a/Wrappers/Python/demos/qualitymetrics.py +++ /dev/null @@ -1,18 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Wed Feb 21 13:34:32 2018 -# quality metrics -@authors: Daniil Kazantsev, Edoardo Pasca -""" -import numpy as np - -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 diff --git a/Wrappers/Python/setup-regularisers.py.in b/Wrappers/Python/setup-regularisers.py.in deleted file mode 100644 index 462edda..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'}, -) - - -@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); - |