From 494a857b830fce5e786dfc058f68bf78d9673ba6 Mon Sep 17 00:00:00 2001 From: dkazanc Date: Tue, 26 Nov 2019 14:51:34 +0000 Subject: PDTV 3D verision --- demos/demo_cpu_regularisers.py | 2 +- demos/demo_cpu_regularisers3D.py | 51 ++++++++++- src/Core/CMakeLists.txt | 2 +- src/Core/regularisers_CPU/FGP_TV_core.c | 34 -------- src/Core/regularisers_CPU/FGP_TV_core.h | 1 - src/Core/regularisers_CPU/FGP_dTV_core.c | 141 +++++++++---------------------- src/Core/regularisers_CPU/FGP_dTV_core.h | 2 - src/Core/regularisers_CPU/PD_TV_core.c | 95 +++++++++++++++++++-- src/Core/regularisers_CPU/PD_TV_core.h | 5 +- src/Core/regularisers_CPU/utils.c | 37 +++++++- src/Core/regularisers_CPU/utils.h | 1 + src/Python/setup-regularisers.py.in | 1 - src/Python/src/cpu_regularisers.pyx | 29 ++++++- test/test_run_test.py | 2 +- 14 files changed, 246 insertions(+), 157 deletions(-) diff --git a/demos/demo_cpu_regularisers.py b/demos/demo_cpu_regularisers.py index d0bbe63..9888743 100644 --- a/demos/demo_cpu_regularisers.py +++ b/demos/demo_cpu_regularisers.py @@ -179,7 +179,7 @@ pars = {'algorithm' : PD_TV, \ 'tolerance_constant':1e-06,\ 'methodTV': 0 ,\ 'nonneg': 1, - 'lipschitz_const' : 12} + 'lipschitz_const' : 6} print ("#############PD TV CPU####################") start_time = timeit.default_timer() diff --git a/demos/demo_cpu_regularisers3D.py b/demos/demo_cpu_regularisers3D.py index fc1e8e6..349300f 100644 --- a/demos/demo_cpu_regularisers3D.py +++ b/demos/demo_cpu_regularisers3D.py @@ -12,7 +12,7 @@ import matplotlib.pyplot as plt import numpy as np import os import timeit -from ccpi.filters.regularisers import ROF_TV, FGP_TV, SB_TV, TGV, LLT_ROF, FGP_dTV, NDF, Diff4th +from ccpi.filters.regularisers import ROF_TV, FGP_TV, PD_TV, SB_TV, TGV, LLT_ROF, FGP_dTV, NDF, Diff4th from ccpi.supp.qualitymetrics import QualityTools ############################################################################### def printParametersToString(pars): @@ -30,8 +30,7 @@ def printParametersToString(pars): return txt ############################################################################### -# filename = os.path.join( "data" ,"lena_gray_512.tif") -filename = "/home/algol/Documents/DEV/CCPi-Regularisation-Toolkit/test/lena_gray_512.tif" +filename = os.path.join( "data" ,"lena_gray_512.tif") # read image Im = plt.imread(filename) @@ -169,7 +168,53 @@ a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, verticalalignment='top', bbox=props) imgplot = plt.imshow(fgp_cpu3D[10,:,:], cmap="gray") plt.title('{}'.format('Recovered volume on the CPU using FGP-TV')) +#%% +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +print ("_______________PD-TV (3D)__________________") +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") + +## plot +fig = plt.figure() +plt.suptitle('Performance of PD-TV regulariser using the CPU') +a=fig.add_subplot(1,2,1) +a.set_title('Noisy Image') +imgplot = plt.imshow(noisyVol[10,:,:],cmap="gray") +# set parameters +pars = {'algorithm' : PD_TV, \ + 'input' : noisyVol,\ + 'regularisation_parameter':0.02, \ + 'number_of_iterations' :1000 ,\ + 'tolerance_constant': 1e-06,\ + 'methodTV': 0 ,\ + 'lipschitz_const' : 12,\ + 'nonneg': 0} + +print ("#############FGP TV GPU####################") +start_time = timeit.default_timer() +(pd_cpu3D,info_vec_cpu) = PD_TV(pars['input'], + pars['regularisation_parameter'], + pars['number_of_iterations'], + pars['tolerance_constant'], + pars['methodTV'], + pars['nonneg'], + pars['lipschitz_const'], 'cpu') + +Qtools = QualityTools(idealVol, pd_cpu3D) +pars['rmse'] = Qtools.rmse() + +txtstr = printParametersToString(pars) +txtstr += "%s = %.3fs" % ('elapsed time',timeit.default_timer() - start_time) +print (txtstr) +a=fig.add_subplot(1,2,2) + +# these are matplotlib.patch.Patch properties +props = dict(boxstyle='round', facecolor='wheat', alpha=0.75) +# place a text box in upper left in axes coords +a.text(0.15, 0.25, txtstr, transform=a.transAxes, fontsize=14, + verticalalignment='top', bbox=props) +imgplot = plt.imshow(pd_cpu3D[10,:,:], cmap="gray") +plt.title('{}'.format('Recovered volume on the CPU using PD-TV')) #%% print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") print ("_______________SB-TV (3D)_________________") diff --git a/src/Core/CMakeLists.txt b/src/Core/CMakeLists.txt index ac7ec91..c1bd7bb 100644 --- a/src/Core/CMakeLists.txt +++ b/src/Core/CMakeLists.txt @@ -66,7 +66,7 @@ message("Adding regularisers as a shared library") add_library(cilreg SHARED ${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/FGP_TV_core.c ${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/SB_TV_core.c - ${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/PD_TV_core.c + ${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/PD_TV_core.c ${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/TGV_core.c ${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/Diffusion_core.c ${CMAKE_CURRENT_SOURCE_DIR}/regularisers_CPU/Diffus4th_order_core.c diff --git a/src/Core/regularisers_CPU/FGP_TV_core.c b/src/Core/regularisers_CPU/FGP_TV_core.c index ff67af2..12f2254 100644 --- a/src/Core/regularisers_CPU/FGP_TV_core.c +++ b/src/Core/regularisers_CPU/FGP_TV_core.c @@ -254,40 +254,6 @@ float Grad_func3D(float *P1, float *P2, float *P3, float *D, float *R1, float *R }}} return 1; } -float Proj_func3D(float *P1, float *P2, float *P3, int methTV, long DimTotal) -{ - float val1, val2, val3, denom, sq_denom; - long i; - if (methTV == 0) { - /* isotropic TV*/ -#pragma omp parallel for shared(P1,P2,P3) private(i,val1,val2,val3,sq_denom) - for(i=0; i 1.0f) { - sq_denom = 1.0f/sqrtf(denom); - P1[i] = P1[i]*sq_denom; - P2[i] = P2[i]*sq_denom; - P3[i] = P3[i]*sq_denom; - } - } - } - else { - /* anisotropic TV*/ -#pragma omp parallel for shared(P1,P2,P3) private(i,val1,val2,val3) - for(i=0; i 3) break; } } - + free(P3); free(P3_prev); free(R3); free(InputRef_z); } if (epsil != 0.0f) free(Output_prev); free(P1); free(P2); free(P1_prev); free(P2_prev); free(R1); free(R2); free(InputRef_x); free(InputRef_y); - + /*adding info into info_vector */ infovector[0] = (float)(ll); /*iterations number (if stopped earlier based on tolerance)*/ infovector[1] = re; /* reached tolerance */ - + return 0; } @@ -185,7 +185,6 @@ float dTV_FGP_CPU_main(float *Input, float *InputRef, float *Output, float *info /********************************************************************/ /***************************2D Functions*****************************/ /********************************************************************/ - float GradNorm_func2D(float *B, float *B_x, float *B_y, float eta, long dimX, long dimY) { long i,j,index; @@ -249,47 +248,17 @@ float Grad_dfunc2D(float *P1, float *P2, float *D, float *R1, float *R2, float * /* boundary conditions */ if (i == dimX-1) val1 = 0.0f; else val1 = D[index] - D[j*dimX + (i+1)]; if (j == dimY-1) val2 = 0.0f; else val2 = D[index] - D[(j+1)*dimX + i]; - + in_prod = val1*B_x[index] + val2*B_y[index]; /* calculate inner product */ val1 = val1 - in_prod*B_x[index]; val2 = val2 - in_prod*B_y[index]; - + P1[index] = R1[index] + multip*val1; P2[index] = R2[index] + multip*val2; - + }} return 1; } -float Proj_dfunc2D(float *P1, float *P2, int methTV, long DimTotal) -{ - float val1, val2, denom, sq_denom; - long i; - if (methTV == 0) { - /* isotropic TV*/ -#pragma omp parallel for shared(P1,P2) private(i,denom,sq_denom) - for(i=0; i 1.0f) { - sq_denom = 1.0f/sqrtf(denom); - P1[i] = P1[i]*sq_denom; - P2[i] = P2[i]*sq_denom; - } - } - } - else { - /* anisotropic TV*/ -#pragma omp parallel for shared(P1,P2) private(i,val1,val2) - for(i=0; i 1.0f) { - sq_denom = 1.0f/sqrtf(denom); - P1[i] = P1[i]*sq_denom; - P2[i] = P2[i]*sq_denom; - P3[i] = P3[i]*sq_denom; - } - } - } - else { - /* anisotropic TV*/ -#pragma omp parallel for shared(P1,P2,P3) private(i,val1,val2,val3) - for(i=0; i 3) break; + } + /*get updated solution*/ + + getX(U, U_old, theta, DimTotal); + } + free(P1); free(P2); free(P3); free(U_old); } /*adding info into info_vector */ infovector[0] = (float)(ll); /*iterations number (if stopped earlier based on tolerance)*/ @@ -134,7 +173,7 @@ float DivProj2D(float *U, float *Input, float *P1, float *P2, long dimX, long di { long i,j,index; float P_v1, P_v2, div_var; - #pragma omp parallel for shared(U,Input,P1,P2) private(i, j, P_v1, P_v2, div_var) + #pragma omp parallel for shared(U,Input,P1,P2) private(index, i, j, P_v1, P_v2, div_var) for(j=0; j 1.0f) { + sq_denom = 1.0f/sqrtf(denom); + P1[i] = P1[i]*sq_denom; + P2[i] = P2[i]*sq_denom; + P3[i] = P3[i]*sq_denom; + } + } + } + else { + /* anisotropic TV*/ +#pragma omp parallel for shared(P1,P2,P3) private(i,val1,val2,val3) + for(i=0; i