diff options
Diffstat (limited to 'Wrappers/Python/src')
| -rw-r--r-- | Wrappers/Python/src/cpu_regularizers.cpp | 756 |
1 files changed, 0 insertions, 756 deletions
diff --git a/Wrappers/Python/src/cpu_regularizers.cpp b/Wrappers/Python/src/cpu_regularizers.cpp deleted file mode 100644 index 3529ebd..0000000 --- a/Wrappers/Python/src/cpu_regularizers.cpp +++ /dev/null @@ -1,756 +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 <iostream> -#include <cmath> - -#include <boost/python.hpp> -#include <boost/python/numpy.hpp> -#include "boost/tuple/tuple.hpp" - -#include "SplitBregman_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 <windows.h> -// 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<float *>(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<float>(); - - 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<float *>(npU.get_data()); - U_old = reinterpret_cast<float *>(npU_old.get_data()); - Dx = reinterpret_cast<float *>(npDx.get_data()); - Dy = reinterpret_cast<float *>(npDy.get_data()); - Bx = reinterpret_cast<float *>(npBx.get_data()); - By = reinterpret_cast<float *>(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<np::ndarray>(npU); - result.append<int>(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<float>(); - - 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<float *>(npU.get_data()); - U_old = reinterpret_cast<float *>(npU_old.get_data()); - Dx = reinterpret_cast<float *>(npDx.get_data()); - Dy = reinterpret_cast<float *>(npDy.get_data()); - Dz = reinterpret_cast<float *>(npDz.get_data()); - Bx = reinterpret_cast<float *>(npBx.get_data()); - By = reinterpret_cast<float *>(npBy.get_data()); - Bz = reinterpret_cast<float *>(npBz.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_seidel3D(U, A, Dx, Dy, Dz, Bx, By, Bz, dimX, dimY, dimZ, lambda, mu); - - if (methTV == 1) updDxDyDz_shrinkAniso3D(U, Dx, Dy, Dz, Bx, By, Bz, dimX, dimY, dimZ, lambda); - else updDxDyDz_shrinkIso3D(U, Dx, Dy, Dz, Bx, By, Bz, dimX, dimY, dimZ, lambda); - - updBxByBz3D(U, Dx, Dy, Dz, Bx, By, Bz, dimX, dimY, dimZ); - - /* calculate norm to terminate earlier */ - re = 0.0f; re1 = 0.0f; - for (j = 0; j<dimX*dimY*dimZ; j++) - { - re += pow(U[j] - U_old[j], 2); - re1 += pow(U[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; - } - /*printf("%f %i %i \n", re, ll, count); */ - re_old = re; - } - //printf("SB iterations stopped at iteration: %i\n", ll); - result.append<np::ndarray>(npU); - result.append<int>(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<float *>(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<float>(); - - - 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<np::ndarray>(npU); - - U = reinterpret_cast<float *>(npU.get_data()); - U_old = reinterpret_cast<float *>(npU_old.get_data()); - D1 = reinterpret_cast<float *>(npD1.get_data()); - D2 = reinterpret_cast<float *>(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<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; - - } /*end of iterations*/ - // printf("HO iterations stopped at iteration: %i\n", ll); - result.append<np::ndarray>(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<float>(); - - - 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<unsigned short>()); - Map = reinterpret_cast<unsigned short *>(npMap.get_data()); - if (switcher != 0) { - //Map = (unsigned short*)mxGetPr(plhs[1] = mxCreateNumericArray(3, dim_array, mxUINT16_CLASS, mxREAL)); - - Map = reinterpret_cast<unsigned short *>(npMap.get_data()); - } - - U = reinterpret_cast<float *>(npU.get_data()); - U_old = reinterpret_cast<float *>(npU_old.get_data()); - D1 = reinterpret_cast<float *>(npD1.get_data()); - D2 = reinterpret_cast<float *>(npD2.get_data()); - D3 = reinterpret_cast<float *>(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<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; - - } /*end of iterations*/ - //printf("HO iterations stopped at iteration: %i\n", ll); - result.append<np::ndarray>(npU); - if (switcher != 0) result.append<np::ndarray>(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<float *>(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<float>(); - - 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<float *>(npB.get_data()); - Ap = reinterpret_cast<float *>(npAp.get_data()); - Bp = reinterpret_cast<float *>(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<np::ndarray>(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<float>(); - - 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<float *>(npB.get_data()); - Ap = reinterpret_cast<float *>(npAp.get_data()); - Bp = reinterpret_cast<float *>(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<np::ndarray>(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<float *>(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<float>(); - - 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<float *>(npU.get_data()); - U_old = reinterpret_cast<float *>(npU_old.get_data()); - P1 = reinterpret_cast<float *>(npP1.get_data()); - P2 = reinterpret_cast<float *>(npP2.get_data()); - Q1 = reinterpret_cast<float *>(npQ1.get_data()); - Q2 = reinterpret_cast<float *>(npQ2.get_data()); - Q3 = reinterpret_cast<float *>(npQ3.get_data()); - V1 = reinterpret_cast<float *>(npV1.get_data()); - V1_old = reinterpret_cast<float *>(npV1_old.get_data()); - V2 = reinterpret_cast<float *>(npV2.get_data()); - V2_old = reinterpret_cast<float *>(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<np::ndarray>(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<uint8_t>(); - np::dtype dt2 = np::dtype::get_builtin<uint16_t>(); - - def("SplitBregman_TV", SplitBregman_TV); - def("LLT_model", LLT_model); - def("PatchBased_Regul", PatchBased_Regul); - def("TGV_PD", TGV_PD); -} |