From cfcc4be4413f65a0b9c4ef197687e3a167eff0e8 Mon Sep 17 00:00:00 2001 From: Daniil Kazantsev Date: Wed, 6 Mar 2019 23:34:55 +0000 Subject: cont1 --- src/Core/regularisers_CPU/FGP_TV_core.c | 16 +++++++------- src/Core/regularisers_CPU/ROF_TV_core.c | 39 +++++++++++++++++++++++++++------ src/Core/regularisers_CPU/ROF_TV_core.h | 10 +++++---- src/Python/ccpi/filters/regularisers.py | 5 +++-- src/Python/src/cpu_regularisers.pyx | 28 +++++++++++++---------- 5 files changed, 66 insertions(+), 32 deletions(-) (limited to 'src') diff --git a/src/Core/regularisers_CPU/FGP_TV_core.c b/src/Core/regularisers_CPU/FGP_TV_core.c index fc50f13..3c3cbd4 100644 --- a/src/Core/regularisers_CPU/FGP_TV_core.c +++ b/src/Core/regularisers_CPU/FGP_TV_core.c @@ -3,8 +3,8 @@ This work is part of the Core Imaging Library developed by Visual Analytics and Imaging System Group of the Science Technology Facilities Council, STFC -Copyright 2017 Daniil Kazantsev -Copyright 2017 Srikanth Nagella, Edoardo Pasca +Copyright 2019 Daniil Kazantsev +Copyright 2019 Srikanth Nagella, Edoardo Pasca Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. @@ -79,6 +79,11 @@ float TV_FGP_CPU_main(float *Input, float *Output, float *infovector, float lamb tkp1 = (1.0f + sqrt(1.0f + 4.0f*tk*tk))*0.5f; Rupd_func2D(P1, P1_prev, P2, P2_prev, R1, R2, tkp1, tk, DimTotal); + /*storing old values*/ + copyIm(P1, P1_prev, (long)(dimX), (long)(dimY), 1l); + copyIm(P2, P2_prev, (long)(dimX), (long)(dimY), 1l); + tk = tkp1; + /* check early stopping criteria */ if (epsil != 0.0f) { for(j=0; j 4) break; copyIm(Output, Output_prev, (long)(dimX), (long)(dimY), 1l); } - - /*storing old values*/ - copyIm(P1, P1_prev, (long)(dimX), (long)(dimY), 1l); - copyIm(P2, P2_prev, (long)(dimX), (long)(dimY), 1l); - tk = tkp1; - } + } if (epsil != 0.0f) free(Output_prev); free(P1); free(P2); free(P1_prev); free(P2_prev); free(R1); free(R2); } diff --git a/src/Core/regularisers_CPU/ROF_TV_core.c b/src/Core/regularisers_CPU/ROF_TV_core.c index 1858442..1b7cf76 100644 --- a/src/Core/regularisers_CPU/ROF_TV_core.c +++ b/src/Core/regularisers_CPU/ROF_TV_core.c @@ -19,7 +19,7 @@ #include "ROF_TV_core.h" -#define EPS 1.0e-12 +#define EPS 1.0e-15 #define MAX(x, y) (((x) > (y)) ? (x) : (y)) #define MIN(x, y) (((x) < (y)) ? (x) : (y)) @@ -37,20 +37,25 @@ int sign(float x) { * 2. lambda - regularization parameter [REQUIRED] * 3. tau - marching step for explicit scheme, ~1 is recommended [REQUIRED] * 4. Number of iterations, for explicit scheme >= 150 is recommended [REQUIRED] + * 5. eplsilon: tolerance constant * * Output: * [1] Regularized image/volume + * [2] Information vector which contains [iteration no., reached tolerance] * * This function is based on the paper by * [1] Rudin, Osher, Fatemi, "Nonlinear Total Variation based noise removal algorithms" */ /* Running iterations of TV-ROF function */ -float TV_ROF_CPU_main(float *Input, float *Output, float lambdaPar, int iterationsNumb, float tau, int dimX, int dimY, int dimZ) +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 *D1, *D2, *D3; + float *D1=NULL, *D2=NULL, *D3=NULL, *Output_prev=NULL; + float re, re1; + re = 0.0f; re1 = 0.0f; + int count = 0; int i; - long DimTotal; + long DimTotal,j; DimTotal = (long)(dimX*dimY*dimZ); D1 = calloc(DimTotal, sizeof(float)); @@ -59,6 +64,7 @@ float TV_ROF_CPU_main(float *Input, float *Output, float lambdaPar, int iteratio /* copy into output */ copyIm(Input, Output, (long)(dimX), (long)(dimY), (long)(dimZ)); + if (epsil != 0.0f) Output_prev = calloc(DimTotal, sizeof(float)); /* start TV iterations */ for(i=0; i < iterationsNumb; i++) { @@ -67,9 +73,28 @@ float TV_ROF_CPU_main(float *Input, float *Output, float lambdaPar, int iteratio 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)); + + /* check early stopping criteria */ + if (epsil != 0.0f) { + for(j=0; j 4) break; + copyIm(Output, Output_prev, (long)(dimX), (long)(dimY), (long)(dimZ)); + } } free(D1);free(D2); free(D3); - return *Output; + if (epsil != 0.0f) free(Output_prev); + + /*adding info into info_vector */ + infovector[0] = (float)(i); /*iterations number (if stopped earlier based on tolerance)*/ + infovector[1] = re; /* reached tolerance */ + + return 0; } /* calculate differences 1 */ @@ -264,7 +289,7 @@ 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*(2.0f*lambda*(dv1 + dv2 + dv3) - (B[index] - A[index])); + B[index] += tau*(lambda*(dv1 + dv2 + dv3) - (B[index] - A[index])); }}} } else { @@ -282,7 +307,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*(2.0f*lambda*(dv1 + dv2) - (B[index] - A[index])); + B[index] += tau*(lambda*(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 4e320e9..d6949fa 100644 --- a/src/Core/regularisers_CPU/ROF_TV_core.h +++ b/src/Core/regularisers_CPU/ROF_TV_core.h @@ -31,11 +31,13 @@ limitations under the License. * Input Parameters: * 1. Noisy image/volume [REQUIRED] * 2. lambda - regularization parameter [REQUIRED] - * 3. Number of iterations, for explicit scheme >= 150 is recommended [REQUIRED] - * 4. tau - marching step for explicit scheme, ~1 is recommended [REQUIRED] + * 3. tau - marching step for explicit scheme, ~1 is recommended [REQUIRED] + * 4. Number of iterations, for explicit scheme >= 150 is recommended [REQUIRED] + * 5. eplsilon: tolerance constant * * Output: * [1] Regularized image/volume + * [2] Information vector which contains [iteration no., reached tolerance] * * This function is based on the paper by * [1] Rudin, Osher, Fatemi, "Nonlinear Total Variation based noise removal algorithms" @@ -46,7 +48,7 @@ limitations under the License. #ifdef __cplusplus extern "C" { #endif -CCPI_EXPORT float TV_ROF_CPU_main(float *Input, float *Output, float lambdaPar, int iterationsNumb, float tau, int dimX, int dimY, int dimZ); +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_kernel(float *D1, float *D2, float *D3, float *B, float *A, float lambda, float tau, long dimX, long dimY, long dimZ); CCPI_EXPORT float D1_func(float *A, float *D1, long dimX, long dimY, long dimZ); @@ -54,4 +56,4 @@ 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); #ifdef __cplusplus } -#endif \ No newline at end of file +#endif diff --git a/src/Python/ccpi/filters/regularisers.py b/src/Python/ccpi/filters/regularisers.py index fb2c999..67f432b 100644 --- a/src/Python/ccpi/filters/regularisers.py +++ b/src/Python/ccpi/filters/regularisers.py @@ -11,12 +11,13 @@ except ImportError: from ccpi.filters.cpu_regularisers import NDF_INPAINT_CPU, NVM_INPAINT_CPU def ROF_TV(inputData, regularisation_parameter, iterations, - time_marching_parameter,device='cpu'): + time_marching_parameter,tolerance_param,device='cpu'): if device == 'cpu': return TV_ROF_CPU(inputData, regularisation_parameter, iterations, - time_marching_parameter) + time_marching_parameter, + tolerance_param) elif device == 'gpu' and gpu_enabled: return TV_ROF_GPU(inputData, regularisation_parameter, diff --git a/src/Python/src/cpu_regularisers.pyx b/src/Python/src/cpu_regularisers.pyx index 49cdf94..aeca141 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 lambdaPar, int iterationsNumb, float tau, int dimX, int dimY, int dimZ); +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_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 lambdaPar, int iterationsNumb, float epsil, int methodTV, int printM, int dimX, int dimY, int dimZ); cdef extern float LLT_ROF_CPU_main(float *Input, float *Output, float lambdaROF, float lambdaLLT, int iterationsNumb, float tau, int dimX, int dimY, int dimZ); @@ -37,32 +37,36 @@ cdef extern float TV_energy3D(float *U, float *U0, float *E_val, float lambdaPar #****************************************************************# #********************** Total-variation ROF *********************# #****************************************************************# -def TV_ROF_CPU(inputData, regularisation_parameter, iterationsNumb, marching_step_parameter): +def TV_ROF_CPU(inputData, regularisation_parameter, iterationsNumb, marching_step_parameter,tolerance_param): if inputData.ndim == 2: - return TV_ROF_2D(inputData, regularisation_parameter, iterationsNumb, marching_step_parameter) + return TV_ROF_2D(inputData, regularisation_parameter, iterationsNumb, marching_step_parameter,tolerance_param) elif inputData.ndim == 3: - return TV_ROF_3D(inputData, regularisation_parameter, iterationsNumb, marching_step_parameter) + 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, int iterationsNumb, - float marching_step_parameter): + float marching_step_parameter, + float tolerance_param): cdef long dims[2] dims[0] = inputData.shape[0] dims[1] = inputData.shape[1] 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], regularisation_parameter, iterationsNumb, marching_step_parameter, 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) - return outputData + return (outputData,infovec) def TV_ROF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, float regularisation_parameter, int iterationsNumb, - float marching_step_parameter): + float marching_step_parameter, + float tolerance_param): cdef long dims[3] dims[0] = inputData.shape[0] dims[1] = inputData.shape[1] @@ -70,11 +74,13 @@ def TV_ROF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 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], regularisation_parameter, iterationsNumb, marching_step_parameter, 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]) - return outputData + return (outputData,infovec) #****************************************************************# #********************** Total-variation FGP *********************# -- cgit v1.2.3