summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/Core/regularisers_CPU/ROF_TV_core.c18
-rw-r--r--src/Core/regularisers_CPU/ROF_TV_core.h4
-rw-r--r--src/Python/src/cpu_regularisers.pyx31
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<dimY; j++) {
for(i=0; i<dimX; i++) {
for(k=0; k<dimZ; k++) {
index = (dimX*dimY)*k + j*dimX+i;
+ lambda_val = *(lambda + index* lambda_is_arr);
/* symmetric boundary conditions (Neuman) */
i1 = i + 1; if (i1 >= 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<dimY; j++) {
for(i=0; i<dimX; i++) {
index = j*dimX+i;
+ lambda_val = *(lambda + index* lambda_is_arr);
/* symmetric boundary conditions (Neuman) */
i1 = i + 1; if (i1 >= 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], &reg[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], &reg[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)
#****************************************************************#
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-08 06:38:58 +0000