cont1

7 files changed, 101 insertions, 62 deletions

diff --git a/demos/SoftwareX_supp/Demo_VolumeDenoise.py b/demos/SoftwareX_supp/Demo_VolumeDenoise.py
index 17cdf4d..6e7ea46 100644
--- a/demos/SoftwareX_supp/Demo_VolumeDenoise.py
+++ b/demos/SoftwareX_supp/Demo_VolumeDenoise.py
@@ -120,19 +120,21 @@ print("ROF-TV (gpu) ____ RMSE: {}, MMSIM: {}, run time: {} sec".format(RMSE_rof,
 #%%
 print ("#############FGP TV CPU####################")
 # set parameters
-pars = {'algorithm':FGP_TV, \
+pars =  {'algorithm' : FGP_TV, \
         'input' : phantom_noise,\
-        'regularisation_parameter':0.04,\
-        'number_of_iterations': 300,\
-        'time_marching_parameter': 0.0025,\
-        'tolerance_constant':1e-05,\
-        }
+        'regularisation_parameter':0.05, \
+        'number_of_iterations' :100 ,\
+        'tolerance_constant':1e-04,\
+        'methodTV': 0 ,\
+        'nonneg': 0}
 
 tic=timeit.default_timer()
-fgp_cpu3D = FGP_TV(pars['input'],
-             pars['regularisation_parameter'],
-             pars['number_of_iterations'],
-             pars['time_marching_parameter'],'cpu')
+(fgp_cpu3D, infoFGP) = FGP_TV(pars['input'], 
+              pars['regularisation_parameter'],
+              pars['number_of_iterations'],
+              pars['tolerance_constant'], 
+              pars['methodTV'],
+              pars['nonneg'],'cpu')
 toc=timeit.default_timer()
 
 Run_time_fgp = toc - tic
@@ -149,19 +151,21 @@ print("FGP-TV (cpu) ____ RMSE: {}, MMSIM: {}, run time: {} sec".format(RMSE_rof,
 #%%
 print ("#############FGP TV GPU####################")
 # set parameters
-pars = {'algorithm':FGP_TV, \
+pars = {'algorithm' : FGP_TV, \
         'input' : phantom_noise,\
-        'regularisation_parameter':0.04,\
-        'number_of_iterations': 300,\
-        'time_marching_parameter': 0.0025,\
-        'tolerance_constant':1e-05,\
-        }
+        'regularisation_parameter':0.05, \
+        'number_of_iterations' :80 ,\
+        'tolerance_constant':1e-04,\
+        'methodTV': 0 ,\
+        'nonneg': 0}
 
 tic=timeit.default_timer()
-fgp_gpu3D = FGP_TV(pars['input'],
-             pars['regularisation_parameter'],
-             pars['number_of_iterations'],
-             pars['time_marching_parameter'],'gpu')
+(fgp_gpu3D)  = FGP_TV(pars['input'], 
+              pars['regularisation_parameter'],
+              pars['number_of_iterations'],
+              pars['tolerance_constant'], 
+              pars['methodTV'],
+              pars['nonneg'],'gpu')
 toc=timeit.default_timer()
 
 Run_time_fgp = toc - tic
diff --git a/demos/demo_cpu_regularisers.py b/demos/demo_cpu_regularisers.py
index 8fa7022..32b97ce 100644
--- a/demos/demo_cpu_regularisers.py
+++ b/demos/demo_cpu_regularisers.py
@@ -32,7 +32,7 @@ def printParametersToString(pars):
 ###############################################################################
 
 #filename = os.path.join( "data" ,"lena_gray_512.tif")
-filename = "/home/kjy41806/Documents/SOFT/CCPi-Regularisation-Toolkit/test/lena_gray_512.tif"
+filename = "/home/algol/Documents/DEV/CCPi-Regularisation-Toolkit/test/lena_gray_512.tif"
 
 # read image
 Im = plt.imread(filename)
@@ -86,15 +86,17 @@ imgplot = plt.imshow(u0,cmap="gray")
 pars = {'algorithm': ROF_TV, \
         'input' : u0,\
         'regularisation_parameter':0.02,\
-        'number_of_iterations': 2000,\
-        'time_marching_parameter': 0.0025
-        }
+        'number_of_iterations': 5000,\
+        'time_marching_parameter': 0.000385,\
+        'tolerance_constant':1e-06}
+
 print ("#############ROF TV CPU####################")
 start_time = timeit.default_timer()
-rof_cpu = ROF_TV(pars['input'],
+(rof_cpu,info_vec) = ROF_TV(pars['input'],
              pars['regularisation_parameter'],
              pars['number_of_iterations'],
-             pars['time_marching_parameter'],'cpu')
+             pars['time_marching_parameter'],
+             pars['tolerance_constant'], 'cpu')
 
 Qtools = QualityTools(Im, rof_cpu)
 pars['rmse'] = Qtools.rmse()
@@ -127,12 +129,11 @@ imgplot = plt.imshow(u0,cmap="gray")
 # set parameters
 pars = {'algorithm' : FGP_TV, \
         'input' : u0,\
-        'regularisation_parameter':0.04, \
+        'regularisation_parameter':0.02, \
         'number_of_iterations' :800 ,\
-        'tolerance_constant':1e-07,\
+        'tolerance_constant':1e-06,\
         'methodTV': 0 ,\
-        'nonneg': 0
-        }
+        'nonneg': 0}
         
 print ("#############FGP TV CPU####################")
 start_time = timeit.default_timer()
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 *********************#