Cmake/Cython fixes to compile ROF-FGP

12 files changed, 51 insertions, 48 deletions

diff --git a/Core/CMakeLists.txt b/Core/CMakeLists.txt
index 3e3f89e..68b7111 100644
--- a/Core/CMakeLists.txt
+++ b/Core/CMakeLists.txt
@@ -129,7 +129,8 @@ if (CUDA_FOUND)
   CUDA_ADD_LIBRARY(cilregcuda SHARED
 	${CMAKE_CURRENT_SOURCE_DIR}/regularizers_GPU/Diffus_HO/Diff4th_GPU_kernel.cu
 	${CMAKE_CURRENT_SOURCE_DIR}/regularizers_GPU/NL_Regul/NLM_GPU_kernel.cu
-	${CMAKE_CURRENT_SOURCE_DIR}/regularizers_GPU/TV_ROF/TV_ROF_GPU.cu
+	${CMAKE_CURRENT_SOURCE_DIR}/regularizers_GPU/TV_ROF/TV_ROF_GPU_core.cu
+	${CMAKE_CURRENT_SOURCE_DIR}/regularizers_GPU/TV_FGP/TV_FGP_GPU_core.cu
   )
   if (UNIX)
     message ("I'd install into ${CMAKE_INSTALL_PREFIX}/lib")
diff --git a/Core/regularizers_CPU/FGP_TV_core.c b/Core/regularizers_CPU/FGP_TV_core.c
index 2f1439d..5365631 100644
--- a/Core/regularizers_CPU/FGP_TV_core.c
+++ b/Core/regularizers_CPU/FGP_TV_core.c
@@ -23,7 +23,7 @@ limitations under the License.
  *
  * Input Parameters:
  * 1. Noisy image/volume 
- * 2. lambda - regularization parameter 
+ * 2. lambdaPar - regularization parameter 
  * 3. Number of iterations
  * 4. eplsilon: tolerance constant 
  * 5. TV-type: methodTV - 'iso' (0) or 'l1' (1)
@@ -37,7 +37,7 @@ limitations under the License.
  * [1] Amir Beck and Marc Teboulle, "Fast Gradient-Based Algorithms for Constrained Total Variation Image Denoising and Deblurring Problems"
  */
  
-float TV_FGP_CPU(float *Input, float *Output, float lambda, int iter, float epsil, int methodTV, int nonneg, int printM, int dimX, int dimY, int dimZ)
+float TV_FGP_CPU_main(float *Input, float *Output, float lambdaPar, int iter, float epsil, int methodTV, int nonneg, int printM, int dimX, int dimY, int dimZ)
 {
 	int ll, j, DimTotal;
 	float re, re1;
@@ -62,13 +62,13 @@ float TV_FGP_CPU(float *Input, float *Output, float lambda, int iter, float epsi
         for(ll=0; ll<iter; ll++) {
             
             /* computing the gradient of the objective function */
-            Obj_func2D(Input, Output, R1, R2, lambda, dimX, dimY);
+            Obj_func2D(Input, Output, R1, R2, lambdaPar, dimX, dimY);
             
             /* apply nonnegativity */
             if (nonneg == 1) for(j=0; j<dimX*dimY; j++) {if (Output[j] < 0.0f) Output[j] = 0.0f;}            
             
             /*Taking a step towards minus of the gradient*/
-            Grad_func2D(P1, P2, Output, R1, R2, lambda, dimX, dimY);
+            Grad_func2D(P1, P2, Output, R1, R2, lambdaPar, dimX, dimY);
             
             /* projection step */
             Proj_func2D(P1, P2, methodTV, dimX, dimY);
@@ -117,13 +117,13 @@ float TV_FGP_CPU(float *Input, float *Output, float lambda, int iter, float epsi
         for(ll=0; ll<iter; ll++) {
             
             /* computing the gradient of the objective function */
-            Obj_func3D(Input, Output, R1, R2, R3, lambda, dimX, dimY, dimZ);
+            Obj_func3D(Input, Output, R1, R2, R3, lambdaPar, dimX, dimY, dimZ);
             
             /* apply nonnegativity */
             if (nonneg == 1) for(j=0; j<dimX*dimY*dimZ; j++) {if (Output[j] < 0.0f) Output[j] = 0.0f;}  
             
             /*Taking a step towards minus of the gradient*/
-            Grad_func3D(P1, P2, P3, Output, R1, R2, R3, lambda, dimX, dimY, dimZ);
+            Grad_func3D(P1, P2, P3, Output, R1, R2, R3, lambdaPar, dimX, dimY, dimZ);
             
             /* projection step */
             Proj_func3D(P1, P2, P3, methodTV, dimX, dimY, dimZ);
diff --git a/Core/regularizers_CPU/FGP_TV_core.h b/Core/regularizers_CPU/FGP_TV_core.h
index 98ceaec..360251d 100644
--- a/Core/regularizers_CPU/FGP_TV_core.h
+++ b/Core/regularizers_CPU/FGP_TV_core.h
@@ -47,7 +47,7 @@ limitations under the License.
 #ifdef __cplusplus
 extern "C" {
 #endif
-CCPI_EXPORT float TV_FGP_CPU(float *Input, float *Output, float lambda, int iter, float epsil, int methodTV, int nonneg, int printM, int dimX, int dimY, int dimZ);
+CCPI_EXPORT float TV_FGP_CPU_main(float *Input, float *Output, float lambdaPar, int iter, float epsil, int methodTV, int nonneg, int printM, int dimX, int dimY, int dimZ);
 
 CCPI_EXPORT float Obj_func2D(float *A, float *D, float *R1, float *R2, float lambda, int dimX, int dimY);
 CCPI_EXPORT float Grad_func2D(float *P1, float *P2, float *D, float *R1, float *R2, float lambda, int dimX, int dimY);
diff --git a/Core/regularizers_CPU/ROF_TV_core.c b/Core/regularizers_CPU/ROF_TV_core.c
index b2c6f00..a058f54 100644
--- a/Core/regularizers_CPU/ROF_TV_core.c
+++ b/Core/regularizers_CPU/ROF_TV_core.c
@@ -19,7 +19,7 @@
 
 #include "ROF_TV_core.h"
 
-#define EPS 1.0e-4
+#define EPS 1.0e-5
 #define MAX(x, y) (((x) > (y)) ? (x) : (y))
 #define MIN(x, y) (((x) < (y)) ? (x) : (y))
 
@@ -46,7 +46,7 @@ int sign(float x) {
  */
 
 /* Running iterations of TV-ROF function */
-float TV_ROF_CPU(float *Input, float *Output, int dimX, int dimY, int dimZ, int iterationsNumb, float tau, float lambda)
+float TV_ROF_CPU_main(float *Input, float *Output, int dimX, int dimY, int dimZ, int iterationsNumb, float tau, float lambda)
 {
     float *D1, *D2, *D3;
     int i, DimTotal;
diff --git a/Core/regularizers_CPU/ROF_TV_core.h b/Core/regularizers_CPU/ROF_TV_core.h
index b32d0d5..63c0419 100644
--- a/Core/regularizers_CPU/ROF_TV_core.h
+++ b/Core/regularizers_CPU/ROF_TV_core.h
@@ -47,7 +47,7 @@ limitations under the License.
 extern "C" {
 #endif
 CCPI_EXPORT float TV_kernel(float *D1, float *D2, float *D3, float *B, float *A, float lambda, float tau, int dimY, int dimX, int dimZ);
-CCPI_EXPORT float TV_ROF_CPU(float *Input, float *Output, int dimX, int dimY, int dimZ, int iterationsNumb, float tau, float lambda);
+CCPI_EXPORT float TV_ROF_CPU_main(float *Input, float *Output, int dimX, int dimY, int dimZ, int iterationsNumb, float tau, float lambda);
 CCPI_EXPORT float D1_func(float *A, float *D1, int dimY, int dimX, int dimZ);
 CCPI_EXPORT float D2_func(float *A, float *D2, int dimY, int dimX, int dimZ);
 CCPI_EXPORT float D3_func(float *A, float *D3, int dimY, int dimX, int dimZ);
diff --git a/Core/regularizers_GPU/TV_FGP/TV_FGP_GPU_core.cu b/Core/regularizers_GPU/TV_FGP/TV_FGP_GPU_core.cu
index 0533a85..cbe9b5e 100755
--- a/Core/regularizers_GPU/TV_FGP/TV_FGP_GPU_core.cu
+++ b/Core/regularizers_GPU/TV_FGP/TV_FGP_GPU_core.cu
@@ -25,7 +25,7 @@ limitations under the License.
  *
  * Input Parameters:
  * 1. Noisy image/volume 
- * 2. lambda - regularization parameter 
+ * 2. lambdaPar - regularization parameter 
  * 3. Number of iterations
  * 4. eplsilon: tolerance constant 
  * 5. TV-type: methodTV - 'iso' (0) or 'l1' (1)
@@ -341,7 +341,7 @@ __global__ void copy_kernel3D(float *Input, float* Output, int N, int M, int Z,
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 
 ////////////MAIN HOST FUNCTION ///////////////
-extern "C" void TV_FGP_GPU(float *Input, float *Output, float lambda, int iter, float epsil, int methodTV, int nonneg, int printM, int dimX, int dimY, int dimZ)
+extern "C" void TV_FGP_GPU(float *Input, float *Output, float lambdaPar, int iter, float epsil, int methodTV, int nonneg, int printM, int dimX, int dimY, int dimZ)
 {
     int deviceCount = -1; // number of devices
     cudaGetDeviceCount(&deviceCount);
@@ -383,13 +383,13 @@ extern "C" void TV_FGP_GPU(float *Input, float *Output, float lambda, int iter,
         cudaMemset(R2, 0, ImSize*sizeof(float));
 
         /********************** Run CUDA 2D kernel here ********************/    
-        multip = (1.0f/(8.0f*lambda));
+        multip = (1.0f/(8.0f*lambdaPar));
     
         /* The main kernel */
         for (i = 0; i < iter; i++) {
         
             /* computing the gradient of the objective function */
-            Obj_func2D_kernel<<<dimGrid,dimBlock>>>(d_input, d_update, R1, R2, dimX, dimY, ImSize, lambda);
+            Obj_func2D_kernel<<<dimGrid,dimBlock>>>(d_input, d_update, R1, R2, dimX, dimY, ImSize, lambdaPar);
             checkCudaErrors( cudaDeviceSynchronize() );
             checkCudaErrors(cudaPeekAtLastError() );
             
@@ -493,13 +493,13 @@ extern "C" void TV_FGP_GPU(float *Input, float *Output, float lambda, int iter,
             cudaMemset(R2, 0, ImSize*sizeof(float));
             cudaMemset(R3, 0, ImSize*sizeof(float));
             /********************** Run CUDA 3D kernel here ********************/    
-            multip = (1.0f/(8.0f*lambda));
+            multip = (1.0f/(8.0f*lambdaPar));
     
             /* The main kernel */
         for (i = 0; i < iter; i++) {
         
             /* computing the gradient of the objective function */
-            Obj_func3D_kernel<<<dimGrid,dimBlock>>>(d_input, d_update, R1, R2, R3, dimX, dimY, dimZ, ImSize, lambda);
+            Obj_func3D_kernel<<<dimGrid,dimBlock>>>(d_input, d_update, R1, R2, R3, dimX, dimY, dimZ, ImSize, lambdaPar);
             checkCudaErrors( cudaDeviceSynchronize() );
             checkCudaErrors(cudaPeekAtLastError() );
         
diff --git a/Core/regularizers_GPU/TV_FGP/TV_FGP_GPU_core.h b/Core/regularizers_GPU/TV_FGP/TV_FGP_GPU_core.h
index 15c7120..2b8fdcf 100755
--- a/Core/regularizers_GPU/TV_FGP/TV_FGP_GPU_core.h
+++ b/Core/regularizers_GPU/TV_FGP/TV_FGP_GPU_core.h
@@ -5,6 +5,6 @@
 #ifndef _TV_FGP_GPU_
 #define _TV_FGP_GPU_
 
-extern "C" void TV_FGP_GPU(float *Input, float *Output, float lambda, int iter, float epsil, int methodTV, int nonneg, int printM, int dimX, int dimY, int dimZ);   
+extern "C" void TV_FGP_GPU(float *Input, float *Output, float lambdaPar, int iter, float epsil, int methodTV, int nonneg, int printM, int dimX, int dimY, int dimZ);   
 
 #endif 
diff --git a/Wrappers/Python/demo/test_cpu_regularizers.py b/Wrappers/Python/demo/test_cpu_regularizers.py
index b08c418..53b8538 100644
--- a/Wrappers/Python/demo/test_cpu_regularizers.py
+++ b/Wrappers/Python/demo/test_cpu_regularizers.py
@@ -11,10 +11,9 @@ import numpy as np
 import os    
 from enum import Enum
 import timeit
-from ccpi.filters.cpu_regularizers_boost import SplitBregman_TV , FGP_TV ,\
-                                                 LLT_model, PatchBased_Regul ,\
+from ccpi.filters.cpu_regularizers_boost import SplitBregman_TV, LLT_model, PatchBased_Regul ,\
                                                  TGV_PD
-from ccpi.filters.cpu_regularizers_cython import ROF_TV
+from ccpi.filters.cpu_regularizers_cython import TV_ROF_CPU, TV_FGP_CPU
 
 ###############################################################################
 #https://stackoverflow.com/questions/13875989/comparing-image-in-url-to-image-in-filesystem-in-python/13884956#13884956
@@ -128,21 +127,25 @@ imgplot = plt.imshow(splitbregman,\
                      )
 
 ###################### FGP_TV #########################################
-# u = FGP_TV(single(u0), 0.05, 100, 1e-04);
+
 start_time = timeit.default_timer()
-pars = {'algorithm' : FGP_TV , \
+pars = {'algorithm' : TV_FGP_CPU , \
         'input' : u0,
         'regularization_parameter':0.05, \
         'number_of_iterations' :200 ,\
-        'tolerance_constant':1e-4,\
-        'TV_penalty': 0
+        'tolerance_constant':1e-5,\
+        'methodTV': 0 ,\
+        'nonneg': 0 ,\
+        'printingOut': 0
 }
 
-out = FGP_TV (pars['input'], 
+out = TV_FGP_CPU (pars['input'], 
               pars['regularization_parameter'],
               pars['number_of_iterations'],
               pars['tolerance_constant'], 
-              pars['TV_penalty'])  
+              pars['methodTV'],
+              pars['nonneg'],
+              pars['printingOut'])  
 
 fgp = out[0]
 rms = rmse(Im, fgp)
@@ -282,13 +285,13 @@ imgplot = plt.imshow(tgv, cmap="gray")
 
 start_time = timeit.default_timer()
 
-pars = {'algorithm': ROF_TV , \
+pars = {'algorithm': TV_ROF_CPU , \
         'input' : u0,\
         'regularization_parameter':0.04,\
         'marching_step': 0.0025,\
         'number_of_iterations': 300
         }
-rof = ROF_TV(pars['input'],
+rof = TV_ROF_CPU(pars['input'],
              pars['number_of_iterations'],
              pars['regularization_parameter'],
              pars['marching_step'] 
diff --git a/Wrappers/Python/src/cpu_regularizers.cpp b/Wrappers/Python/src/cpu_regularizers.cpp
index 43d5d11..b8156d5 100644
--- a/Wrappers/Python/src/cpu_regularizers.cpp
+++ b/Wrappers/Python/src/cpu_regularizers.cpp
@@ -1040,7 +1040,6 @@ BOOST_PYTHON_MODULE(cpu_regularizers_boost)
 	np::dtype dt2 = np::dtype::get_builtin<uint16_t>();
 
 	def("SplitBregman_TV", SplitBregman_TV);
-	def("FGP_TV", FGP_TV);
 	def("LLT_model", LLT_model);
 	def("PatchBased_Regul", PatchBased_Regul);
 	def("TGV_PD", TGV_PD);
diff --git a/Wrappers/Python/src/cpu_regularizers.pyx b/Wrappers/Python/src/cpu_regularizers.pyx
index b8089a8..448da31 100644
--- a/Wrappers/Python/src/cpu_regularizers.pyx
+++ b/Wrappers/Python/src/cpu_regularizers.pyx
@@ -18,8 +18,8 @@ import cython
 import numpy as np
 cimport numpy as np
 
-cdef extern float TV_ROF_CPU(float *Input, float *Output, int dimX, int dimY, int dimZ, int iterationsNumb, float tau, float flambda);
-cdef extern float TV_FGP_CPU(float *Input, float *Output, float lambda, int iter, float epsil, int methodTV, int nonneg, int printM, int dimX, int dimY, int dimZ);
+cdef extern float TV_ROF_CPU_main(float *Input, float *Output, int dimX, int dimY, int dimZ, int iterationsNumb, float tau, float flambda);
+cdef extern float TV_FGP_CPU_main(float *Input, float *Output, float lambdaPar, int iter, float epsil, int methodTV, int nonneg, int printM, int dimX, int dimY, int dimZ);								  
 
 # Can we use the same name here in "def" as the C function? 
 def TV_ROF_CPU(inputData, iterations, regularization_parameter, 
@@ -45,11 +45,10 @@ def TV_ROF_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
             np.zeros([dims[0],dims[1]], dtype='float32')
                    
     #/* Run ROF iterations for 2D data */
-    TV_ROF_CPU(&inputData[0,0], &B[0,0], dims[0], dims[1], 1, iterations,
+    TV_ROF_CPU_main(&inputData[0,0], &B[0,0], dims[0], dims[1], 1, iterations,
                marching_step_parameter, regularization_parameter)
     
-    return B
-        
+    return B        
             
 def TV_ROF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
                      int iterations,
@@ -65,7 +64,7 @@ def TV_ROF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
             np.zeros([dims[0],dims[1],dims[2]], dtype='float32')
            
     #/* Run ROF iterations for 3D data */
-    TV_FGP_CPU(&inputData[0,0,0], &B[0,0,0], dims[0], dims[1], dims[2], iterations,
+    TV_ROF_CPU_main(&inputData[0,0,0], &B[0,0,0], dims[0], dims[1], dims[2], iterations,
              marching_step_parameter, regularization_parameter)
 
     return B    
@@ -88,11 +87,11 @@ def TV_FGP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
     dims[0] = inputData.shape[0]
     dims[1] = inputData.shape[1]
     
-    cdef np.ndarray[np.float32_t, ndim=2, mode="c"] B = \
+    cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \
             np.zeros([dims[0],dims[1]], dtype='float32')
                    
     #/* Run ROF iterations for 2D data */
-    TV_FGP_CPU(&inputData[0,0], &B[0,0], regularization_parameter, 
+    TV_FGP_CPU_main(&inputData[0,0], &outputData[0,0], regularization_parameter, 
                        iterations, 
                        tolerance_param,
                        methodTV,
@@ -100,8 +99,7 @@ def TV_FGP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData,
                        printM,
                        dims[0], dims[1], 1)
     
-    return B
-        
+    return outputData        
             
 def TV_FGP_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 
  float regularization_parameter,
@@ -115,15 +113,16 @@ def TV_FGP_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData,
     dims[1] = inputData.shape[1]
     dims[2] = inputData.shape[2]
     
-    cdef np.ndarray[np.float32_t, ndim=3, mode="c"] B = \
+    cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \
             np.zeros([dims[0],dims[1],dims[2]], dtype='float32')
            
     #/* Run ROF iterations for 3D data */
-    TV_FGP_CPU(&inputData[0,0, 0], &B[0,0, 0], iterations, 
+    TV_FGP_CPU_main(&inputData[0,0,0], &outputData[0,0, 0], regularization_parameter,
+					   iterations, 
                        tolerance_param,
                        methodTV,
                        nonneg,
                        printM,
                        dims[0], dims[1], dims[2])
 
-    return B  
+    return outputData 
diff --git a/Wrappers/Python/src/gpu_regularizers.pyx b/Wrappers/Python/src/gpu_regularizers.pyx
index a14a20d..e99bfa7 100644
--- a/Wrappers/Python/src/gpu_regularizers.pyx
+++ b/Wrappers/Python/src/gpu_regularizers.pyx
@@ -26,7 +26,7 @@ cdef extern void NLM_GPU_kernel(float *A, float* B, float *Eucl_Vec,
                                 int SearchW, int SimilW, 
                                 int SearchW_real, float denh2, float lambdaf);
 cdef extern void TV_ROF_GPU(float* Input, float* Output, int N, int M, int Z, int iter, float tau, float lambdaf);
-cdef extern void TV_FGP_GPU(float *Input, float *Output, float lambda, int iter, float epsil, int methodTV, int nonneg, int printM, int N, int M, int Z);
+cdef extern void TV_FGP_GPU(float *Input, float *Output, float lambdaPar, int iter, float epsil, int methodTV, int nonneg, int printM, int N, int M, int Z);
 
 cdef extern float pad_crop(float *A, float *Ap, 
                            int OldSizeX, int OldSizeY, int OldSizeZ, 
diff --git a/Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py b/Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py
index 6344021..e162afa 100644
--- a/Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py
+++ b/Wrappers/Python/test/test_cpu_vs_gpu_regularizers.py
@@ -12,7 +12,7 @@ import matplotlib.pyplot as plt
 import numpy as np
 import os    
 import timeit
-from ccpi.filters.gpu_regularizers import Diff4thHajiaboli, NML, TV_ROF_GPU
+from ccpi.filters.gpu_regularizers import Diff4thHajiaboli, NML, GPU_ROF_TV
 from ccpi.filters.cpu_regularizers_cython import TV_ROF_CPU
 ###############################################################################
 def printParametersToString(pars):
@@ -56,11 +56,11 @@ imgplot = plt.imshow(u0,cmap="gray")
 
 
 # set parameters
-pars = {'algorithm': ROF_TV , \
+pars = {'algorithm': TV_ROF_CPU , \
         'input' : u0,\
         'regularization_parameter':0.04,\
         'time_marching_parameter': 0.0025,\
-        'number_of_iterations': 600
+        'number_of_iterations': 1200
         }
 print ("#################ROF TV CPU#####################")
 start_time = timeit.default_timer()
@@ -89,13 +89,14 @@ plt.title('{}'.format('CPU results'))
 
 print ("#################ROF TV GPU#####################")
 start_time = timeit.default_timer()
-rof_gpu = TV_ROF_GPU(pars['input'], 
+rof_gpu = GPU_ROF_TV(pars['input'], 
                      pars['number_of_iterations'], 
                      pars['time_marching_parameter'], 
                      pars['regularization_parameter'])
                      
 rms = rmse(Im, rof_gpu)
 pars['rmse'] = rms
+pars['algorithm'] = GPU_ROF_TV
 txtstr = printParametersToString(pars)
 txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time)
 print (txtstr)