summaryrefslogtreecommitdiffstats
path: root/Wrappers/Python
diff options
context:
space:
mode:
Diffstat (limited to 'Wrappers/Python')
-rw-r--r--Wrappers/Python/src/cpu_regularizers.pyx12
-rw-r--r--Wrappers/Python/src/gpu_regularizers.pyx307
2 files changed, 17 insertions, 302 deletions
diff --git a/Wrappers/Python/src/cpu_regularizers.pyx b/Wrappers/Python/src/cpu_regularizers.pyx
index 60e8627..f993e54 100644
--- a/Wrappers/Python/src/cpu_regularizers.pyx
+++ b/Wrappers/Python/src/cpu_regularizers.pyx
@@ -21,6 +21,10 @@ cimport numpy as np
cdef extern float TV_ROF_CPU_main(float *Input, float *Output, float lambdaPar, int iterationsNumb, float tau, int dimX, int dimY, int dimZ);
cdef extern float TV_FGP_CPU_main(float *Input, float *Output, float lambdaPar, int iterationsNumb, float epsil, int methodTV, int nonneg, int printM, int dimX, int dimY, int dimZ);
+
+#****************************************************************#
+#********************** Total-variation ROF *********************#
+#****************************************************************#
def TV_ROF_CPU(inputData, regularization_parameter, iterationsNumb, marching_step_parameter):
if inputData.ndim == 2:
return TV_ROF_2D(inputData, regularization_parameter, iterationsNumb, marching_step_parameter)
@@ -58,8 +62,12 @@ def TV_ROF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
# Run ROF iterations for 3D data
TV_ROF_CPU_main(&inputData[0,0,0], &outputData[0,0,0], regularization_parameter, iterationsNumb, marching_step_parameter, dims[0], dims[1], dims[2])
- return outputData
-
+ return outputData
+
+#****************************************************************#
+#********************** Total-variation FGP *********************#
+#****************************************************************#
+#******** Total-variation Fast-Gradient-Projection (FGP)*********#
def TV_FGP_CPU(inputData, regularization_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, printM):
if inputData.ndim == 2:
return TV_FGP_2D(inputData, regularization_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, printM)
diff --git a/Wrappers/Python/src/gpu_regularizers.pyx b/Wrappers/Python/src/gpu_regularizers.pyx
index f96156a..a44bd1d 100644
--- a/Wrappers/Python/src/gpu_regularizers.pyx
+++ b/Wrappers/Python/src/gpu_regularizers.pyx
@@ -18,59 +18,9 @@ import cython
import numpy as np
cimport numpy as np
-
-cdef extern void Diff4th_GPU_kernel(float* A, float* B, int N, int M, int Z,
- float sigma, int iter, float tau, float lambdaf);
-cdef extern void NLM_GPU_kernel(float *A, float* B, float *Eucl_Vec,
- int N, int M, int Z, int dimension,
- int SearchW, int SimilW,
- int SearchW_real, float denh2, float lambdaf);
cdef extern void TV_ROF_GPU_main(float* Input, float* Output, float lambdaPar, int iter, float tau, int N, int M, int Z);
cdef extern void TV_FGP_GPU_main(float *Input, float *Output, float lambdaPar, int iter, float epsil, int methodTV, int nonneg, int printM, int N, int M, int Z);
-# padding function
-cdef extern float pad_crop(float *A, float *Ap,
- int OldSizeX, int OldSizeY, int OldSizeZ,
- int NewSizeX, int NewSizeY, int NewSizeZ,
- int padXY, int switchpad_crop);
-
-#Diffusion 4th order regularizer
-def Diff4thHajiaboli(inputData,
- edge_preserv_parameter,
- iterations,
- time_marching_parameter,
- regularization_parameter):
- if inputData.ndim == 2:
- return Diff4thHajiaboli2D(inputData,
- edge_preserv_parameter,
- iterations,
- time_marching_parameter,
- regularization_parameter)
- elif inputData.ndim == 3:
- return Diff4thHajiaboli3D(inputData,
- edge_preserv_parameter,
- iterations,
- time_marching_parameter,
- regularization_parameter)
-# patch-based nonlocal regularization
-def NML(inputData,
- SearchW_real,
- SimilW,
- h,
- lambdaf):
- if inputData.ndim == 2:
- return NML2D(inputData,
- SearchW_real,
- SimilW,
- h,
- lambdaf)
- elif inputData.ndim == 3:
- return NML3D(inputData,
- SearchW_real,
- SimilW,
- h,
- lambdaf)
-
# Total-variation Rudin-Osher-Fatemi (ROF)
def TV_ROF_GPU(inputData,
regularization_parameter,
@@ -111,255 +61,10 @@ def TV_FGP_GPU(inputData,
methodTV,
nonneg,
printM)
+
+#****************************************************************#
+#********************** Total-variation ROF *********************#
#****************************************************************#
-def Diff4thHajiaboli2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
- float edge_preserv_parameter,
- int iterations,
- float time_marching_parameter,
- float regularization_parameter):
-
- cdef long dims[2]
- dims[0] = inputData.shape[0]
- dims[1] = inputData.shape[1]
- N = dims[0] + 2;
- M = dims[1] + 2;
-
- #A_L = (float*)mxGetData(mxCreateNumericMatrix(N, M, mxSINGLE_CLASS, mxREAL));
- #B_L = (float*)mxGetData(mxCreateNumericMatrix(N, M, mxSINGLE_CLASS, mxREAL));
- cdef np.ndarray[np.float32_t, ndim=2, mode="c"] A_L = \
- np.zeros([N,M], dtype='float32')
- cdef np.ndarray[np.float32_t, ndim=2, mode="c"] B_L = \
- np.zeros([N,M], dtype='float32')
- cdef np.ndarray[np.float32_t, ndim=2, mode="c"] B = \
- np.zeros([dims[0],dims[1]], dtype='float32')
- #A = inputData
-
- # copy A to the bigger A_L with boundaries
- #pragma omp parallel for shared(A_L, A) private(i,j)
- cdef int i, j;
- for i in range(N):
- for j in range(M):
- if (((i > 0) and (i < N-1)) and ((j > 0) and (j < M-1))):
- #A_L[i*M+j] = A[(i-1)*(dims[1])+(j-1)]
- A_L[i][j] = inputData[i-1][j-1]
-
- # Running CUDA code here
- #Diff4th_GPU_kernel(A_L, B_L, N, M, Z, (float)sigma, iter, (float)tau, lambda);
- Diff4th_GPU_kernel(
- #<float*> A_L.data, <float*> B_L.data,
- &A_L[0,0], &B_L[0,0],
- N, M, 0,
- edge_preserv_parameter,
- iterations ,
- time_marching_parameter,
- regularization_parameter);
- # copy the processed B_L to a smaller B
- for i in range(N):
- for j in range(M):
- if (((i > 0) and (i < N-1)) and ((j > 0) and (j < M-1))):
- B[i-1][j-1] = B_L[i][j]
- ##pragma omp parallel for shared(B_L, B) private(i,j)
- #for (i=0; i < N; i++) {
- # for (j=0; j < M; j++) {
- # if (((i > 0) && (i < N-1)) && ((j > 0) && (j < M-1))) B[(i-1)*(dims[1])+(j-1)] = B_L[i*M+j];
- # }}
-
- return B
-
-def Diff4thHajiaboli3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
- float edge_preserv_parameter,
- int iterations,
- float time_marching_parameter,
- float regularization_parameter):
-
- cdef long dims[3]
- dims[0] = inputData.shape[0]
- dims[1] = inputData.shape[1]
- dims[2] = inputData.shape[2]
- N = dims[0] + 2
- M = dims[1] + 2
- Z = dims[2] + 2
-
-
- #A_L = (float*)mxGetData(mxCreateNumericMatrix(N, M, mxSINGLE_CLASS, mxREAL));
- #B_L = (float*)mxGetData(mxCreateNumericMatrix(N, M, mxSINGLE_CLASS, mxREAL));
- cdef np.ndarray[np.float32_t, ndim=3, mode="c"] A_L = \
- np.zeros([N,M,Z], dtype='float32')
- cdef np.ndarray[np.float32_t, ndim=3, mode="c"] B_L = \
- np.zeros([N,M,Z], dtype='float32')
- cdef np.ndarray[np.float32_t, ndim=3, mode="c"] B = \
- np.zeros([dims[0],dims[1],dims[2]], dtype='float32')
- #A = inputData
-
- # copy A to the bigger A_L with boundaries
- #pragma omp parallel for shared(A_L, A) private(i,j)
- cdef int i, j, k;
- for i in range(N):
- for j in range(M):
- for k in range(Z):
- if (((i > 0) and (i < N-1)) and \
- ((j > 0) and (j < M-1)) and \
- ((k > 0) and (k < Z-1))):
- A_L[i][j][k] = inputData[i-1][j-1][k-1];
-
- # Running CUDA code here
- #Diff4th_GPU_kernel(A_L, B_L, N, M, Z, (float)sigma, iter, (float)tau, lambda);
- Diff4th_GPU_kernel(
- #<float*> A_L.data, <float*> B_L.data,
- &A_L[0,0,0], &B_L[0,0,0],
- N, M, Z,
- edge_preserv_parameter,
- iterations ,
- time_marching_parameter,
- regularization_parameter);
- # copy the processed B_L to a smaller B
- for i in range(N):
- for j in range(M):
- for k in range(Z):
- if (((i > 0) and (i < N-1)) and \
- ((j > 0) and (j < M-1)) and \
- ((k > 0) and (k < Z-1))):
- B[i-1][j-1][k-1] = B_L[i][j][k]
-
-
- return B
-
-
-def 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
-
-# Total-variation ROF
def ROFTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
float regularization_parameter,
int iterations,
@@ -404,8 +109,10 @@ def ROFTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
dims[0], dims[1], dims[2]);
return outputData
-
-# Total-variation Fast-Gradient-Projection (FGP)
+#****************************************************************#
+#********************** Total-variation FGP *********************#
+#****************************************************************#
+#******** Total-variation Fast-Gradient-Projection (FGP)*********#
def FGPTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
float regularization_parameter,
int iterations,
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-27 08:37:55 +0000