From 9c57162fef5822369c09b7b36b3637dae5e397c0 Mon Sep 17 00:00:00 2001 From: dkazanc Date: Tue, 16 Apr 2019 17:21:04 +0100 Subject: succesfull variable regulariser implementation --- src/Core/regularisers_CPU/ROF_TV_core.c | 18 ++++++++++-------- src/Core/regularisers_CPU/ROF_TV_core.h | 4 ++-- src/Python/src/cpu_regularisers.pyx | 31 ++++++++++++++++++++++--------- 3 files changed, 34 insertions(+), 19 deletions(-) diff --git a/src/Core/regularisers_CPU/ROF_TV_core.c b/src/Core/regularisers_CPU/ROF_TV_core.c index 6d23eef..955cabc 100644 --- a/src/Core/regularisers_CPU/ROF_TV_core.c +++ b/src/Core/regularisers_CPU/ROF_TV_core.c @@ -48,7 +48,7 @@ int sign(float x) { */ /* Running iterations of TV-ROF function */ -float TV_ROF_CPU_main(float *Input, float *Output, float *infovector, float lambdaPar, int iterationsNumb, float tau, float epsil, int dimX, int dimY, int dimZ) +float TV_ROF_CPU_main(float *Input, float *Output, float *infovector, float *lambdaPar, int lambda_is_arr, int iterationsNumb, float tau, float epsil, int dimX, int dimY, int dimZ) { float *D1=NULL, *D2=NULL, *D3=NULL, *Output_prev=NULL; float re, re1; @@ -74,7 +74,7 @@ float TV_ROF_CPU_main(float *Input, float *Output, float *infovector, float lamb D1_func(Output, D1, (long)(dimX), (long)(dimY), (long)(dimZ)); D2_func(Output, D2, (long)(dimX), (long)(dimY), (long)(dimZ)); if (dimZ > 1) D3_func(Output, D3, (long)(dimX), (long)(dimY), (long)(dimZ)); - TV_kernel(D1, D2, D3, Output, Input, lambdaPar, tau, (long)(dimX), (long)(dimY), (long)(dimZ)); + TV_kernel(D1, D2, D3, Output, Input, lambdaPar, lambda_is_arr, tau, (long)(dimX), (long)(dimY), (long)(dimZ)); /* check early stopping criteria */ if ((epsil != 0.0f) && (i % 5 == 0)) { @@ -267,17 +267,18 @@ float D3_func(float *A, float *D3, long dimX, long dimY, long dimZ) } /* calculate divergence */ -float TV_kernel(float *D1, float *D2, float *D3, float *B, float *A, float lambda, float tau, long dimX, long dimY, long dimZ) +float TV_kernel(float *D1, float *D2, float *D3, float *B, float *A, float *lambda, int lambda_is_arr, float tau, long dimX, long dimY, long dimZ) { - float dv1, dv2, dv3; + float dv1, dv2, dv3, lambda_val; long index,i,j,k,i1,i2,k1,j1,j2,k2; if (dimZ > 1) { -#pragma omp parallel for shared (D1, D2, D3, B, dimX, dimY, dimZ) private(index, i, j, k, i1, j1, k1, i2, j2, k2, dv1,dv2,dv3) +#pragma omp parallel for shared (D1, D2, D3, B, dimX, dimY, dimZ) private(index, i, j, k, i1, j1, k1, i2, j2, k2, dv1,dv2,dv3,lambda_val) for(j=0; j= dimX) i1 = i-1; i2 = i - 1; if (i2 < 0) i2 = i+1; @@ -291,14 +292,15 @@ float TV_kernel(float *D1, float *D2, float *D3, float *B, float *A, float lambd dv2 = D2[index] - D2[(dimX*dimY)*k + j*dimX+i2]; dv3 = D3[index] - D3[(dimX*dimY)*k2 + j*dimX+i]; - B[index] += tau*(lambda*(dv1 + dv2 + dv3) - (B[index] - A[index])); + B[index] += tau*(lambda_val*(dv1 + dv2 + dv3) - (B[index] - A[index])); }}} } else { -#pragma omp parallel for shared (D1, D2, B, dimX, dimY) private(index, i, j, i1, j1, i2, j2,dv1,dv2) +#pragma omp parallel for shared (D1, D2, B, dimX, dimY) private(index, i, j, i1, j1, i2, j2,dv1,dv2,lambda_val) for(j=0; j= dimX) i1 = i-1; i2 = i - 1; if (i2 < 0) i2 = i+1; @@ -309,7 +311,7 @@ float TV_kernel(float *D1, float *D2, float *D3, float *B, float *A, float lambd dv1 = D1[index] - D1[j2*dimX + i]; dv2 = D2[index] - D2[j*dimX + i2]; - B[index] += tau*(lambda*(dv1 + dv2) - (B[index] - A[index])); + B[index] += tau*(lambda_val*(dv1 + dv2) - (B[index] - A[index])); }} } return *B; diff --git a/src/Core/regularisers_CPU/ROF_TV_core.h b/src/Core/regularisers_CPU/ROF_TV_core.h index d6949fa..8b5e2e4 100644 --- a/src/Core/regularisers_CPU/ROF_TV_core.h +++ b/src/Core/regularisers_CPU/ROF_TV_core.h @@ -48,9 +48,9 @@ limitations under the License. #ifdef __cplusplus extern "C" { #endif -CCPI_EXPORT float TV_ROF_CPU_main(float *Input, float *Output, float *infovector, float lambdaPar, int iterationsNumb, float tau, float epsil, int dimX, int dimY, int dimZ); +CCPI_EXPORT float TV_ROF_CPU_main(float *Input, float *Output, float *infovector, float *lambdaPar, int lambda_is_arr, int iterationsNumb, float tau, float epsil, int dimX, int dimY, int dimZ); -CCPI_EXPORT float TV_kernel(float *D1, float *D2, float *D3, float *B, float *A, float lambda, float tau, long dimX, long dimY, long dimZ); +CCPI_EXPORT float TV_kernel(float *D1, float *D2, float *D3, float *B, float *A, float *lambda, int lambda_is_arr, float tau, long dimX, long dimY, long dimZ); CCPI_EXPORT float D1_func(float *A, float *D1, long dimX, long dimY, long dimZ); CCPI_EXPORT float D2_func(float *A, float *D2, long dimX, long dimY, long dimZ); CCPI_EXPORT float D3_func(float *A, float *D3, long dimX, long dimY, long dimZ); diff --git a/src/Python/src/cpu_regularisers.pyx b/src/Python/src/cpu_regularisers.pyx index a63ecfa..a2c4c32 100644 --- a/src/Python/src/cpu_regularisers.pyx +++ b/src/Python/src/cpu_regularisers.pyx @@ -18,7 +18,7 @@ import cython import numpy as np cimport numpy as np -cdef extern float TV_ROF_CPU_main(float *Input, float *Output, float *infovector, float lambdaPar, int iterationsNumb, float tau, float epsil, int dimX, int dimY, int dimZ); +cdef extern float TV_ROF_CPU_main(float *Input, float *Output, float *infovector, float *lambdaPar, int arrayscal, int iterationsNumb, float tau, float epsil, int dimX, int dimY, int dimZ); cdef extern float TV_FGP_CPU_main(float *Input, float *Output, float *infovector, float lambdaPar, int iterationsNumb, float epsil, int methodTV, int nonneg, int dimX, int dimY, int dimZ); cdef extern float SB_TV_CPU_main(float *Input, float *Output, float *infovector, float mu, int iter, float epsil, int methodTV, int dimX, int dimY, int dimZ); cdef extern float LLT_ROF_CPU_main(float *Input, float *Output, float *infovector, float lambdaROF, float lambdaLLT, int iterationsNumb, float tau, float epsil, int dimX, int dimY, int dimZ); @@ -45,26 +45,33 @@ def TV_ROF_CPU(inputData, regularisation_parameter, iterationsNumb, marching_ste return TV_ROF_3D(inputData, regularisation_parameter, iterationsNumb, marching_step_parameter,tolerance_param) def TV_ROF_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, - float regularisation_parameter, + regularisation_parameter, int iterationsNumb, float marching_step_parameter, float tolerance_param): cdef long dims[2] dims[0] = inputData.shape[0] dims[1] = inputData.shape[1] - + cdef float lambdareg + cdef np.ndarray[np.float32_t, ndim=2, mode="c"] reg cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \ np.zeros([dims[0],dims[1]], dtype='float32') cdef np.ndarray[np.float32_t, ndim=1, mode="c"] infovec = \ np.ones([2], dtype='float32') # Run ROF iterations for 2D data - TV_ROF_CPU_main(&inputData[0,0], &outputData[0,0], &infovec[0], regularisation_parameter, iterationsNumb, marching_step_parameter, tolerance_param, dims[1], dims[0], 1) - + # TV_ROF_CPU_main(&inputData[0,0], &outputData[0,0], &infovec[0], regularisation_parameter, iterationsNumb, marching_step_parameter, tolerance_param, dims[1], dims[0], 1) + # Run ROF iterations for 2D data + if isinstance (regularisation_parameter, np.ndarray): + reg = regularisation_parameter.copy() + TV_ROF_CPU_main(&inputData[0,0], &outputData[0,0], &infovec[0], ®[0,0], 1, iterationsNumb, marching_step_parameter, tolerance_param, dims[1], dims[0], 1) + else: # supposedly this would be a float + lambdareg = regularisation_parameter; + TV_ROF_CPU_main(&inputData[0,0], &outputData[0,0], &infovec[0], &lambdareg, 0, iterationsNumb, marching_step_parameter, tolerance_param, dims[1], dims[0], 1) return (outputData,infovec) def TV_ROF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, - float regularisation_parameter, + regularisation_parameter, int iterationsNumb, float marching_step_parameter, float tolerance_param): @@ -72,15 +79,21 @@ def TV_ROF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, dims[0] = inputData.shape[0] dims[1] = inputData.shape[1] dims[2] = inputData.shape[2] - + cdef float lambdareg + cdef np.ndarray[np.float32_t, ndim=3, mode="c"] reg cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \ np.zeros([dims[0],dims[1],dims[2]], dtype='float32') cdef np.ndarray[np.float32_t, ndim=1, mode="c"] infovec = \ np.ones([2], dtype='float32') # Run ROF iterations for 3D data - TV_ROF_CPU_main(&inputData[0,0,0], &outputData[0,0,0], &infovec[0], regularisation_parameter, iterationsNumb, marching_step_parameter, tolerance_param, dims[2], dims[1], dims[0]) - + #TV_ROF_CPU_main(&inputData[0,0,0], &outputData[0,0,0], &infovec[0], regularisation_parameter, iterationsNumb, marching_step_parameter, tolerance_param, dims[2], dims[1], dims[0]) + if isinstance (regularisation_parameter, np.ndarray): + reg = regularisation_parameter.copy() + TV_ROF_CPU_main(&inputData[0,0,0], &outputData[0,0,0], &infovec[0], ®[0,0,0], 1, iterationsNumb, marching_step_parameter, tolerance_param, dims[2], dims[1], dims[0]) + else: # supposedly this would be a float + lambdareg = regularisation_parameter + TV_ROF_CPU_main(&inputData[0,0,0], &outputData[0,0,0], &infovec[0], &lambdareg, 0, iterationsNumb, marching_step_parameter, tolerance_param, dims[2], dims[1], dims[0]) return (outputData,infovec) #****************************************************************# -- cgit v1.2.3