diff --git a/Demos/Python/DemoFISTA_artifacts2D.py b/Demos/Python/DemoFISTA_artifacts2D.py index a258c3d3f..6c20edb52 100644 --- a/Demos/Python/DemoFISTA_artifacts2D.py +++ b/Demos/Python/DemoFISTA_artifacts2D.py @@ -210,7 +210,7 @@ device='gpu') RecFISTA_Huber_reg = RectoolsIR.FISTA(noisy_zing_stripe, - huber_data_threshold=3.5, + huber_data_threshold=4.0, iterationsFISTA = 20, regularisation = 'ROF_TV', regularisation_parameter = 0.01, @@ -250,7 +250,7 @@ device='gpu') RecFISTA_Huber_RingModel = RectoolsIR.FISTA(noisy_zing_stripe, - huber_data_threshold=3.5, + huber_data_threshold=4.0, ring_model_horiz_size=9, iterationsFISTA = 20, regularisation = 'ROF_TV', diff --git a/Demos/Python/DemoFISTA_artifacts3D.py b/Demos/Python/DemoFISTA_artifacts3D.py index 0c5b1f3ac..9818bc596 100644 --- a/Demos/Python/DemoFISTA_artifacts3D.py +++ b/Demos/Python/DemoFISTA_artifacts3D.py @@ -26,8 +26,8 @@ print ("Building 3D phantom using TomoPhantom software") tic=timeit.default_timer() -model = 9 # select a model number from the library -N_size = 256 # Define phantom dimensions using a scalar value (cubic phantom) +model = 13 # select a model number from the library +N_size = 150 # Define phantom dimensions using a scalar value (cubic phantom) path = os.path.dirname(tomophantom.__file__) path_library3D = os.path.join(path, "Phantom3DLibrary.dat") #This will generate a N_size x N_size x N_size phantom (3D) @@ -63,7 +63,7 @@ projData3D_analyt= TomoP3D.ModelSino(model, N_size, Horiz_det, Vert_det, angles, path_library3D) intens_max = 70 -sliceSel = 150 +sliceSel = 100 plt.figure() plt.subplot(131) plt.imshow(projData3D_analyt[:,sliceSel,:],vmin=0, vmax=intens_max) @@ -105,15 +105,41 @@ #%% # initialise tomobar DIRECT reconstruction class ONCE from tomobar.methodsDIR import RecToolsDIR -RectoolsDIR = RecToolsDIR(DetectorsDimH = Horiz_det, # DetectorsDimH # detector dimension (horizontal) +Rectools = RecToolsDIR(DetectorsDimH = Horiz_det, # DetectorsDimH # detector dimension (horizontal) DetectorsDimV = Vert_det, # DetectorsDimV # detector dimension (vertical) for 3D case only - CenterRotOffset = None, # Center of Rotation (CoR) scalar (for 3D case only) + CenterRotOffset = 0.0, # Center of Rotation (CoR) scalar (for 3D case only) AnglesVec = angles_rad, # array of angles in radians ObjSize = N_size, # a scalar to define reconstructed object dimensions device = 'gpu') #%% +# testing RING +Forwproj = Rectools.FORWPROJ(phantom_tm) +residual = Forwproj - projData3D_norm +#residual = np.swapaxes(residual,0,1) +#%% +from tomobar.supp.addmodules import RING_WEIGHTS +ring_model_horiz_size = 13 +ring_model_vert_size=0 +ring_model_slices_size=0 + +offset_rings = RING_WEIGHTS(residual, ring_model_horiz_size, ring_model_vert_size, ring_model_slices_size) + +intens_max = 70 +sliceSel = 120 +plt.figure() +plt.subplot(131) +plt.imshow(offset_rings[:,sliceSel,:]) +plt.title('2D Projection (erroneous)') +plt.subplot(132) +plt.imshow(offset_rings[sliceSel,:,:]) +plt.title('Sinogram view') +plt.subplot(133) +plt.imshow(offset_rings[:,:,sliceSel]) +plt.title('Tangentogram view') +plt.show() +#%% print ("Reconstruction using FBP from tomobar") -recNumerical= RectoolsDIR.FBP(projData3D_norm) # FBP reconstruction +recNumerical= Rectools.FBP(projData3D_norm) # FBP reconstruction sliceSel = int(0.5*N_size) max_val = 1 @@ -150,7 +176,7 @@ datafidelity='LS',# data fidelity, choose LS, PWLS (wip), GH (wip), Student (wip) nonnegativity='ENABLE', # enable nonnegativity constraint (set to 'ENABLE') OS_number = 10, # the number of subsets, NONE/(or > 1) ~ classical / ordered subsets - tolerance = 1e-08, # tolerance to stop outer iterations earlier + tolerance = 1e-10, # tolerance to stop outer iterations earlier device='gpu') lc = RectoolsIR.powermethod() # calculate Lipschitz constant #%% @@ -165,6 +191,11 @@ regularisation_iterations = 300, lipschitz_const = lc) + +Qtools = QualityTools(phantom_tm, RecFISTA_reg) +RMSE_FISTA_TV = Qtools.rmse() +print("RMSE for FISTA-OS-TV is {}".format(RMSE_FISTA_TV)) + sliceSel = int(0.5*N_size) max_val = 1 plt.figure() @@ -190,6 +221,10 @@ regularisation_iterations = 300, lipschitz_const = lc) +Qtools = QualityTools(phantom_tm, RecFISTA_Huber_TV) +RMSE_FISTA_HUBER_TV = Qtools.rmse() +print("RMSE for FISTA-OS-Huber-TV is {}".format(RMSE_FISTA_HUBER_TV)) + sliceSel = int(0.5*N_size) max_val = 1 plt.figure() @@ -206,18 +241,20 @@ plt.title('3D Huber Rec, sagittal') plt.show() #%% -#%% # Run FISTA reconstrucion algorithm with 3D regularisation RecFISTA_HuberRING_TV = RectoolsIR.FISTA(projData3D_norm, iterationsFISTA = 20, huber_data_threshold = 1.0, ring_model_horiz_size= 7, - ring_model_vert_size = 1, regularisation = 'FGP_TV', regularisation_parameter = 0.00015, regularisation_iterations = 300, lipschitz_const = lc) +Qtools = QualityTools(phantom_tm, RecFISTA_HuberRING_TV) +RMSE_FISTA_HUBER_RING_TV = Qtools.rmse() +print("RMSE for FISTA-OS-Huber-Ring-TV is {}".format(RMSE_FISTA_HUBER_RING_TV)) + sliceSel = int(0.5*N_size) max_val = 1 plt.figure() diff --git a/src/Core/functions_CPU/RingWeights_core.c b/src/Core/functions_CPU/RingWeights_core.c index dbacc7052..c2d9d28ca 100644 --- a/src/Core/functions_CPU/RingWeights_core.c +++ b/src/Core/functions_CPU/RingWeights_core.c @@ -46,28 +46,28 @@ void swap(float *xp, float *yp) } -float RingWeights_main(float *residual, float *weights, int horiz_window_halfsize, int vert_window_halfsize, long anglesDim, long detectorsDim, long slices) +float RingWeights_main(float *residual, float *weights, int horiz_window_halfsize, int vert_window_halfsize, int slices_window_halfsize, long anglesDim, long detectorsDim, long slices) { long i, j, k; int full_window; + full_window = (int)(2*slices_window_halfsize+1)*(2*vert_window_halfsize+1)*(2*horiz_window_halfsize+1); if (slices == 1) { /****************2D INPUT ***************/ - full_window = (2*horiz_window_halfsize+1); #pragma omp parallel for shared(residual, weights) private(j, i) for(j=0; j= 0) && (j1 < detectorsDim)) { - Values_Vec[counter] = residual[i*detectorsDim+j1]; + for(l=-vert_window_halfsize; l < vert_window_halfsize; l++) { + i1 = i + l; + if ((j1 >= 0) && (j1 < detectorsDim) && (i1 >= 0) && (i1 < anglesDim)) { + Values_Vec[counter] = residual[i1*detectorsDim+j1]; } else Values_Vec[counter] = residual[index]; counter ++; - } + }} + /* perform sorting of the vector array */ for (x = 0; x < counter-1; x++) { for (y = 0; y < counter-x-1; y++) { @@ -102,59 +107,39 @@ float RingWeights2D(float *residual, float *weights, int horiz_window_halfsize, } } } - /* include diagonal values in the estimation of the backround */ - int i1; - float *Values_Vec_diag1; - Values_Vec_diag1 = (float*) calloc (full_window, sizeof(float)); - counter = 0; - for(k=-horiz_window_halfsize; k < horiz_window_halfsize; k++) { - j1 = j + k; - i1 = i + k; - if ((j1 >= 0) && (j1 < detectorsDim) && (i1 >= 0) && (i1 < anglesDim)) { - Values_Vec_diag1[counter] = residual[i1*detectorsDim+j1]; - } - else Values_Vec_diag1[counter] = residual[index]; - counter ++; - } - /* perform sorting of the vector array */ - for (x = 0; x < counter-1; x++) { - for (y = 0; y < counter-x-1; y++) { - if (Values_Vec_diag1[y] > Values_Vec_diag1[y+1]) { - swap(&Values_Vec_diag1[y], &Values_Vec_diag1[y+1]); - } - } - } - weights[index] = residual[index] - 0.5f*(Values_Vec[horiz_window_halfsize] + Values_Vec_diag1[horiz_window_halfsize]); + weights[index] = residual[index] - Values_Vec[midval]; free(Values_Vec); - free(Values_Vec_diag1); return *weights; } -float RingWeights3D(float *residual, float *weights, int horiz_window_halfsize, int vert_window_halfsize, int full_window, long anglesDim, long detectorsDim, long slices, long i, long j, long k) +float RingWeights3D(float *residual, float *weights, int horiz_window_halfsize, int vert_window_halfsize, int slices_window_halfsize, int full_window, long anglesDim, long detectorsDim, long slices, long i, long j, long k) { float *Values_Vec; - long k1, j1, l, v, index; + long k1, j1, i1, l, v, m, index; int counter, x, y, midval; midval = (int)(0.5f*full_window) - 1; index = (detectorsDim*anglesDim*k) + i*detectorsDim+j; Values_Vec = (float*) calloc (full_window, sizeof(float)); - - /* intiate the estimation of the backround using strictly horizontal values */ +/* counter = 0; for(l=-horiz_window_halfsize; l < horiz_window_halfsize; l++) { - for(v=-vert_window_halfsize; v < vert_window_halfsize; v++) { - j1 = j + l; - k1 = k + v; - if ((j1 >= 0) && (j1 < detectorsDim) && (k1 >= 0) && (k1 < slices)) { - Values_Vec[counter] = residual[(detectorsDim*anglesDim*k1) + i*detectorsDim+j1]; + j1 = j + l; + for(m=-vert_window_halfsize; m < vert_window_halfsize; m++) { + i1 = i + m; + for(v=-slices_window_halfsize; v < slices_window_halfsize; v++) { + k1 = k + v; + if ((j1 >= 0) && (j1 < detectorsDim) && (k1 >= 0) && (k1 < slices) && (i1 >= 0) && (i1 < anglesDim)) { + Values_Vec[counter] = residual[(detectorsDim*anglesDim*k1) + i1*detectorsDim+j1]; } else Values_Vec[counter] = residual[index]; counter ++; - }} + }}} +*/ /* perform sorting of the vector array */ + /* for (x = 0; x < counter-1; x++) { for (y = 0; y < counter-x-1; y++) { if (Values_Vec[y] > Values_Vec[y+1]) { @@ -163,6 +148,9 @@ float RingWeights3D(float *residual, float *weights, int horiz_window_halfsize, } } weights[index] = residual[index] - Values_Vec[midval]; + */ + j1 = j + 15; + if ((j1 > 0) && (j1 < detectorsDim)) weights[index] = residual[(detectorsDim*anglesDim*k) + i*detectorsDim+j1]; free(Values_Vec); return *weights; diff --git a/src/Core/functions_CPU/RingWeights_core.h b/src/Core/functions_CPU/RingWeights_core.h index 0e85916db..706e38b3a 100644 --- a/src/Core/functions_CPU/RingWeights_core.h +++ b/src/Core/functions_CPU/RingWeights_core.h @@ -39,11 +39,11 @@ #ifdef __cplusplus extern "C" { #endif -float RingWeights_main(float *residual, float *weights, int horiz_window_halfsize, int vert_window_halfsize, long anglesDim, long detectorsDim, long slices); +float RingWeights_main(float *residual, float *weights, int horiz_window_halfsize, int vert_window_halfsize, int slices_window_halfsize, long anglesDim, long detectorsDim, long slices); /************2D functions ***********/ -float RingWeights2D(float *residual, float *weights, int horiz_window_halfsize, int full_window, long anglesDim, long detectorsDim, long i, long j); +float RingWeights2D(float *residual, float *weights, int horiz_window_halfsize, int vert_window_halfsize, int full_window, long anglesDim, long detectorsDim, long i, long j); /************3D functions ***********/ -float RingWeights3D(float *residual, float *weights, int horiz_window_halfsize, int vert_window_halfsize, int full_window, long anglesDim, long detectorsDim, long slices, long i, long j, long k); +float RingWeights3D(float *residual, float *weights, int horiz_window_halfsize, int vert_window_halfsize, int slices_window_halfsize, int full_window, long anglesDim, long detectorsDim, long slices, long i, long j, long k); #ifdef __cplusplus } #endif diff --git a/src/Python/src/cpu_wrappers.pyx b/src/Python/src/cpu_wrappers.pyx index 1fd590809..804062a9b 100644 --- a/src/Python/src/cpu_wrappers.pyx +++ b/src/Python/src/cpu_wrappers.pyx @@ -15,17 +15,18 @@ import cython import numpy as np cimport numpy as np -cdef extern float RingWeights_main(float *residual, float *weights, int horiz_window_halfsize, int vert_window_halfsize, long anglesDim, long detectorsDim, long slices); +cdef extern float RingWeights_main(float *residual, float *weights, int horiz_window_halfsize, int vert_window_halfsize, int slices_window_halfsize, long anglesDim, long detectorsDim, long slices); ############################################################################## -def RING_WEIGHTS(residual, horiz_window_halfsize, vert_window_halfsize): +def RING_WEIGHTS(residual, horiz_window_halfsize, vert_window_halfsize, slices_window_halfsize): if residual.ndim == 2: - return RING_WEIGHTS_2D(residual, horiz_window_halfsize) + return RING_WEIGHTS_2D(residual, horiz_window_halfsize, vert_window_halfsize) elif residual.ndim == 3: - return RING_WEIGHTS_3D(residual, horiz_window_halfsize, vert_window_halfsize) + return RING_WEIGHTS_3D(residual, horiz_window_halfsize, vert_window_halfsize, slices_window_halfsize) def RING_WEIGHTS_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] residual, - int horiz_window_halfsize): + int horiz_window_halfsize, + int vert_window_halfsize): cdef long dims[2] dims[0] = residual.shape[0] @@ -34,13 +35,13 @@ def RING_WEIGHTS_2D(np.ndarray[np.float32_t, ndim=2, mode="c"] residual, cdef np.ndarray[np.float32_t, ndim=2, mode="c"] weights = \ np.zeros([dims[0],dims[1]], dtype='float32') - RingWeights_main(&residual[0,0], &weights[0,0], horiz_window_halfsize, 0, dims[1], dims[0], 1); - + RingWeights_main(&residual[0,0], &weights[0,0], horiz_window_halfsize, vert_window_halfsize, 0, dims[1], dims[0], 1); return weights def RING_WEIGHTS_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] residual, int horiz_window_halfsize, - int vert_window_halfsize): + int vert_window_halfsize, + int slices_window_halfsize): cdef long dims[3] dims[0] = residual.shape[0] @@ -50,6 +51,5 @@ def RING_WEIGHTS_3D(np.ndarray[np.float32_t, ndim=3, mode="c"] residual, cdef np.ndarray[np.float32_t, ndim=3, mode="c"] weights = \ np.zeros([dims[0],dims[1],dims[2]], dtype='float32') - RingWeights_main(&residual[0,0,0], &weights[0,0,0], horiz_window_halfsize, vert_window_halfsize, dims[2], dims[1], dims[0]); - + RingWeights_main(&residual[0,0,0], &weights[0,0,0], horiz_window_halfsize, vert_window_halfsize, slices_window_halfsize, dims[1], dims[2], dims[0]); return weights diff --git a/src/Python/tomobar/methodsDIR.py b/src/Python/tomobar/methodsDIR.py index 0703db253..1732ae501 100644 --- a/src/Python/tomobar/methodsDIR.py +++ b/src/Python/tomobar/methodsDIR.py @@ -77,7 +77,7 @@ def FORWPROJ(self, image): sinogram = Atools.forwproj(image) if (self.geom == '3D'): from tomobar.supp.astraOP import AstraTools3D - Atools = AstraTools3D(self.DetectorsDimH, self.DetectorsDimV, self.AnglesVec, self.ObjSize) # initiate 3D ASTRA class object + Atools = AstraTools3D(self.DetectorsDimH, self.DetectorsDimV, self.AnglesVec, self.CenterRotOffset, self.ObjSize) # initiate 3D ASTRA class object sinogram = Atools.forwproj(image) return sinogram def BACKPROJ(self, sinogram): diff --git a/src/Python/tomobar/methodsIR.py b/src/Python/tomobar/methodsIR.py index 25aa7fd38..a6a93d33a 100644 --- a/src/Python/tomobar/methodsIR.py +++ b/src/Python/tomobar/methodsIR.py @@ -172,7 +172,8 @@ def FISTA(self, lambdaR_L1 = 0.0, # regularisation parameter for GH data model alpha_ring = 50, # GH data model convergence accelerator (use carefully -> can generate artifacts) ring_model_horiz_size = None, # enable a better model to supress ring artifacts, size of the window defines a possible thickness of ring artifacts - ring_model_vert_size = 1, # for 3D case only define a vertical window + ring_model_vert_size = 0, # define a vertical window + ring_model_slices_size = 0, # 3D case only, define a slices window huber_data_threshold = 0.0, # threshold parameter for __Huber__ data fidelity student_data_threshold = 0.0, # threshold parameter for __Students't__ data fidelity initialise = None, # initialise reconstruction with an array @@ -321,7 +322,7 @@ def FISTA(self, multHuber[(np.where(np.abs(res) > huber_data_threshold))] = np.divide(huber_data_threshold, np.abs(res[(np.where(np.abs(res) > huber_data_threshold))])) else: # Apply smarter Huber model to supress ring artifacts - offset_rings = RING_WEIGHTS(res, ring_model_horiz_size, ring_model_vert_size) + offset_rings = RING_WEIGHTS(res, ring_model_horiz_size, ring_model_vert_size, ring_model_slices_size) tempRes = res+offset_rings multHuber[(np.where(np.abs(tempRes) > huber_data_threshold))] = np.divide(huber_data_threshold, np.abs(tempRes[(np.where(np.abs(tempRes) > huber_data_threshold))])) if (self.OS_number > 1): @@ -337,7 +338,7 @@ def FISTA(self, multStudent = np.divide(2.0, student_data_threshold**2 + res**2) else: # Apply smarter Students't model to supress ring artifacts - offset_rings = RING_WEIGHTS(res, ring_model_horiz_size, ring_model_vert_size) + offset_rings = RING_WEIGHTS(res, ring_model_horiz_size, ring_model_vert_size, ring_model_slices_size) tempRes = res+offset_rings multStudent = np.divide(2.0, student_data_threshold**2 + tempRes**2) if (self.OS_number > 1):