diff options
Diffstat (limited to 'src/Python')
| -rw-r--r-- | src/Python/ccpi/filters/regularisers.py | 142 | ||||
| -rw-r--r-- | src/Python/src/cpu_regularisers.pyx | 503 | ||||
| -rw-r--r-- | src/Python/src/gpu_regularisers.pyx | 474 |
3 files changed, 594 insertions, 525 deletions
diff --git a/src/Python/ccpi/filters/regularisers.py b/src/Python/ccpi/filters/regularisers.py index 09b465a..398e11c 100644 --- a/src/Python/ccpi/filters/regularisers.py +++ b/src/Python/ccpi/filters/regularisers.py @@ -7,7 +7,7 @@ try: from ccpi.filters.gpu_regularisers import TV_ROF_GPU, TV_FGP_GPU, TV_SB_GPU, dTV_FGP_GPU, NDF_GPU, Diff4th_GPU, TGV_GPU, LLT_ROF_GPU, PATCHSEL_GPU gpu_enabled = True except ImportError: - gpu_enabled = False + gpu_enabled = False from ccpi.filters.cpu_regularisers import NDF_INPAINT_CPU, NVM_INPAINT_CPU def ROF_TV(inputData, regularisation_parameter, iterations, @@ -15,13 +15,13 @@ def ROF_TV(inputData, regularisation_parameter, iterations, if device == 'cpu': return TV_ROF_CPU(inputData, regularisation_parameter, - iterations, + iterations, time_marching_parameter, tolerance_param) elif device == 'gpu' and gpu_enabled: return TV_ROF_GPU(inputData, regularisation_parameter, - iterations, + iterations, time_marching_parameter, tolerance_param) else: @@ -35,14 +35,14 @@ def FGP_TV(inputData, regularisation_parameter,iterations, if device == 'cpu': return TV_FGP_CPU(inputData, regularisation_parameter, - iterations, + iterations, tolerance_param, methodTV, nonneg) elif device == 'gpu' and gpu_enabled: return TV_FGP_GPU(inputData, regularisation_parameter, - iterations, + iterations, tolerance_param, methodTV, nonneg) @@ -56,13 +56,13 @@ def SB_TV(inputData, regularisation_parameter, iterations, if device == 'cpu': return TV_SB_CPU(inputData, regularisation_parameter, - iterations, + iterations, tolerance_param, methodTV) elif device == 'gpu' and gpu_enabled: return TV_SB_GPU(inputData, regularisation_parameter, - iterations, + iterations, tolerance_param, methodTV) else: @@ -81,91 +81,115 @@ def LLT_ROF(inputData, regularisation_parameterROF, regularisation_parameterLLT, raise ValueError ('GPU is not available') raise ValueError('Unknown device {0}. Expecting gpu or cpu'\ .format(device)) -def FGP_dTV(inputData, refdata, regularisation_parameter, iterations, - tolerance_param, eta_const, methodTV, nonneg, printM, device='cpu'): +def TGV(inputData, regularisation_parameter, alpha1, alpha0, iterations, + LipshitzConst, tolerance_param, device='cpu'): if device == 'cpu': - return dTV_FGP_CPU(inputData, - refdata, - regularisation_parameter, - iterations, - tolerance_param, - eta_const, - methodTV, - nonneg, - printM) + return TGV_CPU(inputData, + regularisation_parameter, + alpha1, + alpha0, + iterations, + LipshitzConst, + tolerance_param) elif device == 'gpu' and gpu_enabled: - return dTV_FGP_GPU(inputData, - refdata, - regularisation_parameter, - iterations, - tolerance_param, - eta_const, - methodTV, - nonneg, - printM) + return TGV_GPU(inputData, + regularisation_parameter, + alpha1, + alpha0, + iterations, + LipshitzConst, + tolerance_param) else: if not gpu_enabled and device == 'gpu': raise ValueError ('GPU is not available') raise ValueError('Unknown device {0}. Expecting gpu or cpu'\ .format(device)) -def TNV(inputData, regularisation_parameter, iterations, tolerance_param): - return TNV_CPU(inputData, - regularisation_parameter, - iterations, - tolerance_param) def NDF(inputData, regularisation_parameter, edge_parameter, iterations, - time_marching_parameter, penalty_type, device='cpu'): + time_marching_parameter, penalty_type, tolerance_param, device='cpu'): if device == 'cpu': return NDF_CPU(inputData, regularisation_parameter, edge_parameter, - iterations, + iterations, time_marching_parameter, - penalty_type) + penalty_type, + tolerance_param) elif device == 'gpu' and gpu_enabled: return NDF_GPU(inputData, regularisation_parameter, edge_parameter, - iterations, + iterations, time_marching_parameter, - penalty_type) + penalty_type, + tolerance_param) else: if not gpu_enabled and device == 'gpu': raise ValueError ('GPU is not available') raise ValueError('Unknown device {0}. Expecting gpu or cpu'\ .format(device)) def Diff4th(inputData, regularisation_parameter, edge_parameter, iterations, - time_marching_parameter, device='cpu'): + time_marching_parameter, tolerance_param, device='cpu'): if device == 'cpu': return Diff4th_CPU(inputData, regularisation_parameter, edge_parameter, - iterations, - time_marching_parameter) + iterations, + time_marching_parameter, + tolerance_param) elif device == 'gpu' and gpu_enabled: return Diff4th_GPU(inputData, regularisation_parameter, edge_parameter, - iterations, - time_marching_parameter) + iterations, + time_marching_parameter, + tolerance_param) else: if not gpu_enabled and device == 'gpu': raise ValueError ('GPU is not available') raise ValueError('Unknown device {0}. Expecting gpu or cpu'\ .format(device)) - +def FGP_dTV(inputData, refdata, regularisation_parameter, iterations, + tolerance_param, eta_const, methodTV, nonneg, device='cpu'): + if device == 'cpu': + return dTV_FGP_CPU(inputData, + refdata, + regularisation_parameter, + iterations, + tolerance_param, + eta_const, + methodTV, + nonneg) + elif device == 'gpu' and gpu_enabled: + return dTV_FGP_GPU(inputData, + refdata, + regularisation_parameter, + iterations, + tolerance_param, + eta_const, + methodTV, + nonneg) + else: + if not gpu_enabled and device == 'gpu': + raise ValueError ('GPU is not available') + raise ValueError('Unknown device {0}. Expecting gpu or cpu'\ + .format(device)) +def TNV(inputData, regularisation_parameter, iterations, tolerance_param): + return TNV_CPU(inputData, + regularisation_parameter, + iterations, + tolerance_param) def PatchSelect(inputData, searchwindow, patchwindow, neighbours, edge_parameter, device='cpu'): if device == 'cpu': return PATCHSEL_CPU(inputData, searchwindow, patchwindow, - neighbours, + neighbours, edge_parameter) elif device == 'gpu' and gpu_enabled: return PATCHSEL_GPU(inputData, searchwindow, patchwindow, - neighbours, + neighbours, edge_parameter) else: if not gpu_enabled and device == 'gpu': @@ -177,36 +201,14 @@ def NLTV(inputData, H_i, H_j, H_k, Weights, regularisation_parameter, iterations return NLTV_CPU(inputData, H_i, H_j, - H_k, + H_k, Weights, regularisation_parameter, iterations) - -def TGV(inputData, regularisation_parameter, alpha1, alpha0, iterations, - LipshitzConst, device='cpu'): - if device == 'cpu': - return TGV_CPU(inputData, - regularisation_parameter, - alpha1, - alpha0, - iterations, - LipshitzConst) - elif device == 'gpu' and gpu_enabled: - return TGV_GPU(inputData, - regularisation_parameter, - alpha1, - alpha0, - iterations, - LipshitzConst) - else: - if not gpu_enabled and device == 'gpu': - raise ValueError ('GPU is not available') - raise ValueError('Unknown device {0}. Expecting gpu or cpu'\ - .format(device)) def NDF_INP(inputData, maskData, regularisation_parameter, edge_parameter, iterations, time_marching_parameter, penalty_type): - return NDF_INPAINT_CPU(inputData, maskData, regularisation_parameter, + return NDF_INPAINT_CPU(inputData, maskData, regularisation_parameter, edge_parameter, iterations, time_marching_parameter, penalty_type) - + def NVM_INP(inputData, maskData, SW_increment, iterations): return NVM_INPAINT_CPU(inputData, maskData, SW_increment, iterations) diff --git a/src/Python/src/cpu_regularisers.pyx b/src/Python/src/cpu_regularisers.pyx index f2276bb..add641b 100644 --- a/src/Python/src/cpu_regularisers.pyx +++ b/src/Python/src/cpu_regularisers.pyx @@ -22,11 +22,11 @@ cdef extern float TV_ROF_CPU_main(float *Input, float *Output, float *infovector 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 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 lambdaPar, float alpha1, float alpha0, int iterationsNumb, float L2, int dimX, int dimY, int dimZ); -cdef extern float Diffusion_CPU_main(float *Input, float *Output, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int penaltytype, int dimX, int dimY, int dimZ); -cdef extern float Diffus4th_CPU_main(float *Input, float *Output, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, 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 dTV_FGP_CPU_main(float *Input, float *InputRef, float *Output, float lambdaPar, int iterationsNumb, float epsil, float eta, int methodTV, int nonneg, int printM, 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, int switchM); 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); @@ -43,7 +43,7 @@ def TV_ROF_CPU(inputData, regularisation_parameter, iterationsNumb, marching_ste 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, +def TV_ROF_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, float regularisation_parameter, int iterationsNumb, float marching_step_parameter, @@ -51,18 +51,18 @@ def TV_ROF_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 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 ROF iterations for 2D data + + # Run ROF iterations for 2D data TV_ROF_CPU_main(&inputData[0,0], &outputData[0,0], &infovec[0], regularisation_parameter, 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, + +def TV_ROF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, float regularisation_parameter, int iterationsNumb, float marching_step_parameter, @@ -71,13 +71,13 @@ def TV_ROF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 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.ones([2], dtype='float32') - - # Run ROF iterations for 3D data + + # 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]) return (outputData,infovec) @@ -92,9 +92,9 @@ def TV_FGP_CPU(inputData, regularisation_parameter, iterationsNumb, tolerance_pa 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, +def TV_FGP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, float regularisation_parameter, - int iterationsNumb, + int iterationsNumb, float tolerance_param, int methodTV, int nonneg): @@ -102,35 +102,35 @@ def TV_FGP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 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, + 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, + +def TV_FGP_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, float regularisation_parameter, - int iterationsNumb, + 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 = \ @@ -138,7 +138,7 @@ def TV_FGP_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, #/* Run FGP-TV iterations for 3D data */ TV_FGP_CPU_main(&inputData[0,0,0], &outputData[0,0,0], &infovec[0], regularisation_parameter, - iterationsNumb, + iterationsNumb, tolerance_param, methodTV, nonneg, @@ -155,54 +155,54 @@ def TV_SB_CPU(inputData, regularisation_parameter, iterationsNumb, tolerance_par 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, +def TV_SB_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, float regularisation_parameter, - int iterationsNumb, + 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, + 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, + +def TV_SB_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, float regularisation_parameter, - int iterationsNumb, + 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, + iterationsNumb, tolerance_param, methodTV, dims[2], dims[1], dims[0]) - return (outputData,infovec) + return (outputData,infovec) #***************************************************************# #******************* ROF - LLT regularisation ******************# #***************************************************************# @@ -212,288 +212,321 @@ def LLT_ROF_CPU(inputData, regularisation_parameterROF, regularisation_parameter 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, +def LLT_ROF_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, float regularisation_parameterROF, float regularisation_parameterLLT, - int iterations, + 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, + 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) + return (outputData,infovec) -def LLT_ROF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, +def LLT_ROF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, float regularisation_parameterROF, float regularisation_parameterLLT, - int iterations, + 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, + LLT_ROF_CPU_main(&inputData[0,0,0], &outputData[0,0,0], &infovec[0], regularisation_parameterROF, regularisation_parameterLLT, iterations, time_marching_parameter, - tolerance_param, + tolerance_param, dims[2], dims[1], dims[0]) return (outputData,infovec) #***************************************************************# #***************** Total Generalised Variation *****************# #***************************************************************# -def TGV_CPU(inputData, regularisation_parameter, alpha1, alpha0, iterations, LipshitzConst): +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) + 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) + 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, +def TGV_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, float regularisation_parameter, float alpha1, float alpha0, - int iterationsNumb, - float LipshitzConst): - + 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], regularisation_parameter, + TGV_main(&inputData[0,0], &outputData[0,0], &infovec[0], regularisation_parameter, alpha1, alpha0, - iterationsNumb, + iterationsNumb, LipshitzConst, + tolerance_param, dims[1],dims[0],1) - return outputData -def TGV_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, + 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): - + 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], regularisation_parameter, + TGV_main(&inputData[0,0,0], &outputData[0,0,0], &infovec[0], regularisation_parameter, alpha1, alpha0, - iterationsNumb, + iterationsNumb, LipshitzConst, + tolerance_param, dims[2], dims[1], dims[0]) - return outputData - + return (outputData,infovec) #****************************************************************# -#**************Directional Total-variation FGP ******************# +#***************Nonlinear (Isotropic) Diffusion******************# #****************************************************************# -#******** Directional TV Fast-Gradient-Projection (FGP)*********# -def dTV_FGP_CPU(inputData, refdata, regularisation_parameter, iterationsNumb, tolerance_param, eta_const, methodTV, nonneg, printM): +def NDF_CPU(inputData, regularisation_parameter, edge_parameter, iterationsNumb,time_marching_parameter, penalty_type,tolerance_param): if inputData.ndim == 2: - return dTV_FGP_2D(inputData, refdata, regularisation_parameter, iterationsNumb, tolerance_param, eta_const, methodTV, nonneg, printM) + return NDF_2D(inputData, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type, tolerance_param) elif inputData.ndim == 3: - return dTV_FGP_3D(inputData, refdata, regularisation_parameter, iterationsNumb, tolerance_param, eta_const, methodTV, nonneg, printM) + return NDF_3D(inputData, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type, tolerance_param) -def dTV_FGP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, - np.ndarray[np.float32_t, ndim=2, mode="c"] refdata, +def NDF_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, float regularisation_parameter, - int iterationsNumb, - float tolerance_param, - float eta_const, - int methodTV, - int nonneg, - int printM): - + 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') - - #/* Run FGP-dTV iterations for 2D data */ - dTV_FGP_CPU_main(&inputData[0,0], &refdata[0,0], &outputData[0,0], regularisation_parameter, - iterationsNumb, - tolerance_param, - eta_const, - methodTV, - nonneg, - printM, - dims[1], dims[0], 1) - - return outputData - -def dTV_FGP_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, - np.ndarray[np.float32_t, ndim=3, mode="c"] refdata, + 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, - int iterationsNumb, - float tolerance_param, - float eta_const, - int methodTV, - int nonneg, - int printM): + 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') - - #/* Run FGP-dTV iterations for 3D data */ - dTV_FGP_CPU_main(&inputData[0,0,0], &refdata[0,0,0], &outputData[0,0,0], regularisation_parameter, - iterationsNumb, - tolerance_param, - eta_const, - methodTV, - nonneg, - printM, - dims[2], dims[1], dims[0]) - return 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) #****************************************************************# -#*********************Total Nuclear Variation********************# +#*************Anisotropic Fourth-Order diffusion*****************# #****************************************************************# -def TNV_CPU(inputData, regularisation_parameter, iterationsNumb, tolerance_param): +def Diff4th_CPU(inputData, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter,tolerance_param): if inputData.ndim == 2: - return + return Diff4th_2D(inputData, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter,tolerance_param) elif inputData.ndim == 3: - return TNV_3D(inputData, regularisation_parameter, iterationsNumb, tolerance_param) + 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] -def TNV_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, + 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') - - # 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 + 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) #****************************************************************# -#***************Nonlinear (Isotropic) Diffusion******************# +#**************Directional Total-variation FGP ******************# #****************************************************************# -def NDF_CPU(inputData, regularisation_parameter, edge_parameter, iterationsNumb,time_marching_parameter, penalty_type): +#******** 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 NDF_2D(inputData, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type) + return dTV_FGP_2D(inputData, refdata, regularisation_parameter, iterationsNumb, tolerance_param, eta_const, methodTV, nonneg) elif inputData.ndim == 3: - return NDF_3D(inputData, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type) + return dTV_FGP_3D(inputData, refdata, regularisation_parameter, iterationsNumb, tolerance_param, eta_const, methodTV, nonneg) -def NDF_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, +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, - float edge_parameter, - int iterationsNumb, - float time_marching_parameter, - int penalty_type): + 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') - - # Run Nonlinear Diffusion iterations for 2D data - Diffusion_CPU_main(&inputData[0,0], &outputData[0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type, dims[1], dims[0], 1) - return outputData - -def NDF_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, + 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, - float edge_parameter, - int iterationsNumb, - float time_marching_parameter, - int penalty_type): + 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') - - # Run Nonlinear Diffusion iterations for 3D data - Diffusion_CPU_main(&inputData[0,0,0], &outputData[0,0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type, dims[2], dims[1], dims[0]) + 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') - return outputData + #/* 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) #****************************************************************# -#*************Anisotropic Fourth-Order diffusion*****************# +#*********************Total Nuclear Variation********************# #****************************************************************# -def Diff4th_CPU(inputData, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter): +def TNV_CPU(inputData, regularisation_parameter, iterationsNumb, tolerance_param): if inputData.ndim == 2: - return Diff4th_2D(inputData, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter) + return elif inputData.ndim == 3: - return Diff4th_3D(inputData, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter) + return TNV_3D(inputData, regularisation_parameter, iterationsNumb, tolerance_param) -def Diff4th_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, +def TNV_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, float regularisation_parameter, - float edge_parameter, - int iterationsNumb, - float time_marching_parameter): - 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') - - # Run Anisotropic Fourth-Order diffusion for 2D data - Diffus4th_CPU_main(&inputData[0,0], &outputData[0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, dims[1], dims[0], 1) - return outputData - -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') - - # Run Anisotropic Fourth-Order diffusion for 3D data - Diffus4th_CPU_main(&inputData[0,0,0], &outputData[0,0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, dims[2], dims[1], dims[0]) + # 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******************# #****************************************************************# @@ -511,14 +544,14 @@ def PatchSel_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 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') @@ -536,16 +569,16 @@ def PatchSel_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 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') @@ -573,10 +606,10 @@ def NLTV_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, 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) return outputData @@ -590,7 +623,7 @@ def NDF_INPAINT_CPU(inputData, maskData, regularisation_parameter, edge_paramete elif inputData.ndim == 3: return NDF_INP_3D(inputData, maskData, regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type) -def NDF_INP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, +def NDF_INP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, np.ndarray[np.uint8_t, ndim=2, mode="c"] maskData, float regularisation_parameter, float edge_parameter, @@ -605,12 +638,12 @@ def NDF_INP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, cdef np.ndarray[np.float32_t, ndim=2, mode="c"] outputData = \ np.zeros([dims[0],dims[1]], dtype='float32') - - # Run Inpaiting by Diffusion iterations for 2D data + + # Run Inpaiting by Diffusion iterations for 2D data Diffusion_Inpaint_CPU_main(&inputData[0,0], &maskData[0,0], &outputData[0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type, dims[1], dims[0], 1) return outputData - -def NDF_INP_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, + +def NDF_INP_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, np.ndarray[np.uint8_t, ndim=3, mode="c"] maskData, float regularisation_parameter, float edge_parameter, @@ -621,11 +654,11 @@ def NDF_INP_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 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 Inpaiting by Diffusion iterations for 3D data + + # Run Inpaiting by Diffusion iterations for 3D data Diffusion_Inpaint_CPU_main(&inputData[0,0,0], &maskData[0,0,0], &outputData[0,0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type, dims[2], dims[1], dims[0]) return outputData @@ -636,27 +669,27 @@ def NVM_INPAINT_CPU(inputData, maskData, SW_increment, iterationsNumb): if inputData.ndim == 2: return NVM_INP_2D(inputData, maskData, SW_increment, iterationsNumb) elif inputData.ndim == 3: - return + return -def NVM_INP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, +def NVM_INP_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, np.ndarray[np.uint8_t, ndim=2, mode="c"] maskData, int SW_increment, int iterationsNumb): 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') - + np.zeros([dims[0],dims[1]], dtype='float32') + cdef np.ndarray[np.uint8_t, ndim=2, mode="c"] maskData_upd = \ np.zeros([dims[0],dims[1]], dtype='uint8') - - # Run Inpaiting by Nonlocal vertical marching method for 2D data - NonlocalMarching_Inpaint_main(&inputData[0,0], &maskData[0,0], &outputData[0,0], + + # Run Inpaiting by Nonlocal vertical marching method for 2D data + NonlocalMarching_Inpaint_main(&inputData[0,0], &maskData[0,0], &outputData[0,0], &maskData_upd[0,0], SW_increment, iterationsNumb, 1, dims[1], dims[0], 1) - + return (outputData, maskData_upd) @@ -669,36 +702,36 @@ def TV_ENERGY(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, +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 + + # 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, + 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]) diff --git a/src/Python/src/gpu_regularisers.pyx b/src/Python/src/gpu_regularisers.pyx index 6fc4135..84ee981 100644 --- a/src/Python/src/gpu_regularisers.pyx +++ b/src/Python/src/gpu_regularisers.pyx @@ -24,68 +24,68 @@ cdef extern int TV_ROF_GPU_main(float* Input, float* Output, float *infovector, cdef extern int TV_FGP_GPU_main(float *Input, float *Output, float *infovector, float lambdaPar, int iter, float epsil, int methodTV, int nonneg, int N, int M, int Z); cdef extern int TV_SB_GPU_main(float *Input, float *Output, float *infovector, float lambdaPar, int iter, float epsil, int methodTV, int N, int M, int Z); cdef extern int LLT_ROF_GPU_main(float *Input, float *Output, float *infovector, float lambdaROF, float lambdaLLT, int iterationsNumb, float tau, float epsil, int N, int M, int Z); -cdef extern int TGV_GPU_main(float *Input, float *Output, float lambdaPar, float alpha1, float alpha0, int iterationsNumb, float L2, int dimX, int dimY, int dimZ); -cdef extern int NonlDiff_GPU_main(float *Input, float *Output, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int penaltytype, int N, int M, int Z); -cdef extern int dTV_FGP_GPU_main(float *Input, float *InputRef, float *Output, float lambdaPar, int iterationsNumb, float epsil, float eta, int methodTV, int nonneg, int printM, int N, int M, int Z); -cdef extern int Diffus4th_GPU_main(float *Input, float *Output, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int N, int M, int Z); +cdef extern int TGV_GPU_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 int NonlDiff_GPU_main(float *Input, float *Output, float *infovector, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, int penaltytype, float epsil, int N, int M, int Z); +cdef extern int Diffus4th_GPU_main(float *Input, float *Output, float *infovector, float lambdaPar, float sigmaPar, int iterationsNumb, float tau, float epsil, int N, int M, int Z); +cdef extern int dTV_FGP_GPU_main(float *Input, float *InputRef, float *Output, float *infovector, float lambdaPar, int iterationsNumb, float epsil, float eta, int methodTV, int nonneg, int N, int M, int Z); cdef extern int PatchSelect_GPU_main(float *Input, unsigned short *H_i, unsigned short *H_j, float *Weights, int N, int M, int SearchWindow, int SimilarWin, int NumNeighb, float h); # Total-variation Rudin-Osher-Fatemi (ROF) def TV_ROF_GPU(inputData, regularisation_parameter, - iterations, + iterations, time_marching_parameter, tolerance_param): if inputData.ndim == 2: - return ROFTV2D(inputData, + return ROFTV2D(inputData, regularisation_parameter, iterations, time_marching_parameter, tolerance_param) elif inputData.ndim == 3: - return ROFTV3D(inputData, + return ROFTV3D(inputData, regularisation_parameter, - iterations, + iterations, time_marching_parameter, tolerance_param) - + # Total-variation Fast-Gradient-Projection (FGP) def TV_FGP_GPU(inputData, regularisation_parameter, - iterations, + iterations, tolerance_param, methodTV, nonneg): if inputData.ndim == 2: return FGPTV2D(inputData, regularisation_parameter, - iterations, + iterations, tolerance_param, methodTV, nonneg) elif inputData.ndim == 3: return FGPTV3D(inputData, regularisation_parameter, - iterations, + iterations, tolerance_param, methodTV, nonneg) # Total-variation Split Bregman (SB) def TV_SB_GPU(inputData, regularisation_parameter, - iterations, + iterations, tolerance_param, methodTV): if inputData.ndim == 2: return SBTV2D(inputData, regularisation_parameter, - iterations, + iterations, tolerance_param, methodTV) elif inputData.ndim == 3: return SBTV3D(inputData, regularisation_parameter, - iterations, + iterations, tolerance_param, methodTV) # LLT-ROF model @@ -95,90 +95,93 @@ def LLT_ROF_GPU(inputData, regularisation_parameterROF, regularisation_parameter elif inputData.ndim == 3: return LLT_ROF_GPU3D(inputData, regularisation_parameterROF, regularisation_parameterLLT, iterations, time_marching_parameter, tolerance_param) # Total Generilised Variation (TGV) -def TGV_GPU(inputData, regularisation_parameter, alpha1, alpha0, iterations, LipshitzConst): +def TGV_GPU(inputData, regularisation_parameter, alpha1, alpha0, iterations, LipshitzConst, tolerance_param): if inputData.ndim == 2: - return TGV2D(inputData, regularisation_parameter, alpha1, alpha0, iterations, LipshitzConst) + return TGV2D(inputData, regularisation_parameter, alpha1, alpha0, iterations, LipshitzConst, tolerance_param) elif inputData.ndim == 3: - return TGV3D(inputData, regularisation_parameter, alpha1, alpha0, iterations, LipshitzConst) + return TGV3D(inputData, regularisation_parameter, alpha1, alpha0, iterations, LipshitzConst, tolerance_param) # Directional Total-variation Fast-Gradient-Projection (FGP) def dTV_FGP_GPU(inputData, refdata, regularisation_parameter, - iterations, + iterations, tolerance_param, eta_const, methodTV, - nonneg, - printM): + nonneg): if inputData.ndim == 2: return FGPdTV2D(inputData, refdata, regularisation_parameter, - iterations, + iterations, tolerance_param, eta_const, methodTV, - nonneg, - printM) + nonneg) elif inputData.ndim == 3: return FGPdTV3D(inputData, refdata, regularisation_parameter, - iterations, + iterations, tolerance_param, eta_const, methodTV, - nonneg, - printM) + nonneg) # Nonlocal Isotropic Diffusion (NDF) def NDF_GPU(inputData, regularisation_parameter, edge_parameter, - iterations, + iterations, time_marching_parameter, - penalty_type): + penalty_type, + tolerance_param): if inputData.ndim == 2: return NDF_GPU_2D(inputData, regularisation_parameter, edge_parameter, - iterations, + iterations, time_marching_parameter, - penalty_type) + penalty_type, + tolerance_param) elif inputData.ndim == 3: return NDF_GPU_3D(inputData, regularisation_parameter, edge_parameter, - iterations, + iterations, time_marching_parameter, - penalty_type) + penalty_type, + tolerance_param) # Anisotropic Fourth-Order diffusion def Diff4th_GPU(inputData, regularisation_parameter, edge_parameter, - iterations, - time_marching_parameter): + iterations, + time_marching_parameter, + tolerance_param): if inputData.ndim == 2: return Diff4th_2D(inputData, regularisation_parameter, edge_parameter, - iterations, - time_marching_parameter) + iterations, + time_marching_parameter, + tolerance_param) elif inputData.ndim == 3: return Diff4th_3D(inputData, regularisation_parameter, edge_parameter, - iterations, - time_marching_parameter) - + iterations, + time_marching_parameter, + tolerance_param) + #****************************************************************# #********************** Total-variation ROF *********************# #****************************************************************# -def ROFTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, +def ROFTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, float regularisation_parameter, - int iterations, + int iterations, float time_marching_parameter, float tolerance_param): - + cdef long dims[2] dims[0] = inputData.shape[0] dims[1] = inputData.shape[1] @@ -187,25 +190,25 @@ def ROFTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, np.zeros([dims[0],dims[1]], dtype='float32') cdef np.ndarray[np.float32_t, ndim=1, mode="c"] infovec = \ np.ones([2], dtype='float32') - + # Running CUDA code here if (TV_ROF_GPU_main( &inputData[0,0], &outputData[0,0], &infovec[0], regularisation_parameter, - iterations, - time_marching_parameter, + iterations, + time_marching_parameter, tolerance_param, dims[1], dims[0], 1)==0): return (outputData,infovec) else: raise ValueError(CUDAErrorMessage); - -def ROFTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, + +def ROFTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, float regularisation_parameter, - int iterations, + int iterations, float time_marching_parameter, float tolerance_param): - + cdef long dims[3] dims[0] = inputData.shape[0] dims[1] = inputData.shape[1] @@ -215,13 +218,13 @@ def ROFTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 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') - - # Running CUDA code here + + # Running CUDA code here if (TV_ROF_GPU_main( &inputData[0,0,0], &outputData[0,0,0], &infovec[0], regularisation_parameter, - iterations, - time_marching_parameter, + iterations, + time_marching_parameter, tolerance_param, dims[2], dims[1], dims[0])==0): return (outputData,infovec) @@ -231,13 +234,13 @@ def ROFTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, #********************** Total-variation FGP *********************# #****************************************************************# #******** Total-variation Fast-Gradient-Projection (FGP)*********# -def FGPTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, +def FGPTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, float regularisation_parameter, - int iterations, + int iterations, float tolerance_param, int methodTV, int nonneg): - + cdef long dims[2] dims[0] = inputData.shape[0] dims[1] = inputData.shape[1] @@ -246,10 +249,10 @@ def FGPTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, np.zeros([dims[0],dims[1]], dtype='float32') cdef np.ndarray[np.float32_t, ndim=1, mode="c"] infovec = \ np.ones([2], dtype='float32') - - # Running CUDA code here + + # Running CUDA code here if (TV_FGP_GPU_main(&inputData[0,0], &outputData[0,0], &infovec[0], - regularisation_parameter, + regularisation_parameter, iterations, tolerance_param, methodTV, @@ -258,14 +261,14 @@ def FGPTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, return (outputData,infovec) else: raise ValueError(CUDAErrorMessage); - -def FGPTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, + +def FGPTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, float regularisation_parameter, - int iterations, + int iterations, float tolerance_param, int methodTV, int nonneg): - + cdef long dims[3] dims[0] = inputData.shape[0] dims[1] = inputData.shape[1] @@ -276,11 +279,11 @@ def FGPTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 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') - - # Running CUDA code here + + # Running CUDA code here if (TV_FGP_GPU_main(&inputData[0,0,0], &outputData[0,0,0], &infovec[0], - regularisation_parameter, - iterations, + regularisation_parameter, + iterations, tolerance_param, methodTV, nonneg, @@ -293,12 +296,12 @@ def FGPTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, #********************** Total-variation SB *********************# #***************************************************************# #*************** Total-variation Split Bregman (SB)*************# -def SBTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, +def SBTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, float regularisation_parameter, - int iterations, + int iterations, float tolerance_param, int methodTV): - + cdef long dims[2] dims[0] = inputData.shape[0] dims[1] = inputData.shape[1] @@ -307,11 +310,11 @@ def SBTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, np.zeros([dims[0],dims[1]], dtype='float32') cdef np.ndarray[np.float32_t, ndim=1, mode="c"] infovec = \ np.ones([2], dtype='float32') - + # Running CUDA code here if (TV_SB_GPU_main(&inputData[0,0], &outputData[0,0],&infovec[0], - regularisation_parameter, - iterations, + regularisation_parameter, + iterations, tolerance_param, methodTV, dims[1], dims[0], 1)==0): @@ -319,13 +322,13 @@ def SBTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, else: raise ValueError(CUDAErrorMessage); - -def SBTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, + +def SBTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, float regularisation_parameter, - int iterations, + int iterations, float tolerance_param, int methodTV): - + cdef long dims[3] dims[0] = inputData.shape[0] dims[1] = inputData.shape[1] @@ -335,11 +338,11 @@ def SBTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 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') - - # Running CUDA code here + + # Running CUDA code here if (TV_SB_GPU_main(&inputData[0,0,0], &outputData[0,0,0],&infovec[0], - regularisation_parameter , - iterations, + regularisation_parameter , + iterations, tolerance_param, methodTV, dims[2], dims[1], dims[0])==0): @@ -352,13 +355,13 @@ def SBTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, #************************ LLT-ROF model ************************# #***************************************************************# #************Joint LLT-ROF model for higher order **************# -def LLT_ROF_GPU2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, +def LLT_ROF_GPU2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, float regularisation_parameterROF, float regularisation_parameterLLT, - int iterations, + int iterations, float time_marching_parameter, float tolerance_param): - + cdef long dims[2] dims[0] = inputData.shape[0] dims[1] = inputData.shape[1] @@ -367,24 +370,24 @@ def LLT_ROF_GPU2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, np.zeros([dims[0],dims[1]], dtype='float32') cdef np.ndarray[np.float32_t, ndim=1, mode="c"] infovec = \ np.ones([2], dtype='float32') - - # Running CUDA code here - if (LLT_ROF_GPU_main(&inputData[0,0], &outputData[0,0],&infovec[0],regularisation_parameterROF, regularisation_parameterLLT, iterations, - time_marching_parameter, + + # Running CUDA code here + if (LLT_ROF_GPU_main(&inputData[0,0], &outputData[0,0],&infovec[0],regularisation_parameterROF, regularisation_parameterLLT, iterations, + time_marching_parameter, tolerance_param, dims[1],dims[0],1)==0): return (outputData,infovec) else: raise ValueError(CUDAErrorMessage); - -def LLT_ROF_GPU3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, + +def LLT_ROF_GPU3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, float regularisation_parameterROF, float regularisation_parameterLLT, - int iterations, + int iterations, float time_marching_parameter, float tolerance_param): - + cdef long dims[3] dims[0] = inputData.shape[0] dims[1] = inputData.shape[1] @@ -394,11 +397,11 @@ def LLT_ROF_GPU3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, 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') - + # Running CUDA code here - if (LLT_ROF_GPU_main(&inputData[0,0,0], &outputData[0,0,0], &infovec[0], regularisation_parameterROF, regularisation_parameterLLT, - iterations, - time_marching_parameter, + if (LLT_ROF_GPU_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])==0): return (outputData,infovec) @@ -409,38 +412,43 @@ def LLT_ROF_GPU3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, #***************************************************************# #***************** Total Generalised Variation *****************# #***************************************************************# -def TGV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, +def TGV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, float regularisation_parameter, float alpha1, float alpha0, - int iterationsNumb, - float LipshitzConst): - + 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.ones([2], dtype='float32') + #/* Run TGV iterations for 2D data */ - if (TGV_GPU_main(&inputData[0,0], &outputData[0,0], regularisation_parameter, + if (TGV_GPU_main(&inputData[0,0], &outputData[0,0], &infovec[0], regularisation_parameter, alpha1, alpha0, - iterationsNumb, + iterationsNumb, LipshitzConst, + tolerance_param, dims[1],dims[0], 1)==0): - return outputData + return (outputData,infovec) else: raise ValueError(CUDAErrorMessage); -def TGV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, +def TGV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, float regularisation_parameter, float alpha1, float alpha0, - int iterationsNumb, - float LipshitzConst): - + int iterationsNumb, + float LipshitzConst, + float tolerance_param): + cdef long dims[3] dims[0] = inputData.shape[0] dims[1] = inputData.shape[1] @@ -448,178 +456,205 @@ def TGV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, cdef np.ndarray[np.float32_t, ndim=3, mode="c"] outputData = \ np.zeros([dims[0],dims[1],dims[2]], dtype='float32') - - # Running CUDA code here + cdef np.ndarray[np.float32_t, ndim=1, mode="c"] infovec = \ + np.ones([2], dtype='float32') + + # Running CUDA code here if (TGV_GPU_main( - &inputData[0,0,0], &outputData[0,0,0], regularisation_parameter, + &inputData[0,0,0], &outputData[0,0,0], &infovec[0], regularisation_parameter, alpha1, alpha0, - iterationsNumb, + iterationsNumb, LipshitzConst, + tolerance_param, dims[2], dims[1], dims[0])==0): - return outputData; + return (outputData,infovec) else: raise ValueError(CUDAErrorMessage); - #****************************************************************# -#**************Directional Total-variation FGP ******************# +#***************Nonlinear (Isotropic) Diffusion******************# #****************************************************************# -#******** Directional TV Fast-Gradient-Projection (FGP)*********# -def FGPdTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, - np.ndarray[np.float32_t, ndim=2, mode="c"] refdata, +def NDF_GPU_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, float regularisation_parameter, - int iterations, - float tolerance_param, - float eta_const, - int methodTV, - int nonneg, - int printM): - + 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') - - # Running CUDA code here - if (dTV_FGP_GPU_main(&inputData[0,0], &refdata[0,0], &outputData[0,0], - regularisation_parameter, - iterations, - tolerance_param, - eta_const, - methodTV, - nonneg, - printM, - dims[1], dims[0], 1)==0): - return 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') + + #rangecheck = penalty_type < 1 and penalty_type > 3 + #if not rangecheck: +# raise ValueError('Choose penalty type as 1 for Huber, 2 - Perona-Malik, 3 - Tukey Biweight') + + # Run Nonlinear Diffusion iterations for 2D data + # Running CUDA code here + if (NonlDiff_GPU_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)==0): + return (outputData,infovec) else: raise ValueError(CUDAErrorMessage); - - -def FGPdTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, - np.ndarray[np.float32_t, ndim=3, mode="c"] refdata, +def NDF_GPU_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, float regularisation_parameter, - int iterations, - float tolerance_param, - float eta_const, - int methodTV, - int nonneg, - int printM): - + 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') - - # Running CUDA code here - if (dTV_FGP_GPU_main(&inputData[0,0,0], &refdata[0,0,0], &outputData[0,0,0], - regularisation_parameter , - iterations, - tolerance_param, - eta_const, - methodTV, - nonneg, - printM, - dims[2], dims[1], dims[0])==0): - return 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 Nonlinear Diffusion iterations for 3D data + # Running CUDA code here + if (NonlDiff_GPU_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])==0): + return (outputData,infovec) else: raise ValueError(CUDAErrorMessage); - #****************************************************************# -#***************Nonlinear (Isotropic) Diffusion******************# +#************Anisotropic Fourth-Order diffusion******************# #****************************************************************# -def NDF_GPU_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, +def Diff4th_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, float regularisation_parameter, float edge_parameter, - int iterationsNumb, + 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') - - #rangecheck = penalty_type < 1 and penalty_type > 3 - #if not rangecheck: -# raise ValueError('Choose penalty type as 1 for Huber, 2 - Perona-Malik, 3 - Tukey Biweight') - - # Run Nonlinear Diffusion iterations for 2D data - # Running CUDA code here - if (NonlDiff_GPU_main(&inputData[0,0], &outputData[0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type, dims[1], dims[0], 1)==0): - return outputData; + cdef np.ndarray[np.float32_t, ndim=1, mode="c"] infovec = \ + np.ones([2], dtype='float32') + + # Run Anisotropic Fourth-Order diffusion for 2D data + # Running CUDA code here + if (Diffus4th_GPU_main(&inputData[0,0], &outputData[0,0], &infovec[0], + regularisation_parameter, edge_parameter, iterationsNumb, + time_marching_parameter, + tolerance_param, + dims[1], dims[0], 1)==0): + return (outputData,infovec) else: raise ValueError(CUDAErrorMessage); - -def NDF_GPU_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, +def Diff4th_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, float regularisation_parameter, float edge_parameter, - int iterationsNumb, + 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') - - # Run Nonlinear Diffusion iterations for 3D data - # Running CUDA code here - if (NonlDiff_GPU_main(&inputData[0,0,0], &outputData[0,0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, penalty_type, dims[2], dims[1], dims[0])==0): - return 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 Anisotropic Fourth-Order diffusion for 3D data + # Running CUDA code here + if (Diffus4th_GPU_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])==0): + return (outputData,infovec) else: raise ValueError(CUDAErrorMessage); - #****************************************************************# -#************Anisotropic Fourth-Order diffusion******************# +#**************Directional Total-variation FGP ******************# #****************************************************************# -def Diff4th_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, +#******** Directional TV Fast-Gradient-Projection (FGP)*********# +def FGPdTV2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, + np.ndarray[np.float32_t, ndim=2, mode="c"] refdata, float regularisation_parameter, - float edge_parameter, - int iterationsNumb, - float time_marching_parameter): + int iterations, + 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') - - # Run Anisotropic Fourth-Order diffusion for 2D data - # Running CUDA code here - if (Diffus4th_GPU_main(&inputData[0,0], &outputData[0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, dims[1], dims[0], 1)==0): - return 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') + + # Running CUDA code here + if (dTV_FGP_GPU_main(&inputData[0,0], &refdata[0,0], &outputData[0,0], &infovec[0], + regularisation_parameter, + iterations, + tolerance_param, + eta_const, + methodTV, + nonneg, + dims[1], dims[0], 1)==0): + return (outputData,infovec) else: raise ValueError(CUDAErrorMessage); - -def Diff4th_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, + +def FGPdTV3D(np.ndarray[np.float32_t, ndim=3, mode="c"] inputData, + np.ndarray[np.float32_t, ndim=3, mode="c"] refdata, float regularisation_parameter, - float edge_parameter, - int iterationsNumb, - float time_marching_parameter): + int iterations, + 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') - - # Run Anisotropic Fourth-Order diffusion for 3D data - # Running CUDA code here - if (Diffus4th_GPU_main(&inputData[0,0,0], &outputData[0,0,0], regularisation_parameter, edge_parameter, iterationsNumb, time_marching_parameter, dims[2], dims[1], dims[0])==0): - return 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') + + # Running CUDA code here + if (dTV_FGP_GPU_main(&inputData[0,0,0], &refdata[0,0,0], &outputData[0,0,0], &infovec[0], + regularisation_parameter , + iterations, + tolerance_param, + eta_const, + methodTV, + nonneg, + dims[2], dims[1], dims[0])==0): + return (outputData,infovec) else: raise ValueError(CUDAErrorMessage); @@ -639,14 +674,14 @@ def PatchSel_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, cdef long dims[3] dims[0] = neighbours dims[1] = inputData.shape[0] - dims[2] = inputData.shape[1] - + 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') @@ -655,4 +690,3 @@ def PatchSel_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] inputData, return H_i, H_j, Weights; else: raise ValueError(CUDAErrorMessage); - |