From 107eb18c28255c4c8dbdf8245ffb85fe6f7535d6 Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <edo.paskino@gmail.com>
Date: Mon, 29 Jan 2018 13:36:11 +0000
Subject: renamed fista_module_gpu to gpu_regularizers.pyx
Cleaned test_cpu_regularizers.py
---
Wrappers/Python/setup.py | 3 +-
Wrappers/Python/src/cpu_regularizers.pyx | 0
Wrappers/Python/src/fista_module_gpu.pyx | 315 ------------------
Wrappers/Python/src/gpu_regularizers.pyx | 315 ++++++++++++++++++
Wrappers/Python/test/test_cpu_regularizers.py | 445 ++++++++++++++++++++++++++
5 files changed, 762 insertions(+), 316 deletions(-)
create mode 100644 Wrappers/Python/src/cpu_regularizers.pyx
delete mode 100644 Wrappers/Python/src/fista_module_gpu.pyx
create mode 100644 Wrappers/Python/src/gpu_regularizers.pyx
create mode 100644 Wrappers/Python/test/test_cpu_regularizers.py
(limited to 'Wrappers')
diff --git a/Wrappers/Python/setup.py b/Wrappers/Python/setup.py
index 951146a..00c93fc 100644
--- a/Wrappers/Python/setup.py
+++ b/Wrappers/Python/setup.py
@@ -60,7 +60,7 @@ setup(
cmdclass = {'build_ext': build_ext},
ext_modules = [Extension("ccpi.filters.gpu_regularizers",
sources=[
- os.path.join("." , "src", "fista_module_gpu.pyx" ),
+ os.path.join("." , "src", "gpu_regularizers.pyx" ),
],
include_dirs=extra_include_dirs,
library_dirs=extra_library_dirs,
@@ -79,6 +79,7 @@ setup(
cmdclass = {'build_ext': build_ext},
ext_modules = [Extension("ccpi.filters.cpu_regularizers",
sources=[os.path.join("." , "src", "fista_module.cpp" ),
+ os.path.join("." , "src", "cpu_regularizers.pyx" )
# os.path.join("@CMAKE_SOURCE_DIR@" , "main_func" , "regularizers_CPU", "FGP_TV_core.c"),
# os.path.join("@CMAKE_SOURCE_DIR@" , "main_func" , "regularizers_CPU", "SplitBregman_TV_core.c"),
# os.path.join("@CMAKE_SOURCE_DIR@" , "main_func" , "regularizers_CPU", "LLT_model_core.c"),
diff --git a/Wrappers/Python/src/cpu_regularizers.pyx b/Wrappers/Python/src/cpu_regularizers.pyx
new file mode 100644
index 0000000..e69de29
diff --git a/Wrappers/Python/src/fista_module_gpu.pyx b/Wrappers/Python/src/fista_module_gpu.pyx
deleted file mode 100644
index 7658e36..0000000
--- a/Wrappers/Python/src/fista_module_gpu.pyx
+++ /dev/null
@@ -1,315 +0,0 @@
-# distutils: language=c++
-"""
-Copyright 2018 CCPi
-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.
-
-Author: Edoardo Pasca
-"""
-
-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 float pad_crop(float *A, float *Ap,
- int OldSizeX, int OldSizeY, int OldSizeZ,
- int NewSizeX, int NewSizeY, int NewSizeZ,
- int padXY, int switchpad_crop);
-
-def Diff4thHajiaboli(inputData,
- regularization_parameter,
- iterations,
- edge_preserving_parameter):
- if inputData.ndim == 2:
- return Diff4thHajiaboli2D(inputData,
- regularization_parameter,
- iterations,
- edge_preserving_parameter)
- elif inputData.ndim == 3:
- return Diff4thHajiaboli3D(inputData,
- regularization_parameter,
- iterations,
- edge_preserving_parameter)
-
-def Diff4thHajiaboli2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
- float regularization_parameter,
- int iterations,
- float edge_preserving_parameter):
-
- cdef long dims[2]
- dims[0] = inputData.shape[0]
- dims[1] = inputData.shape[1]
- N = dims[0] + 2;
- M = dims[1] + 2;
-
- #time step is sufficiently small for an explicit methods
- tau = 0.01
-
- #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_preserving_parameter,
- iterations ,
- tau,
- 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 regularization_parameter,
- int iterations,
- float edge_preserving_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
-
- # Time Step is small for an explicit methods
- tau = 0.0007;
-
- #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_preserving_parameter,
- iterations ,
- tau,
- 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 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
diff --git a/Wrappers/Python/src/gpu_regularizers.pyx b/Wrappers/Python/src/gpu_regularizers.pyx
new file mode 100644
index 0000000..7658e36
--- /dev/null
+++ b/Wrappers/Python/src/gpu_regularizers.pyx
@@ -0,0 +1,315 @@
+# distutils: language=c++
+"""
+Copyright 2018 CCPi
+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.
+
+Author: Edoardo Pasca
+"""
+
+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 float pad_crop(float *A, float *Ap,
+ int OldSizeX, int OldSizeY, int OldSizeZ,
+ int NewSizeX, int NewSizeY, int NewSizeZ,
+ int padXY, int switchpad_crop);
+
+def Diff4thHajiaboli(inputData,
+ regularization_parameter,
+ iterations,
+ edge_preserving_parameter):
+ if inputData.ndim == 2:
+ return Diff4thHajiaboli2D(inputData,
+ regularization_parameter,
+ iterations,
+ edge_preserving_parameter)
+ elif inputData.ndim == 3:
+ return Diff4thHajiaboli3D(inputData,
+ regularization_parameter,
+ iterations,
+ edge_preserving_parameter)
+
+def Diff4thHajiaboli2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
+ float regularization_parameter,
+ int iterations,
+ float edge_preserving_parameter):
+
+ cdef long dims[2]
+ dims[0] = inputData.shape[0]
+ dims[1] = inputData.shape[1]
+ N = dims[0] + 2;
+ M = dims[1] + 2;
+
+ #time step is sufficiently small for an explicit methods
+ tau = 0.01
+
+ #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_preserving_parameter,
+ iterations ,
+ tau,
+ 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 regularization_parameter,
+ int iterations,
+ float edge_preserving_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
+
+ # Time Step is small for an explicit methods
+ tau = 0.0007;
+
+ #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_preserving_parameter,
+ iterations ,
+ tau,
+ 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 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
diff --git a/Wrappers/Python/test/test_cpu_regularizers.py b/Wrappers/Python/test/test_cpu_regularizers.py
new file mode 100644
index 0000000..ac595e3
--- /dev/null
+++ b/Wrappers/Python/test/test_cpu_regularizers.py
@@ -0,0 +1,445 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Aug 4 11:10:05 2017
+
+@author: ofn77899
+"""
+
+#from ccpi.viewer.CILViewer2D import Converter
+#import vtk
+
+import matplotlib.pyplot as plt
+import numpy as np
+import os
+from enum import Enum
+import timeit
+#from PIL import Image
+#from Regularizer import Regularizer
+#from ccpi.filters.Regularizer import Regularizer
+from ccpi.filters.cpu_regularizers import SplitBregman_TV , FGP_TV , LLT_model, \
+ PatchBased_Regul , TGV_PD
+
+###############################################################################
+#https://stackoverflow.com/questions/13875989/comparing-image-in-url-to-image-in-filesystem-in-python/13884956#13884956
+#NRMSE a normalization of the root of the mean squared error
+#NRMSE is simply 1 - [RMSE / (maxval - minval)]. Where maxval is the maximum
+# intensity from the two images being compared, and respectively the same for
+# minval. RMSE is given by the square root of MSE:
+# sqrt[(sum(A - B) ** 2) / |A|],
+# where |A| means the number of elements in A. By doing this, the maximum value
+# given by RMSE is maxval.
+
+def nrmse(im1, im2):
+ a, b = im1.shape
+ rmse = np.sqrt(np.sum((im2 - im1) ** 2) / float(a * b))
+ max_val = max(np.max(im1), np.max(im2))
+ min_val = min(np.min(im1), np.min(im2))
+ return 1 - (rmse / (max_val - min_val))
+###############################################################################
+def printParametersToString(pars):
+ txt = r''
+ for key, value in pars.items():
+ if key== 'algorithm' :
+ txt += "{0} = {1}".format(key, value.__name__)
+ elif key == 'input':
+ txt += "{0} = {1}".format(key, np.shape(value))
+ else:
+ txt += "{0} = {1}".format(key, value)
+ txt += '\n'
+ return txt
+###############################################################################
+#
+# 2D Regularizers
+#
+###############################################################################
+#Example:
+# figure;
+# Im = double(imread('lena_gray_256.tif'))/255; % loading image
+# u0 = Im + .05*randn(size(Im)); u0(u0 < 0) = 0;
+# u = SplitBregman_TV(single(u0), 10, 30, 1e-04);
+
+# assumes the script is launched from the test directory
+filename = os.path.join(".." , ".." , ".." , "data" ,"lena_gray_512.tif")
+#filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-FISTA_reconstruction\data\lena_gray_512.tif"
+#filename = r"/home/ofn77899/Reconstruction/CCPi-FISTA_Reconstruction/data/lena_gray_512.tif"
+#filename = r'/home/algol/Documents/Python/STD_test_images/lena_gray_512.tif'
+
+Im = plt.imread(filename)
+Im = np.asarray(Im, dtype='float32')
+
+perc = 0.05
+u0 = Im + (perc* np.random.normal(size=np.shape(Im)))
+# map the u0 u0->u0>0
+f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1)
+u0 = f(u0).astype('float32')
+
+## plot
+fig = plt.figure()
+
+a=fig.add_subplot(2,3,1)
+a.set_title('noise')
+imgplot = plt.imshow(u0,cmap="gray")
+
+reg_output = []
+##############################################################################
+# Call regularizer
+
+####################### SplitBregman_TV #####################################
+# u = SplitBregman_TV(single(u0), 10, 30, 1e-04);
+
+start_time = timeit.default_timer()
+pars = {'algorithm' : SplitBregman_TV , \
+ 'input' : u0,
+ 'regularization_parameter':10. , \
+'number_of_iterations' :35 ,\
+'tolerance_constant':0.0001 , \
+'TV_penalty': 0
+}
+
+out = SplitBregman_TV (pars['input'], pars['regularization_parameter'],
+ pars['number_of_iterations'],
+ pars['tolerance_constant'],
+ pars['TV_penalty'])
+plotme = out[0]
+txtstr = printParametersToString(pars)
+txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
+print (txtstr)
+
+
+a=fig.add_subplot(2,3,2)
+
+
+# these are matplotlib.patch.Patch properties
+props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+# place a text box in upper left in axes coords
+a.text(0.05, 0.95, txtstr, transform=a.transAxes, fontsize=14,
+ verticalalignment='top', bbox=props)
+imgplot = plt.imshow(plotme,\
+ #cmap="gray"
+ )
+
+###################### FGP_TV #########################################
+# u = FGP_TV(single(u0), 0.05, 100, 1e-04);
+start_time = timeit.default_timer()
+pars = {'algorithm' : FGP_TV , \
+ 'input' : u0,
+ 'regularization_parameter':5e-4, \
+'number_of_iterations' :10 ,\
+'tolerance_constant':0.001,\
+'TV_penalty': 0
+}
+
+out = FGP_TV (pars['input'], pars['regularization_parameter'],
+ pars['number_of_iterations'],
+ pars['tolerance_constant'],
+ pars['TV_penalty'])
+plotme = out[0]
+txtstr = printParametersToString(pars)
+txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
+print (txtstr)
+
+
+a=fig.add_subplot(2,3,3)
+
+# these are matplotlib.patch.Patch properties
+props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+# place a text box in upper left in axes coords
+imgplot = plt.imshow(plotme, \
+ #cmap="gray"
+ )
+# place a text box in upper left in axes coords
+a.text(0.05, 0.95, txtstr, transform=a.transAxes, fontsize=14,
+ verticalalignment='top', bbox=props)
+
+###################### LLT_model #########################################
+# * u0 = Im + .03*randn(size(Im)); % adding noise
+# [Den] = LLT_model(single(u0), 10, 0.1, 1);
+#Den = LLT_model(single(u0), 25, 0.0003, 300, 0.0001, 0);
+#input, regularization_parameter , time_step, number_of_iterations,
+# tolerance_constant, restrictive_Z_smoothing=0
+
+start_time = timeit.default_timer()
+
+pars = {'algorithm': LLT_model , \
+ 'input' : u0,
+ 'regularization_parameter': 25,\
+ 'time_step':0.0003, \
+'number_of_iterations' :300,\
+'tolerance_constant':0.001,\
+'restrictive_Z_smoothing': 0
+}
+out = LLT_model(pars['input'],
+ pars['regularization_parameter'],
+ pars['time_step'] ,
+ pars['number_of_iterations'],
+ pars['tolerance_constant'],
+ pars['restrictive_Z_smoothing'] )
+
+plotme = out[0]
+txtstr = printParametersToString(pars)
+txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
+print (txtstr)
+a=fig.add_subplot(2,3,4)
+
+# these are matplotlib.patch.Patch properties
+props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+# place a text box in upper left in axes coords
+a.text(0.05, 0.95, txtstr, transform=a.transAxes, fontsize=14,
+ verticalalignment='top', bbox=props)
+imgplot = plt.imshow(plotme,\
+ #cmap="gray"
+ )
+
+
+# ###################### PatchBased_Regul #########################################
+# # Quick 2D denoising example in Matlab:
+# # Im = double(imread('lena_gray_256.tif'))/255; % loading image
+# # u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise
+# # ImDen = PB_Regul_CPU(single(u0), 3, 1, 0.08, 0.05);
+start_time = timeit.default_timer()
+
+pars = {'algorithm': PatchBased_Regul , \
+ 'input' : u0,
+ 'regularization_parameter': 0.05,\
+ 'searching_window_ratio':3, \
+'similarity_window_ratio':1,\
+'PB_filtering_parameter': 0.08
+}
+out = PatchBased_Regul(
+ pars['input'], pars['regularization_parameter'],
+ pars['searching_window_ratio'] ,
+ pars['similarity_window_ratio'] ,
+ pars['PB_filtering_parameter'])
+plotme = out[0]
+txtstr = printParametersToString(pars)
+txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
+print (txtstr)
+
+a=fig.add_subplot(2,3,5)
+
+
+# these are matplotlib.patch.Patch properties
+props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+# place a text box in upper left in axes coords
+a.text(0.05, 0.95, txtstr, transform=a.transAxes, fontsize=14,
+ verticalalignment='top', bbox=props)
+imgplot = plt.imshow(plotme #,cmap="gray"
+ )
+
+
+# ###################### TGV_PD #########################################
+# # Quick 2D denoising example in Matlab:
+# # Im = double(imread('lena_gray_256.tif'))/255; % loading image
+# # u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise
+# # u = PrimalDual_TGV(single(u0), 0.02, 1.3, 1, 550);
+
+start_time = timeit.default_timer()
+
+pars = {'algorithm': TGV_PD , \
+ 'input' : u0,\
+ 'regularization_parameter':0.05,\
+ 'first_order_term': 1.3,\
+ 'second_order_term': 1, \
+'number_of_iterations': 550
+}
+out = TGV_PD(pars['input'], pars['regularization_parameter'],
+ pars['first_order_term'] ,
+ pars['second_order_term'] ,
+ pars['number_of_iterations'])
+plotme = out[0]
+txtstr = printParametersToString(pars)
+txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
+print (txtstr)
+a=fig.add_subplot(2,3,6)
+
+# these are matplotlib.patch.Patch properties
+props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+# place a text box in upper left in axes coords
+a.text(0.05, 0.95, txtstr, transform=a.transAxes, fontsize=14,
+ verticalalignment='top', bbox=props)
+imgplot = plt.imshow(plotme #, cmap="gray")
+ )
+
+
+plt.show()
+
+################################################################################
+##
+## 3D Regularizers
+##
+################################################################################
+##Example:
+## figure;
+## Im = double(imread('lena_gray_256.tif'))/255; % loading image
+## u0 = Im + .05*randn(size(Im)); u0(u0 < 0) = 0;
+## u = SplitBregman_TV(single(u0), 10, 30, 1e-04);
+#
+##filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-Reconstruction\python\test\reconstruction_example.mha"
+#filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-Simpleflex\data\head.mha"
+#
+#reader = vtk.vtkMetaImageReader()
+#reader.SetFileName(os.path.normpath(filename))
+#reader.Update()
+##vtk returns 3D images, let's take just the one slice there is as 2D
+#Im = Converter.vtk2numpy(reader.GetOutput())
+#Im = Im.astype('float32')
+##imgplot = plt.imshow(Im)
+#perc = 0.05
+#u0 = Im + (perc* np.random.normal(size=np.shape(Im)))
+## map the u0 u0->u0>0
+#f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1)
+#u0 = f(u0).astype('float32')
+#converter = Converter.numpy2vtkImporter(u0, reader.GetOutput().GetSpacing(),
+# reader.GetOutput().GetOrigin())
+#converter.Update()
+#writer = vtk.vtkMetaImageWriter()
+#writer.SetInputData(converter.GetOutput())
+#writer.SetFileName(r"C:\Users\ofn77899\Documents\GitHub\CCPi-FISTA_reconstruction\data\noisy_head.mha")
+##writer.Write()
+#
+#
+### plot
+#fig3D = plt.figure()
+##a=fig.add_subplot(3,3,1)
+##a.set_title('Original')
+##imgplot = plt.imshow(Im)
+#sliceNo = 32
+#
+#a=fig3D.add_subplot(2,3,1)
+#a.set_title('noise')
+#imgplot = plt.imshow(u0.T[sliceNo])
+#
+#reg_output3d = []
+#
+###############################################################################
+## Call regularizer
+#
+######################## SplitBregman_TV #####################################
+## u = SplitBregman_TV(single(u0), 10, 30, 1e-04);
+#
+##reg = Regularizer(Regularizer.Algorithm.SplitBregman_TV)
+#
+##out = reg(input=u0, regularization_parameter=10., #number_of_iterations=30,
+## #tolerance_constant=1e-4,
+## TV_Penalty=Regularizer.TotalVariationPenalty.l1)
+#
+#out2 = Regularizer.SplitBregman_TV(input=u0, regularization_parameter=10., number_of_iterations=30,
+# tolerance_constant=1e-4,
+# TV_Penalty=Regularizer.TotalVariationPenalty.l1)
+#
+#
+#pars = out2[-2]
+#reg_output3d.append(out2)
+#
+#a=fig3D.add_subplot(2,3,2)
+#
+#
+#textstr = out2[-1]
+#
+#
+## these are matplotlib.patch.Patch properties
+#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+## place a text box in upper left in axes coords
+#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
+# verticalalignment='top', bbox=props)
+#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
+#
+####################### FGP_TV #########################################
+## u = FGP_TV(single(u0), 0.05, 100, 1e-04);
+#out2 = Regularizer.FGP_TV(input=u0, regularization_parameter=0.005,
+# number_of_iterations=200)
+#pars = out2[-2]
+#reg_output3d.append(out2)
+#
+#a=fig3D.add_subplot(2,3,2)
+#
+#
+#textstr = out2[-1]
+#
+#
+## these are matplotlib.patch.Patch properties
+#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+## place a text box in upper left in axes coords
+#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
+# verticalalignment='top', bbox=props)
+#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
+#
+####################### LLT_model #########################################
+## * u0 = Im + .03*randn(size(Im)); % adding noise
+## [Den] = LLT_model(single(u0), 10, 0.1, 1);
+##Den = LLT_model(single(u0), 25, 0.0003, 300, 0.0001, 0);
+##input, regularization_parameter , time_step, number_of_iterations,
+## tolerance_constant, restrictive_Z_smoothing=0
+#out2 = Regularizer.LLT_model(input=u0, regularization_parameter=25,
+# time_step=0.0003,
+# tolerance_constant=0.0001,
+# number_of_iterations=300)
+#pars = out2[-2]
+#reg_output3d.append(out2)
+#
+#a=fig3D.add_subplot(2,3,2)
+#
+#
+#textstr = out2[-1]
+#
+#
+## these are matplotlib.patch.Patch properties
+#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+## place a text box in upper left in axes coords
+#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
+# verticalalignment='top', bbox=props)
+#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
+#
+####################### PatchBased_Regul #########################################
+## Quick 2D denoising example in Matlab:
+## Im = double(imread('lena_gray_256.tif'))/255; % loading image
+## u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise
+## ImDen = PB_Regul_CPU(single(u0), 3, 1, 0.08, 0.05);
+#
+#out2 = Regularizer.PatchBased_Regul(input=u0, regularization_parameter=0.05,
+# searching_window_ratio=3,
+# similarity_window_ratio=1,
+# PB_filtering_parameter=0.08)
+#pars = out2[-2]
+#reg_output3d.append(out2)
+#
+#a=fig3D.add_subplot(2,3,2)
+#
+#
+#textstr = out2[-1]
+#
+#
+## these are matplotlib.patch.Patch properties
+#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+## place a text box in upper left in axes coords
+#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
+# verticalalignment='top', bbox=props)
+#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
+#
+
+###################### TGV_PD #########################################
+# Quick 2D denoising example in Matlab:
+# Im = double(imread('lena_gray_256.tif'))/255; % loading image
+# u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise
+# u = PrimalDual_TGV(single(u0), 0.02, 1.3, 1, 550);
+
+
+#out2 = Regularizer.TGV_PD(input=u0, regularization_parameter=0.05,
+# first_order_term=1.3,
+# second_order_term=1,
+# number_of_iterations=550)
+#pars = out2[-2]
+#reg_output3d.append(out2)
+#
+#a=fig3D.add_subplot(2,3,2)
+#
+#
+#textstr = out2[-1]
+#
+#
+## these are matplotlib.patch.Patch properties
+#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+## place a text box in upper left in axes coords
+#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14,
+# verticalalignment='top', bbox=props)
+#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo])
--
cgit v1.2.3