diff options
Diffstat (limited to 'demos/DemoRD2.m')
| -rw-r--r-- | demos/DemoRD2.m | 110 |
1 files changed, 88 insertions, 22 deletions
diff --git a/demos/DemoRD2.m b/demos/DemoRD2.m index 0db3595..717a55d 100644 --- a/demos/DemoRD2.m +++ b/demos/DemoRD2.m @@ -1,12 +1,12 @@ % Demonstration of tomographic 3D reconstruction from X-ray synchrotron % dataset (dendrites) using various data fidelities -% ! can take up to 15-20 minutes to run for the whole 3D data ! +% ! It is advisable not to run the whole script, it will take lots of time to reconstruct the whole 3D data using many algorithms ! clear all close all %% % % adding paths addpath('../data/'); -addpath('../main_func/'); addpath('../main_func/regularizers_CPU/'); +addpath('../main_func/'); addpath('../main_func/regularizers_CPU/'); addpath('../main_func/regularizers_GPU/NL_Regul/'); addpath('../main_func/regularizers_GPU/Diffus_HO/'); addpath('../supp/'); load('DendrRawData.mat') % load raw data of 3D dendritic set @@ -30,7 +30,7 @@ Weights3D = single(data_raw3D); % weights for PW model clear data_raw3D %% % set projection/reconstruction geometry here -Z_slices = 3; +Z_slices = 5; det_row_count = Z_slices; proj_geom = astra_create_proj_geom('parallel3d', 1, 1, det_row_count, size_det, angles_rad); vol_geom = astra_create_vol_geom(recon_size,recon_size,Z_slices); @@ -38,70 +38,136 @@ vol_geom = astra_create_vol_geom(recon_size,recon_size,Z_slices); fprintf('%s\n', 'Reconstruction using FBP...'); FBP = iradon(Sino3D(:,:,10), angles,recon_size); figure; imshow(FBP , [0, 3]); title ('FBP reconstruction'); + +%--------FISTA_REC modular reconstruction alogrithms--------- %% -fprintf('%s\n', 'Reconstruction using FISTA-PWLS without regularization...'); +fprintf('%s\n', 'Reconstruction using FISTA-OS-PWLS without regularization...'); clear params params.proj_geom = proj_geom; % pass geometry to the function params.vol_geom = vol_geom; params.sino = Sino3D; -params.iterFISTA = 35; +params.iterFISTA = 12; params.weights = Weights3D; +params.subsets = 16; % the number of ordered subsets params.show = 1; params.maxvalplot = 2.5; params.slice = 2; -tic; [X_fista, output] = FISTA_REC(params); toc; -figure; imshow(X_fista(:,:,params.slice) , [0, 2.5]); title ('FISTA-PWLS reconstruction'); +tic; [X_fista, outputFISTA] = FISTA_REC(params); toc; +figure; imshow(X_fista(:,:,params.slice) , [0, 2.5]); title ('FISTA-OS-PWLS reconstruction'); %% -fprintf('%s\n', 'Reconstruction using FISTA-PWLS-TV...'); +fprintf('%s\n', 'Reconstruction using FISTA-OS-PWLS-TV...'); clear params params.proj_geom = proj_geom; % pass geometry to the function params.vol_geom = vol_geom; params.sino = Sino3D; -params.iterFISTA = 40; +params.iterFISTA = 12; params.Regul_Lambda_FGPTV = 0.005; % TV regularization parameter for FGP-TV params.weights = Weights3D; +params.subsets = 16; % the number of ordered subsets params.show = 1; params.maxvalplot = 2.5; params.slice = 2; -tic; [X_fista_TV] = FISTA_REC(params); toc; -figure; imshow(X_fista_TV(:,:,params.slice) , [0, 2.5]); title ('FISTA-PWLS-TV reconstruction'); -%% +tic; [X_fista_TV, outputTV] = FISTA_REC(params); toc; +figure; imshow(X_fista_TV(:,:,params.slice) , [0, 2.5]); title ('FISTA-OS-PWLS-TV reconstruction'); %% -fprintf('%s\n', 'Reconstruction using FISTA-GH-TV...'); +fprintf('%s\n', 'Reconstruction using FISTA-OS-GH-TV...'); clear params params.proj_geom = proj_geom; % pass geometry to the function params.vol_geom = vol_geom; params.sino = Sino3D; -params.iterFISTA = 40; +params.iterFISTA = 12; params.Regul_Lambda_FGPTV = 0.005; % TV regularization parameter for FGP-TV params.Ring_LambdaR_L1 = 0.002; % Soft-Thresh L1 ring variable parameter params.Ring_Alpha = 21; % to boost ring removal procedure params.weights = Weights3D; +params.subsets = 16; % the number of ordered subsets params.show = 1; params.maxvalplot = 2.5; params.slice = 2; -tic; [X_fista_GH_TV] = FISTA_REC(params); toc; -figure; imshow(X_fista_GH_TV(:,:,params.slice) , [0, 2.5]); title ('FISTA-GH-TV reconstruction'); +tic; [X_fista_GH_TV, outputGHTV] = FISTA_REC(params); toc; +figure; imshow(X_fista_GH_TV(:,:,params.slice) , [0, 2.5]); title ('FISTA-OS-GH-TV reconstruction'); %% -%% -fprintf('%s\n', 'Reconstruction using FISTA-GH-TV-LLT...'); +fprintf('%s\n', 'Reconstruction using FISTA-OS-GH-TV-LLT...'); clear params params.proj_geom = proj_geom; % pass geometry to the function params.vol_geom = vol_geom; params.sino = Sino3D; -params.iterFISTA = 5; +params.iterFISTA = 12; params.Regul_Lambda_FGPTV = 0.005; % TV regularization parameter for FGP-TV -params.Regul_LambdaHO = 250; % regularization parameter for LLT problem +params.Regul_LambdaLLT = 100; % regularization parameter for LLT problem params.Regul_tauLLT = 0.0005; % time-step parameter for the explicit scheme params.Ring_LambdaR_L1 = 0.002; % Soft-Thresh L1 ring variable parameter params.Ring_Alpha = 21; % to boost ring removal procedure params.weights = Weights3D; +params.subsets = 16; % the number of ordered subsets params.show = 1; params.maxvalplot = 2.5; params.slice = 2; -tic; [X_fista_GH_TVLLT] = FISTA_REC(params); toc; -figure; imshow(X_fista_GH_TVLLT(:,:,params.slice) , [0, 2.5]); title ('FISTA-GH-TV-LLT reconstruction'); +tic; [X_fista_GH_TVLLT, outputGH_TVLLT] = FISTA_REC(params); toc; +figure; imshow(X_fista_GH_TVLLT(:,:,params.slice) , [0, 2.5]); title ('FISTA-OS-GH-TV-LLT reconstruction'); + +%% +fprintf('%s\n', 'Reconstruction using FISTA-OS-GH-HigherOrderDiffusion...'); +% !GPU version! +clear params +params.proj_geom = proj_geom; % pass geometry to the function +params.vol_geom = vol_geom; +params.sino = Sino3D(:,:,1:5); +params.iterFISTA = 25; +params.Regul_LambdaDiffHO = 2; % DiffHO regularization parameter +params.Regul_DiffHO_EdgePar = 0.05; % threshold parameter +params.Regul_Iterations = 150; +params.Ring_LambdaR_L1 = 0.002; % Soft-Thresh L1 ring variable parameter +params.Ring_Alpha = 21; % to boost ring removal procedure +params.weights = Weights3D(:,:,1:5); +params.subsets = 16; % the number of ordered subsets +params.show = 1; +params.maxvalplot = 2.5; params.slice = 1; + +tic; [X_fista_GH_HO, outputHO] = FISTA_REC(params); toc; +figure; imshow(X_fista_GH_HO(:,:,params.slice) , [0, 2.5]); title ('FISTA-OS-HigherOrderDiffusion reconstruction'); + +%% +fprintf('%s\n', 'Reconstruction using FISTA-PB...'); +% !GPU version! +clear params +params.proj_geom = proj_geom; % pass geometry to the function +params.vol_geom = vol_geom; +params.sino = Sino3D(:,:,1); +params.iterFISTA = 25; +params.Regul_LambdaPatchBased_GPU = 3; % PB regularization parameter +params.Regul_PB_h = 0.04; % threhsold parameter +params.Regul_PB_SearchW = 3; +params.Regul_PB_SimilW = 1; +params.Ring_LambdaR_L1 = 0.002; % Soft-Thresh L1 ring variable parameter +params.Ring_Alpha = 21; % to boost ring removal procedure +params.weights = Weights3D(:,:,1); +params.show = 1; +params.maxvalplot = 2.5; params.slice = 1; + +tic; [X_fista_GH_PB, outputPB] = FISTA_REC(params); toc; +figure; imshow(X_fista_GH_PB(:,:,params.slice) , [0, 2.5]); title ('FISTA-OS-PB reconstruction'); %% +fprintf('%s\n', 'Reconstruction using FISTA-OS-GH-TGV...'); +% still testing... +clear params +params.proj_geom = proj_geom; % pass geometry to the function +params.vol_geom = vol_geom; +params.sino = Sino3D; +params.iterFISTA = 12; +params.Regul_LambdaTGV = 0.5; % TGV regularization parameter +params.Regul_Iterations = 5; +params.Ring_LambdaR_L1 = 0.002; % Soft-Thresh L1 ring variable parameter +params.Ring_Alpha = 21; % to boost ring removal procedure +params.weights = Weights3D; +params.subsets = 16; % the number of ordered subsets +params.show = 1; +params.maxvalplot = 2.5; params.slice = 1; + +tic; [X_fista_GH_TGV, outputTGV] = FISTA_REC(params); toc; +figure; imshow(X_fista_GH_TGV(:,:,params.slice) , [0, 2.5]); title ('FISTA-OS-GH-TGV reconstruction'); + + %% % fprintf('%s\n', 'Reconstruction using FISTA-Student-TV...'); % clear params |