From e6ab844cd82080dab3a5f257fc15f4c1a20b498c Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Fri, 26 Jan 2018 16:28:38 +0000 Subject: cython wrapper builds without errors TODO: test --- Wrappers/Python/setup.py | 1 - Wrappers/Python/src/fista_module_gpu.pyx | 322 ++++++++++++++++--------------- 2 files changed, 167 insertions(+), 156 deletions(-) diff --git a/Wrappers/Python/setup.py b/Wrappers/Python/setup.py index c535a34..951146a 100644 --- a/Wrappers/Python/setup.py +++ b/Wrappers/Python/setup.py @@ -61,7 +61,6 @@ setup( ext_modules = [Extension("ccpi.filters.gpu_regularizers", sources=[ os.path.join("." , "src", "fista_module_gpu.pyx" ), - #os.path.join("." , "src", "multiply.pyx" ) ], include_dirs=extra_include_dirs, library_dirs=extra_library_dirs, diff --git a/Wrappers/Python/src/fista_module_gpu.pyx b/Wrappers/Python/src/fista_module_gpu.pyx index f18181b..7658e36 100644 --- a/Wrappers/Python/src/fista_module_gpu.pyx +++ b/Wrappers/Python/src/fista_module_gpu.pyx @@ -15,7 +15,6 @@ Author: Edoardo Pasca """ import cython - import numpy as np cimport numpy as np @@ -62,9 +61,12 @@ def Diff4thHajiaboli2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, #A_L = (float*)mxGetData(mxCreateNumericMatrix(N, M, mxSINGLE_CLASS, mxREAL)); #B_L = (float*)mxGetData(mxCreateNumericMatrix(N, M, mxSINGLE_CLASS, mxREAL)); - A_L = np.zeros((N,M), dtype=np.float) - B_L = np.zeros((N,M), dtype=np.float) - B = np.zeros((dims[0],dims[1]), dtype=np.float) + cdef np.ndarray[np.float32_t, ndim=2, mode="c"] A_L = \ + np.zeros([N,M], dtype='float32') + cdef np.ndarray[np.float32_t, ndim=2, mode="c"] B_L = \ + np.zeros([N,M], dtype='float32') + cdef np.ndarray[np.float32_t, ndim=2, mode="c"] B = \ + np.zeros([dims[0],dims[1]], dtype='float32') #A = inputData # copy A to the bigger A_L with boundaries @@ -85,7 +87,7 @@ def Diff4thHajiaboli2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, edge_preserving_parameter, iterations , tau, - regularization_parameter) + regularization_parameter); # copy the processed B_L to a smaller B for i in range(N): for j in range(M): @@ -117,9 +119,12 @@ def Diff4thHajiaboli3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, #A_L = (float*)mxGetData(mxCreateNumericMatrix(N, M, mxSINGLE_CLASS, mxREAL)); #B_L = (float*)mxGetData(mxCreateNumericMatrix(N, M, mxSINGLE_CLASS, mxREAL)); - A_L = np.zeros((N,M,Z), dtype=np.float) - B_L = np.zeros((N,M,Z), dtype=np.float) - B = np.zeros((dims[0],dims[1],dims[2]), dtype=np.float) + cdef np.ndarray[np.float32_t, ndim=3, mode="c"] A_L = \ + np.zeros([N,M,Z], dtype='float32') + cdef np.ndarray[np.float32_t, ndim=3, mode="c"] B_L = \ + np.zeros([N,M,Z], dtype='float32') + cdef np.ndarray[np.float32_t, ndim=3, mode="c"] B = \ + np.zeros([dims[0],dims[1],dims[2]], dtype='float32') #A = inputData # copy A to the bigger A_L with boundaries @@ -142,7 +147,7 @@ def Diff4thHajiaboli3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, edge_preserving_parameter, iterations , tau, - regularization_parameter) + regularization_parameter); # copy the processed B_L to a smaller B for i in range(N): for j in range(M): @@ -155,149 +160,156 @@ def Diff4thHajiaboli3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, return B - -#def NML(inputData, -# regularization_parameter, -# iterations, -# edge_preserving_parameter): -# if inputData.ndim == 2: -# return NML2D(inputData, -# regularization_parameter, -# iterations, -# edge_preserving_parameter) -# elif inputData.ndim == 3: -# return NML3D(inputData, -# regularization_parameter, -# iterations, -# edge_preserving_parameter) -# -# #SearchW_real = (int) mxGetScalar(prhs[1]); /* the searching window ratio */ -# #SimilW = (int) mxGetScalar(prhs[2]); /* the similarity window ratio */ -# #h = (float) mxGetScalar(prhs[3]); /* parameter for the PB filtering function */ -# #lambda = (float) mxGetScalar(prhs[4]); -# -#def NML2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, -# SearchW_real, -# SimilW, -# h, -# lambdaf): -# N = inputData.shape[0] -# M = inputData.shape[1] -# Z = 0 -# numdims = inputData.ndim -# -# if h < 0: -# raise ValueError('Parameter for the PB filtering function must be > 0') -# -# SearchW = SearchW_real + 2*SimilW; -# -# SearchW_full = 2*SearchW + 1; #/* the full searching window size */ -# SimilW_full = 2*SimilW + 1; #/* the full similarity window size */ -# h2 = h*h; -# -# padXY = SearchW + 2*SimilW; #/* padding sizes */ -# newsizeX = N + 2*(padXY); #/* the X size of the padded array */ -# newsizeY = M + 2*(padXY); #/* the Y size of the padded array */ -# #newsizeZ = Z + 2*(padXY); #/* the Z size of the padded array */ -# -# #output -# B = np.zeros((N,M), dtype=np.float ) -# #/*allocating memory for the padded arrays */ -# -# Ap = np.zeros((newsizeX, newsizeY), dtype=np.float) -# Bp = np.zeros((newsizeX, newsizeY), dtype=np.float) -# Eucl_Vec = np.zeros((SimilW_full*SimilW_full), dtype=np.float) -# -# #/*Gaussian kernel */ -# cdef int count, i_n, j_n; -# cdef float val; -# count = 0 -# for i_n in range (-SimilW, SimilW +1): -# for j_n in range(-SimilW, SimilW +1): -# val = (float)(i_n*i_n + j_n*j_n)/(2*SimilW*SimilW) -# Eucl_Vec[count] = np.exp(-val) -# count = count + 1 -# -# #/*Perform padding of image A to the size of [newsizeX * newsizeY] */ -# switchpad_crop = 0; # /*padding*/ -# pad_crop(&inputData[0,0], &Ap[0,0], M, N, 0, newsizeY, newsizeX, 0, padXY, -# switchpad_crop); -# -# #/* Do PB regularization with the padded array */ -# NLM_GPU_kernel(&Ap[0,0], &Bp[0,0], &Eucl_Vec[0,0], newsizeY, newsizeX, 0, -# numdims, SearchW, SimilW, SearchW_real, -# h2, lambdaf); -# -# switchpad_crop = 1; #/*cropping*/ -# pad_crop(&Bp[0,0], &B[0,0], M, N, 0, newsizeX, newsizeY, 0, padXY, -# switchpad_crop) -# -# return B -# -#def NML3D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, -# SearchW_real, -# SimilW, -# h, -# lambdaf): -# N = inputData.shape[0] -# M = inputData.shape[1] -# Z = inputData.shape[2] -# numdims = inputData.ndim -# -# if h < 0: -# raise ValueError('Parameter for the PB filtering function must be > 0') -# -# SearchW = SearchW_real + 2*SimilW; -# -# SearchW_full = 2*SearchW + 1; #/* the full searching window size */ -# SimilW_full = 2*SimilW + 1; #/* the full similarity window size */ -# h2 = h*h; -# -# padXY = SearchW + 2*SimilW; #/* padding sizes */ -# newsizeX = N + 2*(padXY); #/* the X size of the padded array */ -# newsizeY = M + 2*(padXY); #/* the Y size of the padded array */ -# newsizeZ = Z + 2*(padXY); #/* the Z size of the padded array */ -# -# #output -# B = np.zeros((N,M,Z), dtype=np.float ) -# #/*allocating memory for the padded arrays */ -# -# Ap = np.zeros((newsizeX, newsizeY, newsizeZ), dtype=np.float) -# Bp = np.zeros((newsizeX, newsizeY, newsizeZ), dtype=np.float) -# Eucl_Vec = np.zeros((SimilW_full*SimilW_full*SimilW_full), dtype=np.float) -# -# #/*Gaussian kernel */ -# cdef int count, i_n, j_n, k_n; -# cdef float val; -# count = 0 -# for i_n in range (-SimilW, SimilW +1): -# for j_n in range(-SimilW, SimilW +1): -# for k_n in range(-SimilW, SimilW+1): -# val = (i_n*i_n + j_n*j_n + k_n*k_n)/(2*SimilW*SimilW*SimilW) -# Eucl_Vec[count] = np.exp(-val) -# count = count + 1 -# -# #/*Perform padding of image A to the size of [newsizeX * newsizeY] */ -# switchpad_crop = 0; # /*padding*/ -# pad_crop(&inputData[0,0,0], &Ap[0,0,0], -# M, N, Z, -# newsizeY, newsizeX, newsizeZ, -# padXY, -# switchpad_crop); -# -# #/* Do PB regularization with the padded array */ -# NLM_GPU_kernel(&Ap[0,0,0], &Bp[0,0,0], &Eucl_Vec[0,0,0], -# newsizeY, newsizeX, newsizeZ, -# numdims, SearchW, SimilW, SearchW_real, -# h2, lambdaf); -# -# switchpad_crop = 1; #/*cropping*/ -# pad_crop(&Bp[0,0,0], &B[0,0,0], -# M, N, Z, -# newsizeX, newsizeY, newsizeZ, -# padXY, -# switchpad_crop) -# -# return B -# -# \ No newline at end of file +def NML(inputData, + regularization_parameter, + iterations, + edge_preserving_parameter): + if inputData.ndim == 2: + return NML2D(inputData, + regularization_parameter, + iterations, + edge_preserving_parameter) + elif inputData.ndim == 3: + return NML3D(inputData, + regularization_parameter, + iterations, + edge_preserving_parameter) + + #SearchW_real = (int) mxGetScalar(prhs[1]); /* the searching window ratio */ + #SimilW = (int) mxGetScalar(prhs[2]); /* the similarity window ratio */ + #h = (float) mxGetScalar(prhs[3]); /* parameter for the PB filtering function */ + #lambda = (float) mxGetScalar(prhs[4]); + +def NML2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, + SearchW_real, + SimilW, + h, + lambdaf): + N = inputData.shape[0] + M = inputData.shape[1] + Z = 0 + numdims = inputData.ndim + + if h < 0: + raise ValueError('Parameter for the PB filtering function must be > 0') + + SearchW = SearchW_real + 2*SimilW; + + SearchW_full = 2*SearchW + 1; #/* the full searching window size */ + SimilW_full = 2*SimilW + 1; #/* the full similarity window size */ + h2 = h*h; + + padXY = SearchW + 2*SimilW; #/* padding sizes */ + newsizeX = N + 2*(padXY); #/* the X size of the padded array */ + newsizeY = M + 2*(padXY); #/* the Y size of the padded array */ + #newsizeZ = Z + 2*(padXY); #/* the Z size of the padded array */ + + #output + cdef np.ndarray[np.float32_t, ndim=2, mode="c"] B = \ + np.zeros([N,M], dtype='float32') + #/*allocating memory for the padded arrays */ + + cdef np.ndarray[np.float32_t, ndim=2, mode="c"] Ap = \ + np.zeros([newsizeX, newsizeY], dtype='float32') + cdef np.ndarray[np.float32_t, ndim=2, mode="c"] Bp = \ + np.zeros([newsizeX, newsizeY], dtype='float32') + cdef np.ndarray[np.float32_t, ndim=1, mode="c"] Eucl_Vec = \ + np.zeros([SimilW_full*SimilW_full], dtype='float32') + + #/*Gaussian kernel */ + cdef int count, i_n, j_n; + cdef float val; + count = 0 + for i_n in range (-SimilW, SimilW +1): + for j_n in range(-SimilW, SimilW +1): + val = (float)(i_n*i_n + j_n*j_n)/(2*SimilW*SimilW) + Eucl_Vec[count] = np.exp(-val) + count = count + 1 + + #/*Perform padding of image A to the size of [newsizeX * newsizeY] */ + switchpad_crop = 0; # /*padding*/ + pad_crop(&inputData[0,0], &Ap[0,0], M, N, 0, newsizeY, newsizeX, 0, padXY, + switchpad_crop); + + #/* Do PB regularization with the padded array */ + NLM_GPU_kernel(&Ap[0,0], &Bp[0,0], &Eucl_Vec[0], newsizeY, newsizeX, 0, + numdims, SearchW, SimilW, SearchW_real, + h2, lambdaf); + + switchpad_crop = 1; #/*cropping*/ + pad_crop(&Bp[0,0], &B[0,0], M, N, 0, newsizeX, newsizeY, 0, padXY, + switchpad_crop) + + return B + +def NML3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, + SearchW_real, + SimilW, + h, + lambdaf): + N = inputData.shape[0] + M = inputData.shape[1] + Z = inputData.shape[2] + numdims = inputData.ndim + + if h < 0: + raise ValueError('Parameter for the PB filtering function must be > 0') + + SearchW = SearchW_real + 2*SimilW; + + SearchW_full = 2*SearchW + 1; #/* the full searching window size */ + SimilW_full = 2*SimilW + 1; #/* the full similarity window size */ + h2 = h*h; + + padXY = SearchW + 2*SimilW; #/* padding sizes */ + newsizeX = N + 2*(padXY); #/* the X size of the padded array */ + newsizeY = M + 2*(padXY); #/* the Y size of the padded array */ + newsizeZ = Z + 2*(padXY); #/* the Z size of the padded array */ + + #output + cdef np.ndarray[np.float32_t, ndim=3, mode="c"] B = \ + np.zeros([N,M,Z], dtype='float32') + #/*allocating memory for the padded arrays */ + + cdef np.ndarray[np.float32_t, ndim=3, mode="c"] Ap = \ + np.zeros([newsizeX, newsizeY, newsizeZ], dtype='float32') + cdef np.ndarray[np.float32_t, ndim=3, mode="c"] Bp = \ + np.zeros([newsizeX, newsizeY, newsizeZ], dtype='float32') + cdef np.ndarray[np.float32_t, ndim=1, mode="c"] Eucl_Vec = \ + np.zeros([SimilW_full*SimilW_full*SimilW_full], + dtype='float32') + + + #/*Gaussian kernel */ + cdef int count, i_n, j_n, k_n; + cdef float val; + count = 0 + for i_n in range (-SimilW, SimilW +1): + for j_n in range(-SimilW, SimilW +1): + for k_n in range(-SimilW, SimilW+1): + val = (i_n*i_n + j_n*j_n + k_n*k_n)/(2*SimilW*SimilW*SimilW) + Eucl_Vec[count] = np.exp(-val) + count = count + 1 + + #/*Perform padding of image A to the size of [newsizeX * newsizeY] */ + switchpad_crop = 0; # /*padding*/ + pad_crop(&inputData[0,0,0], &Ap[0,0,0], + M, N, Z, + newsizeY, newsizeX, newsizeZ, + padXY, + switchpad_crop); + + #/* Do PB regularization with the padded array */ + NLM_GPU_kernel(&Ap[0,0,0], &Bp[0,0,0], &Eucl_Vec[0], + newsizeY, newsizeX, newsizeZ, + numdims, SearchW, SimilW, SearchW_real, + h2, lambdaf); + + switchpad_crop = 1; #/*cropping*/ + pad_crop(&Bp[0,0,0], &B[0,0,0], + M, N, Z, + newsizeX, newsizeY, newsizeZ, + padXY, + switchpad_crop) + + return B -- cgit v1.2.3