From 09efb291c7edc57cff7300d16adc627fce13bcc7 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Tue, 13 Feb 2018 16:36:26 +0000 Subject: removed and renamed a few things --- Wrappers/Python/FindAnacondaEnvironment.cmake | 154 ---- Wrappers/Python/fista_module.cpp | 1047 ------------------------- Wrappers/Python/setup-regularizers.py.in | 2 +- Wrappers/Python/src/cpu_regularizers.cpp | 1047 +++++++++++++++++++++++++ Wrappers/Python/src/fista_module.cpp | 1047 ------------------------- Wrappers/Python/src/multiply.pyx | 49 -- 6 files changed, 1048 insertions(+), 2298 deletions(-) delete mode 100644 Wrappers/Python/FindAnacondaEnvironment.cmake delete mode 100644 Wrappers/Python/fista_module.cpp create mode 100644 Wrappers/Python/src/cpu_regularizers.cpp delete mode 100644 Wrappers/Python/src/fista_module.cpp delete mode 100644 Wrappers/Python/src/multiply.pyx (limited to 'Wrappers') diff --git a/Wrappers/Python/FindAnacondaEnvironment.cmake b/Wrappers/Python/FindAnacondaEnvironment.cmake deleted file mode 100644 index 6475128..0000000 --- a/Wrappers/Python/FindAnacondaEnvironment.cmake +++ /dev/null @@ -1,154 +0,0 @@ -# Copyright 2017 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. - -# #.rst: -# FindAnacondaEnvironment -# -------------- -# -# Find Python executable and library for a specific Anaconda environment -# -# This module finds the Python interpreter for a specific Anaconda enviroment, -# if installed and determines where the include files and libraries are. -# This code sets the following variables: -# -# :: -# PYTHONINTERP_FOUND - if the Python interpret has been found -# PYTHON_EXECUTABLE - the Python interpret found -# PYTHON_LIBRARY - path to the python library -# PYTHON_INCLUDE_PATH - path to where Python.h is found (deprecated) -# PYTHON_INCLUDE_DIRS - path to where Python.h is found -# PYTHONLIBS_VERSION_STRING - version of the Python libs found (since CMake 2.8.8) -# PYTHON_VERSION_MAJOR - major Python version -# PYTHON_VERSION_MINOR - minor Python version -# PYTHON_VERSION_PATCH - patch Python version - - - -function (findPythonForAnacondaEnvironment env) - if (WIN32) - file(TO_CMAKE_PATH ${env}/python.exe PYTHON_EXECUTABLE) - elseif (UNIX) - file(TO_CMAKE_PATH ${env}/bin/python PYTHON_EXECUTABLE) - endif() - - - message("findPythonForAnacondaEnvironment Found Python Executable" ${PYTHON_EXECUTABLE}) - ####### FROM FindPythonInterpr ######## - # determine python version string - if(PYTHON_EXECUTABLE) - execute_process(COMMAND "${PYTHON_EXECUTABLE}" -c - "import sys; sys.stdout.write(';'.join([str(x) for x in sys.version_info[:3]]))" - OUTPUT_VARIABLE _VERSION - RESULT_VARIABLE _PYTHON_VERSION_RESULT - ERROR_QUIET) - if(NOT _PYTHON_VERSION_RESULT) - string(REPLACE ";" "." _PYTHON_VERSION_STRING "${_VERSION}") - list(GET _VERSION 0 _PYTHON_VERSION_MAJOR) - list(GET _VERSION 1 _PYTHON_VERSION_MINOR) - list(GET _VERSION 2 _PYTHON_VERSION_PATCH) - if(PYTHON_VERSION_PATCH EQUAL 0) - # it's called "Python 2.7", not "2.7.0" - string(REGEX REPLACE "\\.0$" "" _PYTHON_VERSION_STRING "${PYTHON_VERSION_STRING}") - endif() - else() - # sys.version predates sys.version_info, so use that - execute_process(COMMAND "${PYTHON_EXECUTABLE}" -c "import sys; sys.stdout.write(sys.version)" - OUTPUT_VARIABLE _VERSION - RESULT_VARIABLE _PYTHON_VERSION_RESULT - ERROR_QUIET) - if(NOT _PYTHON_VERSION_RESULT) - string(REGEX REPLACE " .*" "" _PYTHON_VERSION_STRING "${_VERSION}") - string(REGEX REPLACE "^([0-9]+)\\.[0-9]+.*" "\\1" _PYTHON_VERSION_MAJOR "${PYTHON_VERSION_STRING}") - string(REGEX REPLACE "^[0-9]+\\.([0-9])+.*" "\\1" _PYTHON_VERSION_MINOR "${PYTHON_VERSION_STRING}") - if(PYTHON_VERSION_STRING MATCHES "^[0-9]+\\.[0-9]+\\.([0-9]+)") - set(PYTHON_VERSION_PATCH "${CMAKE_MATCH_1}") - else() - set(PYTHON_VERSION_PATCH "0") - endif() - else() - # sys.version was first documented for Python 1.5, so assume - # this is older. - set(PYTHON_VERSION_STRING "1.4" PARENT_SCOPE) - set(PYTHON_VERSION_MAJOR "1" PARENT_SCOPE) - set(PYTHON_VERSION_MINOR "4" PARENT_SCOPE) - set(PYTHON_VERSION_PATCH "0" PARENT_SCOPE) - endif() - endif() - unset(_PYTHON_VERSION_RESULT) - unset(_VERSION) - endif() - ############################################### - - set (PYTHON_EXECUTABLE ${PYTHON_EXECUTABLE} PARENT_SCOPE) - set (PYTHONINTERP_FOUND "ON" PARENT_SCOPE) - set (PYTHON_VERSION_STRING ${_PYTHON_VERSION_STRING} PARENT_SCOPE) - set (PYTHON_VERSION_MAJOR ${_PYTHON_VERSION_MAJOR} PARENT_SCOPE) - set (PYTHON_VERSION_MINOR ${_PYTHON_VERSION_MINOR} PARENT_SCOPE) - set (PYTHON_VERSION_PATCH ${_PYTHON_VERSION_PATCH} PARENT_SCOPE) - message("My version found " ${PYTHON_VERSION_STRING}) - ## find conda executable - if (WIN32) - set (CONDA_EXECUTABLE ${env}/Script/conda PARENT_SCOPE) - elseif(UNIX) - set (CONDA_EXECUTABLE ${env}/bin/conda PARENT_SCOPE) - endif() -endfunction() - - - -set(Python_ADDITIONAL_VERSIONS 3.5) - -find_package(PythonInterp) -if (PYTHONINTERP_FOUND) - - message("Found interpret " ${PYTHON_EXECUTABLE}) - message("Python Library " ${PYTHON_LIBRARY}) - message("Python Include Dir " ${PYTHON_INCLUDE_DIR}) - message("Python Include Path " ${PYTHON_INCLUDE_PATH}) - - foreach(pv ${PYTHON_VERSION_STRING}) - message("Found interpret " ${pv}) - endforeach() -endif() - - - -find_package(PythonLibs) -if (PYTHONLIB_FOUND) - message("Found PythonLibs PYTHON_LIBRARIES " ${PYTHON_LIBRARIES}) - message("Found PythonLibs PYTHON_INCLUDE_PATH " ${PYTHON_INCLUDE_PATH}) - message("Found PythonLibs PYTHON_INCLUDE_DIRS " ${PYTHON_INCLUDE_DIRS}) - message("Found PythonLibs PYTHONLIBS_VERSION_STRING " ${PYTHONLIBS_VERSION_STRING} ) -else() - message("No PythonLibs Found") -endif() - - - - -function(findPythonPackagesPath) - execute_process(COMMAND ${PYTHON_EXECUTABLE} -c "from distutils.sysconfig import *; print (get_python_lib())" - RESULT_VARIABLE PYTHON_CVPY_PROCESS - OUTPUT_VARIABLE PYTHON_STD_PACKAGES_PATH - OUTPUT_STRIP_TRAILING_WHITESPACE) - #message("STD_PACKAGES " ${PYTHON_STD_PACKAGES_PATH}) - if("${PYTHON_STD_PACKAGES_PATH}" MATCHES "site-packages") - set(_PYTHON_PACKAGES_PATH "python${PYTHON_VERSION_MAJOR_MINOR}/site-packages") - endif() - - SET(PYTHON_PACKAGES_PATH "${PYTHON_STD_PACKAGES_PATH}" PARENT_SCOPE) - -endfunction() - - diff --git a/Wrappers/Python/fista_module.cpp b/Wrappers/Python/fista_module.cpp deleted file mode 100644 index 3876cad..0000000 --- a/Wrappers/Python/fista_module.cpp +++ /dev/null @@ -1,1047 +0,0 @@ -/* -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/setup-regularizers.py.in b/Wrappers/Python/setup-regularizers.py.in index 3f4d2d7..a125261 100644 --- a/Wrappers/Python/setup-regularizers.py.in +++ b/Wrappers/Python/setup-regularizers.py.in @@ -52,7 +52,7 @@ setup( version=cil_version, cmdclass = {'build_ext': build_ext}, ext_modules = [Extension("ccpi.filters.cpu_regularizers_boost", - sources=[os.path.join("." , "src", "fista_module.cpp" )], + sources=[os.path.join("." , "src", "cpu_regularizers.cpp" )], include_dirs=extra_include_dirs, library_dirs=extra_library_dirs, extra_compile_args=extra_compile_args, diff --git a/Wrappers/Python/src/cpu_regularizers.cpp b/Wrappers/Python/src/cpu_regularizers.cpp new file mode 100644 index 0000000..e311570 --- /dev/null +++ b/Wrappers/Python/src/cpu_regularizers.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 4) { + 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_boost) +{ + np::initialize(); + + //To specify that this module is a package + bp::object package = bp::scope(); + package.attr("__path__") = "cpu_regularizers_boost"; + + 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.cpp b/Wrappers/Python/src/fista_module.cpp deleted file mode 100644 index e311570..0000000 --- a/Wrappers/Python/src/fista_module.cpp +++ /dev/null @@ -1,1047 +0,0 @@ -/* -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 4) { - 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_boost) -{ - np::initialize(); - - //To specify that this module is a package - bp::object package = bp::scope(); - package.attr("__path__") = "cpu_regularizers_boost"; - - 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/multiply.pyx b/Wrappers/Python/src/multiply.pyx deleted file mode 100644 index 65df1c6..0000000 --- a/Wrappers/Python/src/multiply.pyx +++ /dev/null @@ -1,49 +0,0 @@ -""" -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