From 953a1b6ab10ccaba48a0bd9031885794b9489063 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Thu, 25 Jan 2018 12:08:21 +0000 Subject: compiles the GPU kernels --- Core/CMakeLists.txt | 34 ++++++++++++++++++++++++++++++++-- 1 file changed, 32 insertions(+), 2 deletions(-) diff --git a/Core/CMakeLists.txt b/Core/CMakeLists.txt index 97f8a54..070b45c 100644 --- a/Core/CMakeLists.txt +++ b/Core/CMakeLists.txt @@ -35,7 +35,8 @@ if (Boost_FOUND) endif() ## Build the regularizers package as a library -message("Creating Iterative Reconstruction CGLS as shared library") +message("Creating Regularizers as shared library") + message("CMAKE_CXX_FLAGS ${CMAKE_CXX_FLAGS}") set(CMAKE_BUILD_TYPE "Release") @@ -85,6 +86,9 @@ include_directories(cilreg PUBLIC # EXPORT_FILE_NAME CCPiCore_Exports.h # STATIC_DEFINE cilreg_BUILT_AS_STATIC #) + +## Install + if (UNIX) message ("I'd install into ${CMAKE_INSTALL_PREFIX}/lib") install(TARGETS cilreg @@ -100,7 +104,33 @@ message ("I'd install into ${CMAKE_INSTALL_PREFIX} lib bin") ) endif() +# GPU Regularizers + +find_package(CUDA) +if (CUDA_FOUND) + CUDA_ADD_LIBRARY(cilregcuda SHARED + ${CMAKE_CURRENT_SOURCE_DIR}/regularizers_GPU/Diffus_HO/Diff4th_GPU_kernel.cu + ${CMAKE_CURRENT_SOURCE_DIR}/regularizers_GPU/NL_Regul/NLM_GPU_kernel.cu + ) + if (UNIX) + message ("I'd install into ${CMAKE_INSTALL_PREFIX}/lib") + install(TARGETS cilregcuda + LIBRARY DESTINATION lib + CONFIGURATIONS ${CMAKE_BUILD_TYPE} + ) + elseif(WIN32) + message ("I'd install into ${CMAKE_INSTALL_PREFIX} lib bin") + install(TARGETS cilregcuda + RUNTIME DESTINATION bin + ARCHIVE DESTINATION lib + CONFIGURATIONS ${CMAKE_BUILD_TYPE} + ) + endif() +else() + message("CUDA NOT FOUND") +endif() + #add_executable(regularizer_test ${CMAKE_CURRENT_SOURCE_DIR}/test/test_regularizer.cpp) -#target_link_libraries (regularizer_test LINK_PUBLIC regularizers_lib) \ No newline at end of file +#target_link_libraries (regularizer_test LINK_PUBLIC regularizers_lib) -- cgit v1.2.3 From 6cdaf1cfd380b3c40908d6b82e24afb99e005c71 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Thu, 25 Jan 2018 15:39:24 +0000 Subject: added info --- .../regularizers_GPU/Diffus_HO/Diff4thHajiaboli_GPU.cpp | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/Wrappers/Matlab/mex_compile/regularizers_GPU/Diffus_HO/Diff4thHajiaboli_GPU.cpp b/Wrappers/Matlab/mex_compile/regularizers_GPU/Diffus_HO/Diff4thHajiaboli_GPU.cpp index 5a8c7c0..7fd6077 100644 --- a/Wrappers/Matlab/mex_compile/regularizers_GPU/Diffus_HO/Diff4thHajiaboli_GPU.cpp +++ b/Wrappers/Matlab/mex_compile/regularizers_GPU/Diffus_HO/Diff4thHajiaboli_GPU.cpp @@ -25,6 +25,17 @@ * Linux/Matlab compilation: * compile in terminal: nvcc -Xcompiler -fPIC -shared -o Diff4th_GPU_kernel.o Diff4th_GPU_kernel.cu * then compile in Matlab: mex -I/usr/local/cuda-7.5/include -L/usr/local/cuda-7.5/lib64 -lcudart Diff4thHajiaboli_GPU.cpp Diff4th_GPU_kernel.o + + void mexFunction( +int nlhs, mxArray *plhs[], +int nrhs, const mxArray *prhs[]) +where: +prhs Array of pointers to the INPUT mxArrays +nrhs int number of INPUT mxArrays + +nlhs Array of pointers to the OUTPUT mxArrays +plhs int number of OUTPUT mxArrays + */ void mexFunction( -- cgit v1.2.3 From 886d74aa7bddf2cf5972ab6516ace2dcb764e844 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Thu, 25 Jan 2018 15:41:36 +0000 Subject: added cython wrapper for gpu regularizers --- Wrappers/Python/setup.py | 22 +- Wrappers/Python/src/fista_module.cpp | 1047 ++++++++++++++++++++++++++++++ Wrappers/Python/src/fista_module_gpu.pyx | 154 +++++ Wrappers/Python/src/multiply.pyx | 49 ++ 4 files changed, 1271 insertions(+), 1 deletion(-) create mode 100644 Wrappers/Python/src/fista_module.cpp create mode 100644 Wrappers/Python/src/fista_module_gpu.pyx create mode 100644 Wrappers/Python/src/multiply.pyx diff --git a/Wrappers/Python/setup.py b/Wrappers/Python/setup.py index d2129b0..c535a34 100644 --- a/Wrappers/Python/setup.py +++ b/Wrappers/Python/setup.py @@ -53,13 +53,33 @@ else: else: extra_libraries += ['boost_python', 'boost_numpy','gomp'] +setup( + name='ccpi', + description='CCPi Core Imaging Library - Image Regularizers', + version=cil_version, + cmdclass = {'build_ext': build_ext}, + ext_modules = [Extension("ccpi.filters.gpu_regularizers", + sources=[ + os.path.join("." , "src", "fista_module_gpu.pyx" ), + #os.path.join("." , "src", "multiply.pyx" ) + ], + include_dirs=extra_include_dirs, + library_dirs=extra_library_dirs, + extra_compile_args=extra_compile_args, + libraries=extra_libraries ), + + ], + zip_safe = False, + packages = {'ccpi','ccpi.filters'}, +) + setup( name='ccpi', description='CCPi Core Imaging Library - Image Regularizers', version=cil_version, cmdclass = {'build_ext': build_ext}, ext_modules = [Extension("ccpi.filters.cpu_regularizers", - sources=[os.path.join("." , "fista_module.cpp" ), + sources=[os.path.join("." , "src", "fista_module.cpp" ), # os.path.join("@CMAKE_SOURCE_DIR@" , "main_func" , "regularizers_CPU", "FGP_TV_core.c"), # os.path.join("@CMAKE_SOURCE_DIR@" , "main_func" , "regularizers_CPU", "SplitBregman_TV_core.c"), # os.path.join("@CMAKE_SOURCE_DIR@" , "main_func" , "regularizers_CPU", "LLT_model_core.c"), diff --git a/Wrappers/Python/src/fista_module.cpp b/Wrappers/Python/src/fista_module.cpp new file mode 100644 index 0000000..3876cad --- /dev/null +++ b/Wrappers/Python/src/fista_module.cpp @@ -0,0 +1,1047 @@ +/* +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. +*/ + +#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION + +#include +#include + +#include +#include +#include "boost/tuple/tuple.hpp" + +#include "SplitBregman_TV_core.h" +#include "FGP_TV_core.h" +#include "LLT_model_core.h" +#include "PatchBased_Regul_core.h" +#include "TGV_PD_core.h" +#include "utils.h" + + + +#if defined(_WIN32) || defined(_WIN32) || defined(__WIN32__) || defined(_WIN64) +#include +// this trick only if compiler is MSVC +__if_not_exists(uint8_t) { typedef __int8 uint8_t; } +__if_not_exists(uint16_t) { typedef __int8 uint16_t; } +#endif + +namespace bp = boost::python; +namespace np = boost::python::numpy; + +/*! in the Matlab implementation this is called as +void mexFunction( +int nlhs, mxArray *plhs[], +int nrhs, const mxArray *prhs[]) +where: +prhs Array of pointers to the INPUT mxArrays +nrhs int number of INPUT mxArrays + +nlhs Array of pointers to the OUTPUT mxArrays +plhs int number of OUTPUT mxArrays + +*********************************************************** + +*********************************************************** +double mxGetScalar(const mxArray *pm); +args: pm Pointer to an mxArray; cannot be a cell mxArray, a structure mxArray, or an empty mxArray. +Returns: Pointer to the value of the first real (nonimaginary) element of the mxArray. In C, mxGetScalar returns a double. +*********************************************************** +char *mxArrayToString(const mxArray *array_ptr); +args: array_ptr Pointer to mxCHAR array. +Returns: C-style string. Returns NULL on failure. Possible reasons for failure include out of memory and specifying an array that is not an mxCHAR array. +Description: Call mxArrayToString to copy the character data of an mxCHAR array into a C-style string. +*********************************************************** +mxClassID mxGetClassID(const mxArray *pm); +args: pm Pointer to an mxArray +Returns: Numeric identifier of the class (category) of the mxArray that pm points to.For user-defined types, +mxGetClassId returns a unique value identifying the class of the array contents. +Use mxIsClass to determine whether an array is of a specific user-defined type. + +mxClassID Value MATLAB Type MEX Type C Primitive Type +mxINT8_CLASS int8 int8_T char, byte +mxUINT8_CLASS uint8 uint8_T unsigned char, byte +mxINT16_CLASS int16 int16_T short +mxUINT16_CLASS uint16 uint16_T unsigned short +mxINT32_CLASS int32 int32_T int +mxUINT32_CLASS uint32 uint32_T unsigned int +mxINT64_CLASS int64 int64_T long long +mxUINT64_CLASS uint64 uint64_T unsigned long long +mxSINGLE_CLASS single float float +mxDOUBLE_CLASS double double double + +**************************************************************** +double *mxGetPr(const mxArray *pm); +args: pm Pointer to an mxArray of type double +Returns: Pointer to the first element of the real data. Returns NULL in C (0 in Fortran) if there is no real data. +**************************************************************** +mxArray *mxCreateNumericArray(mwSize ndim, const mwSize *dims, +mxClassID classid, mxComplexity ComplexFlag); +args: ndimNumber of dimensions. If you specify a value for ndim that is less than 2, mxCreateNumericArray automatically sets the number of dimensions to 2. +dims Dimensions array. Each element in the dimensions array contains the size of the array in that dimension. +For example, in C, setting dims[0] to 5 and dims[1] to 7 establishes a 5-by-7 mxArray. Usually there are ndim elements in the dims array. +classid Identifier for the class of the array, which determines the way the numerical data is represented in memory. +For example, specifying mxINT16_CLASS in C causes each piece of numerical data in the mxArray to be represented as a 16-bit signed integer. +ComplexFlag If the mxArray you are creating is to contain imaginary data, set ComplexFlag to mxCOMPLEX in C (1 in Fortran). Otherwise, set ComplexFlag to mxREAL in C (0 in Fortran). +Returns: Pointer to the created mxArray, if successful. If unsuccessful in a standalone (non-MEX file) application, returns NULL in C (0 in Fortran). +If unsuccessful in a MEX file, the MEX file terminates and returns control to the MATLAB prompt. The function is unsuccessful when there is not +enough free heap space to create the mxArray. +*/ + + + +bp::list SplitBregman_TV(np::ndarray input, double d_mu, int iter, double d_epsil, int methTV) { + + // the result is in the following list + bp::list result; + + int number_of_dims, dimX, dimY, dimZ, ll, j, count; + //const int *dim_array; + float *A, *U = NULL, *U_old = NULL, *Dx = NULL, *Dy = NULL, *Dz = NULL, *Bx = NULL, *By = NULL, *Bz = NULL, lambda, mu, epsil, re, re1, re_old; + + //number_of_dims = mxGetNumberOfDimensions(prhs[0]); + //dim_array = mxGetDimensions(prhs[0]); + + number_of_dims = input.get_nd(); + int dim_array[3]; + + dim_array[0] = input.shape(0); + dim_array[1] = input.shape(1); + if (number_of_dims == 2) { + dim_array[2] = -1; + } + else { + dim_array[2] = input.shape(2); + } + + // Parameter handling is be done in Python + ///*Handling Matlab input data*/ + //if ((nrhs < 2) || (nrhs > 5)) mexErrMsgTxt("At least 2 parameters is required: Image(2D/3D), Regularization parameter. The full list of parameters: Image(2D/3D), Regularization parameter, iterations number, tolerance, penalty type ('iso' or 'l1')"); + + ///*Handling Matlab input data*/ + //A = (float *)mxGetData(prhs[0]); /*noisy image (2D/3D) */ + A = reinterpret_cast(input.get_data()); + + //mu = (float)mxGetScalar(prhs[1]); /* regularization parameter */ + mu = (float)d_mu; + + //iter = 35; /* default iterations number */ + + //epsil = 0.0001; /* default tolerance constant */ + epsil = (float)d_epsil; + //methTV = 0; /* default isotropic TV penalty */ + //if ((nrhs == 3) || (nrhs == 4) || (nrhs == 5)) iter = (int)mxGetScalar(prhs[2]); /* iterations number */ + //if ((nrhs == 4) || (nrhs == 5)) epsil = (float)mxGetScalar(prhs[3]); /* tolerance constant */ + //if (nrhs == 5) { + // 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 (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) { mexErrMsgTxt("The input image must be in a single precision"); } + + lambda = 2.0f*mu; + count = 1; + re_old = 0.0f; + /*Handling Matlab output data*/ + dimY = dim_array[0]; dimX = dim_array[1]; dimZ = dim_array[2]; + + if (number_of_dims == 2) { + dimZ = 1; /*2D case*/ + //U = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); + //U_old = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); + //Dx = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); + //Dy = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); + //Bx = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); + //By = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); + bp::tuple shape = bp::make_tuple(dim_array[0], dim_array[1]); + np::dtype dtype = np::dtype::get_builtin(); + + np::ndarray npU = np::zeros(shape, dtype); + np::ndarray npU_old = np::zeros(shape, dtype); + np::ndarray npDx = np::zeros(shape, dtype); + np::ndarray npDy = np::zeros(shape, dtype); + np::ndarray npBx = np::zeros(shape, dtype); + np::ndarray npBy = np::zeros(shape, dtype); + + U = reinterpret_cast(npU.get_data()); + U_old = reinterpret_cast(npU_old.get_data()); + Dx = reinterpret_cast(npDx.get_data()); + Dy = reinterpret_cast(npDy.get_data()); + Bx = reinterpret_cast(npBx.get_data()); + By = reinterpret_cast(npBy.get_data()); + + + + copyIm(A, U, dimX, dimY, dimZ); /*initialize */ + + /* begin outer SB iterations */ + for (ll = 0; ll < iter; ll++) { + + /*storing old values*/ + copyIm(U, U_old, dimX, dimY, dimZ); + + /*GS iteration */ + gauss_seidel2D(U, A, Dx, Dy, Bx, By, dimX, dimY, lambda, mu); + + if (methTV == 1) updDxDy_shrinkAniso2D(U, Dx, Dy, Bx, By, dimX, dimY, lambda); + else updDxDy_shrinkIso2D(U, Dx, Dy, Bx, By, dimX, dimY, lambda); + + updBxBy2D(U, Dx, Dy, Bx, By, dimX, dimY); + + /* calculate norm to terminate earlier */ + re = 0.0f; re1 = 0.0f; + for (j = 0; j < dimX*dimY*dimZ; j++) + { + re += pow(U_old[j] - U[j], 2); + re1 += pow(U_old[j], 2); + } + re = sqrt(re) / sqrt(re1); + if (re < epsil) count++; + if (count > 4) break; + + /* check that the residual norm is decreasing */ + if (ll > 2) { + if (re > re_old) break; + } + re_old = re; + /*printf("%f %i %i \n", re, ll, count); */ + + /*copyIm(U_old, U, dimX, dimY, dimZ); */ + + } + //printf("SB iterations stopped at iteration: %i\n", ll); + result.append(npU); + result.append(ll); + } + if (number_of_dims == 3) { + /*U = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); + U_old = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); + Dx = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); + Dy = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); + Dz = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); + Bx = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); + By = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); + Bz = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));*/ + bp::tuple shape = bp::make_tuple(dim_array[0], dim_array[1], dim_array[2]); + np::dtype dtype = np::dtype::get_builtin(); + + np::ndarray npU = np::zeros(shape, dtype); + np::ndarray npU_old = np::zeros(shape, dtype); + np::ndarray npDx = np::zeros(shape, dtype); + np::ndarray npDy = np::zeros(shape, dtype); + np::ndarray npDz = np::zeros(shape, dtype); + np::ndarray npBx = np::zeros(shape, dtype); + np::ndarray npBy = np::zeros(shape, dtype); + np::ndarray npBz = np::zeros(shape, dtype); + + U = reinterpret_cast(npU.get_data()); + U_old = reinterpret_cast(npU_old.get_data()); + Dx = reinterpret_cast(npDx.get_data()); + Dy = reinterpret_cast(npDy.get_data()); + Dz = reinterpret_cast(npDz.get_data()); + Bx = reinterpret_cast(npBx.get_data()); + By = reinterpret_cast(npBy.get_data()); + Bz = reinterpret_cast(npBz.get_data()); + + copyIm(A, U, dimX, dimY, dimZ); /*initialize */ + + /* begin outer SB iterations */ + for (ll = 0; ll 4) break; + + /* check that the residual norm is decreasing */ + if (ll > 2) { + if (re > re_old) break; + } + /*printf("%f %i %i \n", re, ll, count); */ + re_old = re; + } + //printf("SB iterations stopped at iteration: %i\n", ll); + result.append(npU); + result.append(ll); + } + return result; + + } + + + +bp::list FGP_TV(np::ndarray input, double d_mu, int iter, double d_epsil, int methTV) { + + // the result is in the following list + bp::list result; + + int number_of_dims, dimX, dimY, dimZ, ll, j, count; + float *A, *D = NULL, *D_old = NULL, *P1 = NULL, *P2 = NULL, *P3 = NULL, *P1_old = NULL, *P2_old = NULL, *P3_old = NULL, *R1 = NULL, *R2 = NULL, *R3 = NULL; + float lambda, tk, tkp1, re, re1, re_old, epsil, funcval; + + //number_of_dims = mxGetNumberOfDimensions(prhs[0]); + //dim_array = mxGetDimensions(prhs[0]); + + number_of_dims = input.get_nd(); + int dim_array[3]; + + dim_array[0] = input.shape(0); + dim_array[1] = input.shape(1); + if (number_of_dims == 2) { + dim_array[2] = -1; + } + else { + dim_array[2] = input.shape(2); + } + // Parameter handling is be done in Python + ///*Handling Matlab input data*/ + //if ((nrhs < 2) || (nrhs > 5)) mexErrMsgTxt("At least 2 parameters is required: Image(2D/3D), Regularization parameter. The full list of parameters: Image(2D/3D), Regularization parameter, iterations number, tolerance, penalty type ('iso' or 'l1')"); + + ///*Handling Matlab input data*/ + //A = (float *)mxGetData(prhs[0]); /*noisy image (2D/3D) */ + A = reinterpret_cast(input.get_data()); + + //mu = (float)mxGetScalar(prhs[1]); /* regularization parameter */ + lambda = (float)d_mu; + + //iter = 35; /* default iterations number */ + + //epsil = 0.0001; /* default tolerance constant */ + epsil = (float)d_epsil; + //methTV = 0; /* default isotropic TV penalty */ + //if ((nrhs == 3) || (nrhs == 4) || (nrhs == 5)) iter = (int)mxGetScalar(prhs[2]); /* iterations number */ + //if ((nrhs == 4) || (nrhs == 5)) epsil = (float)mxGetScalar(prhs[3]); /* tolerance constant */ + //if (nrhs == 5) { + // 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 (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) { mexErrMsgTxt("The input image must be in a single precision"); } + + //plhs[1] = mxCreateNumericMatrix(1, 1, mxSINGLE_CLASS, mxREAL); + bp::tuple shape1 = bp::make_tuple(dim_array[0], dim_array[1]); + np::dtype dtype = np::dtype::get_builtin(); + np::ndarray out1 = np::zeros(shape1, dtype); + + //float *funcvalA = (float *)mxGetData(plhs[1]); + float * funcvalA = reinterpret_cast(out1.get_data()); + //if (mxGetClassID(prhs[0]) != 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]; + + tk = 1.0f; + tkp1 = 1.0f; + count = 1; + re_old = 0.0f; + + if (number_of_dims == 2) { + dimZ = 1; /*2D case*/ + /*D = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); + D_old = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); + P1 = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); + P2 = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); + P1_old = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); + P2_old = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); + R1 = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); + R2 = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));*/ + + bp::tuple shape = bp::make_tuple(dim_array[0], dim_array[1]); + np::dtype dtype = np::dtype::get_builtin(); + + + np::ndarray npD = np::zeros(shape, dtype); + np::ndarray npD_old = np::zeros(shape, dtype); + np::ndarray npP1 = np::zeros(shape, dtype); + np::ndarray npP2 = np::zeros(shape, dtype); + np::ndarray npP1_old = np::zeros(shape, dtype); + np::ndarray npP2_old = np::zeros(shape, dtype); + np::ndarray npR1 = np::zeros(shape, dtype); + np::ndarray npR2 = np::zeros(shape, dtype); + + D = reinterpret_cast(npD.get_data()); + D_old = reinterpret_cast(npD_old.get_data()); + P1 = reinterpret_cast(npP1.get_data()); + P2 = reinterpret_cast(npP2.get_data()); + P1_old = reinterpret_cast(npP1_old.get_data()); + P2_old = reinterpret_cast(npP2_old.get_data()); + R1 = reinterpret_cast(npR1.get_data()); + R2 = reinterpret_cast(npR2.get_data()); + + /* begin iterations */ + for (ll = 0; ll 3) { + Obj_func2D(A, D, P1, P2, lambda, dimX, dimY); + funcval = 0.0f; + for (j = 0; j 2) { + if (re > re_old) { + Obj_func2D(A, D, P1, P2, lambda, dimX, dimY); + funcval = 0.0f; + for (j = 0; j(npD); + result.append(out1); + result.append(ll); + } + if (number_of_dims == 3) { + /*D = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); + D_old = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); + P1 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); + P2 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); + P3 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); + P1_old = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); + P2_old = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); + P3_old = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); + R1 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); + R2 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); + R3 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL));*/ + bp::tuple shape = bp::make_tuple(dim_array[0], dim_array[1], dim_array[2]); + np::dtype dtype = np::dtype::get_builtin(); + + np::ndarray npD = np::zeros(shape, dtype); + np::ndarray npD_old = np::zeros(shape, dtype); + np::ndarray npP1 = np::zeros(shape, dtype); + np::ndarray npP2 = np::zeros(shape, dtype); + np::ndarray npP3 = np::zeros(shape, dtype); + np::ndarray npP1_old = np::zeros(shape, dtype); + np::ndarray npP2_old = np::zeros(shape, dtype); + np::ndarray npP3_old = np::zeros(shape, dtype); + np::ndarray npR1 = np::zeros(shape, dtype); + np::ndarray npR2 = np::zeros(shape, dtype); + np::ndarray npR3 = np::zeros(shape, dtype); + + D = reinterpret_cast(npD.get_data()); + D_old = reinterpret_cast(npD_old.get_data()); + P1 = reinterpret_cast(npP1.get_data()); + P2 = reinterpret_cast(npP2.get_data()); + P3 = reinterpret_cast(npP3.get_data()); + P1_old = reinterpret_cast(npP1_old.get_data()); + P2_old = reinterpret_cast(npP2_old.get_data()); + P3_old = reinterpret_cast(npP3_old.get_data()); + R1 = reinterpret_cast(npR1.get_data()); + R2 = reinterpret_cast(npR2.get_data()); + R3 = reinterpret_cast(npR3.get_data()); + /* begin iterations */ + for (ll = 0; ll 3) { + Obj_func3D(A, D, P1, P2, P3, lambda, dimX, dimY, dimZ); + funcval = 0.0f; + for (j = 0; j 2) { + if (re > re_old) { + Obj_func3D(A, D, P1, P2, P3, lambda, dimX, dimY, dimZ); + funcval = 0.0f; + for (j = 0; j(npD); + result.append(out1); + result.append(ll); + } + + return result; +} + +bp::list LLT_model(np::ndarray input, double d_lambda, double d_tau, int iter, double d_epsil, int switcher) { + // the result is in the following list + bp::list result; + + int number_of_dims, dimX, dimY, dimZ, ll, j, count; + //const int *dim_array; + float *U0, *U = NULL, *U_old = NULL, *D1 = NULL, *D2 = NULL, *D3 = NULL, lambda, tau, re, re1, epsil, re_old; + unsigned short *Map = NULL; + + number_of_dims = input.get_nd(); + int dim_array[3]; + + dim_array[0] = input.shape(0); + dim_array[1] = input.shape(1); + if (number_of_dims == 2) { + dim_array[2] = -1; + } + else { + dim_array[2] = input.shape(2); + } + + ///*Handling Matlab input data*/ + //U0 = (float *)mxGetData(prhs[0]); /*origanal noise image/volume*/ + //if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) { mexErrMsgTxt("The input in single precision is required"); } + //lambda = (float)mxGetScalar(prhs[1]); /*regularization parameter*/ + //tau = (float)mxGetScalar(prhs[2]); /* time-step */ + //iter = (int)mxGetScalar(prhs[3]); /*iterations number*/ + //epsil = (float)mxGetScalar(prhs[4]); /* tolerance constant */ + //switcher = (int)mxGetScalar(prhs[5]); /*switch on (1) restrictive smoothing in Z dimension*/ + + U0 = reinterpret_cast(input.get_data()); + lambda = (float)d_lambda; + tau = (float)d_tau; + // iter is passed as parameter + epsil = (float)d_epsil; + // switcher is passed as parameter + /*Handling Matlab output data*/ + dimX = dim_array[0]; dimY = dim_array[1]; dimZ = 1; + + if (number_of_dims == 2) { + /*2D case*/ + /*U = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); + U_old = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); + D1 = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); + D2 = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));*/ + + bp::tuple shape = bp::make_tuple(dim_array[0], dim_array[1]); + np::dtype dtype = np::dtype::get_builtin(); + + + np::ndarray npU = np::zeros(shape, dtype); + np::ndarray npU_old = np::zeros(shape, dtype); + np::ndarray npD1 = np::zeros(shape, dtype); + np::ndarray npD2 = np::zeros(shape, dtype); + + //result.append(npU); + + U = reinterpret_cast(npU.get_data()); + U_old = reinterpret_cast(npU_old.get_data()); + D1 = reinterpret_cast(npD1.get_data()); + D2 = reinterpret_cast(npD2.get_data()); + + /*Copy U0 to U*/ + copyIm(U0, U, dimX, dimY, dimZ); + + count = 1; + re_old = 0.0f; + + for (ll = 0; ll < iter; ll++) { + copyIm(U, U_old, dimX, dimY, dimZ); + + /*estimate inner derrivatives */ + der2D(U, D1, D2, dimX, dimY, dimZ); + /* calculate div^2 and update */ + div_upd2D(U0, U, D1, D2, dimX, dimY, dimZ, lambda, tau); + + /* calculate norm to terminate earlier */ + re = 0.0f; re1 = 0.0f; + for (j = 0; j 4) break; + + /* check that the residual norm is decreasing */ + if (ll > 2) { + if (re > re_old) break; + } + re_old = re; + + } /*end of iterations*/ + // printf("HO iterations stopped at iteration: %i\n", ll); + result.append(npU); + + } + else if (number_of_dims == 3) { + /*3D case*/ + dimZ = dim_array[2]; + /*U = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); + U_old = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); + D1 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); + D2 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); + D3 = (float*)mxGetPr(mxCreateNumericArray(3, dim_array, mxSINGLE_CLASS, mxREAL)); + if (switcher != 0) { + Map = (unsigned short*)mxGetPr(plhs[1] = mxCreateNumericArray(3, dim_array, mxUINT16_CLASS, mxREAL)); + }*/ + bp::tuple shape = bp::make_tuple(dim_array[0], dim_array[1], dim_array[2]); + np::dtype dtype = np::dtype::get_builtin(); + + + np::ndarray npU = np::zeros(shape, dtype); + np::ndarray npU_old = np::zeros(shape, dtype); + np::ndarray npD1 = np::zeros(shape, dtype); + np::ndarray npD2 = np::zeros(shape, dtype); + np::ndarray npD3 = np::zeros(shape, dtype); + np::ndarray npMap = np::zeros(shape, np::dtype::get_builtin()); + Map = reinterpret_cast(npMap.get_data()); + if (switcher != 0) { + //Map = (unsigned short*)mxGetPr(plhs[1] = mxCreateNumericArray(3, dim_array, mxUINT16_CLASS, mxREAL)); + + Map = reinterpret_cast(npMap.get_data()); + } + + U = reinterpret_cast(npU.get_data()); + U_old = reinterpret_cast(npU_old.get_data()); + D1 = reinterpret_cast(npD1.get_data()); + D2 = reinterpret_cast(npD2.get_data()); + D3 = reinterpret_cast(npD2.get_data()); + + /*Copy U0 to U*/ + copyIm(U0, U, dimX, dimY, dimZ); + + count = 1; + re_old = 0.0f; + + + if (switcher == 1) { + /* apply restrictive smoothing */ + calcMap(U, Map, dimX, dimY, dimZ); + /*clear outliers */ + cleanMap(Map, dimX, dimY, dimZ); + } + for (ll = 0; ll < iter; ll++) { + + copyIm(U, U_old, dimX, dimY, dimZ); + + /*estimate inner derrivatives */ + der3D(U, D1, D2, D3, dimX, dimY, dimZ); + /* calculate div^2 and update */ + div_upd3D(U0, U, D1, D2, D3, Map, switcher, dimX, dimY, dimZ, lambda, tau); + + /* calculate norm to terminate earlier */ + re = 0.0f; re1 = 0.0f; + for (j = 0; j 4) break; + + /* check that the residual norm is decreasing */ + if (ll > 2) { + if (re > re_old) break; + } + re_old = re; + + } /*end of iterations*/ + //printf("HO iterations stopped at iteration: %i\n", ll); + result.append(npU); + if (switcher != 0) result.append(npMap); + + } + return result; +} + + +bp::list PatchBased_Regul(np::ndarray input, double d_lambda, int SearchW_real, int SimilW, double d_h) { + // the result is in the following list + bp::list result; + + int N, M, Z, numdims, SearchW, /*SimilW, SearchW_real,*/ padXY, newsizeX, newsizeY, newsizeZ, switchpad_crop; + //const int *dims; + float *A, *B = NULL, *Ap = NULL, *Bp = NULL, h, lambda; + + numdims = input.get_nd(); + int dims[3]; + + dims[0] = input.shape(0); + dims[1] = input.shape(1); + if (numdims == 2) { + dims[2] = -1; + } + else { + dims[2] = input.shape(2); + } + /*numdims = mxGetNumberOfDimensions(prhs[0]); + dims = mxGetDimensions(prhs[0]);*/ + + N = dims[0]; + M = dims[1]; + Z = dims[2]; + + //if ((numdims < 2) || (numdims > 3)) { mexErrMsgTxt("The input should be 2D image or 3D volume"); } + //if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) { mexErrMsgTxt("The input in single precision is required"); } + + //if (nrhs != 5) mexErrMsgTxt("Five inputs reqired: Image(2D,3D), SearchW, SimilW, Threshold, Regularization parameter"); + + ///*Handling inputs*/ + //A = (float *)mxGetData(prhs[0]); /* the image to regularize/filter */ + A = reinterpret_cast(input.get_data()); + //SearchW_real = (int)mxGetScalar(prhs[1]); /* the searching window ratio */ + //SimilW = (int)mxGetScalar(prhs[2]); /* the similarity window ratio */ + //h = (float)mxGetScalar(prhs[3]); /* parameter for the PB filtering function */ + //lambda = (float)mxGetScalar(prhs[4]); /* regularization parameter */ + + //if (h <= 0) mexErrMsgTxt("Parmeter for the PB penalty function should be > 0"); + //if (lambda <= 0) mexErrMsgTxt(" Regularization parmeter should be > 0"); + + lambda = (float)d_lambda; + h = (float)d_h; + SearchW = SearchW_real + 2 * SimilW; + + /* SearchW_full = 2*SearchW + 1; */ /* the full searching window size */ + /* SimilW_full = 2*SimilW + 1; */ /* the full similarity window size */ + + + padXY = SearchW + 2 * SimilW; /* padding sizes */ + newsizeX = N + 2 * (padXY); /* the X size of the padded array */ + newsizeY = M + 2 * (padXY); /* the Y size of the padded array */ + newsizeZ = Z + 2 * (padXY); /* the Z size of the padded array */ + int N_dims[] = { newsizeX, newsizeY, newsizeZ }; + /******************************2D case ****************************/ + if (numdims == 2) { + ///*Handling output*/ + //B = (float*)mxGetData(plhs[0] = mxCreateNumericMatrix(N, M, mxSINGLE_CLASS, mxREAL)); + ///*allocating memory for the padded arrays */ + //Ap = (float*)mxGetData(mxCreateNumericMatrix(newsizeX, newsizeY, mxSINGLE_CLASS, mxREAL)); + //Bp = (float*)mxGetData(mxCreateNumericMatrix(newsizeX, newsizeY, mxSINGLE_CLASS, mxREAL)); + ///**************************************************************************/ + + bp::tuple shape = bp::make_tuple(N, M); + np::dtype dtype = np::dtype::get_builtin(); + + np::ndarray npB = np::zeros(shape, dtype); + + shape = bp::make_tuple(newsizeX, newsizeY); + np::ndarray npAp = np::zeros(shape, dtype); + np::ndarray npBp = np::zeros(shape, dtype); + B = reinterpret_cast(npB.get_data()); + Ap = reinterpret_cast(npAp.get_data()); + Bp = reinterpret_cast(npBp.get_data()); + + /*Perform padding of image A to the size of [newsizeX * newsizeY] */ + switchpad_crop = 0; /*padding*/ + pad_crop(A, Ap, M, N, 0, newsizeY, newsizeX, 0, padXY, switchpad_crop); + + /* Do PB regularization with the padded array */ + PB_FUNC2D(Ap, Bp, newsizeY, newsizeX, padXY, SearchW, SimilW, (float)h, (float)lambda); + + switchpad_crop = 1; /*cropping*/ + pad_crop(Bp, B, M, N, 0, newsizeY, newsizeX, 0, padXY, switchpad_crop); + + result.append(npB); + } + else + { + /******************************3D case ****************************/ + ///*Handling output*/ + //B = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, dims, mxSINGLE_CLASS, mxREAL)); + ///*allocating memory for the padded arrays */ + //Ap = (float*)mxGetPr(mxCreateNumericArray(3, N_dims, mxSINGLE_CLASS, mxREAL)); + //Bp = (float*)mxGetPr(mxCreateNumericArray(3, N_dims, mxSINGLE_CLASS, mxREAL)); + /**************************************************************************/ + bp::tuple shape = bp::make_tuple(dims[0], dims[1], dims[2]); + bp::tuple shape_AB = bp::make_tuple(N_dims[0], N_dims[1], N_dims[2]); + np::dtype dtype = np::dtype::get_builtin(); + + np::ndarray npB = np::zeros(shape, dtype); + np::ndarray npAp = np::zeros(shape_AB, dtype); + np::ndarray npBp = np::zeros(shape_AB, dtype); + B = reinterpret_cast(npB.get_data()); + Ap = reinterpret_cast(npAp.get_data()); + Bp = reinterpret_cast(npBp.get_data()); + /*Perform padding of image A to the size of [newsizeX * newsizeY * newsizeZ] */ + switchpad_crop = 0; /*padding*/ + pad_crop(A, Ap, M, N, Z, newsizeY, newsizeX, newsizeZ, padXY, switchpad_crop); + + /* Do PB regularization with the padded array */ + PB_FUNC3D(Ap, Bp, newsizeY, newsizeX, newsizeZ, padXY, SearchW, SimilW, (float)h, (float)lambda); + + switchpad_crop = 1; /*cropping*/ + pad_crop(Bp, B, M, N, Z, newsizeY, newsizeX, newsizeZ, padXY, switchpad_crop); + + result.append(npB); + } /*end else ndims*/ + + return result; +} + +bp::list TGV_PD(np::ndarray input, double d_lambda, double d_alpha1, double d_alpha0, int iter) { + // the result is in the following list + bp::list result; + int number_of_dims, /*iter,*/ dimX, dimY, dimZ, ll; + //const int *dim_array; + float *A, *U, *U_old, *P1, *P2, *Q1, *Q2, *Q3, *V1, *V1_old, *V2, *V2_old, lambda, L2, tau, sigma, alpha1, alpha0; + + //number_of_dims = mxGetNumberOfDimensions(prhs[0]); + //dim_array = mxGetDimensions(prhs[0]); + number_of_dims = input.get_nd(); + int dim_array[3]; + + dim_array[0] = input.shape(0); + dim_array[1] = input.shape(1); + if (number_of_dims == 2) { + dim_array[2] = -1; + } + else { + dim_array[2] = input.shape(2); + } + /*Handling Matlab input data*/ + //A = (float *)mxGetData(prhs[0]); /*origanal noise image/volume*/ + //if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS) { mexErrMsgTxt("The input in single precision is required"); } + + A = reinterpret_cast(input.get_data()); + + //lambda = (float)mxGetScalar(prhs[1]); /*regularization parameter*/ + //alpha1 = (float)mxGetScalar(prhs[2]); /*first-order term*/ + //alpha0 = (float)mxGetScalar(prhs[3]); /*second-order term*/ + //iter = (int)mxGetScalar(prhs[4]); /*iterations number*/ + //if (nrhs != 5) mexErrMsgTxt("Five input parameters is reqired: Image(2D/3D), Regularization parameter, alpha1, alpha0, Iterations"); + lambda = (float)d_lambda; + alpha1 = (float)d_alpha1; + alpha0 = (float)d_alpha0; + + /*Handling Matlab output data*/ + dimX = dim_array[0]; dimY = dim_array[1]; + + if (number_of_dims == 2) { + /*2D case*/ + dimZ = 1; + bp::tuple shape = bp::make_tuple(dim_array[0], dim_array[1]); + np::dtype dtype = np::dtype::get_builtin(); + + np::ndarray npU = np::zeros(shape, dtype); + np::ndarray npP1 = np::zeros(shape, dtype); + np::ndarray npP2 = np::zeros(shape, dtype); + np::ndarray npQ1 = np::zeros(shape, dtype); + np::ndarray npQ2 = np::zeros(shape, dtype); + np::ndarray npQ3 = np::zeros(shape, dtype); + np::ndarray npV1 = np::zeros(shape, dtype); + np::ndarray npV1_old = np::zeros(shape, dtype); + np::ndarray npV2 = np::zeros(shape, dtype); + np::ndarray npV2_old = np::zeros(shape, dtype); + np::ndarray npU_old = np::zeros(shape, dtype); + + U = reinterpret_cast(npU.get_data()); + U_old = reinterpret_cast(npU_old.get_data()); + P1 = reinterpret_cast(npP1.get_data()); + P2 = reinterpret_cast(npP2.get_data()); + Q1 = reinterpret_cast(npQ1.get_data()); + Q2 = reinterpret_cast(npQ2.get_data()); + Q3 = reinterpret_cast(npQ3.get_data()); + V1 = reinterpret_cast(npV1.get_data()); + V1_old = reinterpret_cast(npV1_old.get_data()); + V2 = reinterpret_cast(npV2.get_data()); + V2_old = reinterpret_cast(npV2_old.get_data()); + //U = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); + + /*dual variables*/ + /*P1 = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); + P2 = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); + + Q1 = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); + Q2 = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); + Q3 = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); + + U_old = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); + + V1 = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); + V1_old = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); + V2 = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL)); + V2_old = (float*)mxGetPr(mxCreateNumericArray(2, dim_array, mxSINGLE_CLASS, mxREAL));*/ + /*printf("%i \n", i);*/ + L2 = 12.0; /*Lipshitz constant*/ + tau = 1.0 / pow(L2, 0.5); + sigma = 1.0 / pow(L2, 0.5); + + /*Copy A to U*/ + copyIm(A, U, dimX, dimY, dimZ); + /* Here primal-dual iterations begin for 2D */ + for (ll = 0; ll < iter; ll++) { + + /* Calculate Dual Variable P */ + DualP_2D(U, V1, V2, P1, P2, dimX, dimY, dimZ, sigma); + + /*Projection onto convex set for P*/ + ProjP_2D(P1, P2, dimX, dimY, dimZ, alpha1); + + /* Calculate Dual Variable Q */ + DualQ_2D(V1, V2, Q1, Q2, Q3, dimX, dimY, dimZ, sigma); + + /*Projection onto convex set for Q*/ + ProjQ_2D(Q1, Q2, Q3, dimX, dimY, dimZ, alpha0); + + /*saving U into U_old*/ + copyIm(U, U_old, dimX, dimY, dimZ); + + /*adjoint operation -> divergence and projection of P*/ + DivProjP_2D(U, A, P1, P2, dimX, dimY, dimZ, lambda, tau); + + /*get updated solution U*/ + newU(U, U_old, dimX, dimY, dimZ); + + /*saving V into V_old*/ + copyIm(V1, V1_old, dimX, dimY, dimZ); + copyIm(V2, V2_old, dimX, dimY, dimZ); + + /* upd V*/ + UpdV_2D(V1, V2, P1, P2, Q1, Q2, Q3, dimX, dimY, dimZ, tau); + + /*get new V*/ + newU(V1, V1_old, dimX, dimY, dimZ); + newU(V2, V2_old, dimX, dimY, dimZ); + } /*end of iterations*/ + + result.append(npU); + } + + + + + return result; +} + +BOOST_PYTHON_MODULE(cpu_regularizers) +{ + np::initialize(); + + //To specify that this module is a package + bp::object package = bp::scope(); + package.attr("__path__") = "cpu_regularizers"; + + np::dtype dt1 = np::dtype::get_builtin(); + np::dtype dt2 = np::dtype::get_builtin(); + + def("SplitBregman_TV", SplitBregman_TV); + def("FGP_TV", FGP_TV); + def("LLT_model", LLT_model); + def("PatchBased_Regul", PatchBased_Regul); + def("TGV_PD", TGV_PD); +} diff --git a/Wrappers/Python/src/fista_module_gpu.pyx b/Wrappers/Python/src/fista_module_gpu.pyx new file mode 100644 index 0000000..9d5b15a --- /dev/null +++ b/Wrappers/Python/src/fista_module_gpu.pyx @@ -0,0 +1,154 @@ +# distutils: language=c++ +""" +Copyright 2018 CCPi +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + http://www.apache.org/licenses/LICENSE-2.0 +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. + +Author: Edoardo Pasca +""" + +import cython + +import numpy as np +cimport numpy as np + + +cdef extern void Diff4th_GPU_kernel(float* A, float* B, int N, int M, int Z, + float sigma, int iter, float tau, float lambdaf); +cdef extern void NLM_GPU_kernel(float *A, float* B, float *Eucl_Vec, + int N, int M, int Z, int dimension, + int SearchW, int SimilW, + int SearchW_real, float denh2, float lambdaf); + +def Diff4thHajiaboli(inputData, + regularization_parameter, + iterations, + edge_preserving_parameter): + if inputData.ndims == 2: + return Diff4thHajiaboli2D(inputData, + regularization_parameter, + iterations, + edge_preserving_parameter) + elif inputData.ndims == 3: + return Diff4thHajiaboli3D(inputData, + regularization_parameter, + iterations, + edge_preserving_parameter) + +def Diff4thHajiaboli2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, + float regularization_parameter, + int iterations, + float edge_preserving_parameter): + + cdef long dims[2] + dims[0] = inputData.shape[0] + dims[1] = inputData.shape[1] + N = dims[0] + 2; + M = dims[1] + 2; + + #time step is sufficiently small for an explicit methods + tau = 0.01 + + #A_L = (float*)mxGetData(mxCreateNumericMatrix(N, M, mxSINGLE_CLASS, mxREAL)); + #B_L = (float*)mxGetData(mxCreateNumericMatrix(N, M, mxSINGLE_CLASS, mxREAL)); + A_L = np.zeros((N,M), dtype=np.float) + B_L = np.zeros((N,M), dtype=np.float) + B = np.zeros((dims[0],dims[1]), dtype=np.float) + #A = inputData + + # copy A to the bigger A_L with boundaries + #pragma omp parallel for shared(A_L, A) private(i,j) + cdef int i, j; + for i in range(N): + for j in range(M): + if (((i > 0) and (i < N-1)) and ((j > 0) and (j < M-1))): + #A_L[i*M+j] = A[(i-1)*(dims[1])+(j-1)] + A_L[i][j] = inputData[i-1][j-1] + + # Running CUDA code here + #Diff4th_GPU_kernel(A_L, B_L, N, M, Z, (float)sigma, iter, (float)tau, lambda); +# Diff4th_GPU_kernel( +# # A_L.data, B_L.data, +# &A_L[0,0], &B_L[0,0], +# N, M, 0, +# edge_preserving_parameter, +# iterations , +# tau, +# regularization_parameter) + # copy the processed B_L to a smaller B + for i in range(N): + for j in range(M): + if (((i > 0) and (i < N-1)) and ((j > 0) and (j < M-1))): + B[i-1][j-1] = B_L[i][j] + ##pragma omp parallel for shared(B_L, B) private(i,j) + #for (i=0; i < N; i++) { + # for (j=0; j < M; j++) { + # if (((i > 0) && (i < N-1)) && ((j > 0) && (j < M-1))) B[(i-1)*(dims[1])+(j-1)] = B_L[i*M+j]; + # }} + + return B + +def Diff4thHajiaboli3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, + float regularization_parameter, + int iterations, + float edge_preserving_parameter): + + cdef long dims[3] + dims[0] = inputData.shape[0] + dims[1] = inputData.shape[1] + dims[2] = inputData.shape[2] + N = dims[0] + 2 + M = dims[1] + 2 + Z = dims[2] + 2 + + # Time Step is small for an explicit methods + tau = 0.0007; + + #A_L = (float*)mxGetData(mxCreateNumericMatrix(N, M, mxSINGLE_CLASS, mxREAL)); + #B_L = (float*)mxGetData(mxCreateNumericMatrix(N, M, mxSINGLE_CLASS, mxREAL)); + A_L = np.zeros((N,M,Z), dtype=np.float) + B_L = np.zeros((N,M,Z), dtype=np.float) + B = np.zeros((dims[0],dims[1],dims[2]), dtype=np.float) + #A = inputData + + # copy A to the bigger A_L with boundaries + #pragma omp parallel for shared(A_L, A) private(i,j) + cdef int i, j, k; + for i in range(N): + for j in range(M): + for k in range(Z): + if (((i > 0) and (i < N-1)) and \ + ((j > 0) and (j < M-1)) and \ + ((k > 0) and (k < Z-1))): + A_L[i][j][k] = inputData[i-1][j-1][k-1]; + + # Running CUDA code here + #Diff4th_GPU_kernel(A_L, B_L, N, M, Z, (float)sigma, iter, (float)tau, lambda); +# Diff4th_GPU_kernel( +# # A_L.data, B_L.data, +# &A_L[0,0,0], &B_L[0,0,0], +# N, M, Z, +# edge_preserving_parameter, +# iterations , +# tau, +# regularization_parameter) + # copy the processed B_L to a smaller B + for i in range(N): + for j in range(M): + for k in range(Z): + if (((i > 0) and (i < N-1)) and \ + ((j > 0) and (j < M-1)) and \ + ((k > 0) and (k < Z-1))): + B[i-1][j-1][k-1] = B_L[i][j][k] + + + return B + + diff --git a/Wrappers/Python/src/multiply.pyx b/Wrappers/Python/src/multiply.pyx new file mode 100644 index 0000000..65df1c6 --- /dev/null +++ b/Wrappers/Python/src/multiply.pyx @@ -0,0 +1,49 @@ +""" +multiply.pyx + +simple cython test of accessing a numpy array's data + +the C function: c_multiply multiplies all the values in a 2-d array by a scalar, in place. + +""" + +import cython + +# import both numpy and the Cython declarations for numpy +import numpy as np +cimport numpy as np + +# declare the interface to the C code +cdef extern void c_multiply (double* array, double value, int m, int n) + +@cython.boundscheck(False) +@cython.wraparound(False) +def multiply(np.ndarray[double, ndim=2, mode="c"] input not None, double value): + """ + multiply (arr, value) + + Takes a numpy arry as input, and multiplies each elemetn by value, in place + + param: array -- a 2-d numpy array of np.float64 + param: value -- a number that will be multiplied by each element in the array + + """ + cdef int m, n + + m, n = input.shape[0], input.shape[1] + + c_multiply (&input[0,0], value, m, n) + + return None + +def multiply2(np.ndarray[double, ndim=2, mode="c"] input not None, double value): + """ + this method works fine, but is not as future-proof the nupy API might change, etc. + """ + cdef int m, n + + m, n = input.shape[0], input.shape[1] + + c_multiply ( input.data, value, m, n) + + return None \ No newline at end of file -- cgit v1.2.3 From 9c341f9187d5f973e3d14425392a306e796524b8 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Thu, 25 Jan 2018 15:52:51 +0000 Subject: ndarray has ndim not ndims --- Wrappers/Python/src/fista_module_gpu.pyx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Wrappers/Python/src/fista_module_gpu.pyx b/Wrappers/Python/src/fista_module_gpu.pyx index 9d5b15a..da86c0a 100644 --- a/Wrappers/Python/src/fista_module_gpu.pyx +++ b/Wrappers/Python/src/fista_module_gpu.pyx @@ -31,12 +31,12 @@ def Diff4thHajiaboli(inputData, regularization_parameter, iterations, edge_preserving_parameter): - if inputData.ndims == 2: + if inputData.ndim == 2: return Diff4thHajiaboli2D(inputData, regularization_parameter, iterations, edge_preserving_parameter) - elif inputData.ndims == 3: + elif inputData.ndim == 3: return Diff4thHajiaboli3D(inputData, regularization_parameter, iterations, -- cgit v1.2.3 From b0af00c358ababaa47afdc581734018b2faf4f0f Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Thu, 25 Jan 2018 21:53:55 +0000 Subject: WIP --- Core/regularizers_GPU/NL_Regul/NLM_GPU_kernel.cu | 26 ++++++++ Core/regularizers_GPU/NL_Regul/NLM_GPU_kernel.h | 2 +- Wrappers/Python/src/fista_module_gpu.pyx | 75 +++++++++++++++++++----- 3 files changed, 86 insertions(+), 17 deletions(-) diff --git a/Core/regularizers_GPU/NL_Regul/NLM_GPU_kernel.cu b/Core/regularizers_GPU/NL_Regul/NLM_GPU_kernel.cu index 1089539..0f18b41 100644 --- a/Core/regularizers_GPU/NL_Regul/NLM_GPU_kernel.cu +++ b/Core/regularizers_GPU/NL_Regul/NLM_GPU_kernel.cu @@ -237,3 +237,29 @@ extern "C" void NLM_GPU_kernel(float *A, float* B, float *Eucl_Vec, int N, int M checkCudaErrors( cudaMemcpy(B,Bd,N*M*Z*sizeof(float),cudaMemcpyDeviceToHost) ); cudaFree(Ad); cudaFree(Bd); cudaFree(Eucl_Vec_d); } + +float pad_crop(float *A, float *Ap, int OldSizeX, int OldSizeY, int OldSizeZ, int NewSizeX, int NewSizeY, int NewSizeZ, int padXY, int switchpad_crop) +{ + /* padding-cropping function */ + int i,j,k; + if (NewSizeZ > 1) { + for (i=0; i < NewSizeX; i++) { + for (j=0; j < NewSizeY; j++) { + for (k=0; k < NewSizeZ; k++) { + if (((i >= padXY) && (i < NewSizeX-padXY)) && ((j >= padXY) && (j < NewSizeY-padXY)) && ((k >= padXY) && (k < NewSizeZ-padXY))) { + if (switchpad_crop == 0) Ap[NewSizeX*NewSizeY*k + i*NewSizeY+j] = A[OldSizeX*OldSizeY*(k - padXY) + (i-padXY)*(OldSizeY)+(j-padXY)]; + else Ap[OldSizeX*OldSizeY*(k - padXY) + (i-padXY)*(OldSizeY)+(j-padXY)] = A[NewSizeX*NewSizeY*k + i*NewSizeY+j]; + } + }}} + } + else { + for (i=0; i < NewSizeX; i++) { + for (j=0; j < NewSizeY; j++) { + if (((i >= padXY) && (i < NewSizeX-padXY)) && ((j >= padXY) && (j < NewSizeY-padXY))) { + if (switchpad_crop == 0) Ap[i*NewSizeY+j] = A[(i-padXY)*(OldSizeY)+(j-padXY)]; + else Ap[(i-padXY)*(OldSizeY)+(j-padXY)] = A[i*NewSizeY+j]; + } + }} + } + return *Ap; +} \ No newline at end of file diff --git a/Core/regularizers_GPU/NL_Regul/NLM_GPU_kernel.h b/Core/regularizers_GPU/NL_Regul/NLM_GPU_kernel.h index f370d0d..3c2bbc5 100644 --- a/Core/regularizers_GPU/NL_Regul/NLM_GPU_kernel.h +++ b/Core/regularizers_GPU/NL_Regul/NLM_GPU_kernel.h @@ -3,5 +3,5 @@ #include "CCPiDefines.h" extern "C" CCPI_EXPORT void NLM_GPU_kernel(float *A, float* B, float *Eucl_Vec, int N, int M, int Z, int dimension, int SearchW, int SimilW, int SearchW_real, float denh2, float lambda); - +extern "C" CCPI_EXPORT float pad_crop(float *A, float *Ap, int OldSizeX, int OldSizeY, int OldSizeZ, int NewSizeX, int NewSizeY, int NewSizeZ, int padXY, int switchpad_crop); #endif diff --git a/Wrappers/Python/src/fista_module_gpu.pyx b/Wrappers/Python/src/fista_module_gpu.pyx index da86c0a..41cf4a6 100644 --- a/Wrappers/Python/src/fista_module_gpu.pyx +++ b/Wrappers/Python/src/fista_module_gpu.pyx @@ -74,14 +74,14 @@ def Diff4thHajiaboli2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, # Running CUDA code here #Diff4th_GPU_kernel(A_L, B_L, N, M, Z, (float)sigma, iter, (float)tau, lambda); -# Diff4th_GPU_kernel( -# # A_L.data, B_L.data, -# &A_L[0,0], &B_L[0,0], -# N, M, 0, -# edge_preserving_parameter, -# iterations , -# tau, -# regularization_parameter) + Diff4th_GPU_kernel( + # A_L.data, B_L.data, + &A_L[0,0], &B_L[0,0], + N, M, 0, + edge_preserving_parameter, + iterations , + tau, + regularization_parameter) # copy the processed B_L to a smaller B for i in range(N): for j in range(M): @@ -131,14 +131,14 @@ def Diff4thHajiaboli3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, # Running CUDA code here #Diff4th_GPU_kernel(A_L, B_L, N, M, Z, (float)sigma, iter, (float)tau, lambda); -# Diff4th_GPU_kernel( -# # A_L.data, B_L.data, -# &A_L[0,0,0], &B_L[0,0,0], -# N, M, Z, -# edge_preserving_parameter, -# iterations , -# tau, -# regularization_parameter) + Diff4th_GPU_kernel( + # A_L.data, B_L.data, + &A_L[0,0,0], &B_L[0,0,0], + N, M, Z, + edge_preserving_parameter, + iterations , + tau, + regularization_parameter) # copy the processed B_L to a smaller B for i in range(N): for j in range(M): @@ -152,3 +152,46 @@ def Diff4thHajiaboli3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, return B +def NML(inputData, + regularization_parameter, + iterations, + edge_preserving_parameter): + if inputData.ndim == 2: + return NML2D(inputData, + regularization_parameter, + iterations, + edge_preserving_parameter) + elif inputData.ndim == 3: + return NML3D(inputData, + regularization_parameter, + iterations, + edge_preserving_parameter) + + #SearchW_real = (int) mxGetScalar(prhs[1]); /* the searching window ratio */ + #SimilW = (int) mxGetScalar(prhs[2]); /* the similarity window ratio */ + #h = (float) mxGetScalar(prhs[3]); /* parameter for the PB filtering function */ + #lambda = (float) mxGetScalar(prhs[4]); + +def NML2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, + SearchW_real, + SimilW, + h, + lambdaf): + N, M = inputData.shape + if h < 0: + raise ValueError('Parameter for the PB filtering function must be > 0') + + SearchW = SearchW_real + 2*SimilW; + + SearchW_full = 2*SearchW + 1; #/* the full searching window size */ + SimilW_full = 2*SimilW + 1; #/* the full similarity window size */ + h2 = h*h; + + padXY = SearchW + 2*SimilW; #/* padding sizes */ + newsizeX = N + 2*(padXY); #/* the X size of the padded array */ + newsizeY = M + 2*(padXY); #/* the Y size of the padded array */ + newsizeZ = Z + 2*(padXY); #/* the Z size of the padded array */ + + B = np.zeros((N,M), dtype=np.float ) + + \ No newline at end of file -- cgit v1.2.3 From 15d24fd2b0cb9ee62ff83afe88b61e5b8ac63cae Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Fri, 26 Jan 2018 10:25:17 +0000 Subject: added NML but commented out --- Wrappers/Python/src/fista_module_gpu.pyx | 192 ++++++++++++++++++++++++------- 1 file changed, 149 insertions(+), 43 deletions(-) diff --git a/Wrappers/Python/src/fista_module_gpu.pyx b/Wrappers/Python/src/fista_module_gpu.pyx index 41cf4a6..f18181b 100644 --- a/Wrappers/Python/src/fista_module_gpu.pyx +++ b/Wrappers/Python/src/fista_module_gpu.pyx @@ -26,6 +26,10 @@ cdef extern void NLM_GPU_kernel(float *A, float* B, float *Eucl_Vec, int N, int M, int Z, int dimension, int SearchW, int SimilW, int SearchW_real, float denh2, float lambdaf); +cdef extern float pad_crop(float *A, float *Ap, + int OldSizeX, int OldSizeY, int OldSizeZ, + int NewSizeX, int NewSizeY, int NewSizeZ, + int padXY, int switchpad_crop); def Diff4thHajiaboli(inputData, regularization_parameter, @@ -152,46 +156,148 @@ def Diff4thHajiaboli3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, return B -def NML(inputData, - regularization_parameter, - iterations, - edge_preserving_parameter): - if inputData.ndim == 2: - return NML2D(inputData, - regularization_parameter, - iterations, - edge_preserving_parameter) - elif inputData.ndim == 3: - return NML3D(inputData, - regularization_parameter, - iterations, - edge_preserving_parameter) - - #SearchW_real = (int) mxGetScalar(prhs[1]); /* the searching window ratio */ - #SimilW = (int) mxGetScalar(prhs[2]); /* the similarity window ratio */ - #h = (float) mxGetScalar(prhs[3]); /* parameter for the PB filtering function */ - #lambda = (float) mxGetScalar(prhs[4]); - -def NML2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, - SearchW_real, - SimilW, - h, - lambdaf): - N, M = inputData.shape - if h < 0: - raise ValueError('Parameter for the PB filtering function must be > 0') - - SearchW = SearchW_real + 2*SimilW; - - SearchW_full = 2*SearchW + 1; #/* the full searching window size */ - SimilW_full = 2*SimilW + 1; #/* the full similarity window size */ - h2 = h*h; - - padXY = SearchW + 2*SimilW; #/* padding sizes */ - newsizeX = N + 2*(padXY); #/* the X size of the padded array */ - newsizeY = M + 2*(padXY); #/* the Y size of the padded array */ - newsizeZ = Z + 2*(padXY); #/* the Z size of the padded array */ - - B = np.zeros((N,M), dtype=np.float ) - - \ No newline at end of file +#def NML(inputData, +# regularization_parameter, +# iterations, +# edge_preserving_parameter): +# if inputData.ndim == 2: +# return NML2D(inputData, +# regularization_parameter, +# iterations, +# edge_preserving_parameter) +# elif inputData.ndim == 3: +# return NML3D(inputData, +# regularization_parameter, +# iterations, +# edge_preserving_parameter) +# +# #SearchW_real = (int) mxGetScalar(prhs[1]); /* the searching window ratio */ +# #SimilW = (int) mxGetScalar(prhs[2]); /* the similarity window ratio */ +# #h = (float) mxGetScalar(prhs[3]); /* parameter for the PB filtering function */ +# #lambda = (float) mxGetScalar(prhs[4]); +# +#def NML2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, +# SearchW_real, +# SimilW, +# h, +# lambdaf): +# N = inputData.shape[0] +# M = inputData.shape[1] +# Z = 0 +# numdims = inputData.ndim +# +# if h < 0: +# raise ValueError('Parameter for the PB filtering function must be > 0') +# +# SearchW = SearchW_real + 2*SimilW; +# +# SearchW_full = 2*SearchW + 1; #/* the full searching window size */ +# SimilW_full = 2*SimilW + 1; #/* the full similarity window size */ +# h2 = h*h; +# +# padXY = SearchW + 2*SimilW; #/* padding sizes */ +# newsizeX = N + 2*(padXY); #/* the X size of the padded array */ +# newsizeY = M + 2*(padXY); #/* the Y size of the padded array */ +# #newsizeZ = Z + 2*(padXY); #/* the Z size of the padded array */ +# +# #output +# B = np.zeros((N,M), dtype=np.float ) +# #/*allocating memory for the padded arrays */ +# +# Ap = np.zeros((newsizeX, newsizeY), dtype=np.float) +# Bp = np.zeros((newsizeX, newsizeY), dtype=np.float) +# Eucl_Vec = np.zeros((SimilW_full*SimilW_full), dtype=np.float) +# +# #/*Gaussian kernel */ +# cdef int count, i_n, j_n; +# cdef float val; +# count = 0 +# for i_n in range (-SimilW, SimilW +1): +# for j_n in range(-SimilW, SimilW +1): +# val = (float)(i_n*i_n + j_n*j_n)/(2*SimilW*SimilW) +# Eucl_Vec[count] = np.exp(-val) +# count = count + 1 +# +# #/*Perform padding of image A to the size of [newsizeX * newsizeY] */ +# switchpad_crop = 0; # /*padding*/ +# pad_crop(&inputData[0,0], &Ap[0,0], M, N, 0, newsizeY, newsizeX, 0, padXY, +# switchpad_crop); +# +# #/* Do PB regularization with the padded array */ +# NLM_GPU_kernel(&Ap[0,0], &Bp[0,0], &Eucl_Vec[0,0], newsizeY, newsizeX, 0, +# numdims, SearchW, SimilW, SearchW_real, +# h2, lambdaf); +# +# switchpad_crop = 1; #/*cropping*/ +# pad_crop(&Bp[0,0], &B[0,0], M, N, 0, newsizeX, newsizeY, 0, padXY, +# switchpad_crop) +# +# return B +# +#def NML3D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, +# SearchW_real, +# SimilW, +# h, +# lambdaf): +# N = inputData.shape[0] +# M = inputData.shape[1] +# Z = inputData.shape[2] +# numdims = inputData.ndim +# +# if h < 0: +# raise ValueError('Parameter for the PB filtering function must be > 0') +# +# SearchW = SearchW_real + 2*SimilW; +# +# SearchW_full = 2*SearchW + 1; #/* the full searching window size */ +# SimilW_full = 2*SimilW + 1; #/* the full similarity window size */ +# h2 = h*h; +# +# padXY = SearchW + 2*SimilW; #/* padding sizes */ +# newsizeX = N + 2*(padXY); #/* the X size of the padded array */ +# newsizeY = M + 2*(padXY); #/* the Y size of the padded array */ +# newsizeZ = Z + 2*(padXY); #/* the Z size of the padded array */ +# +# #output +# B = np.zeros((N,M,Z), dtype=np.float ) +# #/*allocating memory for the padded arrays */ +# +# Ap = np.zeros((newsizeX, newsizeY, newsizeZ), dtype=np.float) +# Bp = np.zeros((newsizeX, newsizeY, newsizeZ), dtype=np.float) +# Eucl_Vec = np.zeros((SimilW_full*SimilW_full*SimilW_full), dtype=np.float) +# +# #/*Gaussian kernel */ +# cdef int count, i_n, j_n, k_n; +# cdef float val; +# count = 0 +# for i_n in range (-SimilW, SimilW +1): +# for j_n in range(-SimilW, SimilW +1): +# for k_n in range(-SimilW, SimilW+1): +# val = (i_n*i_n + j_n*j_n + k_n*k_n)/(2*SimilW*SimilW*SimilW) +# Eucl_Vec[count] = np.exp(-val) +# count = count + 1 +# +# #/*Perform padding of image A to the size of [newsizeX * newsizeY] */ +# switchpad_crop = 0; # /*padding*/ +# pad_crop(&inputData[0,0,0], &Ap[0,0,0], +# M, N, Z, +# newsizeY, newsizeX, newsizeZ, +# padXY, +# switchpad_crop); +# +# #/* Do PB regularization with the padded array */ +# NLM_GPU_kernel(&Ap[0,0,0], &Bp[0,0,0], &Eucl_Vec[0,0,0], +# newsizeY, newsizeX, newsizeZ, +# numdims, SearchW, SimilW, SearchW_real, +# h2, lambdaf); +# +# switchpad_crop = 1; #/*cropping*/ +# pad_crop(&Bp[0,0,0], &B[0,0,0], +# M, N, Z, +# newsizeX, newsizeY, newsizeZ, +# padXY, +# switchpad_crop) +# +# return B +# +# \ No newline at end of file -- cgit v1.2.3 From 9a56fd879fe0f604c1040116463e112eb7e2a8de Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Fri, 26 Jan 2018 16:28:38 +0000 Subject: cython wrapper builds without errors TODO: test --- Wrappers/Python/setup.py | 1 - Wrappers/Python/src/fista_module_gpu.pyx | 322 ++++++++++++++++--------------- 2 files changed, 167 insertions(+), 156 deletions(-) diff --git a/Wrappers/Python/setup.py b/Wrappers/Python/setup.py index c535a34..951146a 100644 --- a/Wrappers/Python/setup.py +++ b/Wrappers/Python/setup.py @@ -61,7 +61,6 @@ setup( ext_modules = [Extension("ccpi.filters.gpu_regularizers", sources=[ os.path.join("." , "src", "fista_module_gpu.pyx" ), - #os.path.join("." , "src", "multiply.pyx" ) ], include_dirs=extra_include_dirs, library_dirs=extra_library_dirs, diff --git a/Wrappers/Python/src/fista_module_gpu.pyx b/Wrappers/Python/src/fista_module_gpu.pyx index f18181b..7658e36 100644 --- a/Wrappers/Python/src/fista_module_gpu.pyx +++ b/Wrappers/Python/src/fista_module_gpu.pyx @@ -15,7 +15,6 @@ Author: Edoardo Pasca """ import cython - import numpy as np cimport numpy as np @@ -62,9 +61,12 @@ def Diff4thHajiaboli2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, #A_L = (float*)mxGetData(mxCreateNumericMatrix(N, M, mxSINGLE_CLASS, mxREAL)); #B_L = (float*)mxGetData(mxCreateNumericMatrix(N, M, mxSINGLE_CLASS, mxREAL)); - A_L = np.zeros((N,M), dtype=np.float) - B_L = np.zeros((N,M), dtype=np.float) - B = np.zeros((dims[0],dims[1]), dtype=np.float) + cdef np.ndarray[np.float32_t, ndim=2, mode="c"] A_L = \ + np.zeros([N,M], dtype='float32') + cdef np.ndarray[np.float32_t, ndim=2, mode="c"] B_L = \ + np.zeros([N,M], dtype='float32') + cdef np.ndarray[np.float32_t, ndim=2, mode="c"] B = \ + np.zeros([dims[0],dims[1]], dtype='float32') #A = inputData # copy A to the bigger A_L with boundaries @@ -85,7 +87,7 @@ def Diff4thHajiaboli2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, edge_preserving_parameter, iterations , tau, - regularization_parameter) + regularization_parameter); # copy the processed B_L to a smaller B for i in range(N): for j in range(M): @@ -117,9 +119,12 @@ def Diff4thHajiaboli3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, #A_L = (float*)mxGetData(mxCreateNumericMatrix(N, M, mxSINGLE_CLASS, mxREAL)); #B_L = (float*)mxGetData(mxCreateNumericMatrix(N, M, mxSINGLE_CLASS, mxREAL)); - A_L = np.zeros((N,M,Z), dtype=np.float) - B_L = np.zeros((N,M,Z), dtype=np.float) - B = np.zeros((dims[0],dims[1],dims[2]), dtype=np.float) + cdef np.ndarray[np.float32_t, ndim=3, mode="c"] A_L = \ + np.zeros([N,M,Z], dtype='float32') + cdef np.ndarray[np.float32_t, ndim=3, mode="c"] B_L = \ + np.zeros([N,M,Z], dtype='float32') + cdef np.ndarray[np.float32_t, ndim=3, mode="c"] B = \ + np.zeros([dims[0],dims[1],dims[2]], dtype='float32') #A = inputData # copy A to the bigger A_L with boundaries @@ -142,7 +147,7 @@ def Diff4thHajiaboli3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, edge_preserving_parameter, iterations , tau, - regularization_parameter) + regularization_parameter); # copy the processed B_L to a smaller B for i in range(N): for j in range(M): @@ -155,149 +160,156 @@ def Diff4thHajiaboli3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, return B - -#def NML(inputData, -# regularization_parameter, -# iterations, -# edge_preserving_parameter): -# if inputData.ndim == 2: -# return NML2D(inputData, -# regularization_parameter, -# iterations, -# edge_preserving_parameter) -# elif inputData.ndim == 3: -# return NML3D(inputData, -# regularization_parameter, -# iterations, -# edge_preserving_parameter) -# -# #SearchW_real = (int) mxGetScalar(prhs[1]); /* the searching window ratio */ -# #SimilW = (int) mxGetScalar(prhs[2]); /* the similarity window ratio */ -# #h = (float) mxGetScalar(prhs[3]); /* parameter for the PB filtering function */ -# #lambda = (float) mxGetScalar(prhs[4]); -# -#def NML2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, -# SearchW_real, -# SimilW, -# h, -# lambdaf): -# N = inputData.shape[0] -# M = inputData.shape[1] -# Z = 0 -# numdims = inputData.ndim -# -# if h < 0: -# raise ValueError('Parameter for the PB filtering function must be > 0') -# -# SearchW = SearchW_real + 2*SimilW; -# -# SearchW_full = 2*SearchW + 1; #/* the full searching window size */ -# SimilW_full = 2*SimilW + 1; #/* the full similarity window size */ -# h2 = h*h; -# -# padXY = SearchW + 2*SimilW; #/* padding sizes */ -# newsizeX = N + 2*(padXY); #/* the X size of the padded array */ -# newsizeY = M + 2*(padXY); #/* the Y size of the padded array */ -# #newsizeZ = Z + 2*(padXY); #/* the Z size of the padded array */ -# -# #output -# B = np.zeros((N,M), dtype=np.float ) -# #/*allocating memory for the padded arrays */ -# -# Ap = np.zeros((newsizeX, newsizeY), dtype=np.float) -# Bp = np.zeros((newsizeX, newsizeY), dtype=np.float) -# Eucl_Vec = np.zeros((SimilW_full*SimilW_full), dtype=np.float) -# -# #/*Gaussian kernel */ -# cdef int count, i_n, j_n; -# cdef float val; -# count = 0 -# for i_n in range (-SimilW, SimilW +1): -# for j_n in range(-SimilW, SimilW +1): -# val = (float)(i_n*i_n + j_n*j_n)/(2*SimilW*SimilW) -# Eucl_Vec[count] = np.exp(-val) -# count = count + 1 -# -# #/*Perform padding of image A to the size of [newsizeX * newsizeY] */ -# switchpad_crop = 0; # /*padding*/ -# pad_crop(&inputData[0,0], &Ap[0,0], M, N, 0, newsizeY, newsizeX, 0, padXY, -# switchpad_crop); -# -# #/* Do PB regularization with the padded array */ -# NLM_GPU_kernel(&Ap[0,0], &Bp[0,0], &Eucl_Vec[0,0], newsizeY, newsizeX, 0, -# numdims, SearchW, SimilW, SearchW_real, -# h2, lambdaf); -# -# switchpad_crop = 1; #/*cropping*/ -# pad_crop(&Bp[0,0], &B[0,0], M, N, 0, newsizeX, newsizeY, 0, padXY, -# switchpad_crop) -# -# return B -# -#def NML3D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, -# SearchW_real, -# SimilW, -# h, -# lambdaf): -# N = inputData.shape[0] -# M = inputData.shape[1] -# Z = inputData.shape[2] -# numdims = inputData.ndim -# -# if h < 0: -# raise ValueError('Parameter for the PB filtering function must be > 0') -# -# SearchW = SearchW_real + 2*SimilW; -# -# SearchW_full = 2*SearchW + 1; #/* the full searching window size */ -# SimilW_full = 2*SimilW + 1; #/* the full similarity window size */ -# h2 = h*h; -# -# padXY = SearchW + 2*SimilW; #/* padding sizes */ -# newsizeX = N + 2*(padXY); #/* the X size of the padded array */ -# newsizeY = M + 2*(padXY); #/* the Y size of the padded array */ -# newsizeZ = Z + 2*(padXY); #/* the Z size of the padded array */ -# -# #output -# B = np.zeros((N,M,Z), dtype=np.float ) -# #/*allocating memory for the padded arrays */ -# -# Ap = np.zeros((newsizeX, newsizeY, newsizeZ), dtype=np.float) -# Bp = np.zeros((newsizeX, newsizeY, newsizeZ), dtype=np.float) -# Eucl_Vec = np.zeros((SimilW_full*SimilW_full*SimilW_full), dtype=np.float) -# -# #/*Gaussian kernel */ -# cdef int count, i_n, j_n, k_n; -# cdef float val; -# count = 0 -# for i_n in range (-SimilW, SimilW +1): -# for j_n in range(-SimilW, SimilW +1): -# for k_n in range(-SimilW, SimilW+1): -# val = (i_n*i_n + j_n*j_n + k_n*k_n)/(2*SimilW*SimilW*SimilW) -# Eucl_Vec[count] = np.exp(-val) -# count = count + 1 -# -# #/*Perform padding of image A to the size of [newsizeX * newsizeY] */ -# switchpad_crop = 0; # /*padding*/ -# pad_crop(&inputData[0,0,0], &Ap[0,0,0], -# M, N, Z, -# newsizeY, newsizeX, newsizeZ, -# padXY, -# switchpad_crop); -# -# #/* Do PB regularization with the padded array */ -# NLM_GPU_kernel(&Ap[0,0,0], &Bp[0,0,0], &Eucl_Vec[0,0,0], -# newsizeY, newsizeX, newsizeZ, -# numdims, SearchW, SimilW, SearchW_real, -# h2, lambdaf); -# -# switchpad_crop = 1; #/*cropping*/ -# pad_crop(&Bp[0,0,0], &B[0,0,0], -# M, N, Z, -# newsizeX, newsizeY, newsizeZ, -# padXY, -# switchpad_crop) -# -# return B -# -# \ No newline at end of file +def NML(inputData, + regularization_parameter, + iterations, + edge_preserving_parameter): + if inputData.ndim == 2: + return NML2D(inputData, + regularization_parameter, + iterations, + edge_preserving_parameter) + elif inputData.ndim == 3: + return NML3D(inputData, + regularization_parameter, + iterations, + edge_preserving_parameter) + + #SearchW_real = (int) mxGetScalar(prhs[1]); /* the searching window ratio */ + #SimilW = (int) mxGetScalar(prhs[2]); /* the similarity window ratio */ + #h = (float) mxGetScalar(prhs[3]); /* parameter for the PB filtering function */ + #lambda = (float) mxGetScalar(prhs[4]); + +def NML2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, + SearchW_real, + SimilW, + h, + lambdaf): + N = inputData.shape[0] + M = inputData.shape[1] + Z = 0 + numdims = inputData.ndim + + if h < 0: + raise ValueError('Parameter for the PB filtering function must be > 0') + + SearchW = SearchW_real + 2*SimilW; + + SearchW_full = 2*SearchW + 1; #/* the full searching window size */ + SimilW_full = 2*SimilW + 1; #/* the full similarity window size */ + h2 = h*h; + + padXY = SearchW + 2*SimilW; #/* padding sizes */ + newsizeX = N + 2*(padXY); #/* the X size of the padded array */ + newsizeY = M + 2*(padXY); #/* the Y size of the padded array */ + #newsizeZ = Z + 2*(padXY); #/* the Z size of the padded array */ + + #output + cdef np.ndarray[np.float32_t, ndim=2, mode="c"] B = \ + np.zeros([N,M], dtype='float32') + #/*allocating memory for the padded arrays */ + + cdef np.ndarray[np.float32_t, ndim=2, mode="c"] Ap = \ + np.zeros([newsizeX, newsizeY], dtype='float32') + cdef np.ndarray[np.float32_t, ndim=2, mode="c"] Bp = \ + np.zeros([newsizeX, newsizeY], dtype='float32') + cdef np.ndarray[np.float32_t, ndim=1, mode="c"] Eucl_Vec = \ + np.zeros([SimilW_full*SimilW_full], dtype='float32') + + #/*Gaussian kernel */ + cdef int count, i_n, j_n; + cdef float val; + count = 0 + for i_n in range (-SimilW, SimilW +1): + for j_n in range(-SimilW, SimilW +1): + val = (float)(i_n*i_n + j_n*j_n)/(2*SimilW*SimilW) + Eucl_Vec[count] = np.exp(-val) + count = count + 1 + + #/*Perform padding of image A to the size of [newsizeX * newsizeY] */ + switchpad_crop = 0; # /*padding*/ + pad_crop(&inputData[0,0], &Ap[0,0], M, N, 0, newsizeY, newsizeX, 0, padXY, + switchpad_crop); + + #/* Do PB regularization with the padded array */ + NLM_GPU_kernel(&Ap[0,0], &Bp[0,0], &Eucl_Vec[0], newsizeY, newsizeX, 0, + numdims, SearchW, SimilW, SearchW_real, + h2, lambdaf); + + switchpad_crop = 1; #/*cropping*/ + pad_crop(&Bp[0,0], &B[0,0], M, N, 0, newsizeX, newsizeY, 0, padXY, + switchpad_crop) + + return B + +def NML3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, + SearchW_real, + SimilW, + h, + lambdaf): + N = inputData.shape[0] + M = inputData.shape[1] + Z = inputData.shape[2] + numdims = inputData.ndim + + if h < 0: + raise ValueError('Parameter for the PB filtering function must be > 0') + + SearchW = SearchW_real + 2*SimilW; + + SearchW_full = 2*SearchW + 1; #/* the full searching window size */ + SimilW_full = 2*SimilW + 1; #/* the full similarity window size */ + h2 = h*h; + + padXY = SearchW + 2*SimilW; #/* padding sizes */ + newsizeX = N + 2*(padXY); #/* the X size of the padded array */ + newsizeY = M + 2*(padXY); #/* the Y size of the padded array */ + newsizeZ = Z + 2*(padXY); #/* the Z size of the padded array */ + + #output + cdef np.ndarray[np.float32_t, ndim=3, mode="c"] B = \ + np.zeros([N,M,Z], dtype='float32') + #/*allocating memory for the padded arrays */ + + cdef np.ndarray[np.float32_t, ndim=3, mode="c"] Ap = \ + np.zeros([newsizeX, newsizeY, newsizeZ], dtype='float32') + cdef np.ndarray[np.float32_t, ndim=3, mode="c"] Bp = \ + np.zeros([newsizeX, newsizeY, newsizeZ], dtype='float32') + cdef np.ndarray[np.float32_t, ndim=1, mode="c"] Eucl_Vec = \ + np.zeros([SimilW_full*SimilW_full*SimilW_full], + dtype='float32') + + + #/*Gaussian kernel */ + cdef int count, i_n, j_n, k_n; + cdef float val; + count = 0 + for i_n in range (-SimilW, SimilW +1): + for j_n in range(-SimilW, SimilW +1): + for k_n in range(-SimilW, SimilW+1): + val = (i_n*i_n + j_n*j_n + k_n*k_n)/(2*SimilW*SimilW*SimilW) + Eucl_Vec[count] = np.exp(-val) + count = count + 1 + + #/*Perform padding of image A to the size of [newsizeX * newsizeY] */ + switchpad_crop = 0; # /*padding*/ + pad_crop(&inputData[0,0,0], &Ap[0,0,0], + M, N, Z, + newsizeY, newsizeX, newsizeZ, + padXY, + switchpad_crop); + + #/* Do PB regularization with the padded array */ + NLM_GPU_kernel(&Ap[0,0,0], &Bp[0,0,0], &Eucl_Vec[0], + newsizeY, newsizeX, newsizeZ, + numdims, SearchW, SimilW, SearchW_real, + h2, lambdaf); + + switchpad_crop = 1; #/*cropping*/ + pad_crop(&Bp[0,0,0], &B[0,0,0], + M, N, Z, + newsizeX, newsizeY, newsizeZ, + padXY, + switchpad_crop) + + return B -- cgit v1.2.3 From 95106cb6620e04c70a3d41d2db84322f9a58afb3 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Mon, 29 Jan 2018 13:36:11 +0000 Subject: renamed fista_module_gpu to gpu_regularizers.pyx Cleaned test_cpu_regularizers.py --- Wrappers/Python/setup.py | 3 +- Wrappers/Python/src/cpu_regularizers.pyx | 0 Wrappers/Python/src/fista_module_gpu.pyx | 315 ------------------ Wrappers/Python/src/gpu_regularizers.pyx | 315 ++++++++++++++++++ Wrappers/Python/test/test_cpu_regularizers.py | 445 ++++++++++++++++++++++++++ 5 files changed, 762 insertions(+), 316 deletions(-) create mode 100644 Wrappers/Python/src/cpu_regularizers.pyx delete mode 100644 Wrappers/Python/src/fista_module_gpu.pyx create mode 100644 Wrappers/Python/src/gpu_regularizers.pyx create mode 100644 Wrappers/Python/test/test_cpu_regularizers.py diff --git a/Wrappers/Python/setup.py b/Wrappers/Python/setup.py index 951146a..00c93fc 100644 --- a/Wrappers/Python/setup.py +++ b/Wrappers/Python/setup.py @@ -60,7 +60,7 @@ setup( cmdclass = {'build_ext': build_ext}, ext_modules = [Extension("ccpi.filters.gpu_regularizers", sources=[ - os.path.join("." , "src", "fista_module_gpu.pyx" ), + os.path.join("." , "src", "gpu_regularizers.pyx" ), ], include_dirs=extra_include_dirs, library_dirs=extra_library_dirs, @@ -79,6 +79,7 @@ setup( cmdclass = {'build_ext': build_ext}, ext_modules = [Extension("ccpi.filters.cpu_regularizers", sources=[os.path.join("." , "src", "fista_module.cpp" ), + os.path.join("." , "src", "cpu_regularizers.pyx" ) # os.path.join("@CMAKE_SOURCE_DIR@" , "main_func" , "regularizers_CPU", "FGP_TV_core.c"), # os.path.join("@CMAKE_SOURCE_DIR@" , "main_func" , "regularizers_CPU", "SplitBregman_TV_core.c"), # os.path.join("@CMAKE_SOURCE_DIR@" , "main_func" , "regularizers_CPU", "LLT_model_core.c"), diff --git a/Wrappers/Python/src/cpu_regularizers.pyx b/Wrappers/Python/src/cpu_regularizers.pyx new file mode 100644 index 0000000..e69de29 diff --git a/Wrappers/Python/src/fista_module_gpu.pyx b/Wrappers/Python/src/fista_module_gpu.pyx deleted file mode 100644 index 7658e36..0000000 --- a/Wrappers/Python/src/fista_module_gpu.pyx +++ /dev/null @@ -1,315 +0,0 @@ -# distutils: language=c++ -""" -Copyright 2018 CCPi -Licensed under the Apache License, Version 2.0 (the "License"); -you may not use this file except in compliance with the License. -You may obtain a copy of the License at - http://www.apache.org/licenses/LICENSE-2.0 -Unless required by applicable law or agreed to in writing, software -distributed under the License is distributed on an "AS IS" BASIS, -WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -See the License for the specific language governing permissions and -limitations under the License. - -Author: Edoardo Pasca -""" - -import cython -import numpy as np -cimport numpy as np - - -cdef extern void Diff4th_GPU_kernel(float* A, float* B, int N, int M, int Z, - float sigma, int iter, float tau, float lambdaf); -cdef extern void NLM_GPU_kernel(float *A, float* B, float *Eucl_Vec, - int N, int M, int Z, int dimension, - int SearchW, int SimilW, - int SearchW_real, float denh2, float lambdaf); -cdef extern float pad_crop(float *A, float *Ap, - int OldSizeX, int OldSizeY, int OldSizeZ, - int NewSizeX, int NewSizeY, int NewSizeZ, - int padXY, int switchpad_crop); - -def Diff4thHajiaboli(inputData, - regularization_parameter, - iterations, - edge_preserving_parameter): - if inputData.ndim == 2: - return Diff4thHajiaboli2D(inputData, - regularization_parameter, - iterations, - edge_preserving_parameter) - elif inputData.ndim == 3: - return Diff4thHajiaboli3D(inputData, - regularization_parameter, - iterations, - edge_preserving_parameter) - -def Diff4thHajiaboli2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, - float regularization_parameter, - int iterations, - float edge_preserving_parameter): - - cdef long dims[2] - dims[0] = inputData.shape[0] - dims[1] = inputData.shape[1] - N = dims[0] + 2; - M = dims[1] + 2; - - #time step is sufficiently small for an explicit methods - tau = 0.01 - - #A_L = (float*)mxGetData(mxCreateNumericMatrix(N, M, mxSINGLE_CLASS, mxREAL)); - #B_L = (float*)mxGetData(mxCreateNumericMatrix(N, M, mxSINGLE_CLASS, mxREAL)); - cdef np.ndarray[np.float32_t, ndim=2, mode="c"] A_L = \ - np.zeros([N,M], dtype='float32') - cdef np.ndarray[np.float32_t, ndim=2, mode="c"] B_L = \ - np.zeros([N,M], dtype='float32') - cdef np.ndarray[np.float32_t, ndim=2, mode="c"] B = \ - np.zeros([dims[0],dims[1]], dtype='float32') - #A = inputData - - # copy A to the bigger A_L with boundaries - #pragma omp parallel for shared(A_L, A) private(i,j) - cdef int i, j; - for i in range(N): - for j in range(M): - if (((i > 0) and (i < N-1)) and ((j > 0) and (j < M-1))): - #A_L[i*M+j] = A[(i-1)*(dims[1])+(j-1)] - A_L[i][j] = inputData[i-1][j-1] - - # Running CUDA code here - #Diff4th_GPU_kernel(A_L, B_L, N, M, Z, (float)sigma, iter, (float)tau, lambda); - Diff4th_GPU_kernel( - # A_L.data, B_L.data, - &A_L[0,0], &B_L[0,0], - N, M, 0, - edge_preserving_parameter, - iterations , - tau, - regularization_parameter); - # copy the processed B_L to a smaller B - for i in range(N): - for j in range(M): - if (((i > 0) and (i < N-1)) and ((j > 0) and (j < M-1))): - B[i-1][j-1] = B_L[i][j] - ##pragma omp parallel for shared(B_L, B) private(i,j) - #for (i=0; i < N; i++) { - # for (j=0; j < M; j++) { - # if (((i > 0) && (i < N-1)) && ((j > 0) && (j < M-1))) B[(i-1)*(dims[1])+(j-1)] = B_L[i*M+j]; - # }} - - return B - -def Diff4thHajiaboli3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, - float regularization_parameter, - int iterations, - float edge_preserving_parameter): - - cdef long dims[3] - dims[0] = inputData.shape[0] - dims[1] = inputData.shape[1] - dims[2] = inputData.shape[2] - N = dims[0] + 2 - M = dims[1] + 2 - Z = dims[2] + 2 - - # Time Step is small for an explicit methods - tau = 0.0007; - - #A_L = (float*)mxGetData(mxCreateNumericMatrix(N, M, mxSINGLE_CLASS, mxREAL)); - #B_L = (float*)mxGetData(mxCreateNumericMatrix(N, M, mxSINGLE_CLASS, mxREAL)); - cdef np.ndarray[np.float32_t, ndim=3, mode="c"] A_L = \ - np.zeros([N,M,Z], dtype='float32') - cdef np.ndarray[np.float32_t, ndim=3, mode="c"] B_L = \ - np.zeros([N,M,Z], dtype='float32') - cdef np.ndarray[np.float32_t, ndim=3, mode="c"] B = \ - np.zeros([dims[0],dims[1],dims[2]], dtype='float32') - #A = inputData - - # copy A to the bigger A_L with boundaries - #pragma omp parallel for shared(A_L, A) private(i,j) - cdef int i, j, k; - for i in range(N): - for j in range(M): - for k in range(Z): - if (((i > 0) and (i < N-1)) and \ - ((j > 0) and (j < M-1)) and \ - ((k > 0) and (k < Z-1))): - A_L[i][j][k] = inputData[i-1][j-1][k-1]; - - # Running CUDA code here - #Diff4th_GPU_kernel(A_L, B_L, N, M, Z, (float)sigma, iter, (float)tau, lambda); - Diff4th_GPU_kernel( - # A_L.data, B_L.data, - &A_L[0,0,0], &B_L[0,0,0], - N, M, Z, - edge_preserving_parameter, - iterations , - tau, - regularization_parameter); - # copy the processed B_L to a smaller B - for i in range(N): - for j in range(M): - for k in range(Z): - if (((i > 0) and (i < N-1)) and \ - ((j > 0) and (j < M-1)) and \ - ((k > 0) and (k < Z-1))): - B[i-1][j-1][k-1] = B_L[i][j][k] - - - return B - -def NML(inputData, - regularization_parameter, - iterations, - edge_preserving_parameter): - if inputData.ndim == 2: - return NML2D(inputData, - regularization_parameter, - iterations, - edge_preserving_parameter) - elif inputData.ndim == 3: - return NML3D(inputData, - regularization_parameter, - iterations, - edge_preserving_parameter) - - #SearchW_real = (int) mxGetScalar(prhs[1]); /* the searching window ratio */ - #SimilW = (int) mxGetScalar(prhs[2]); /* the similarity window ratio */ - #h = (float) mxGetScalar(prhs[3]); /* parameter for the PB filtering function */ - #lambda = (float) mxGetScalar(prhs[4]); - -def NML2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, - SearchW_real, - SimilW, - h, - lambdaf): - N = inputData.shape[0] - M = inputData.shape[1] - Z = 0 - numdims = inputData.ndim - - if h < 0: - raise ValueError('Parameter for the PB filtering function must be > 0') - - SearchW = SearchW_real + 2*SimilW; - - SearchW_full = 2*SearchW + 1; #/* the full searching window size */ - SimilW_full = 2*SimilW + 1; #/* the full similarity window size */ - h2 = h*h; - - padXY = SearchW + 2*SimilW; #/* padding sizes */ - newsizeX = N + 2*(padXY); #/* the X size of the padded array */ - newsizeY = M + 2*(padXY); #/* the Y size of the padded array */ - #newsizeZ = Z + 2*(padXY); #/* the Z size of the padded array */ - - #output - cdef np.ndarray[np.float32_t, ndim=2, mode="c"] B = \ - np.zeros([N,M], dtype='float32') - #/*allocating memory for the padded arrays */ - - cdef np.ndarray[np.float32_t, ndim=2, mode="c"] Ap = \ - np.zeros([newsizeX, newsizeY], dtype='float32') - cdef np.ndarray[np.float32_t, ndim=2, mode="c"] Bp = \ - np.zeros([newsizeX, newsizeY], dtype='float32') - cdef np.ndarray[np.float32_t, ndim=1, mode="c"] Eucl_Vec = \ - np.zeros([SimilW_full*SimilW_full], dtype='float32') - - #/*Gaussian kernel */ - cdef int count, i_n, j_n; - cdef float val; - count = 0 - for i_n in range (-SimilW, SimilW +1): - for j_n in range(-SimilW, SimilW +1): - val = (float)(i_n*i_n + j_n*j_n)/(2*SimilW*SimilW) - Eucl_Vec[count] = np.exp(-val) - count = count + 1 - - #/*Perform padding of image A to the size of [newsizeX * newsizeY] */ - switchpad_crop = 0; # /*padding*/ - pad_crop(&inputData[0,0], &Ap[0,0], M, N, 0, newsizeY, newsizeX, 0, padXY, - switchpad_crop); - - #/* Do PB regularization with the padded array */ - NLM_GPU_kernel(&Ap[0,0], &Bp[0,0], &Eucl_Vec[0], newsizeY, newsizeX, 0, - numdims, SearchW, SimilW, SearchW_real, - h2, lambdaf); - - switchpad_crop = 1; #/*cropping*/ - pad_crop(&Bp[0,0], &B[0,0], M, N, 0, newsizeX, newsizeY, 0, padXY, - switchpad_crop) - - return B - -def NML3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, - SearchW_real, - SimilW, - h, - lambdaf): - N = inputData.shape[0] - M = inputData.shape[1] - Z = inputData.shape[2] - numdims = inputData.ndim - - if h < 0: - raise ValueError('Parameter for the PB filtering function must be > 0') - - SearchW = SearchW_real + 2*SimilW; - - SearchW_full = 2*SearchW + 1; #/* the full searching window size */ - SimilW_full = 2*SimilW + 1; #/* the full similarity window size */ - h2 = h*h; - - padXY = SearchW + 2*SimilW; #/* padding sizes */ - newsizeX = N + 2*(padXY); #/* the X size of the padded array */ - newsizeY = M + 2*(padXY); #/* the Y size of the padded array */ - newsizeZ = Z + 2*(padXY); #/* the Z size of the padded array */ - - #output - cdef np.ndarray[np.float32_t, ndim=3, mode="c"] B = \ - np.zeros([N,M,Z], dtype='float32') - #/*allocating memory for the padded arrays */ - - cdef np.ndarray[np.float32_t, ndim=3, mode="c"] Ap = \ - np.zeros([newsizeX, newsizeY, newsizeZ], dtype='float32') - cdef np.ndarray[np.float32_t, ndim=3, mode="c"] Bp = \ - np.zeros([newsizeX, newsizeY, newsizeZ], dtype='float32') - cdef np.ndarray[np.float32_t, ndim=1, mode="c"] Eucl_Vec = \ - np.zeros([SimilW_full*SimilW_full*SimilW_full], - dtype='float32') - - - #/*Gaussian kernel */ - cdef int count, i_n, j_n, k_n; - cdef float val; - count = 0 - for i_n in range (-SimilW, SimilW +1): - for j_n in range(-SimilW, SimilW +1): - for k_n in range(-SimilW, SimilW+1): - val = (i_n*i_n + j_n*j_n + k_n*k_n)/(2*SimilW*SimilW*SimilW) - Eucl_Vec[count] = np.exp(-val) - count = count + 1 - - #/*Perform padding of image A to the size of [newsizeX * newsizeY] */ - switchpad_crop = 0; # /*padding*/ - pad_crop(&inputData[0,0,0], &Ap[0,0,0], - M, N, Z, - newsizeY, newsizeX, newsizeZ, - padXY, - switchpad_crop); - - #/* Do PB regularization with the padded array */ - NLM_GPU_kernel(&Ap[0,0,0], &Bp[0,0,0], &Eucl_Vec[0], - newsizeY, newsizeX, newsizeZ, - numdims, SearchW, SimilW, SearchW_real, - h2, lambdaf); - - switchpad_crop = 1; #/*cropping*/ - pad_crop(&Bp[0,0,0], &B[0,0,0], - M, N, Z, - newsizeX, newsizeY, newsizeZ, - padXY, - switchpad_crop) - - return B diff --git a/Wrappers/Python/src/gpu_regularizers.pyx b/Wrappers/Python/src/gpu_regularizers.pyx new file mode 100644 index 0000000..7658e36 --- /dev/null +++ b/Wrappers/Python/src/gpu_regularizers.pyx @@ -0,0 +1,315 @@ +# distutils: language=c++ +""" +Copyright 2018 CCPi +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + http://www.apache.org/licenses/LICENSE-2.0 +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. + +Author: Edoardo Pasca +""" + +import cython +import numpy as np +cimport numpy as np + + +cdef extern void Diff4th_GPU_kernel(float* A, float* B, int N, int M, int Z, + float sigma, int iter, float tau, float lambdaf); +cdef extern void NLM_GPU_kernel(float *A, float* B, float *Eucl_Vec, + int N, int M, int Z, int dimension, + int SearchW, int SimilW, + int SearchW_real, float denh2, float lambdaf); +cdef extern float pad_crop(float *A, float *Ap, + int OldSizeX, int OldSizeY, int OldSizeZ, + int NewSizeX, int NewSizeY, int NewSizeZ, + int padXY, int switchpad_crop); + +def Diff4thHajiaboli(inputData, + regularization_parameter, + iterations, + edge_preserving_parameter): + if inputData.ndim == 2: + return Diff4thHajiaboli2D(inputData, + regularization_parameter, + iterations, + edge_preserving_parameter) + elif inputData.ndim == 3: + return Diff4thHajiaboli3D(inputData, + regularization_parameter, + iterations, + edge_preserving_parameter) + +def Diff4thHajiaboli2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, + float regularization_parameter, + int iterations, + float edge_preserving_parameter): + + cdef long dims[2] + dims[0] = inputData.shape[0] + dims[1] = inputData.shape[1] + N = dims[0] + 2; + M = dims[1] + 2; + + #time step is sufficiently small for an explicit methods + tau = 0.01 + + #A_L = (float*)mxGetData(mxCreateNumericMatrix(N, M, mxSINGLE_CLASS, mxREAL)); + #B_L = (float*)mxGetData(mxCreateNumericMatrix(N, M, mxSINGLE_CLASS, mxREAL)); + cdef np.ndarray[np.float32_t, ndim=2, mode="c"] A_L = \ + np.zeros([N,M], dtype='float32') + cdef np.ndarray[np.float32_t, ndim=2, mode="c"] B_L = \ + np.zeros([N,M], dtype='float32') + cdef np.ndarray[np.float32_t, ndim=2, mode="c"] B = \ + np.zeros([dims[0],dims[1]], dtype='float32') + #A = inputData + + # copy A to the bigger A_L with boundaries + #pragma omp parallel for shared(A_L, A) private(i,j) + cdef int i, j; + for i in range(N): + for j in range(M): + if (((i > 0) and (i < N-1)) and ((j > 0) and (j < M-1))): + #A_L[i*M+j] = A[(i-1)*(dims[1])+(j-1)] + A_L[i][j] = inputData[i-1][j-1] + + # Running CUDA code here + #Diff4th_GPU_kernel(A_L, B_L, N, M, Z, (float)sigma, iter, (float)tau, lambda); + Diff4th_GPU_kernel( + # A_L.data, B_L.data, + &A_L[0,0], &B_L[0,0], + N, M, 0, + edge_preserving_parameter, + iterations , + tau, + regularization_parameter); + # copy the processed B_L to a smaller B + for i in range(N): + for j in range(M): + if (((i > 0) and (i < N-1)) and ((j > 0) and (j < M-1))): + B[i-1][j-1] = B_L[i][j] + ##pragma omp parallel for shared(B_L, B) private(i,j) + #for (i=0; i < N; i++) { + # for (j=0; j < M; j++) { + # if (((i > 0) && (i < N-1)) && ((j > 0) && (j < M-1))) B[(i-1)*(dims[1])+(j-1)] = B_L[i*M+j]; + # }} + + return B + +def Diff4thHajiaboli3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, + float regularization_parameter, + int iterations, + float edge_preserving_parameter): + + cdef long dims[3] + dims[0] = inputData.shape[0] + dims[1] = inputData.shape[1] + dims[2] = inputData.shape[2] + N = dims[0] + 2 + M = dims[1] + 2 + Z = dims[2] + 2 + + # Time Step is small for an explicit methods + tau = 0.0007; + + #A_L = (float*)mxGetData(mxCreateNumericMatrix(N, M, mxSINGLE_CLASS, mxREAL)); + #B_L = (float*)mxGetData(mxCreateNumericMatrix(N, M, mxSINGLE_CLASS, mxREAL)); + cdef np.ndarray[np.float32_t, ndim=3, mode="c"] A_L = \ + np.zeros([N,M,Z], dtype='float32') + cdef np.ndarray[np.float32_t, ndim=3, mode="c"] B_L = \ + np.zeros([N,M,Z], dtype='float32') + cdef np.ndarray[np.float32_t, ndim=3, mode="c"] B = \ + np.zeros([dims[0],dims[1],dims[2]], dtype='float32') + #A = inputData + + # copy A to the bigger A_L with boundaries + #pragma omp parallel for shared(A_L, A) private(i,j) + cdef int i, j, k; + for i in range(N): + for j in range(M): + for k in range(Z): + if (((i > 0) and (i < N-1)) and \ + ((j > 0) and (j < M-1)) and \ + ((k > 0) and (k < Z-1))): + A_L[i][j][k] = inputData[i-1][j-1][k-1]; + + # Running CUDA code here + #Diff4th_GPU_kernel(A_L, B_L, N, M, Z, (float)sigma, iter, (float)tau, lambda); + Diff4th_GPU_kernel( + # A_L.data, B_L.data, + &A_L[0,0,0], &B_L[0,0,0], + N, M, Z, + edge_preserving_parameter, + iterations , + tau, + regularization_parameter); + # copy the processed B_L to a smaller B + for i in range(N): + for j in range(M): + for k in range(Z): + if (((i > 0) and (i < N-1)) and \ + ((j > 0) and (j < M-1)) and \ + ((k > 0) and (k < Z-1))): + B[i-1][j-1][k-1] = B_L[i][j][k] + + + return B + +def NML(inputData, + regularization_parameter, + iterations, + edge_preserving_parameter): + if inputData.ndim == 2: + return NML2D(inputData, + regularization_parameter, + iterations, + edge_preserving_parameter) + elif inputData.ndim == 3: + return NML3D(inputData, + regularization_parameter, + iterations, + edge_preserving_parameter) + + #SearchW_real = (int) mxGetScalar(prhs[1]); /* the searching window ratio */ + #SimilW = (int) mxGetScalar(prhs[2]); /* the similarity window ratio */ + #h = (float) mxGetScalar(prhs[3]); /* parameter for the PB filtering function */ + #lambda = (float) mxGetScalar(prhs[4]); + +def NML2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, + SearchW_real, + SimilW, + h, + lambdaf): + N = inputData.shape[0] + M = inputData.shape[1] + Z = 0 + numdims = inputData.ndim + + if h < 0: + raise ValueError('Parameter for the PB filtering function must be > 0') + + SearchW = SearchW_real + 2*SimilW; + + SearchW_full = 2*SearchW + 1; #/* the full searching window size */ + SimilW_full = 2*SimilW + 1; #/* the full similarity window size */ + h2 = h*h; + + padXY = SearchW + 2*SimilW; #/* padding sizes */ + newsizeX = N + 2*(padXY); #/* the X size of the padded array */ + newsizeY = M + 2*(padXY); #/* the Y size of the padded array */ + #newsizeZ = Z + 2*(padXY); #/* the Z size of the padded array */ + + #output + cdef np.ndarray[np.float32_t, ndim=2, mode="c"] B = \ + np.zeros([N,M], dtype='float32') + #/*allocating memory for the padded arrays */ + + cdef np.ndarray[np.float32_t, ndim=2, mode="c"] Ap = \ + np.zeros([newsizeX, newsizeY], dtype='float32') + cdef np.ndarray[np.float32_t, ndim=2, mode="c"] Bp = \ + np.zeros([newsizeX, newsizeY], dtype='float32') + cdef np.ndarray[np.float32_t, ndim=1, mode="c"] Eucl_Vec = \ + np.zeros([SimilW_full*SimilW_full], dtype='float32') + + #/*Gaussian kernel */ + cdef int count, i_n, j_n; + cdef float val; + count = 0 + for i_n in range (-SimilW, SimilW +1): + for j_n in range(-SimilW, SimilW +1): + val = (float)(i_n*i_n + j_n*j_n)/(2*SimilW*SimilW) + Eucl_Vec[count] = np.exp(-val) + count = count + 1 + + #/*Perform padding of image A to the size of [newsizeX * newsizeY] */ + switchpad_crop = 0; # /*padding*/ + pad_crop(&inputData[0,0], &Ap[0,0], M, N, 0, newsizeY, newsizeX, 0, padXY, + switchpad_crop); + + #/* Do PB regularization with the padded array */ + NLM_GPU_kernel(&Ap[0,0], &Bp[0,0], &Eucl_Vec[0], newsizeY, newsizeX, 0, + numdims, SearchW, SimilW, SearchW_real, + h2, lambdaf); + + switchpad_crop = 1; #/*cropping*/ + pad_crop(&Bp[0,0], &B[0,0], M, N, 0, newsizeX, newsizeY, 0, padXY, + switchpad_crop) + + return B + +def NML3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, + SearchW_real, + SimilW, + h, + lambdaf): + N = inputData.shape[0] + M = inputData.shape[1] + Z = inputData.shape[2] + numdims = inputData.ndim + + if h < 0: + raise ValueError('Parameter for the PB filtering function must be > 0') + + SearchW = SearchW_real + 2*SimilW; + + SearchW_full = 2*SearchW + 1; #/* the full searching window size */ + SimilW_full = 2*SimilW + 1; #/* the full similarity window size */ + h2 = h*h; + + padXY = SearchW + 2*SimilW; #/* padding sizes */ + newsizeX = N + 2*(padXY); #/* the X size of the padded array */ + newsizeY = M + 2*(padXY); #/* the Y size of the padded array */ + newsizeZ = Z + 2*(padXY); #/* the Z size of the padded array */ + + #output + cdef np.ndarray[np.float32_t, ndim=3, mode="c"] B = \ + np.zeros([N,M,Z], dtype='float32') + #/*allocating memory for the padded arrays */ + + cdef np.ndarray[np.float32_t, ndim=3, mode="c"] Ap = \ + np.zeros([newsizeX, newsizeY, newsizeZ], dtype='float32') + cdef np.ndarray[np.float32_t, ndim=3, mode="c"] Bp = \ + np.zeros([newsizeX, newsizeY, newsizeZ], dtype='float32') + cdef np.ndarray[np.float32_t, ndim=1, mode="c"] Eucl_Vec = \ + np.zeros([SimilW_full*SimilW_full*SimilW_full], + dtype='float32') + + + #/*Gaussian kernel */ + cdef int count, i_n, j_n, k_n; + cdef float val; + count = 0 + for i_n in range (-SimilW, SimilW +1): + for j_n in range(-SimilW, SimilW +1): + for k_n in range(-SimilW, SimilW+1): + val = (i_n*i_n + j_n*j_n + k_n*k_n)/(2*SimilW*SimilW*SimilW) + Eucl_Vec[count] = np.exp(-val) + count = count + 1 + + #/*Perform padding of image A to the size of [newsizeX * newsizeY] */ + switchpad_crop = 0; # /*padding*/ + pad_crop(&inputData[0,0,0], &Ap[0,0,0], + M, N, Z, + newsizeY, newsizeX, newsizeZ, + padXY, + switchpad_crop); + + #/* Do PB regularization with the padded array */ + NLM_GPU_kernel(&Ap[0,0,0], &Bp[0,0,0], &Eucl_Vec[0], + newsizeY, newsizeX, newsizeZ, + numdims, SearchW, SimilW, SearchW_real, + h2, lambdaf); + + switchpad_crop = 1; #/*cropping*/ + pad_crop(&Bp[0,0,0], &B[0,0,0], + M, N, Z, + newsizeX, newsizeY, newsizeZ, + padXY, + switchpad_crop) + + return B diff --git a/Wrappers/Python/test/test_cpu_regularizers.py b/Wrappers/Python/test/test_cpu_regularizers.py new file mode 100644 index 0000000..ac595e3 --- /dev/null +++ b/Wrappers/Python/test/test_cpu_regularizers.py @@ -0,0 +1,445 @@ +# -*- coding: utf-8 -*- +""" +Created on Fri Aug 4 11:10:05 2017 + +@author: ofn77899 +""" + +#from ccpi.viewer.CILViewer2D import Converter +#import vtk + +import matplotlib.pyplot as plt +import numpy as np +import os +from enum import Enum +import timeit +#from PIL import Image +#from Regularizer import Regularizer +#from ccpi.filters.Regularizer import Regularizer +from ccpi.filters.cpu_regularizers import SplitBregman_TV , FGP_TV , LLT_model, \ + PatchBased_Regul , TGV_PD + +############################################################################### +#https://stackoverflow.com/questions/13875989/comparing-image-in-url-to-image-in-filesystem-in-python/13884956#13884956 +#NRMSE a normalization of the root of the mean squared error +#NRMSE is simply 1 - [RMSE / (maxval - minval)]. Where maxval is the maximum +# intensity from the two images being compared, and respectively the same for +# minval. RMSE is given by the square root of MSE: +# sqrt[(sum(A - B) ** 2) / |A|], +# where |A| means the number of elements in A. By doing this, the maximum value +# given by RMSE is maxval. + +def nrmse(im1, im2): + a, b = im1.shape + rmse = np.sqrt(np.sum((im2 - im1) ** 2) / float(a * b)) + max_val = max(np.max(im1), np.max(im2)) + min_val = min(np.min(im1), np.min(im2)) + return 1 - (rmse / (max_val - min_val)) +############################################################################### +def printParametersToString(pars): + txt = r'' + for key, value in pars.items(): + if key== 'algorithm' : + txt += "{0} = {1}".format(key, value.__name__) + elif key == 'input': + txt += "{0} = {1}".format(key, np.shape(value)) + else: + txt += "{0} = {1}".format(key, value) + txt += '\n' + return txt +############################################################################### +# +# 2D Regularizers +# +############################################################################### +#Example: +# figure; +# Im = double(imread('lena_gray_256.tif'))/255; % loading image +# u0 = Im + .05*randn(size(Im)); u0(u0 < 0) = 0; +# u = SplitBregman_TV(single(u0), 10, 30, 1e-04); + +# assumes the script is launched from the test directory +filename = os.path.join(".." , ".." , ".." , "data" ,"lena_gray_512.tif") +#filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-FISTA_reconstruction\data\lena_gray_512.tif" +#filename = r"/home/ofn77899/Reconstruction/CCPi-FISTA_Reconstruction/data/lena_gray_512.tif" +#filename = r'/home/algol/Documents/Python/STD_test_images/lena_gray_512.tif' + +Im = plt.imread(filename) +Im = np.asarray(Im, dtype='float32') + +perc = 0.05 +u0 = Im + (perc* np.random.normal(size=np.shape(Im))) +# map the u0 u0->u0>0 +f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1) +u0 = f(u0).astype('float32') + +## plot +fig = plt.figure() + +a=fig.add_subplot(2,3,1) +a.set_title('noise') +imgplot = plt.imshow(u0,cmap="gray") + +reg_output = [] +############################################################################## +# Call regularizer + +####################### SplitBregman_TV ##################################### +# u = SplitBregman_TV(single(u0), 10, 30, 1e-04); + +start_time = timeit.default_timer() +pars = {'algorithm' : SplitBregman_TV , \ + 'input' : u0, + 'regularization_parameter':10. , \ +'number_of_iterations' :35 ,\ +'tolerance_constant':0.0001 , \ +'TV_penalty': 0 +} + +out = SplitBregman_TV (pars['input'], pars['regularization_parameter'], + pars['number_of_iterations'], + pars['tolerance_constant'], + pars['TV_penalty']) +plotme = out[0] +txtstr = printParametersToString(pars) +txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) +print (txtstr) + + +a=fig.add_subplot(2,3,2) + + +# these are matplotlib.patch.Patch properties +props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) +# place a text box in upper left in axes coords +a.text(0.05, 0.95, txtstr, transform=a.transAxes, fontsize=14, + verticalalignment='top', bbox=props) +imgplot = plt.imshow(plotme,\ + #cmap="gray" + ) + +###################### FGP_TV ######################################### +# u = FGP_TV(single(u0), 0.05, 100, 1e-04); +start_time = timeit.default_timer() +pars = {'algorithm' : FGP_TV , \ + 'input' : u0, + 'regularization_parameter':5e-4, \ +'number_of_iterations' :10 ,\ +'tolerance_constant':0.001,\ +'TV_penalty': 0 +} + +out = FGP_TV (pars['input'], pars['regularization_parameter'], + pars['number_of_iterations'], + pars['tolerance_constant'], + pars['TV_penalty']) +plotme = out[0] +txtstr = printParametersToString(pars) +txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) +print (txtstr) + + +a=fig.add_subplot(2,3,3) + +# these are matplotlib.patch.Patch properties +props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) +# place a text box in upper left in axes coords +imgplot = plt.imshow(plotme, \ + #cmap="gray" + ) +# place a text box in upper left in axes coords +a.text(0.05, 0.95, txtstr, transform=a.transAxes, fontsize=14, + verticalalignment='top', bbox=props) + +###################### LLT_model ######################################### +# * u0 = Im + .03*randn(size(Im)); % adding noise +# [Den] = LLT_model(single(u0), 10, 0.1, 1); +#Den = LLT_model(single(u0), 25, 0.0003, 300, 0.0001, 0); +#input, regularization_parameter , time_step, number_of_iterations, +# tolerance_constant, restrictive_Z_smoothing=0 + +start_time = timeit.default_timer() + +pars = {'algorithm': LLT_model , \ + 'input' : u0, + 'regularization_parameter': 25,\ + 'time_step':0.0003, \ +'number_of_iterations' :300,\ +'tolerance_constant':0.001,\ +'restrictive_Z_smoothing': 0 +} +out = LLT_model(pars['input'], + pars['regularization_parameter'], + pars['time_step'] , + pars['number_of_iterations'], + pars['tolerance_constant'], + pars['restrictive_Z_smoothing'] ) + +plotme = out[0] +txtstr = printParametersToString(pars) +txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) +print (txtstr) +a=fig.add_subplot(2,3,4) + +# these are matplotlib.patch.Patch properties +props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) +# place a text box in upper left in axes coords +a.text(0.05, 0.95, txtstr, transform=a.transAxes, fontsize=14, + verticalalignment='top', bbox=props) +imgplot = plt.imshow(plotme,\ + #cmap="gray" + ) + + +# ###################### PatchBased_Regul ######################################### +# # Quick 2D denoising example in Matlab: +# # Im = double(imread('lena_gray_256.tif'))/255; % loading image +# # u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise +# # ImDen = PB_Regul_CPU(single(u0), 3, 1, 0.08, 0.05); +start_time = timeit.default_timer() + +pars = {'algorithm': PatchBased_Regul , \ + 'input' : u0, + 'regularization_parameter': 0.05,\ + 'searching_window_ratio':3, \ +'similarity_window_ratio':1,\ +'PB_filtering_parameter': 0.08 +} +out = PatchBased_Regul( + pars['input'], pars['regularization_parameter'], + pars['searching_window_ratio'] , + pars['similarity_window_ratio'] , + pars['PB_filtering_parameter']) +plotme = out[0] +txtstr = printParametersToString(pars) +txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) +print (txtstr) + +a=fig.add_subplot(2,3,5) + + +# these are matplotlib.patch.Patch properties +props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) +# place a text box in upper left in axes coords +a.text(0.05, 0.95, txtstr, transform=a.transAxes, fontsize=14, + verticalalignment='top', bbox=props) +imgplot = plt.imshow(plotme #,cmap="gray" + ) + + +# ###################### TGV_PD ######################################### +# # Quick 2D denoising example in Matlab: +# # Im = double(imread('lena_gray_256.tif'))/255; % loading image +# # u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise +# # u = PrimalDual_TGV(single(u0), 0.02, 1.3, 1, 550); + +start_time = timeit.default_timer() + +pars = {'algorithm': TGV_PD , \ + 'input' : u0,\ + 'regularization_parameter':0.05,\ + 'first_order_term': 1.3,\ + 'second_order_term': 1, \ +'number_of_iterations': 550 +} +out = TGV_PD(pars['input'], pars['regularization_parameter'], + pars['first_order_term'] , + pars['second_order_term'] , + pars['number_of_iterations']) +plotme = out[0] +txtstr = printParametersToString(pars) +txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) +print (txtstr) +a=fig.add_subplot(2,3,6) + +# these are matplotlib.patch.Patch properties +props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) +# place a text box in upper left in axes coords +a.text(0.05, 0.95, txtstr, transform=a.transAxes, fontsize=14, + verticalalignment='top', bbox=props) +imgplot = plt.imshow(plotme #, cmap="gray") + ) + + +plt.show() + +################################################################################ +## +## 3D Regularizers +## +################################################################################ +##Example: +## figure; +## Im = double(imread('lena_gray_256.tif'))/255; % loading image +## u0 = Im + .05*randn(size(Im)); u0(u0 < 0) = 0; +## u = SplitBregman_TV(single(u0), 10, 30, 1e-04); +# +##filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-Reconstruction\python\test\reconstruction_example.mha" +#filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-Simpleflex\data\head.mha" +# +#reader = vtk.vtkMetaImageReader() +#reader.SetFileName(os.path.normpath(filename)) +#reader.Update() +##vtk returns 3D images, let's take just the one slice there is as 2D +#Im = Converter.vtk2numpy(reader.GetOutput()) +#Im = Im.astype('float32') +##imgplot = plt.imshow(Im) +#perc = 0.05 +#u0 = Im + (perc* np.random.normal(size=np.shape(Im))) +## map the u0 u0->u0>0 +#f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1) +#u0 = f(u0).astype('float32') +#converter = Converter.numpy2vtkImporter(u0, reader.GetOutput().GetSpacing(), +# reader.GetOutput().GetOrigin()) +#converter.Update() +#writer = vtk.vtkMetaImageWriter() +#writer.SetInputData(converter.GetOutput()) +#writer.SetFileName(r"C:\Users\ofn77899\Documents\GitHub\CCPi-FISTA_reconstruction\data\noisy_head.mha") +##writer.Write() +# +# +### plot +#fig3D = plt.figure() +##a=fig.add_subplot(3,3,1) +##a.set_title('Original') +##imgplot = plt.imshow(Im) +#sliceNo = 32 +# +#a=fig3D.add_subplot(2,3,1) +#a.set_title('noise') +#imgplot = plt.imshow(u0.T[sliceNo]) +# +#reg_output3d = [] +# +############################################################################### +## Call regularizer +# +######################## SplitBregman_TV ##################################### +## u = SplitBregman_TV(single(u0), 10, 30, 1e-04); +# +##reg = Regularizer(Regularizer.Algorithm.SplitBregman_TV) +# +##out = reg(input=u0, regularization_parameter=10., #number_of_iterations=30, +## #tolerance_constant=1e-4, +## TV_Penalty=Regularizer.TotalVariationPenalty.l1) +# +#out2 = Regularizer.SplitBregman_TV(input=u0, regularization_parameter=10., number_of_iterations=30, +# tolerance_constant=1e-4, +# TV_Penalty=Regularizer.TotalVariationPenalty.l1) +# +# +#pars = out2[-2] +#reg_output3d.append(out2) +# +#a=fig3D.add_subplot(2,3,2) +# +# +#textstr = out2[-1] +# +# +## these are matplotlib.patch.Patch properties +#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) +## place a text box in upper left in axes coords +#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14, +# verticalalignment='top', bbox=props) +#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo]) +# +####################### FGP_TV ######################################### +## u = FGP_TV(single(u0), 0.05, 100, 1e-04); +#out2 = Regularizer.FGP_TV(input=u0, regularization_parameter=0.005, +# number_of_iterations=200) +#pars = out2[-2] +#reg_output3d.append(out2) +# +#a=fig3D.add_subplot(2,3,2) +# +# +#textstr = out2[-1] +# +# +## these are matplotlib.patch.Patch properties +#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) +## place a text box in upper left in axes coords +#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14, +# verticalalignment='top', bbox=props) +#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo]) +# +####################### LLT_model ######################################### +## * u0 = Im + .03*randn(size(Im)); % adding noise +## [Den] = LLT_model(single(u0), 10, 0.1, 1); +##Den = LLT_model(single(u0), 25, 0.0003, 300, 0.0001, 0); +##input, regularization_parameter , time_step, number_of_iterations, +## tolerance_constant, restrictive_Z_smoothing=0 +#out2 = Regularizer.LLT_model(input=u0, regularization_parameter=25, +# time_step=0.0003, +# tolerance_constant=0.0001, +# number_of_iterations=300) +#pars = out2[-2] +#reg_output3d.append(out2) +# +#a=fig3D.add_subplot(2,3,2) +# +# +#textstr = out2[-1] +# +# +## these are matplotlib.patch.Patch properties +#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) +## place a text box in upper left in axes coords +#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14, +# verticalalignment='top', bbox=props) +#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo]) +# +####################### PatchBased_Regul ######################################### +## Quick 2D denoising example in Matlab: +## Im = double(imread('lena_gray_256.tif'))/255; % loading image +## u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise +## ImDen = PB_Regul_CPU(single(u0), 3, 1, 0.08, 0.05); +# +#out2 = Regularizer.PatchBased_Regul(input=u0, regularization_parameter=0.05, +# searching_window_ratio=3, +# similarity_window_ratio=1, +# PB_filtering_parameter=0.08) +#pars = out2[-2] +#reg_output3d.append(out2) +# +#a=fig3D.add_subplot(2,3,2) +# +# +#textstr = out2[-1] +# +# +## these are matplotlib.patch.Patch properties +#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) +## place a text box in upper left in axes coords +#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14, +# verticalalignment='top', bbox=props) +#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo]) +# + +###################### TGV_PD ######################################### +# Quick 2D denoising example in Matlab: +# Im = double(imread('lena_gray_256.tif'))/255; % loading image +# u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise +# u = PrimalDual_TGV(single(u0), 0.02, 1.3, 1, 550); + + +#out2 = Regularizer.TGV_PD(input=u0, regularization_parameter=0.05, +# first_order_term=1.3, +# second_order_term=1, +# number_of_iterations=550) +#pars = out2[-2] +#reg_output3d.append(out2) +# +#a=fig3D.add_subplot(2,3,2) +# +# +#textstr = out2[-1] +# +# +## these are matplotlib.patch.Patch properties +#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) +## place a text box in upper left in axes coords +#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14, +# verticalalignment='top', bbox=props) +#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo]) -- cgit v1.2.3 From fd7236f721017a6e8c082f01545e6b5aaf4c2cac Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Mon, 29 Jan 2018 14:43:46 +0000 Subject: minor changes to test --- Wrappers/Python/test/test_cpu_regularizers.py | 98 +++++++++++++-------------- 1 file changed, 47 insertions(+), 51 deletions(-) diff --git a/Wrappers/Python/test/test_cpu_regularizers.py b/Wrappers/Python/test/test_cpu_regularizers.py index ac595e3..6c97875 100644 --- a/Wrappers/Python/test/test_cpu_regularizers.py +++ b/Wrappers/Python/test/test_cpu_regularizers.py @@ -5,17 +5,12 @@ Created on Fri Aug 4 11:10:05 2017 @author: ofn77899 """ -#from ccpi.viewer.CILViewer2D import Converter -#import vtk import matplotlib.pyplot as plt import numpy as np import os from enum import Enum import timeit -#from PIL import Image -#from Regularizer import Regularizer -#from ccpi.filters.Regularizer import Regularizer from ccpi.filters.cpu_regularizers import SplitBregman_TV , FGP_TV , LLT_model, \ PatchBased_Regul , TGV_PD @@ -67,8 +62,10 @@ filename = os.path.join(".." , ".." , ".." , "data" ,"lena_gray_512.tif") Im = plt.imread(filename) Im = np.asarray(Im, dtype='float32') -perc = 0.05 -u0 = Im + (perc* np.random.normal(size=np.shape(Im))) +perc = 0.15 +u0 = Im + np.random.normal(loc = Im , + scale = perc * Im , + size = np.shape(Im)) # map the u0 u0->u0>0 f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1) u0 = f(u0).astype('float32') @@ -78,7 +75,8 @@ fig = plt.figure() a=fig.add_subplot(2,3,1) a.set_title('noise') -imgplot = plt.imshow(u0,cmap="gray") +imgplot = plt.imshow(u0#,cmap="gray" + ) reg_output = [] ############################################################################## @@ -100,7 +98,7 @@ out = SplitBregman_TV (pars['input'], pars['regularization_parameter'], pars['number_of_iterations'], pars['tolerance_constant'], pars['TV_penalty']) -plotme = out[0] +splitbregman = out[0] txtstr = printParametersToString(pars) txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) print (txtstr) @@ -114,7 +112,7 @@ props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) # place a text box in upper left in axes coords a.text(0.05, 0.95, txtstr, transform=a.transAxes, fontsize=14, verticalalignment='top', bbox=props) -imgplot = plt.imshow(plotme,\ +imgplot = plt.imshow(splitbregman,\ #cmap="gray" ) @@ -124,16 +122,18 @@ start_time = timeit.default_timer() pars = {'algorithm' : FGP_TV , \ 'input' : u0, 'regularization_parameter':5e-4, \ -'number_of_iterations' :10 ,\ -'tolerance_constant':0.001,\ -'TV_penalty': 0 + 'number_of_iterations' :10 ,\ + 'tolerance_constant':0.001,\ + 'TV_penalty': 0 } -out = FGP_TV (pars['input'], pars['regularization_parameter'], - pars['number_of_iterations'], - pars['tolerance_constant'], - pars['TV_penalty']) -plotme = out[0] +out = FGP_TV (pars['input'], + pars['regularization_parameter'], + pars['number_of_iterations'], + pars['tolerance_constant'], + pars['TV_penalty']) + +fgp = out[0] txtstr = printParametersToString(pars) txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) print (txtstr) @@ -144,7 +144,7 @@ a=fig.add_subplot(2,3,3) # these are matplotlib.patch.Patch properties props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) # place a text box in upper left in axes coords -imgplot = plt.imshow(plotme, \ +imgplot = plt.imshow(fgp, \ #cmap="gray" ) # place a text box in upper left in axes coords @@ -152,11 +152,6 @@ a.text(0.05, 0.95, txtstr, transform=a.transAxes, fontsize=14, verticalalignment='top', bbox=props) ###################### LLT_model ######################################### -# * u0 = Im + .03*randn(size(Im)); % adding noise -# [Den] = LLT_model(single(u0), 10, 0.1, 1); -#Den = LLT_model(single(u0), 25, 0.0003, 300, 0.0001, 0); -#input, regularization_parameter , time_step, number_of_iterations, -# tolerance_constant, restrictive_Z_smoothing=0 start_time = timeit.default_timer() @@ -164,18 +159,18 @@ pars = {'algorithm': LLT_model , \ 'input' : u0, 'regularization_parameter': 25,\ 'time_step':0.0003, \ -'number_of_iterations' :300,\ -'tolerance_constant':0.001,\ -'restrictive_Z_smoothing': 0 + 'number_of_iterations' :300,\ + 'tolerance_constant':0.001,\ + 'restrictive_Z_smoothing': 0 } out = LLT_model(pars['input'], - pars['regularization_parameter'], - pars['time_step'] , - pars['number_of_iterations'], - pars['tolerance_constant'], - pars['restrictive_Z_smoothing'] ) + pars['regularization_parameter'], + pars['time_step'] , + pars['number_of_iterations'], + pars['tolerance_constant'], + pars['restrictive_Z_smoothing'] ) -plotme = out[0] +llt = out[0] txtstr = printParametersToString(pars) txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) print (txtstr) @@ -186,7 +181,7 @@ props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) # place a text box in upper left in axes coords a.text(0.05, 0.95, txtstr, transform=a.transAxes, fontsize=14, verticalalignment='top', bbox=props) -imgplot = plt.imshow(plotme,\ +imgplot = plt.imshow(llt,\ #cmap="gray" ) @@ -202,15 +197,15 @@ pars = {'algorithm': PatchBased_Regul , \ 'input' : u0, 'regularization_parameter': 0.05,\ 'searching_window_ratio':3, \ -'similarity_window_ratio':1,\ -'PB_filtering_parameter': 0.08 + 'similarity_window_ratio':1,\ + 'PB_filtering_parameter': 0.08 } -out = PatchBased_Regul( - pars['input'], pars['regularization_parameter'], - pars['searching_window_ratio'] , - pars['similarity_window_ratio'] , - pars['PB_filtering_parameter']) -plotme = out[0] +out = PatchBased_Regul(pars['input'], + pars['regularization_parameter'], + pars['searching_window_ratio'] , + pars['similarity_window_ratio'] , + pars['PB_filtering_parameter']) +pbr = out[0] txtstr = printParametersToString(pars) txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) print (txtstr) @@ -223,7 +218,7 @@ props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) # place a text box in upper left in axes coords a.text(0.05, 0.95, txtstr, transform=a.transAxes, fontsize=14, verticalalignment='top', bbox=props) -imgplot = plt.imshow(plotme #,cmap="gray" +imgplot = plt.imshow(pbr #,cmap="gray" ) @@ -240,13 +235,14 @@ pars = {'algorithm': TGV_PD , \ 'regularization_parameter':0.05,\ 'first_order_term': 1.3,\ 'second_order_term': 1, \ -'number_of_iterations': 550 -} -out = TGV_PD(pars['input'], pars['regularization_parameter'], - pars['first_order_term'] , - pars['second_order_term'] , - pars['number_of_iterations']) -plotme = out[0] + 'number_of_iterations': 550 + } +out = TGV_PD(pars['input'], + pars['regularization_parameter'], + pars['first_order_term'] , + pars['second_order_term'] , + pars['number_of_iterations']) +tgv = out[0] txtstr = printParametersToString(pars) txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) print (txtstr) @@ -257,7 +253,7 @@ props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) # place a text box in upper left in axes coords a.text(0.05, 0.95, txtstr, transform=a.transAxes, fontsize=14, verticalalignment='top', bbox=props) -imgplot = plt.imshow(plotme #, cmap="gray") +imgplot = plt.imshow(tgv #, cmap="gray") ) -- cgit v1.2.3 From 25c43ebd60c141bacf3ff9fd9c5e2f357bc5e4b6 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Mon, 29 Jan 2018 15:46:04 +0000 Subject: finds boost libraries during wrappers build (win) finds and set the correct name for the boost libraries. Builds GPU wrapper if CUDA is found. --- Wrappers/Python/conda-recipe/bld.bat | 7 +- Wrappers/Python/conda-recipe/build.sh | 6 +- Wrappers/Python/conda-recipe/meta.yaml | 1 + Wrappers/Python/setup-regularizers.py.in | 84 ++++++++++++++++++++++ Wrappers/Python/setup.py | 100 -------------------------- Wrappers/Python/src/cpu_regularizers.pyx | 19 +++++ Wrappers/Python/src/fista_module.cpp | 4 +- Wrappers/Python/test/test_cpu_regularizers.py | 5 +- 8 files changed, 118 insertions(+), 108 deletions(-) create mode 100644 Wrappers/Python/setup-regularizers.py.in delete mode 100644 Wrappers/Python/setup.py diff --git a/Wrappers/Python/conda-recipe/bld.bat b/Wrappers/Python/conda-recipe/bld.bat index fe3ddae..850905c 100644 --- a/Wrappers/Python/conda-recipe/bld.bat +++ b/Wrappers/Python/conda-recipe/bld.bat @@ -8,7 +8,10 @@ ROBOCOPY /E "%RECIPE_DIR%\..\.." "%SRC_DIR%\ccpi" ROBOCOPY /E "%RECIPE_DIR%\..\..\..\Core" "%SRC_DIR%\Core" cd %SRC_DIR%\ccpi\Python -%PYTHON% setup.py build_ext +:: issue cmake to create setup.py +cmake . + +%PYTHON% setup-regularizers.py build_ext if errorlevel 1 exit 1 -%PYTHON% setup.py install +%PYTHON% setup-regularizers.py install if errorlevel 1 exit 1 diff --git a/Wrappers/Python/conda-recipe/build.sh b/Wrappers/Python/conda-recipe/build.sh index aaf9a69..9ea4161 100644 --- a/Wrappers/Python/conda-recipe/build.sh +++ b/Wrappers/Python/conda-recipe/build.sh @@ -11,7 +11,9 @@ cd $SRC_DIR/ccpi/Python echo "$SRC_DIR/ccpi/Python" -$PYTHON setup.py build_ext -$PYTHON setup.py install +cmake . + +$PYTHON setup-regularizers.py build_ext +$PYTHON setup-regularizers.py install diff --git a/Wrappers/Python/conda-recipe/meta.yaml b/Wrappers/Python/conda-recipe/meta.yaml index c451b37..8b58738 100644 --- a/Wrappers/Python/conda-recipe/meta.yaml +++ b/Wrappers/Python/conda-recipe/meta.yaml @@ -21,6 +21,7 @@ requirements: - cil_regularizer - vc 14 # [win and py35] - vc 9 # [win and py27] + - cmake run: - python ==2.7 # [py27] diff --git a/Wrappers/Python/setup-regularizers.py.in b/Wrappers/Python/setup-regularizers.py.in new file mode 100644 index 0000000..3f4d2d7 --- /dev/null +++ b/Wrappers/Python/setup-regularizers.py.in @@ -0,0 +1,84 @@ +#!/usr/bin/env python + +import setuptools +from distutils.core import setup +from distutils.extension import Extension +from Cython.Distutils import build_ext + +import os +import sys +import numpy +import platform + +cil_version=os.environ['CIL_VERSION'] +if cil_version == '': + print("Please set the environmental variable CIL_VERSION") + sys.exit(1) + +library_include_path = "" +library_lib_path = "" +try: + library_include_path = os.environ['LIBRARY_INC'] + library_lib_path = os.environ['LIBRARY_LIB'] +except: + library_include_path = os.environ['PREFIX']+'/include' + pass + +extra_include_dirs = [numpy.get_include(), library_include_path] +#extra_library_dirs = [os.path.join(library_include_path, "..", "lib")] +extra_compile_args = [] +extra_library_dirs = [] +extra_compile_args = [] +extra_link_args = [] +extra_libraries = ['cilreg'] + +extra_include_dirs += [os.path.join(".." , ".." , "Core"), + os.path.join(".." , ".." , "Core", "regularizers_CPU"), + os.path.join(".." , ".." , "Core", "regularizers_GPU") , + "."] + +if platform.system() == 'Windows': + extra_compile_args[0:] = ['/DWIN32','/EHsc','/DBOOST_ALL_NO_LIB' , '/openmp' ] +else: + extra_compile_args = ['-fopenmp','-O2', '-funsigned-char', '-Wall', '-std=c++0x'] + extra_libraries += [@EXTRA_OMP_LIB@] + +extra_libraries += ["@BOOST_PYTHON_LIB@", "@BOOST_NUMPY_LIB@"] + + +setup( + name='ccpi', + description='CCPi Core Imaging Library - Image Regularizers', + version=cil_version, + cmdclass = {'build_ext': build_ext}, + ext_modules = [Extension("ccpi.filters.cpu_regularizers_boost", + sources=[os.path.join("." , "src", "fista_module.cpp" )], + include_dirs=extra_include_dirs, + library_dirs=extra_library_dirs, + extra_compile_args=extra_compile_args, + libraries=extra_libraries ), + + ], + zip_safe = False, + packages = {'ccpi','ccpi.filters'}, +) + +setup( + name='ccpi', + description='CCPi Core Imaging Library - Image Regularizers', + version=cil_version, + cmdclass = {'build_ext': build_ext}, + ext_modules = [Extension("ccpi.filters.cpu_regularizers_cython", + sources=[os.path.join("." , "src", "cpu_regularizers.pyx" ) ], + include_dirs=extra_include_dirs, + library_dirs=extra_library_dirs, + extra_compile_args=extra_compile_args, + libraries=extra_libraries ), + + ], + zip_safe = False, + packages = {'ccpi','ccpi.filters'}, +) + + +@SETUP_GPU_WRAPPERS@ \ No newline at end of file diff --git a/Wrappers/Python/setup.py b/Wrappers/Python/setup.py deleted file mode 100644 index 00c93fc..0000000 --- a/Wrappers/Python/setup.py +++ /dev/null @@ -1,100 +0,0 @@ -#!/usr/bin/env python - -import setuptools -from distutils.core import setup -from distutils.extension import Extension -from Cython.Distutils import build_ext - -import os -import sys -import numpy -import platform - -cil_version=os.environ['CIL_VERSION'] -if cil_version == '': - print("Please set the environmental variable CIL_VERSION") - sys.exit(1) - -library_include_path = "" -library_lib_path = "" -try: - library_include_path = os.environ['LIBRARY_INC'] - library_lib_path = os.environ['LIBRARY_LIB'] -except: - library_include_path = os.environ['PREFIX']+'/include' - pass - -extra_include_dirs = [numpy.get_include(), library_include_path] -#extra_library_dirs = [os.path.join(library_include_path, "..", "lib")] -extra_compile_args = [] -extra_library_dirs = [] -extra_compile_args = [] -extra_link_args = [] -extra_libraries = ['cilreg'] - -extra_include_dirs += [os.path.join(".." , ".." , "Core"), - os.path.join(".." , ".." , "Core", "regularizers_CPU"), - os.path.join(".." , ".." , "Core", "regularizers_GPU") , - "."] - -if platform.system() == 'Windows': - - - extra_compile_args[0:] = ['/DWIN32','/EHsc','/DBOOST_ALL_NO_LIB' , '/openmp' ] - - if sys.version_info.major == 3 : - extra_libraries += ['boost_python3-vc140-mt-1_64', 'boost_numpy3-vc140-mt-1_64'] - else: - extra_libraries += ['boost_python-vc90-mt-1_64', 'boost_numpy-vc90-mt-1_64'] -else: - extra_compile_args = ['-fopenmp','-O2', '-funsigned-char', '-Wall', '-std=c++0x'] - if sys.version_info.major == 3: - extra_libraries += ['boost_python3', 'boost_numpy3','gomp'] - else: - extra_libraries += ['boost_python', 'boost_numpy','gomp'] - -setup( - name='ccpi', - description='CCPi Core Imaging Library - Image Regularizers', - version=cil_version, - cmdclass = {'build_ext': build_ext}, - ext_modules = [Extension("ccpi.filters.gpu_regularizers", - sources=[ - os.path.join("." , "src", "gpu_regularizers.pyx" ), - ], - include_dirs=extra_include_dirs, - library_dirs=extra_library_dirs, - extra_compile_args=extra_compile_args, - libraries=extra_libraries ), - - ], - zip_safe = False, - packages = {'ccpi','ccpi.filters'}, -) - -setup( - name='ccpi', - description='CCPi Core Imaging Library - Image Regularizers', - version=cil_version, - cmdclass = {'build_ext': build_ext}, - ext_modules = [Extension("ccpi.filters.cpu_regularizers", - sources=[os.path.join("." , "src", "fista_module.cpp" ), - os.path.join("." , "src", "cpu_regularizers.pyx" ) - # os.path.join("@CMAKE_SOURCE_DIR@" , "main_func" , "regularizers_CPU", "FGP_TV_core.c"), - # os.path.join("@CMAKE_SOURCE_DIR@" , "main_func" , "regularizers_CPU", "SplitBregman_TV_core.c"), - # os.path.join("@CMAKE_SOURCE_DIR@" , "main_func" , "regularizers_CPU", "LLT_model_core.c"), - # os.path.join("@CMAKE_SOURCE_DIR@" , "main_func" , "regularizers_CPU", "PatchBased_Regul_core.c"), - # os.path.join("@CMAKE_SOURCE_DIR@" , "main_func" , "regularizers_CPU", "TGV_PD_core.c"), - # os.path.join("@CMAKE_SOURCE_DIR@" , "main_func" , "regularizers_CPU", "utils.c") - ], - include_dirs=extra_include_dirs, - library_dirs=extra_library_dirs, - extra_compile_args=extra_compile_args, - libraries=extra_libraries ), - - ], - zip_safe = False, - packages = {'ccpi','ccpi.filters'}, -) - - diff --git a/Wrappers/Python/src/cpu_regularizers.pyx b/Wrappers/Python/src/cpu_regularizers.pyx index e69de29..a8f8c8f 100644 --- a/Wrappers/Python/src/cpu_regularizers.pyx +++ b/Wrappers/Python/src/cpu_regularizers.pyx @@ -0,0 +1,19 @@ +# distutils: language=c++ +""" +Copyright 2018 CCPi +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + http://www.apache.org/licenses/LICENSE-2.0 +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. + +Author: Edoardo Pasca +""" + +import cython +import numpy as np +cimport numpy as np diff --git a/Wrappers/Python/src/fista_module.cpp b/Wrappers/Python/src/fista_module.cpp index 3876cad..cef3ecc 100644 --- a/Wrappers/Python/src/fista_module.cpp +++ b/Wrappers/Python/src/fista_module.cpp @@ -1028,13 +1028,13 @@ bp::list TGV_PD(np::ndarray input, double d_lambda, double d_alpha1, double d_al return result; } -BOOST_PYTHON_MODULE(cpu_regularizers) +BOOST_PYTHON_MODULE(cpu_regularizers_boost) { np::initialize(); //To specify that this module is a package bp::object package = bp::scope(); - package.attr("__path__") = "cpu_regularizers"; + package.attr("__path__") = "cpu_regularizers_boost"; np::dtype dt1 = np::dtype::get_builtin(); np::dtype dt2 = np::dtype::get_builtin(); diff --git a/Wrappers/Python/test/test_cpu_regularizers.py b/Wrappers/Python/test/test_cpu_regularizers.py index 6c97875..9713baa 100644 --- a/Wrappers/Python/test/test_cpu_regularizers.py +++ b/Wrappers/Python/test/test_cpu_regularizers.py @@ -11,8 +11,9 @@ import numpy as np import os from enum import Enum import timeit -from ccpi.filters.cpu_regularizers import SplitBregman_TV , FGP_TV , LLT_model, \ - PatchBased_Regul , TGV_PD +from ccpi.filters.cpu_regularizers_boost import SplitBregman_TV , FGP_TV ,\ + LLT_model, PatchBased_Regul ,\ + TGV_PD ############################################################################### #https://stackoverflow.com/questions/13875989/comparing-image-in-url-to-image-in-filesystem-in-python/13884956#13884956 -- cgit v1.2.3 From e3611de2c894db60296afcfd6afac09774a564e8 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Mon, 29 Jan 2018 15:50:01 +0000 Subject: added CMakeLists.txt for Python Wrapper --- Wrappers/Python/CMakeLists.txt | 110 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 110 insertions(+) create mode 100644 Wrappers/Python/CMakeLists.txt diff --git a/Wrappers/Python/CMakeLists.txt b/Wrappers/Python/CMakeLists.txt new file mode 100644 index 0000000..e4c847c --- /dev/null +++ b/Wrappers/Python/CMakeLists.txt @@ -0,0 +1,110 @@ +# Copyright 2018 Edoardo Pasca +cmake_minimum_required (VERSION 3.0) + +project(RegularizerPython) +#https://stackoverflow.com/questions/13298504/using-cmake-with-setup-py + +# The version number. + +set (CIL_VERSION $ENV{CIL_VERSION} CACHE INTERNAL "Core Imaging Library version" FORCE) + +# conda orchestrated build +message("CIL_VERSION ${CIL_VERSION}") +#include (GenerateExportHeader) + +find_package(PythonInterp REQUIRED) +if (PYTHONINTERP_FOUND) + message ("Current Python " ${PYTHON_VERSION_STRING} " found " ${PYTHON_EXECUTABLE}) + if (PYTHON_VERSION_MAJOR EQUAL "3") + set (BOOST_PYTHON "python3") + set (BOOST_NUMPY "numpy3") + else() + set (BOOST_PYTHON "python") + set (BOOST_NUMPY "numpy") + endif() +endif() + +find_package(Boost REQUIRED + COMPONENTS ${BOOST_PYTHON} ${BOOST_NUMPY}) + +if (Boost_FOUND) + message("Boost version " ${Boost_VERSION}) + message("Boost include dir " ${Boost_INCLUDE_DIRS}) + message("Boost library dir " ${Boost_LIBRARY_DIRS}) + message("Boost libraries " ${Boost_LIBRARIES}) +endif() + +## Build the regularizers package as a library +message("Creating Regularizers as shared library") + +message("CMAKE_CXX_FLAGS ${CMAKE_CXX_FLAGS}") + +set(CMAKE_BUILD_TYPE "Release") + +if(WIN32) + set (FLAGS "/DWIN32 /EHsc /DBOOST_ALL_NO_LIB /openmp /DCCPiCore_EXPORTS") + set (CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} /NODEFAULTLIB:MSVCRT.lib") + + set (EXTRA_LIBRARIES + ${Boost_LIBRARIES} + #"tiff" + ) + + message("library lib: ${LIBRARY_LIB}") + +elseif(UNIX) + set (FLAGS "-fopenmp -O2 -funsigned-char -Wall -Wl,--no-undefined -DCCPiReconstructionIterative_EXPORTS -std=c++0x") + set (EXTRA_LIBRARIES + ${Boost_LIBRARIES} + #"tiff" + "gomp" + ) +endif() + + +#set (BOOST_PYTHON_LIB ${Boost_${BOOST_PYTHON}_LIBRARY}) +#set (BOOST_NUMPY_LIB ${Boost_${BOOST_NUMPY}_LIBRARY}) + +list(GET Boost_LIBRARIES 0 place ) +get_filename_component(BOOST_PYTHON_LIB ${place} NAME_WE ) +list(GET Boost_LIBRARIES 1 place ) +get_filename_component(BOOST_NUMPY_LIB ${place} NAME_WE ) + +message ("found " ${BOOST_PYTHON_LIB}) +message ("found " ${BOOST_NUMPY_LIB}) + + + +# GPU Regularizers + +find_package(CUDA) +if (CUDA_FOUND) + message("CUDA FOUND") + set (SETUP_GPU_WRAPPERS "setup( \n\ + name='ccpi', \n\ + description='CCPi Core Imaging Library - Image Regularizers GPU',\n\ + version=cil_version,\n\ + cmdclass = {'build_ext': build_ext},\n\ + ext_modules = [Extension('ccpi.filters.gpu_regularizers',\n\ + sources=[ \n\ + os.path.join('.' , 'src', 'gpu_regularizers.pyx' ),\n\ + ],\n\ + include_dirs=extra_include_dirs, \n\ + library_dirs=extra_library_dirs, \n\ + extra_compile_args=extra_compile_args, \n\ + libraries=extra_libraries ), \n\ + ],\n\ + zip_safe = False, \n\ + packages = {'ccpi','ccpi.filters'},\n\ +)") +else() + message("CUDA NOT FOUND") + set(SETUP_GPU_WRAPPERS "#CUDA NOT FOUND") +endif() + +configure_file("setup-regularizers.py.in" "setup-regularizers.py") + + +#add_executable(regularizer_test ${CMAKE_CURRENT_SOURCE_DIR}/test/test_regularizer.cpp) + +#target_link_libraries (regularizer_test LINK_PUBLIC regularizers_lib) -- cgit v1.2.3 From 382bcab9e036418df88f35e3332aa44fdef6669e Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Mon, 29 Jan 2018 16:02:13 +0000 Subject: removed initial lib from library name (unix) --- Wrappers/Python/CMakeLists.txt | 23 ++++++++++++++++------- 1 file changed, 16 insertions(+), 7 deletions(-) diff --git a/Wrappers/Python/CMakeLists.txt b/Wrappers/Python/CMakeLists.txt index e4c847c..ab4f400 100644 --- a/Wrappers/Python/CMakeLists.txt +++ b/Wrappers/Python/CMakeLists.txt @@ -62,13 +62,22 @@ elseif(UNIX) endif() -#set (BOOST_PYTHON_LIB ${Boost_${BOOST_PYTHON}_LIBRARY}) -#set (BOOST_NUMPY_LIB ${Boost_${BOOST_NUMPY}_LIBRARY}) - -list(GET Boost_LIBRARIES 0 place ) -get_filename_component(BOOST_PYTHON_LIB ${place} NAME_WE ) -list(GET Boost_LIBRARIES 1 place ) -get_filename_component(BOOST_NUMPY_LIB ${place} NAME_WE ) +if (WIN32) + list(GET Boost_LIBRARIES 0 place ) + get_filename_component(BOOST_PYTHON_LIB ${place} NAME_WE ) + list(GET Boost_LIBRARIES 1 place ) + get_filename_component(BOOST_NUMPY_LIB ${place} NAME_WE ) +else() + # on linux the library looks like libboost_numpy3.so: + # we need to get rid of the lib at the beginning + list(GET Boost_LIBRARIES 0 place ) + get_filename_component(place2 ${place} NAME_WE ) + string(REGEX REPLACE "^lib.*" "\\1" $BOOST_PYTHON_LIB "${place2}") + + list(GET Boost_LIBRARIES 1 place ) + get_filename_component(place2 ${place} NAME_WE ) + string(REGEX REPLACE "^lib.*" "\\1" $BOOST_NUMPY_LIB "${place2}") +endif() message ("found " ${BOOST_PYTHON_LIB}) message ("found " ${BOOST_NUMPY_LIB}) -- cgit v1.2.3 From a5979934c10aee7beeae49f254c57adbbb2fa13c Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Mon, 29 Jan 2018 16:14:16 +0000 Subject: fixed Wrapper Build for linux --- Wrappers/Python/CMakeLists.txt | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Wrappers/Python/CMakeLists.txt b/Wrappers/Python/CMakeLists.txt index ab4f400..3f95660 100644 --- a/Wrappers/Python/CMakeLists.txt +++ b/Wrappers/Python/CMakeLists.txt @@ -72,11 +72,11 @@ else() # we need to get rid of the lib at the beginning list(GET Boost_LIBRARIES 0 place ) get_filename_component(place2 ${place} NAME_WE ) - string(REGEX REPLACE "^lib.*" "\\1" $BOOST_PYTHON_LIB "${place2}") + string(REGEX REPLACE "^lib(.*)" "\\1" BOOST_PYTHON_LIB "${place2}") list(GET Boost_LIBRARIES 1 place ) get_filename_component(place2 ${place} NAME_WE ) - string(REGEX REPLACE "^lib.*" "\\1" $BOOST_NUMPY_LIB "${place2}") + string(REGEX REPLACE "^lib(.*)" "\\1" BOOST_NUMPY_LIB "${place2}") endif() message ("found " ${BOOST_PYTHON_LIB}) -- cgit v1.2.3 From 56a94d0e2bf32a5951482dc9dc3c4c4c652755fb Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Mon, 29 Jan 2018 16:32:25 +0000 Subject: added demo directory --- Wrappers/Python/demo/test_cpu_regularizers.py | 442 ++++++++++++++++++++++++++ 1 file changed, 442 insertions(+) create mode 100644 Wrappers/Python/demo/test_cpu_regularizers.py diff --git a/Wrappers/Python/demo/test_cpu_regularizers.py b/Wrappers/Python/demo/test_cpu_regularizers.py new file mode 100644 index 0000000..9713baa --- /dev/null +++ b/Wrappers/Python/demo/test_cpu_regularizers.py @@ -0,0 +1,442 @@ +# -*- coding: utf-8 -*- +""" +Created on Fri Aug 4 11:10:05 2017 + +@author: ofn77899 +""" + + +import matplotlib.pyplot as plt +import numpy as np +import os +from enum import Enum +import timeit +from ccpi.filters.cpu_regularizers_boost import SplitBregman_TV , FGP_TV ,\ + LLT_model, PatchBased_Regul ,\ + TGV_PD + +############################################################################### +#https://stackoverflow.com/questions/13875989/comparing-image-in-url-to-image-in-filesystem-in-python/13884956#13884956 +#NRMSE a normalization of the root of the mean squared error +#NRMSE is simply 1 - [RMSE / (maxval - minval)]. Where maxval is the maximum +# intensity from the two images being compared, and respectively the same for +# minval. RMSE is given by the square root of MSE: +# sqrt[(sum(A - B) ** 2) / |A|], +# where |A| means the number of elements in A. By doing this, the maximum value +# given by RMSE is maxval. + +def nrmse(im1, im2): + a, b = im1.shape + rmse = np.sqrt(np.sum((im2 - im1) ** 2) / float(a * b)) + max_val = max(np.max(im1), np.max(im2)) + min_val = min(np.min(im1), np.min(im2)) + return 1 - (rmse / (max_val - min_val)) +############################################################################### +def printParametersToString(pars): + txt = r'' + for key, value in pars.items(): + if key== 'algorithm' : + txt += "{0} = {1}".format(key, value.__name__) + elif key == 'input': + txt += "{0} = {1}".format(key, np.shape(value)) + else: + txt += "{0} = {1}".format(key, value) + txt += '\n' + return txt +############################################################################### +# +# 2D Regularizers +# +############################################################################### +#Example: +# figure; +# Im = double(imread('lena_gray_256.tif'))/255; % loading image +# u0 = Im + .05*randn(size(Im)); u0(u0 < 0) = 0; +# u = SplitBregman_TV(single(u0), 10, 30, 1e-04); + +# assumes the script is launched from the test directory +filename = os.path.join(".." , ".." , ".." , "data" ,"lena_gray_512.tif") +#filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-FISTA_reconstruction\data\lena_gray_512.tif" +#filename = r"/home/ofn77899/Reconstruction/CCPi-FISTA_Reconstruction/data/lena_gray_512.tif" +#filename = r'/home/algol/Documents/Python/STD_test_images/lena_gray_512.tif' + +Im = plt.imread(filename) +Im = np.asarray(Im, dtype='float32') + +perc = 0.15 +u0 = Im + np.random.normal(loc = Im , + scale = perc * Im , + size = np.shape(Im)) +# map the u0 u0->u0>0 +f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1) +u0 = f(u0).astype('float32') + +## plot +fig = plt.figure() + +a=fig.add_subplot(2,3,1) +a.set_title('noise') +imgplot = plt.imshow(u0#,cmap="gray" + ) + +reg_output = [] +############################################################################## +# Call regularizer + +####################### SplitBregman_TV ##################################### +# u = SplitBregman_TV(single(u0), 10, 30, 1e-04); + +start_time = timeit.default_timer() +pars = {'algorithm' : SplitBregman_TV , \ + 'input' : u0, + 'regularization_parameter':10. , \ +'number_of_iterations' :35 ,\ +'tolerance_constant':0.0001 , \ +'TV_penalty': 0 +} + +out = SplitBregman_TV (pars['input'], pars['regularization_parameter'], + pars['number_of_iterations'], + pars['tolerance_constant'], + pars['TV_penalty']) +splitbregman = out[0] +txtstr = printParametersToString(pars) +txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) +print (txtstr) + + +a=fig.add_subplot(2,3,2) + + +# these are matplotlib.patch.Patch properties +props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) +# place a text box in upper left in axes coords +a.text(0.05, 0.95, txtstr, transform=a.transAxes, fontsize=14, + verticalalignment='top', bbox=props) +imgplot = plt.imshow(splitbregman,\ + #cmap="gray" + ) + +###################### FGP_TV ######################################### +# u = FGP_TV(single(u0), 0.05, 100, 1e-04); +start_time = timeit.default_timer() +pars = {'algorithm' : FGP_TV , \ + 'input' : u0, + 'regularization_parameter':5e-4, \ + 'number_of_iterations' :10 ,\ + 'tolerance_constant':0.001,\ + 'TV_penalty': 0 +} + +out = FGP_TV (pars['input'], + pars['regularization_parameter'], + pars['number_of_iterations'], + pars['tolerance_constant'], + pars['TV_penalty']) + +fgp = out[0] +txtstr = printParametersToString(pars) +txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) +print (txtstr) + + +a=fig.add_subplot(2,3,3) + +# these are matplotlib.patch.Patch properties +props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) +# place a text box in upper left in axes coords +imgplot = plt.imshow(fgp, \ + #cmap="gray" + ) +# place a text box in upper left in axes coords +a.text(0.05, 0.95, txtstr, transform=a.transAxes, fontsize=14, + verticalalignment='top', bbox=props) + +###################### LLT_model ######################################### + +start_time = timeit.default_timer() + +pars = {'algorithm': LLT_model , \ + 'input' : u0, + 'regularization_parameter': 25,\ + 'time_step':0.0003, \ + 'number_of_iterations' :300,\ + 'tolerance_constant':0.001,\ + 'restrictive_Z_smoothing': 0 +} +out = LLT_model(pars['input'], + pars['regularization_parameter'], + pars['time_step'] , + pars['number_of_iterations'], + pars['tolerance_constant'], + pars['restrictive_Z_smoothing'] ) + +llt = out[0] +txtstr = printParametersToString(pars) +txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) +print (txtstr) +a=fig.add_subplot(2,3,4) + +# these are matplotlib.patch.Patch properties +props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) +# place a text box in upper left in axes coords +a.text(0.05, 0.95, txtstr, transform=a.transAxes, fontsize=14, + verticalalignment='top', bbox=props) +imgplot = plt.imshow(llt,\ + #cmap="gray" + ) + + +# ###################### PatchBased_Regul ######################################### +# # Quick 2D denoising example in Matlab: +# # Im = double(imread('lena_gray_256.tif'))/255; % loading image +# # u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise +# # ImDen = PB_Regul_CPU(single(u0), 3, 1, 0.08, 0.05); +start_time = timeit.default_timer() + +pars = {'algorithm': PatchBased_Regul , \ + 'input' : u0, + 'regularization_parameter': 0.05,\ + 'searching_window_ratio':3, \ + 'similarity_window_ratio':1,\ + 'PB_filtering_parameter': 0.08 +} +out = PatchBased_Regul(pars['input'], + pars['regularization_parameter'], + pars['searching_window_ratio'] , + pars['similarity_window_ratio'] , + pars['PB_filtering_parameter']) +pbr = out[0] +txtstr = printParametersToString(pars) +txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) +print (txtstr) + +a=fig.add_subplot(2,3,5) + + +# these are matplotlib.patch.Patch properties +props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) +# place a text box in upper left in axes coords +a.text(0.05, 0.95, txtstr, transform=a.transAxes, fontsize=14, + verticalalignment='top', bbox=props) +imgplot = plt.imshow(pbr #,cmap="gray" + ) + + +# ###################### TGV_PD ######################################### +# # Quick 2D denoising example in Matlab: +# # Im = double(imread('lena_gray_256.tif'))/255; % loading image +# # u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise +# # u = PrimalDual_TGV(single(u0), 0.02, 1.3, 1, 550); + +start_time = timeit.default_timer() + +pars = {'algorithm': TGV_PD , \ + 'input' : u0,\ + 'regularization_parameter':0.05,\ + 'first_order_term': 1.3,\ + 'second_order_term': 1, \ + 'number_of_iterations': 550 + } +out = TGV_PD(pars['input'], + pars['regularization_parameter'], + pars['first_order_term'] , + pars['second_order_term'] , + pars['number_of_iterations']) +tgv = out[0] +txtstr = printParametersToString(pars) +txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) +print (txtstr) +a=fig.add_subplot(2,3,6) + +# these are matplotlib.patch.Patch properties +props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) +# place a text box in upper left in axes coords +a.text(0.05, 0.95, txtstr, transform=a.transAxes, fontsize=14, + verticalalignment='top', bbox=props) +imgplot = plt.imshow(tgv #, cmap="gray") + ) + + +plt.show() + +################################################################################ +## +## 3D Regularizers +## +################################################################################ +##Example: +## figure; +## Im = double(imread('lena_gray_256.tif'))/255; % loading image +## u0 = Im + .05*randn(size(Im)); u0(u0 < 0) = 0; +## u = SplitBregman_TV(single(u0), 10, 30, 1e-04); +# +##filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-Reconstruction\python\test\reconstruction_example.mha" +#filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-Simpleflex\data\head.mha" +# +#reader = vtk.vtkMetaImageReader() +#reader.SetFileName(os.path.normpath(filename)) +#reader.Update() +##vtk returns 3D images, let's take just the one slice there is as 2D +#Im = Converter.vtk2numpy(reader.GetOutput()) +#Im = Im.astype('float32') +##imgplot = plt.imshow(Im) +#perc = 0.05 +#u0 = Im + (perc* np.random.normal(size=np.shape(Im))) +## map the u0 u0->u0>0 +#f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1) +#u0 = f(u0).astype('float32') +#converter = Converter.numpy2vtkImporter(u0, reader.GetOutput().GetSpacing(), +# reader.GetOutput().GetOrigin()) +#converter.Update() +#writer = vtk.vtkMetaImageWriter() +#writer.SetInputData(converter.GetOutput()) +#writer.SetFileName(r"C:\Users\ofn77899\Documents\GitHub\CCPi-FISTA_reconstruction\data\noisy_head.mha") +##writer.Write() +# +# +### plot +#fig3D = plt.figure() +##a=fig.add_subplot(3,3,1) +##a.set_title('Original') +##imgplot = plt.imshow(Im) +#sliceNo = 32 +# +#a=fig3D.add_subplot(2,3,1) +#a.set_title('noise') +#imgplot = plt.imshow(u0.T[sliceNo]) +# +#reg_output3d = [] +# +############################################################################### +## Call regularizer +# +######################## SplitBregman_TV ##################################### +## u = SplitBregman_TV(single(u0), 10, 30, 1e-04); +# +##reg = Regularizer(Regularizer.Algorithm.SplitBregman_TV) +# +##out = reg(input=u0, regularization_parameter=10., #number_of_iterations=30, +## #tolerance_constant=1e-4, +## TV_Penalty=Regularizer.TotalVariationPenalty.l1) +# +#out2 = Regularizer.SplitBregman_TV(input=u0, regularization_parameter=10., number_of_iterations=30, +# tolerance_constant=1e-4, +# TV_Penalty=Regularizer.TotalVariationPenalty.l1) +# +# +#pars = out2[-2] +#reg_output3d.append(out2) +# +#a=fig3D.add_subplot(2,3,2) +# +# +#textstr = out2[-1] +# +# +## these are matplotlib.patch.Patch properties +#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) +## place a text box in upper left in axes coords +#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14, +# verticalalignment='top', bbox=props) +#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo]) +# +####################### FGP_TV ######################################### +## u = FGP_TV(single(u0), 0.05, 100, 1e-04); +#out2 = Regularizer.FGP_TV(input=u0, regularization_parameter=0.005, +# number_of_iterations=200) +#pars = out2[-2] +#reg_output3d.append(out2) +# +#a=fig3D.add_subplot(2,3,2) +# +# +#textstr = out2[-1] +# +# +## these are matplotlib.patch.Patch properties +#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) +## place a text box in upper left in axes coords +#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14, +# verticalalignment='top', bbox=props) +#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo]) +# +####################### LLT_model ######################################### +## * u0 = Im + .03*randn(size(Im)); % adding noise +## [Den] = LLT_model(single(u0), 10, 0.1, 1); +##Den = LLT_model(single(u0), 25, 0.0003, 300, 0.0001, 0); +##input, regularization_parameter , time_step, number_of_iterations, +## tolerance_constant, restrictive_Z_smoothing=0 +#out2 = Regularizer.LLT_model(input=u0, regularization_parameter=25, +# time_step=0.0003, +# tolerance_constant=0.0001, +# number_of_iterations=300) +#pars = out2[-2] +#reg_output3d.append(out2) +# +#a=fig3D.add_subplot(2,3,2) +# +# +#textstr = out2[-1] +# +# +## these are matplotlib.patch.Patch properties +#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) +## place a text box in upper left in axes coords +#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14, +# verticalalignment='top', bbox=props) +#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo]) +# +####################### PatchBased_Regul ######################################### +## Quick 2D denoising example in Matlab: +## Im = double(imread('lena_gray_256.tif'))/255; % loading image +## u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise +## ImDen = PB_Regul_CPU(single(u0), 3, 1, 0.08, 0.05); +# +#out2 = Regularizer.PatchBased_Regul(input=u0, regularization_parameter=0.05, +# searching_window_ratio=3, +# similarity_window_ratio=1, +# PB_filtering_parameter=0.08) +#pars = out2[-2] +#reg_output3d.append(out2) +# +#a=fig3D.add_subplot(2,3,2) +# +# +#textstr = out2[-1] +# +# +## these are matplotlib.patch.Patch properties +#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) +## place a text box in upper left in axes coords +#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14, +# verticalalignment='top', bbox=props) +#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo]) +# + +###################### TGV_PD ######################################### +# Quick 2D denoising example in Matlab: +# Im = double(imread('lena_gray_256.tif'))/255; % loading image +# u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise +# u = PrimalDual_TGV(single(u0), 0.02, 1.3, 1, 550); + + +#out2 = Regularizer.TGV_PD(input=u0, regularization_parameter=0.05, +# first_order_term=1.3, +# second_order_term=1, +# number_of_iterations=550) +#pars = out2[-2] +#reg_output3d.append(out2) +# +#a=fig3D.add_subplot(2,3,2) +# +# +#textstr = out2[-1] +# +# +## these are matplotlib.patch.Patch properties +#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) +## place a text box in upper left in axes coords +#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14, +# verticalalignment='top', bbox=props) +#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo]) -- cgit v1.2.3 From 0b1ee0c14e8962edbd496652ebcddc09274fad62 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Tue, 30 Jan 2018 11:36:59 +0000 Subject: fix GPU regularizer --- Wrappers/Python/CMakeLists.txt | 3 ++- Wrappers/Python/src/gpu_regularizers.pyx | 39 ++++++++++++++++---------------- 2 files changed, 21 insertions(+), 21 deletions(-) diff --git a/Wrappers/Python/CMakeLists.txt b/Wrappers/Python/CMakeLists.txt index 3f95660..16a2137 100644 --- a/Wrappers/Python/CMakeLists.txt +++ b/Wrappers/Python/CMakeLists.txt @@ -89,7 +89,8 @@ message ("found " ${BOOST_NUMPY_LIB}) find_package(CUDA) if (CUDA_FOUND) message("CUDA FOUND") - set (SETUP_GPU_WRAPPERS "setup( \n\ + set (SETUP_GPU_WRAPPERS "extra_libraries += ['cilregcuda']\n\ +setup( \n\ name='ccpi', \n\ description='CCPi Core Imaging Library - Image Regularizers GPU',\n\ version=cil_version,\n\ diff --git a/Wrappers/Python/src/gpu_regularizers.pyx b/Wrappers/Python/src/gpu_regularizers.pyx index 7658e36..dc625c3 100644 --- a/Wrappers/Python/src/gpu_regularizers.pyx +++ b/Wrappers/Python/src/gpu_regularizers.pyx @@ -44,7 +44,25 @@ def Diff4thHajiaboli(inputData, regularization_parameter, iterations, edge_preserving_parameter) - + +def NML(inputData, + SearchW_real, + SimilW, + h, + lambdaf): + if inputData.ndim == 2: + return NML2D(inputData, + SearchW_real, + SimilW, + h, + lambdaf) + elif inputData.ndim == 3: + return NML3D(inputData, + SearchW_real, + SimilW, + h, + lambdaf) + def Diff4thHajiaboli2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, float regularization_parameter, int iterations, @@ -160,25 +178,6 @@ def Diff4thHajiaboli3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, return B -def NML(inputData, - regularization_parameter, - iterations, - edge_preserving_parameter): - if inputData.ndim == 2: - return NML2D(inputData, - regularization_parameter, - iterations, - edge_preserving_parameter) - elif inputData.ndim == 3: - return NML3D(inputData, - regularization_parameter, - iterations, - edge_preserving_parameter) - - #SearchW_real = (int) mxGetScalar(prhs[1]); /* the searching window ratio */ - #SimilW = (int) mxGetScalar(prhs[2]); /* the similarity window ratio */ - #h = (float) mxGetScalar(prhs[3]); /* parameter for the PB filtering function */ - #lambda = (float) mxGetScalar(prhs[4]); def NML2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, SearchW_real, -- cgit v1.2.3