bugs in cython files

5 files changed, 47 insertions, 49 deletions

diff --git a/Core/regularisers_CPU/TNV_core.c b/Core/regularisers_CPU/TNV_core.c
index 8e61869..c51d6cd 100755
--- a/Core/regularisers_CPU/TNV_core.c
+++ b/Core/regularisers_CPU/TNV_core.c
@@ -251,18 +251,19 @@ float gradient(float *u_upd, float *gradx_upd, float *grady_upd, int dimX, int d
     // Compute discrete gradient using forward differences
 #pragma omp parallel for shared(gradx_upd,grady_upd,u_upd) private(i, j, k, l)
     for(k = 0; k < dimZ; k++)   {
-        for(j = 0; j < dimX; j++)   {
-            l = j * dimY;
-            for(i = 0; i < dimY; i++)   {
+        for(j = 0; j < dimY; j++)   {
+            l = j * dimX;           
+            for(i = 0; i < dimX; i++)   {
                 // Derivatives in the x-direction
-                if(i != dimY-1)
+                if(i != dimX-1)
                     gradx_upd[(dimX*dimY)*k + i+l] = u_upd[(dimX*dimY)*k + i+1+l] - u_upd[(dimX*dimY)*k + i+l];
                 else
                     gradx_upd[(dimX*dimY)*k + i+l] = 0.0f;
                 
                 // Derivatives in the y-direction
-                if(j != dimX-1)
-                    grady_upd[(dimX*dimY)*k + i+l] = u_upd[(dimX*dimY)*k + i+dimY+l] -u_upd[(dimX*dimY)*k + i+l];
+                if(j != dimY-1)
+                    //grady_upd[(dimX*dimY)*k + i+l] = u_upd[(dimX*dimY)*k + i+dimY+l] -u_upd[(dimX*dimY)*k + i+l];
+                    grady_upd[(dimX*dimY)*k + i+l] = u_upd[(dimX*dimY)*k + i+(j+1)*dimX] -u_upd[(dimX*dimY)*k + i+l];
                 else
                     grady_upd[(dimX*dimY)*k + i+l] = 0.0f;
             }}}
@@ -301,8 +302,8 @@ float proxF(float *gx, float *gy, float *vx, float *vy, float sigma, int p, int
             
             for(k = 0; k < dimZ; k++)
             {
-                valuex = vx[(dimX*dimY)*k + (i)*dimY + (j)];
-                valuey = vy[(dimX*dimY)*k + (i)*dimY + (j)];
+                valuex = vx[(dimX*dimY)*k + (j)*dimX + (i)];
+                valuey = vy[(dimX*dimY)*k + (j)*dimX + (i)];
                 
                 M1 += (valuex * valuex);
                 M2 += (valuex * valuey);
@@ -412,9 +413,9 @@ float proxF(float *gx, float *gy, float *vx, float *vy, float sigma, int p, int
             
             for(k = 0; k < dimZ; k++)
             {
-                gx[(dimX*dimY)*k + (i)*dimY + (j)] = vx[(dimX*dimY)*k + (i)*dimY + (j)] * t1 + vy[(dimX*dimY)*k + (i)*dimY + (j)] * t2;
-                gy[(dimX*dimY)*k + (i)*dimY + (j)] = vx[(dimX*dimY)*k + (i)*dimY + (j)] * t2 + vy[(dimX*dimY)*k + (i)*dimY + (j)] * t3;
-            }
+                gx[(dimX*dimY)*k + j*dimX + i] = vx[(dimX*dimY)*k + j*dimX + i] * t1 + vy[(dimX*dimY)*k + j*dimX + i] * t2;
+                gy[(dimX*dimY)*k + j*dimX + i] = vx[(dimX*dimY)*k + j*dimX + i] * t2 + vy[(dimX*dimY)*k + j*dimX + i] * t3;
+            }           
             
             // Delete allocated memory
             free(proj);
@@ -428,20 +429,21 @@ float divergence(float *qx_upd, float *qy_upd, float *div_upd, int dimX, int dim
     int i, j, k, l;
 #pragma omp parallel for shared(qx_upd,qy_upd,div_upd) private(i, j, k, l)
     for(k = 0; k < dimZ; k++)   {
-        for(j = 0; j < dimX; j++)   {
-            l = j * dimY;
-            for(i = 0; i < dimY; i++)   {
-                if(i != dimY-1)
+        for(j = 0; j < dimY; j++)   {
+            l = j * dimX;            
+            for(i = 0; i < dimX; i++)   {
+                if(i != dimX-1)
                 {
                     // ux[k][i+l] = u[k][i+1+l] - u[k][i+l]
                     div_upd[(dimX*dimY)*k + i+1+l] -= qx_upd[(dimX*dimY)*k + i+l];
                     div_upd[(dimX*dimY)*k + i+l] += qx_upd[(dimX*dimY)*k + i+l];
                 }
                 
-                if(j != dimX-1)
+                if(j != dimY-1)
                 {
                     // uy[k][i+l] = u[k][i+width+l] - u[k][i+l]
-                    div_upd[(dimX*dimY)*k + i+dimY+l] -= qy_upd[(dimX*dimY)*k + i+l];
+                    //div_upd[(dimX*dimY)*k + i+dimY+l] -= qy_upd[(dimX*dimY)*k + i+l];
+                    div_upd[(dimX*dimY)*k + i+(j+1)*dimX] -= qy_upd[(dimX*dimY)*k + i+l];                    
                     div_upd[(dimX*dimY)*k + i+l] += qy_upd[(dimX*dimY)*k + i+l];
                 }
             }
diff --git a/Wrappers/Python/demos/demo_cpu_inpainters.py b/Wrappers/Python/demos/demo_cpu_inpainters.py
index 348d235..7f452c1 100644
--- a/Wrappers/Python/demos/demo_cpu_inpainters.py
+++ b/Wrappers/Python/demos/demo_cpu_inpainters.py
@@ -72,7 +72,7 @@ pars = {'algorithm' : NDF_INP, \
         'maskData' : mask,\
         'regularisation_parameter':5000,\
         'edge_parameter':0,\
-        'number_of_iterations' :1000 ,\
+        'number_of_iterations' :5000 ,\
         'time_marching_parameter':0.000075,\
         'penalty_type':0
         }
diff --git a/Wrappers/Python/demos/demo_cpu_regularisers.py b/Wrappers/Python/demos/demo_cpu_regularisers.py
index f803870..986e3e9 100644
--- a/Wrappers/Python/demos/demo_cpu_regularisers.py
+++ b/Wrappers/Python/demos/demo_cpu_regularisers.py
@@ -44,29 +44,30 @@ u0 = Im + np.random.normal(loc = 0 ,
 u_ref = Im + np.random.normal(loc = 0 ,
                                   scale = 0.01 * Im , 
                                   size = np.shape(Im))
- 
+(N,M) = np.shape(u0)
 # map the u0 u0->u0>0
 # f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1)
 u0 = u0.astype('float32')
 u_ref = u_ref.astype('float32')
 
 # change dims to check that modules work with non-squared images
-(N,M) = np.shape(u0)
-u_ref2 = np.zeros([N,M-100],dtype='float32')
-u_ref2[:,0:M-100] = u_ref[:,0:M-100]
+"""
+M = M-100
+u_ref2 = np.zeros([N,M],dtype='float32')
+u_ref2[:,0:M] = u_ref[:,0:M]
 u_ref = u_ref2
 del u_ref2
 
-u02 = np.zeros([N,M-100],dtype='float32')
-u02[:,0:M-100] = u0[:,0:M-100]
+u02 = np.zeros([N,M],dtype='float32')
+u02[:,0:M] = u0[:,0:M]
 u0 = u02
 del u02
 
-Im2 = np.zeros([N,M-100],dtype='float32')
-Im2[:,0:M-100] = Im[:,0:M-100]
+Im2 = np.zeros([N,M],dtype='float32')
+Im2[:,0:M] = Im[:,0:M]
 Im = Im2
 del Im2
-
+"""
 #%%
 print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
 print ("_______________ROF-TV (2D)_________________")
@@ -305,7 +306,6 @@ a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
 imgplot = plt.imshow(fgp_dtv_cpu, cmap="gray")
 plt.title('{}'.format('CPU results'))
 
-
 print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
 print ("__________Total nuclear Variation__________")
 print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
@@ -318,9 +318,8 @@ a.set_title('Noisy Image')
 imgplot = plt.imshow(u0,cmap="gray")
 
 channelsNo = 5
-N = 512
-noisyVol = np.zeros((channelsNo,N,N),dtype='float32')
-idealVol = np.zeros((channelsNo,N,N),dtype='float32')
+noisyVol = np.zeros((channelsNo,N,M),dtype='float32')
+idealVol = np.zeros((channelsNo,N,M),dtype='float32')
 
 for i in range (channelsNo):
     noisyVol[i,:,:] = Im + np.random.normal(loc = 0 , scale = perc * Im , size = np.shape(Im))
@@ -361,25 +360,19 @@ plt.title('{}'.format('CPU results'))
 # Uncomment to test 3D regularisation performance 
 #%%
 """
-N = 512
 slices = 20
-
-filename = os.path.join(".." , ".." , ".." , "data" ,"lena_gray_512.tif")
-Im = plt.imread(filename)
-Im = np.asarray(Im, dtype='float32')
-
-Im = Im/255
 perc = 0.05
 
-noisyVol = np.zeros((slices,N,N),dtype='float32')
-noisyRef = np.zeros((slices,N,N),dtype='float32')
-idealVol = np.zeros((slices,N,N),dtype='float32')
+noisyVol = np.zeros((slices,N,M),dtype='float32')
+noisyRef = np.zeros((slices,N,M),dtype='float32')
+idealVol = np.zeros((slices,N,M),dtype='float32')
 
 for i in range (slices):
     noisyVol[i,:,:] = Im + np.random.normal(loc = 0 , scale = perc * Im , size = np.shape(Im))
     noisyRef[i,:,:] = Im + np.random.normal(loc = 0 , scale = 0.01 * Im , size = np.shape(Im))
     idealVol[i,:,:] = Im
 
+
 print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
 print ("_______________ROF-TV (3D)_________________")
 print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
@@ -420,6 +413,7 @@ a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14,
 imgplot = plt.imshow(rof_cpu3D[10,:,:], cmap="gray")
 plt.title('{}'.format('Recovered volume on the CPU using ROF-TV'))
 
+
 print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
 print ("_______________FGP-TV (3D)__________________")
 print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
diff --git a/Wrappers/Python/demos/demo_gpu_regularisers.py b/Wrappers/Python/demos/demo_gpu_regularisers.py
index dfdceee..f3ed50c 100644
--- a/Wrappers/Python/demos/demo_gpu_regularisers.py
+++ b/Wrappers/Python/demos/demo_gpu_regularisers.py
@@ -44,26 +44,28 @@ u0 = Im + np.random.normal(loc = 0 ,
 u_ref = Im + np.random.normal(loc = 0 ,
                                   scale = 0.01 * Im , 
                                   size = np.shape(Im))
+(N,M) = np.shape(u0)
 # map the u0 u0->u0>0
 # f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1)
 u0 = u0.astype('float32')
 u_ref = u_ref.astype('float32')
-
-(N,M) = np.shape(u0)
-u_ref2 = np.zeros([N,M-100],dtype='float32')
-u_ref2[:,0:M-100] = u_ref[:,0:M-100]
+"""
+M = M-100
+u_ref2 = np.zeros([N,M],dtype='float32')
+u_ref2[:,0:M] = u_ref[:,0:M]
 u_ref = u_ref2
 del u_ref2
 
-u02 = np.zeros([N,M-100],dtype='float32')
-u02[:,0:M-100] = u0[:,0:M-100]
+u02 = np.zeros([N,M],dtype='float32')
+u02[:,0:M] = u0[:,0:M]
 u0 = u02
 del u02
 
-Im2 = np.zeros([N,M-100],dtype='float32')
-Im2[:,0:M-100] = Im[:,0:M-100]
+Im2 = np.zeros([N,M],dtype='float32')
+Im2[:,0:M] = Im[:,0:M]
 Im = Im2
 del Im2
+"""
 
 print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
 print ("____________ROF-TV regulariser_____________")
diff --git a/data/SinoInpaint.mat b/data/SinoInpaint.mat
index 58174a9..d748fb4 100644
--- a/data/SinoInpaint.mat
+++ b/data/SinoInpaint.mat