cython wrapper builds without errors

TODO: test

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