From 985fee04ac1abef2aaa69f282ae6c207e438b4af Mon Sep 17 00:00:00 2001 From: algol Date: Wed, 2 May 2018 11:01:57 +0100 Subject: bugs in cython files --- Core/regularisers_CPU/TNV_core.c | 36 +++++++++++----------- Wrappers/Python/demos/demo_cpu_inpainters.py | 2 +- Wrappers/Python/demos/demo_cpu_regularisers.py | 40 +++++++++++-------------- Wrappers/Python/demos/demo_gpu_regularisers.py | 18 ++++++----- data/SinoInpaint.mat | Bin 3584100 -> 3335061 bytes 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 Binary files a/data/SinoInpaint.mat and b/data/SinoInpaint.mat differ -- cgit v1.2.3