succesfull variable regulariser implementation

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)
 
 #****************************************************************#