diff --git a/demos/demo_cpu_regularisers3D.py b/demos/demo_cpu_regularisers3D.py index eb379ab..b3ddf84 100644 --- a/demos/demo_cpu_regularisers3D.py +++ b/demos/demo_cpu_regularisers3D.py @@ -31,7 +31,8 @@ def printParametersToString(pars): return txt ############################################################################### -os.chdir(os.path.join("..", "demos")) +# os.chdir(os.path.join("..", "demos")) +os.chdir("C:/Users/ofn77899/Dev/CCPi-Regularisation-Toolkit/demos") filename = os.path.join( "data" ,"lena_gray_512.tif") # read image @@ -70,7 +71,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') @@ -90,14 +91,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} @@ -108,7 +110,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() @@ -126,6 +128,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)__________________") @@ -141,8 +144,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} @@ -171,6 +174,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)__________________") diff --git a/src/Python/src/cpu_regularisers.pyx b/src/Python/src/cpu_regularisers.pyx deleted file mode 100644 index 655a508..0000000 --- a/src/Python/src/cpu_regularisers.pyx +++ /dev/null @@ -1,734 +0,0 @@ -# distutils: language=c++ -""" -Copyright 2018 CCPi -Licensed under the Apache License, Version 2.0 (the "License"); -you may not use this file except in compliance with the License. -You may obtain a copy of the License at - http://www.apache.org/licenses/LICENSE-2.0 -Unless required by applicable law or agreed to in writing, software -distributed under the License is distributed on an "AS IS" BASIS, -WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -See the License for the specific language governing permissions and -limitations under the License. - -Author: Edoardo Pasca, Daniil Kazantsev -""" - -import cython -import numpy as np -cimport numpy as np - -cdef extern float TV_ROF_CPU_main(float *Input, float *Output, float *infovector, float *lambdaPar, int lambda_is_arr, 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 PDTV_CPU_main(float *Input, float *U, float *infovector, float lambdaPar, int iterationsNumb, float epsil, float lipschitz_const, int methodTV, int nonneg, int dimX, int dimY, int dimZ); -cdef extern float SB_TV_CPU_main(float *Input, float *Output, float *infovector, float mu, int iter, float epsil, int methodTV, int dimX, int dimY, int dimZ); -cdef extern float LLT_ROF_CPU_main(float *Input, float *Output, float *infovector, float lambdaROF, float lambdaLLT, int iterationsNumb, float tau, float epsil, int dimX, int dimY, int dimZ); -cdef extern float TGV_main(float *Input, float *Output, float *infovector, float lambdaPar, float alpha1, float alpha0, int iterationsNumb, float L2, float epsil, int dimX, int dimY, int dimZ); -cdef extern float Diffusion_CPU_main(float *Input, float *Output, float *infovector, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int penaltytype, float epsil, int dimX, int dimY, int dimZ); -cdef extern float Diffus4th_CPU_main(float *Input, float *Output, float *infovector, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, float epsil, int dimX, int dimY, int dimZ); -cdef extern float dTV_FGP_CPU_main(float *Input, float *InputRef, float *Output, float *infovector, float lambdaPar, int iterationsNumb, float epsil, float eta, int methodTV, int nonneg, int dimX, int dimY, int dimZ); -cdef extern float TNV_CPU_main(float *Input, float *u, float lambdaPar, int maxIter, float tol, int dimX, int dimY, int dimZ); -cdef extern float PatchSelect_CPU_main(float *Input, unsigned short *H_i, unsigned short *H_j, unsigned short *H_k, float *Weights, int dimX, int dimY, int dimZ, int SearchWindow, int SimilarWin, int NumNeighb, float h); -cdef extern float Nonlocal_TV_CPU_main(float *A_orig, float *Output, unsigned short *H_i, unsigned short *H_j, unsigned short *H_k, float *Weights, int dimX, int dimY, int dimZ, int NumNeighb, float lambdaReg, int IterNumb, int switchM); - -cdef extern float TV_energy2D(float *U, float *U0, float *E_val, float lambdaPar, int type, int dimX, int dimY); -cdef extern float TV_energy3D(float *U, float *U0, float *E_val, float lambdaPar, int type, int dimX, int dimY, int dimZ); -#****************************************************************# -#********************** Total-variation ROF *********************# -#****************************************************************# -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,tolerance_param) - elif inputData.ndim == 3: - 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, - regularisation_parameter, - int iterationsNumb, - float marching_step_parameter, - float tolerance_param): - cdef long dims[2] - dims[0] = inputData.shape[0] - dims[1] = inputData.shape[1] - cdef float lambdareg - cdef np.ndarray[np.float32_t, ndim=2, mode="c"] reg - 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') - - if isinstance (regularisation_parameter, np.ndarray): - reg = regularisation_parameter.copy() - TV_ROF_CPU_main(&inputData[0,0], &outputData[0,0], &infovec[0], ®[0,0], 1, iterationsNumb, marching_step_parameter, tolerance_param, dims[1], dims[0], 1) - else: # supposedly this would be a float - lambdareg = regularisation_parameter; - TV_ROF_CPU_main(&inputData[0,0], &outputData[0,0], &infovec[0], &lambdareg, 0, iterationsNumb, marching_step_parameter, tolerance_param, dims[1], dims[0], 1) - return (outputData,infovec) - -def TV_ROF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, - regularisation_parameter, - int iterationsNumb, - float marching_step_parameter, - float tolerance_param): - cdef long dims[3] - dims[0] = inputData.shape[0] - dims[1] = inputData.shape[1] - dims[2] = inputData.shape[2] - cdef float lambdareg - cdef np.ndarray[np.float32_t, ndim=3, mode="c"] reg - 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], &infovec[0], regularisation_parameter, iterationsNumb, marching_step_parameter, tolerance_param, dims[2], dims[1], dims[0]) - if isinstance (regularisation_parameter, np.ndarray): - reg = regularisation_parameter.copy() - TV_ROF_CPU_main(&inputData[0,0,0], &outputData[0,0,0], &infovec[0], ®[0,0,0], 1, iterationsNumb, marching_step_parameter, tolerance_param, dims[2], dims[1], dims[0]) - else: # supposedly this would be a float - lambdareg = regularisation_parameter - TV_ROF_CPU_main(&inputData[0,0,0], &outputData[0,0,0], &infovec[0], &lambdareg, 0, iterationsNumb, marching_step_parameter, tolerance_param, dims[2], dims[1], dims[0]) - return (outputData,infovec) - -#****************************************************************# -#********************** Total-variation FGP *********************# -#****************************************************************# -#******** Total-variation Fast-Gradient-Projection (FGP)*********# -def TV_FGP_CPU(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg): - if inputData.ndim == 2: - return TV_FGP_2D(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg) - elif inputData.ndim == 3: - return TV_FGP_3D(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg) - -def TV_FGP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, - float regularisation_parameter, - int iterationsNumb, - float tolerance_param, - int methodTV, - int nonneg): - - 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 FGP-TV iterations for 2D data */ - TV_FGP_CPU_main(&inputData[0,0], &outputData[0,0], &infovec[0], regularisation_parameter, - iterationsNumb, - tolerance_param, - methodTV, - nonneg, - dims[1],dims[0],1) - - return (outputData,infovec) - -def TV_FGP_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, - float regularisation_parameter, - int iterationsNumb, - float tolerance_param, - int methodTV, - int nonneg): - - cdef long dims[3] - dims[0] = inputData.shape[0] - dims[1] = inputData.shape[1] - dims[2] = inputData.shape[2] - - 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.zeros([2], dtype='float32') - - #/* Run FGP-TV iterations for 3D data */ - TV_FGP_CPU_main(&inputData[0,0,0], &outputData[0,0,0], &infovec[0], regularisation_parameter, - iterationsNumb, - tolerance_param, - methodTV, - nonneg, - dims[2], dims[1], dims[0]) - return (outputData,infovec) - -#****************************************************************# -#****************** Total-variation Primal-dual *****************# -#****************************************************************# -def TV_PD_CPU(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, lipschitz_const): - if inputData.ndim == 2: - return TV_PD_2D(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, lipschitz_const) - elif inputData.ndim == 3: - return TV_PD_3D(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV, nonneg, lipschitz_const) - -def TV_PD_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, - float regularisation_parameter, - int iterationsNumb, - float tolerance_param, - int methodTV, - int nonneg, - float lipschitz_const): - - 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 FGP-TV iterations for 2D data */ - PDTV_CPU_main(&inputData[0,0], &outputData[0,0], &infovec[0], regularisation_parameter, - iterationsNumb, - tolerance_param, - lipschitz_const, - methodTV, - nonneg, - dims[1],dims[0], 1) - return (outputData,infovec) - -def TV_PD_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, - float regularisation_parameter, - int iterationsNumb, - float tolerance_param, - int methodTV, - int nonneg, - float lipschitz_const): - - cdef long dims[3] - dims[0] = inputData.shape[0] - dims[1] = inputData.shape[1] - dims[2] = inputData.shape[2] - - 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.zeros([2], dtype='float32') - - #/* Run FGP-TV iterations for 3D data */ - PDTV_CPU_main(&inputData[0,0,0], &outputData[0,0,0], &infovec[0], regularisation_parameter, - iterationsNumb, - tolerance_param, - lipschitz_const, - methodTV, - nonneg, - dims[2], dims[1], dims[0]) - return (outputData,infovec) - -#***************************************************************# -#********************** Total-variation SB *********************# -#***************************************************************# -#*************** Total-variation Split Bregman (SB)*************# -def TV_SB_CPU(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV): - if inputData.ndim == 2: - return TV_SB_2D(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV) - elif inputData.ndim == 3: - return TV_SB_3D(inputData, regularisation_parameter, iterationsNumb, tolerance_param, methodTV) - -def TV_SB_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, - float regularisation_parameter, - int iterationsNumb, - float tolerance_param, - int methodTV): - - 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.zeros([2], dtype='float32') - - #/* Run SB-TV iterations for 2D data */ - SB_TV_CPU_main(&inputData[0,0], &outputData[0,0], &infovec[0], - regularisation_parameter, - iterationsNumb, - tolerance_param, - methodTV, - dims[1],dims[0], 1) - - return (outputData,infovec) - -def TV_SB_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, - float regularisation_parameter, - int iterationsNumb, - float tolerance_param, - int methodTV): - cdef long dims[3] - dims[0] = inputData.shape[0] - dims[1] = inputData.shape[1] - dims[2] = inputData.shape[2] - - 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.zeros([2], dtype='float32') - - #/* Run SB-TV iterations for 3D data */ - SB_TV_CPU_main(&inputData[0,0,0], &outputData[0,0,0], &infovec[0], - regularisation_parameter, - iterationsNumb, - tolerance_param, - methodTV, - dims[2], dims[1], dims[0]) - return (outputData,infovec) -#***************************************************************# -#******************* ROF - LLT regularisation ******************# -#***************************************************************# -def LLT_ROF_CPU(inputData, regularisation_parameterROF, regularisation_parameterLLT, iterations, time_marching_parameter, tolerance_param): - if inputData.ndim == 2: - return LLT_ROF_2D(inputData, regularisation_parameterROF, regularisation_parameterLLT, iterations, time_marching_parameter, tolerance_param) - elif inputData.ndim == 3: - return LLT_ROF_3D(inputData, regularisation_parameterROF, regularisation_parameterLLT, iterations, time_marching_parameter, tolerance_param) - -def LLT_ROF_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, - float regularisation_parameterROF, - float regularisation_parameterLLT, - int iterations, - float time_marching_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.zeros([2], dtype='float32') - - #/* Run ROF-LLT iterations for 2D data */ - LLT_ROF_CPU_main(&inputData[0,0], &outputData[0,0], &infovec[0], regularisation_parameterROF, regularisation_parameterLLT, iterations, time_marching_parameter, - tolerance_param, - dims[1],dims[0],1) - return (outputData,infovec) - -def LLT_ROF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, - float regularisation_parameterROF, - float regularisation_parameterLLT, - int iterations, - float time_marching_parameter, - float tolerance_param): - - cdef long dims[3] - dims[0] = inputData.shape[0] - dims[1] = inputData.shape[1] - dims[2] = inputData.shape[2] - - 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.zeros([2], dtype='float32') - - #/* Run ROF-LLT iterations for 3D data */ - LLT_ROF_CPU_main(&inputData[0,0,0], &outputData[0,0,0], &infovec[0], regularisation_parameterROF, regularisation_parameterLLT, iterations, - time_marching_parameter, - tolerance_param, - dims[2], dims[1], dims[0]) - return (outputData,infovec) -#***************************************************************# -#***************** Total Generalised Variation *****************# -#***************************************************************# -def TGV_CPU(inputData, regularisation_parameter, alpha1, alpha0, iterations, LipshitzConst, tolerance_param): - if inputData.ndim == 2: - return TGV_2D(inputData, regularisation_parameter, alpha1, alpha0, - iterations, LipshitzConst, tolerance_param) - elif inputData.ndim == 3: - return TGV_3D(inputData, regularisation_parameter, alpha1, alpha0, - iterations, LipshitzConst, tolerance_param) - -def TGV_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, - float regularisation_parameter, - float alpha1, - float alpha0, - int iterationsNumb, - float LipshitzConst, - 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.zeros([2], dtype='float32') - - #/* Run TGV iterations for 2D data */ - TGV_main(&inputData[0,0], &outputData[0,0], &infovec[0], regularisation_parameter, - alpha1, - alpha0, - iterationsNumb, - LipshitzConst, - tolerance_param, - dims[1],dims[0],1) - return (outputData,infovec) -def TGV_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, - float regularisation_parameter, - float alpha1, - float alpha0, - int iterationsNumb, - float LipshitzConst, - float tolerance_param): - - cdef long dims[3] - dims[0] = inputData.shape[0] - dims[1] = inputData.shape[1] - dims[2] = inputData.shape[2] - - 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.zeros([2], dtype='float32') - - #/* Run TGV iterations for 3D data */ - TGV_main(&inputData[0,0,0], &outputData[0,0,0], &infovec[0], regularisation_parameter, - alpha1, - alpha0, - iterationsNumb, - LipshitzConst, - tolerance_param, - dims[2], dims[1], dims[0]) - return (outputData,infovec) - -#****************************************************************# -#***************Nonlinear (Isotropic) Diffusion******************# -#****************************************************************# -def NDF_CPU(inputData, regularisation_parameter, edge_parameter, iterationsNumb,time_marching_parameter, penalty_type,tolerance_param): - if inputData.ndim == 2: - return NDF_2D(inputData, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type, tolerance_param) - elif inputData.ndim == 3: - return NDF_3D(inputData, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type, tolerance_param) - -def NDF_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, - float regularisation_parameter, - float edge_parameter, - int iterationsNumb, - float time_marching_parameter, - int penalty_type, - 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.zeros([2], dtype='float32') - - # Run Nonlinear Diffusion iterations for 2D data - Diffusion_CPU_main(&inputData[0,0], &outputData[0,0], &infovec[0], - regularisation_parameter, edge_parameter, iterationsNumb, - time_marching_parameter, penalty_type, - tolerance_param, - dims[1], dims[0], 1) - return (outputData,infovec) - -def NDF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, - float regularisation_parameter, - float edge_parameter, - int iterationsNumb, - float time_marching_parameter, - int penalty_type, - float tolerance_param): - cdef long dims[3] - dims[0] = inputData.shape[0] - dims[1] = inputData.shape[1] - dims[2] = inputData.shape[2] - - 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.zeros([2], dtype='float32') - - # Run Nonlinear Diffusion iterations for 3D data - Diffusion_CPU_main(&inputData[0,0,0], &outputData[0,0,0], &infovec[0], - regularisation_parameter, edge_parameter, iterationsNumb, - time_marching_parameter, penalty_type, - tolerance_param, - dims[2], dims[1], dims[0]) - return (outputData,infovec) - -#****************************************************************# -#*************Anisotropic Fourth-Order diffusion*****************# -#****************************************************************# -def Diff4th_CPU(inputData, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter,tolerance_param): - if inputData.ndim == 2: - return Diff4th_2D(inputData, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter,tolerance_param) - elif inputData.ndim == 3: - return Diff4th_3D(inputData, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter,tolerance_param) - -def Diff4th_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, - float regularisation_parameter, - float edge_parameter, - int iterationsNumb, - float time_marching_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.zeros([2], dtype='float32') - - # Run Anisotropic Fourth-Order diffusion for 2D data - Diffus4th_CPU_main(&inputData[0,0], &outputData[0,0], &infovec[0], - regularisation_parameter, - edge_parameter, iterationsNumb, - time_marching_parameter, - tolerance_param, - dims[1], dims[0], 1) - return (outputData,infovec) - -def Diff4th_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, - float regularisation_parameter, - float edge_parameter, - int iterationsNumb, - float time_marching_parameter, - float tolerance_param): - cdef long dims[3] - dims[0] = inputData.shape[0] - dims[1] = inputData.shape[1] - dims[2] = inputData.shape[2] - - 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.zeros([2], dtype='float32') - - # Run Anisotropic Fourth-Order diffusion for 3D data - Diffus4th_CPU_main(&inputData[0,0,0], &outputData[0,0,0], &infovec[0], - regularisation_parameter, edge_parameter, - iterationsNumb, time_marching_parameter, - tolerance_param, - dims[2], dims[1], dims[0]) - return (outputData,infovec) -#****************************************************************# -#**************Directional Total-variation FGP ******************# -#****************************************************************# -#******** Directional TV Fast-Gradient-Projection (FGP)*********# -def dTV_FGP_CPU(inputData, refdata, regularisation_parameter, iterationsNumb, tolerance_param, eta_const, methodTV, nonneg): - if inputData.ndim == 2: - return dTV_FGP_2D(inputData, refdata, regularisation_parameter, iterationsNumb, tolerance_param, eta_const, methodTV, nonneg) - elif inputData.ndim == 3: - return dTV_FGP_3D(inputData, refdata, regularisation_parameter, iterationsNumb, tolerance_param, eta_const, methodTV, nonneg) - -def dTV_FGP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, - np.ndarray[np.float32_t, ndim=2, mode="c"] refdata, - float regularisation_parameter, - int iterationsNumb, - float tolerance_param, - float eta_const, - int methodTV, - int nonneg): - - 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.zeros([2], dtype='float32') - - #/* Run FGP-dTV iterations for 2D data */ - dTV_FGP_CPU_main(&inputData[0,0], &refdata[0,0], &outputData[0,0], &infovec[0], - regularisation_parameter, - iterationsNumb, - tolerance_param, - eta_const, - methodTV, - nonneg, - dims[1], dims[0], 1) - return (outputData,infovec) - -def dTV_FGP_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, - np.ndarray[np.float32_t, ndim=3, mode="c"] refdata, - float regularisation_parameter, - int iterationsNumb, - float tolerance_param, - float eta_const, - int methodTV, - int nonneg): - cdef long dims[3] - dims[0] = inputData.shape[0] - dims[1] = inputData.shape[1] - dims[2] = inputData.shape[2] - - 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.zeros([2], dtype='float32') - - #/* Run FGP-dTV iterations for 3D data */ - dTV_FGP_CPU_main(&inputData[0,0,0], &refdata[0,0,0], &outputData[0,0,0], &infovec[0], - regularisation_parameter, - iterationsNumb, - tolerance_param, - eta_const, - methodTV, - nonneg, - dims[2], dims[1], dims[0]) - return (outputData,infovec) - -#****************************************************************# -#*********************Total Nuclear Variation********************# -#****************************************************************# -def TNV_CPU(inputData, regularisation_parameter, iterationsNumb, tolerance_param): - if inputData.ndim == 2: - return - elif inputData.ndim == 3: - return TNV_3D(inputData, regularisation_parameter, iterationsNumb, tolerance_param) - -def TNV_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, - float regularisation_parameter, - int iterationsNumb, - float tolerance_param): - cdef long dims[3] - dims[0] = inputData.shape[0] - dims[1] = inputData.shape[1] - dims[2] = inputData.shape[2] - - cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \ - np.zeros([dims[0],dims[1],dims[2]], dtype='float32') - - # Run TNV iterations for 3D (X,Y,Channels) data - TNV_CPU_main(&inputData[0,0,0], &outputData[0,0,0], regularisation_parameter, iterationsNumb, tolerance_param, dims[2], dims[1], dims[0]) - return outputData -#****************************************************************# -#***************Patch-based weights calculation******************# -#****************************************************************# -def PATCHSEL_CPU(inputData, searchwindow, patchwindow, neighbours, edge_parameter): - if inputData.ndim == 2: - return PatchSel_2D(inputData, searchwindow, patchwindow, neighbours, edge_parameter) - elif inputData.ndim == 3: - return 1 -def PatchSel_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, - int searchwindow, - int patchwindow, - int neighbours, - float edge_parameter): - cdef long dims[3] - dims[0] = neighbours - dims[1] = inputData.shape[0] - dims[2] = inputData.shape[1] - - - cdef np.ndarray[np.float32_t, ndim=3, mode="c"] Weights = \ - np.zeros([dims[0], dims[1],dims[2]], dtype='float32') - - cdef np.ndarray[np.uint16_t, ndim=3, mode="c"] H_i = \ - np.zeros([dims[0], dims[1],dims[2]], dtype='uint16') - - cdef np.ndarray[np.uint16_t, ndim=3, mode="c"] H_j = \ - np.zeros([dims[0], dims[1],dims[2]], dtype='uint16') - - # Run patch-based weight selection function - PatchSelect_CPU_main(&inputData[0,0], &H_j[0,0,0], &H_i[0,0,0], &H_i[0,0,0], &Weights[0,0,0], dims[2], dims[1], 0, searchwindow, patchwindow, neighbours, edge_parameter) - return H_i, H_j, Weights -""" -def PatchSel_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, - int searchwindow, - int patchwindow, - int neighbours, - float edge_parameter): - cdef long dims[4] - dims[0] = inputData.shape[0] - dims[1] = inputData.shape[1] - dims[2] = inputData.shape[2] - dims[3] = neighbours - - cdef np.ndarray[np.float32_t, ndim=4, mode="c"] Weights = \ - np.zeros([dims[3],dims[0],dims[1],dims[2]], dtype='float32') - - cdef np.ndarray[np.uint16_t, ndim=4, mode="c"] H_i = \ - np.zeros([dims[3],dims[0],dims[1],dims[2]], dtype='uint16') - - cdef np.ndarray[np.uint16_t, ndim=4, mode="c"] H_j = \ - np.zeros([dims[3],dims[0],dims[1],dims[2]], dtype='uint16') - - cdef np.ndarray[np.uint16_t, ndim=4, mode="c"] H_k = \ - np.zeros([dims[3],dims[0],dims[1],dims[2]], dtype='uint16') - - # Run patch-based weight selection function - PatchSelect_CPU_main(&inputData[0,0,0], &H_i[0,0,0,0], &H_j[0,0,0,0], &H_k[0,0,0,0], &Weights[0,0,0,0], dims[2], dims[1], dims[0], searchwindow, patchwindow, neighbours, edge_parameter, 1) - return H_i, H_j, H_k, Weights -""" - -#****************************************************************# -#***************Non-local Total Variation******************# -#****************************************************************# -def NLTV_CPU(inputData, H_i, H_j, H_k, Weights, regularisation_parameter, iterations): - if inputData.ndim == 2: - return NLTV_2D(inputData, H_i, H_j, Weights, regularisation_parameter, iterations) - elif inputData.ndim == 3: - return 1 -def NLTV_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, - np.ndarray[np.uint16_t, ndim=3, mode="c"] H_i, - np.ndarray[np.uint16_t, ndim=3, mode="c"] H_j, - np.ndarray[np.float32_t, ndim=3, mode="c"] Weights, - float regularisation_parameter, - int iterations): - - cdef long dims[2] - dims[0] = inputData.shape[0] - dims[1] = inputData.shape[1] - neighbours = H_i.shape[0] - - cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \ - np.zeros([dims[0],dims[1]], dtype='float32') - - # Run nonlocal TV regularisation - Nonlocal_TV_CPU_main(&inputData[0,0], &outputData[0,0], &H_i[0,0,0], &H_j[0,0,0], &H_i[0,0,0], &Weights[0,0,0], dims[1], dims[0], 0, neighbours, regularisation_parameter, iterations, 1) - return outputData - -#****************************************************************# -#***************Calculation of TV-energy functional**************# -#****************************************************************# -def TV_ENERGY(inputData, inputData0, regularisation_parameter, typeFunctional): - if inputData.ndim == 2: - return TV_ENERGY_2D(inputData, inputData0, regularisation_parameter, typeFunctional) - elif inputData.ndim == 3: - return TV_ENERGY_3D(inputData, inputData0, regularisation_parameter, typeFunctional) - -def TV_ENERGY_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, - np.ndarray[np.float32_t, ndim=2, mode="c"] inputData0, - float regularisation_parameter, - int typeFunctional): - - cdef long dims[2] - dims[0] = inputData.shape[0] - dims[1] = inputData.shape[1] - - cdef np.ndarray[np.float32_t, ndim=1, mode="c"] outputData = \ - np.zeros([1], dtype='float32') - - # run function - TV_energy2D(&inputData[0,0], &inputData0[0,0], &outputData[0], regularisation_parameter, typeFunctional, dims[1], dims[0]) - - return outputData - -def TV_ENERGY_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, - np.ndarray[np.float32_t, ndim=3, mode="c"] inputData0, - float regularisation_parameter, - int typeFunctional): - - cdef long dims[3] - dims[0] = inputData.shape[0] - dims[1] = inputData.shape[1] - dims[2] = inputData.shape[2] - - cdef np.ndarray[np.float32_t, ndim=1, mode="c"] outputData = \ - np.zeros([1], dtype='float32') - - # Run function - TV_energy3D(&inputData[0,0,0], &inputData0[0,0,0], &outputData[0], regularisation_parameter, typeFunctional, dims[2], dims[1], dims[0]) - - return outputData