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