From cbeb8a1174498751e38a3de8cd6fe55abae20192 Mon Sep 17 00:00:00 2001 From: Daniil Kazantsev Date: Thu, 3 Aug 2017 00:10:53 +0100 Subject: some cleaning inside C functions, updated mex_compile file and the work on ordered-subset FISTA started --- demos/DemoRD2.m | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) (limited to 'demos') diff --git a/demos/DemoRD2.m b/demos/DemoRD2.m index 0db3595..293c1cd 100644 --- a/demos/DemoRD2.m +++ b/demos/DemoRD2.m @@ -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 = 20; 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); @@ -44,7 +44,7 @@ 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 = 1; params.weights = Weights3D; params.show = 1; params.maxvalplot = 2.5; params.slice = 2; -- cgit v1.2.3 From d7281b95f70dd3707f82640972a52f4155c2fde5 Mon Sep 17 00:00:00 2001 From: Daniil Kazantsev Date: Thu, 3 Aug 2017 15:25:50 +0100 Subject: Ordered Subset FISTA added and demos updated in DemoRD2 --- demos/DemoRD2.m | 86 +++++++++++++++++++++++++++++++++++++++++++-------------- 1 file changed, 65 insertions(+), 21 deletions(-) (limited to 'demos') diff --git a/demos/DemoRD2.m b/demos/DemoRD2.m index 293c1cd..f177e26 100644 --- a/demos/DemoRD2.m +++ b/demos/DemoRD2.m @@ -1,6 +1,6 @@ % 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 %% @@ -30,7 +30,7 @@ Weights3D = single(data_raw3D); % weights for PW model clear data_raw3D %% % set projection/reconstruction geometry here -Z_slices = 20; +Z_slices = 1; 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,114 @@ 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 = 1; +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-PB...'); +% % very-slow on CPU +% clear params +% params.proj_geom = proj_geom; % pass geometry to the function +% params.vol_geom = vol_geom; +% params.sino = Sino3D(:,:,1); +% params.iterFISTA = 12; +% params.Regul_LambdaPatchBased = 1; % PB regularization parameter +% params.Regul_PB_h = 0.1; % threhsold parameter +% 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.subsets = 16; % the number of ordered subsets +% 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 -- cgit v1.2.3 From 62ab6cd46c3f1c189328c8d41899db7444c7ac29 Mon Sep 17 00:00:00 2001 From: Daniil Kazantsev Date: Mon, 11 Sep 2017 09:36:13 +0100 Subject: 2 new GPU regularizers, FGP objective fixed, FISTA_REC updated --- demos/DemoRD2.m | 62 ++++++++++++++++++++++++++++++++++++++------------------- 1 file changed, 42 insertions(+), 20 deletions(-) (limited to 'demos') diff --git a/demos/DemoRD2.m b/demos/DemoRD2.m index f177e26..717a55d 100644 --- a/demos/DemoRD2.m +++ b/demos/DemoRD2.m @@ -6,7 +6,7 @@ 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 = 1; +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); @@ -107,25 +107,46 @@ 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-PB...'); -% % very-slow on CPU -% clear params -% params.proj_geom = proj_geom; % pass geometry to the function -% params.vol_geom = vol_geom; -% params.sino = Sino3D(:,:,1); -% params.iterFISTA = 12; -% params.Regul_LambdaPatchBased = 1; % PB regularization parameter -% params.Regul_PB_h = 0.1; % threhsold parameter -% 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.subsets = 16; % the number of ordered subsets -% 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-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... @@ -146,6 +167,7 @@ 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 -- cgit v1.2.3 From a203949c84484fe2641e39451f033d20d445b1f3 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Wed, 23 Aug 2017 16:51:18 +0100 Subject: export/import data from hdf5 Added file to export the data from DemoRD2.m to HDF5 to pass it to Python. Added file to import the data from DemoRD2.m from HDF5. --- demos/exportDemoRD2Data.m | 35 +++++++++++++++++++++++++++++++++++ 1 file changed, 35 insertions(+) create mode 100644 demos/exportDemoRD2Data.m (limited to 'demos') diff --git a/demos/exportDemoRD2Data.m b/demos/exportDemoRD2Data.m new file mode 100644 index 0000000..028353b --- /dev/null +++ b/demos/exportDemoRD2Data.m @@ -0,0 +1,35 @@ +clear all +close all +%% +% % adding paths +addpath('../data/'); +addpath('../main_func/'); addpath('../main_func/regularizers_CPU/'); +addpath('../supp/'); + +load('DendrRawData.mat') % load raw data of 3D dendritic set +angles_rad = angles*(pi/180); % conversion to radians +size_det = size(data_raw3D,1); % detectors dim +angSize = size(data_raw3D, 2); % angles dim +slices_tot = size(data_raw3D, 3); % no of slices +recon_size = 950; % reconstruction size + +Sino3D = zeros(size_det, angSize, slices_tot, 'single'); % log-corrected sino +% normalizing the data +for jj = 1:slices_tot + sino = data_raw3D(:,:,jj); + for ii = 1:angSize + Sino3D(:,ii,jj) = log((flats_ar(:,jj)-darks_ar(:,jj))./(single(sino(:,ii)) - darks_ar(:,jj))); + end +end + +Sino3D = Sino3D.*1000; +Weights3D = single(data_raw3D); % weights for PW model +clear data_raw3D + +hdf5write('DendrData.h5', '/Weights3D', Weights3D) +hdf5write('DendrData.h5', '/Sino3D', Sino3D, 'WriteMode', 'append') +hdf5write('DendrData.h5', '/angles_rad', angles_rad, 'WriteMode', 'append') +hdf5write('DendrData.h5', '/size_det', size_det, 'WriteMode', 'append') +hdf5write('DendrData.h5', '/angSize', angSize, 'WriteMode', 'append') +hdf5write('DendrData.h5', '/slices_tot', slices_tot, 'WriteMode', 'append') +hdf5write('DendrData.h5', '/recon_size', recon_size, 'WriteMode', 'append') \ No newline at end of file -- cgit v1.2.3 From a32469e9d96c6039b0a5e8e1b22f577e87484c40 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Wed, 11 Oct 2017 15:11:56 +0100 Subject: comment out clear of data_raw3D --- demos/exportDemoRD2Data.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'demos') diff --git a/demos/exportDemoRD2Data.m b/demos/exportDemoRD2Data.m index 028353b..0ca036b 100644 --- a/demos/exportDemoRD2Data.m +++ b/demos/exportDemoRD2Data.m @@ -24,7 +24,7 @@ end Sino3D = Sino3D.*1000; Weights3D = single(data_raw3D); % weights for PW model -clear data_raw3D +%clear data_raw3D hdf5write('DendrData.h5', '/Weights3D', Weights3D) hdf5write('DendrData.h5', '/Sino3D', Sino3D, 'WriteMode', 'append') -- cgit v1.2.3