Skip to content

Commit

Permalink
Merge pull request #205 from TomographicImaging/fix_3D_regularisers
Browse files Browse the repository at this point in the history 
ctypes: pass the dimension parameters in the correct order
  • Loading branch information
@casperdcl
casperdcl authored Jun 6, 2024
2 parents b0e1888 + 6f6de90 commit 3c56d78
Show file tree
Hide file tree
Showing 8 changed files with 826 additions and 2,409 deletions.
21 changes: 13 additions & 8 deletions demos/demo_cpu_regularisers3D.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -70,7 +70,7 @@ def printParametersToString(pars):
Im = Im2
del Im2
"""
slices = 15
slices = 20

noisyVol = np.zeros((slices,N,M),dtype='float32')
noisyRef = np.zeros((slices,N,M),dtype='float32')
Expand All @@ -90,14 +90,15 @@ def printParametersToString(pars):
fig = plt.figure()
plt.suptitle('Performance of ROF-TV regulariser using the CPU')
a=fig.add_subplot(1,2,1)
a.set_title('Noisy 15th slice of a volume')
imgplot = plt.imshow(noisyVol[10,:,:],cmap="gray")
slice_num = 10
a.set_title(f'Noisy slice {slice_num} of a volume')
imgplot = plt.imshow(noisyVol[slice_num,:,:],cmap="gray")

# set parameters
pars = {'algorithm': ROF_TV, \
'input' : noisyVol,\
'regularisation_parameter':0.02,\
'number_of_iterations': 7000,\
'regularisation_parameter':0.02 * 100,\
'number_of_iterations': 200,\
'time_marching_parameter': 0.0007,\
'tolerance_constant':1e-06}

Expand All @@ -108,7 +109,7 @@ def printParametersToString(pars):
pars['regularisation_parameter'],
pars['number_of_iterations'],
pars['time_marching_parameter'],
pars['tolerance_constant'], device='cpu', infovector=info_vec_cpu)
pars['tolerance_constant'], device='gpu', infovector=info_vec_cpu)

Qtools = QualityTools(idealVol, rof_cpu3D)
pars['rmse'] = Qtools.rmse()
Expand All @@ -126,6 +127,7 @@ def printParametersToString(pars):
imgplot = plt.imshow(rof_cpu3D[10,:,:], cmap="gray")
plt.title('{}'.format('Recovered volume on the CPU using ROF-TV'))

print (f"information vector: {info_vec_cpu}")
#%%
print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
print ("_______________FGP-TV (3D)__________________")
Expand All @@ -141,8 +143,8 @@ def printParametersToString(pars):
# set parameters
pars = {'algorithm' : FGP_TV, \
'input' : noisyVol,\
'regularisation_parameter':0.02, \
'number_of_iterations' :1000 ,\
'regularisation_parameter': 0.08,# 0.02 * 10000,
'number_of_iterations' :2000 ,\
'tolerance_constant':1e-06,\
'methodTV': 0 ,\
'nonneg': 0}
Expand Down Expand Up @@ -171,6 +173,9 @@ def printParametersToString(pars):
verticalalignment='top', bbox=props)
imgplot = plt.imshow(fgp_cpu3D[10,:,:], cmap="gray")
plt.title('{}'.format('Recovered volume on the CPU using FGP-TV'))

#%%
imgplot = plt.imshow(fgp_cpu3D[:,1,:], cmap="gray")
#%%
print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
print ("_______________PD-TV (3D)__________________")
Expand Down
34 changes: 17 additions & 17 deletions src/Python/ccpi/filters/TV.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,7 +34,7 @@ def TV_ROF_CPU(inputData, regularisation_parameter, iterationsNumb,
infovector = np.zeros((2,), dtype='float32')
infovector_p = infovector.ctypes.data_as(ctypes.POINTER(ctypes.c_float))

dims = list(inputData.shape)
dims = list(inputData.shape)[::-1]
if inputData.ndim == 2:
dims.append(1)

Expand Down Expand Up @@ -79,7 +79,7 @@ def TV_FGP_CPU(inputData, lambdaPar, iterationsNumb, epsil, methodTV, nonneg, ou
infovector = np.zeros((2,), dtype='float32')
infovector_p = infovector.ctypes.data_as(ctypes.POINTER(ctypes.c_float))

dims = list(inputData.shape)
dims = list(inputData.shape)[::-1]
if inputData.ndim == 2:
dims.append(1)

Expand Down Expand Up @@ -122,7 +122,7 @@ def PDTV_CPU(inputData, lambdaPar, iterationsNumb, epsil, lipschitz_const, metho
infovector = np.zeros((2,), dtype='float32')
infovector_p = infovector.ctypes.data_as(ctypes.POINTER(ctypes.c_float))

dims = list(inputData.shape)
dims = list(inputData.shape)[::-1]
if inputData.ndim == 2:
dims.append(1)

Expand Down Expand Up @@ -161,7 +161,7 @@ def SB_TV_CPU(inputData, mu, iter, epsil, methodTV, out=None, infovector=None):
infovector = np.zeros((2,), dtype='float32')
infovector_p = infovector.ctypes.data_as(ctypes.POINTER(ctypes.c_float))

dims = list(inputData.shape)
dims = list(inputData.shape)[::-1]
if inputData.ndim == 2:
dims.append(1)

Expand Down Expand Up @@ -199,7 +199,7 @@ def LLT_ROF_CPU(inputData, lambdaROF, lambdaLLT, iterationsNumb, tau, epsil, out
infovector = np.zeros((2,), dtype='float32')
infovector_p = infovector.ctypes.data_as(ctypes.POINTER(ctypes.c_float))

dims = list(inputData.shape)
dims = list(inputData.shape)[::-1]
if inputData.ndim == 2:
dims.append(1)

Expand Down Expand Up @@ -239,7 +239,7 @@ def TGV_CPU(inputData, lambdaPar, alpha1, alpha0, iterationsNumb, L2, epsil, out
infovector = np.zeros((2,), dtype='float32')
infovector_p = infovector.ctypes.data_as(ctypes.POINTER(ctypes.c_float))

dims = list(inputData.shape)
dims = list(inputData.shape)[::-1]
if inputData.ndim == 2:
dims.append(1)

Expand Down Expand Up @@ -281,7 +281,7 @@ def dTV_FGP_CPU(inputData, inputRef, lambdaPar, iterationsNumb, epsil, eta, meth
infovector = np.zeros((2,), dtype='float32')
infovector_p = infovector.ctypes.data_as(ctypes.POINTER(ctypes.c_float))

dims = list(inputData.shape)
dims = list(inputData.shape)[::-1]
if inputData.ndim == 2:
dims.append(1)

Expand Down Expand Up @@ -312,7 +312,7 @@ def TNV(inputData, lambdaPar, maxIter, tol, out=None):
out = np.zeros_like(inputData)
out_p = out.ctypes.data_as(ctypes.POINTER(ctypes.c_float))

dims = list(inputData.shape)
dims = list(inputData.shape)[::-1]
if inputData.ndim != 3:
raise ValueError('Input data must be 2D + 1 channel')

Expand Down Expand Up @@ -439,7 +439,7 @@ def TV_ENERGY(U, U0, lambdaPar, type, E_val=None):
ctypes.c_int, # dimY (int)
]
cilreg.TV_energy2D.restype = ctypes.c_float # return value is float
result = cilreg.TV_energy2D(u_p, u0_p, e_val_p, lambdaPar, type, dims[0], dims[1])
result = cilreg.TV_energy2D(u_p, u0_p, e_val_p, lambdaPar, type, dims[1], dims[0])
elif U.ndim == 3:
# float TV_energy3D(float *U, float *U0, float *E_val, float lambdaPar, int type, int dimX, int dimY, int dimZ);
cilreg.TV_energy3D.argtypes = [
Expand All @@ -455,7 +455,7 @@ def TV_ENERGY(U, U0, lambdaPar, type, E_val=None):
cilreg.TV_energy3D.restype = ctypes.c_float # return value is float

# float TV_energy3D(float *U, float *U0, float *E_val, float lambdaPar, int type, int dimX, int dimY, int dimZ);
result = cilreg.TV_energy3D(u_p, u0_p, e_val_p, lambdaPar, type, dims[0], dims[1], dims[2])
result = cilreg.TV_energy3D(u_p, u0_p, e_val_p, lambdaPar, type, dims[2], dims[1], dims[0])
else:
raise ValueError(f"TV_ENERGY: Only 2D and 3D data are supported. Got {U.ndim}")
return E_val
Expand Down Expand Up @@ -503,7 +503,7 @@ def TV_ROF_GPU(inputData, regularisation_parameter, iterationsNumb,
infovector = np.zeros((2,), dtype='float32')
infovector_p = infovector.ctypes.data_as(ctypes.POINTER(ctypes.c_float))

dims = list(inputData.shape)
dims = list(inputData.shape)[::-1]
if inputData.ndim == 2:
dims.append(1)

Expand Down Expand Up @@ -551,7 +551,7 @@ def TV_FGP_GPU(inputData, lambdaPar, iterationsNumb, epsil, methodTV, nonneg, gp
infovector = np.zeros((2,), dtype='float32')
infovector_p = infovector.ctypes.data_as(ctypes.POINTER(ctypes.c_float))

dims = list(inputData.shape)
dims = list(inputData.shape)[::-1]
if inputData.ndim == 2:
dims.append(1)
if gpu_device == 'gpu':
Expand Down Expand Up @@ -594,7 +594,7 @@ def PDTV_GPU(Input, lambdaPar, iter, epsil, lipschitz_const, methodTV, nonneg, g
infovector = np.zeros_like(Input)
infovector_p = infovector.ctypes.data_as(ctypes.POINTER(ctypes.c_float))

dims = list(Input.shape)
dims = list(Input.shape)[::-1]
if Input.ndim == 2:
dims.append(1)
# int TV_PD_GPU_main(float *Input, float *Output, float *infovector, float lambdaPar, int iter, float epsil,
Expand Down Expand Up @@ -635,7 +635,7 @@ def SB_TV_GPU(inputData, lambdaPar, iterationsNumb, epsil, methodTV,
infovector = np.zeros((2,), dtype='float32')
infovector_p = infovector.ctypes.data_as(ctypes.POINTER(ctypes.c_float))

dims = list(inputData.shape)
dims = list(inputData.shape)[::-1]
if inputData.ndim == 2:
dims.append(1)

Expand Down Expand Up @@ -676,7 +676,7 @@ def LLT_ROF_GPU(inputData, lambdaROF, lambdaLLT, iterationsNumb, tau, epsil, gpu
infovector = np.zeros((2,), dtype='float32')
infovector_p = infovector.ctypes.data_as(ctypes.POINTER(ctypes.c_float))

dims = list(inputData.shape)
dims = list(inputData.shape)[::-1]
if inputData.ndim == 2:
dims.append(1)

Expand Down Expand Up @@ -719,7 +719,7 @@ def TGV_GPU(inputData, lambdaPar, alpha1, alpha0, iterationsNumb, L2, epsil,
infovector = np.zeros((2,), dtype='float32')
infovector_p = infovector.ctypes.data_as(ctypes.POINTER(ctypes.c_float))

dims = list(inputData.shape)
dims = list(inputData.shape)[::-1]
if inputData.ndim == 2:
dims.append(1)
dims = dims[::-1]
Expand Down Expand Up @@ -765,7 +765,7 @@ def dTV_FGP_GPU(inputData, inputRef, lambdaPar, iterationsNumb, epsil, eta,
infovector = np.zeros((2,), dtype='float32')
infovector_p = infovector.ctypes.data_as(ctypes.POINTER(ctypes.c_float))

dims = list(inputData.shape)
dims = list(inputData.shape)[::-1]
if inputData.ndim == 2:
dims.append(1)

Expand Down
8 changes: 4 additions & 4 deletions src/Python/ccpi/filters/diffusion.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,7 @@ def NDF_CPU(inputData, lambdaPar, sigmaPar, iterationsNumb, tau, penaltytype, ep
infovector = np.zeros((2,), dtype='float32')
infovector_p = infovector.ctypes.data_as(ctypes.POINTER(ctypes.c_float))

dims = list(inputData.shape)
dims = list(inputData.shape)[::-1]
if inputData.ndim == 2:
dims.append(1)

Expand Down Expand Up @@ -74,7 +74,7 @@ def Diffus4th_CPU(inputData, lambdaPar, sigmaPar, iterationsNumb, tau, epsil, ou
infovector = np.zeros((2,), dtype='float32')
infovector_p = infovector.ctypes.data_as(ctypes.POINTER(ctypes.c_float))

dims = list(inputData.shape)
dims = list(inputData.shape)[::-1]
if inputData.ndim == 2:
dims.append(1)

Expand Down Expand Up @@ -117,7 +117,7 @@ def NDF_GPU(inputData, lambdaPar, sigmaPar, iterationsNumb, tau,
infovector = np.zeros((2,), dtype='float32')
infovector_p = infovector.ctypes.data_as(ctypes.POINTER(ctypes.c_float))

dims = list(inputData.shape)
dims = list(inputData.shape)[::-1]
if inputData.ndim == 2:
dims.append(1)

Expand Down Expand Up @@ -159,7 +159,7 @@ def Diffus4th_GPU(inputData, lambdaPar, sigmaPar, iterationsNumb, tau, epsil, gp
infovector = np.zeros((2,), dtype='float32')
infovector_p = infovector.ctypes.data_as(ctypes.POINTER(ctypes.c_float))

dims = list(inputData.shape)
dims = list(inputData.shape)[::-1]
if inputData.ndim == 2:
dims.append(1)

Expand Down
Loading

0 comments on commit 3c56d78

Please sign in to comment.