From cfcc4be4413f65a0b9c4ef197687e3a167eff0e8 Mon Sep 17 00:00:00 2001
From: Daniil Kazantsev <dkazanc@hotmail.com>
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<DimTotal; j++)
@@ -91,12 +96,7 @@ float TV_FGP_CPU_main(float *Input, float *Output, float *infovector, float lamb
if (count > 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<DimTotal; j++)
+ {
+ re += pow(Output[j] - Output_prev[j],2);
+ re1 += pow(Output[j],2);
+ }
+ re = sqrt(re)/sqrt(re1);
+ if (re < epsil) count++;
+ if (count > 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