From 952322eb1394f18cedfbb51488e76ed52728d98d Mon Sep 17 00:00:00 2001 From: Daniil Kazantsev Date: Thu, 19 Oct 2017 22:16:43 +0100 Subject: GH-OS fixed, demos updated, add_artifacts function --- demos/DemoRD1.m | 100 ------------------- demos/DemoRD2.m | 188 ------------------------------------ demos/Demo_Phantom3D_Cone.m | 5 +- demos/Demo_Phantom3D_Parallel.m | 67 ++++++++++--- demos/Demo_RealData3D_Parallel.m | 186 +++++++++++++++++++++++++++++++++++ main_func/FISTA_REC.m | 4 + main_func/regularizers_CPU/FGP_TV.c | 6 +- supp/sino_add_artifacts.m | 33 +++++++ 8 files changed, 284 insertions(+), 305 deletions(-) delete mode 100644 demos/DemoRD1.m delete mode 100644 demos/DemoRD2.m create mode 100644 demos/Demo_RealData3D_Parallel.m create mode 100644 supp/sino_add_artifacts.m diff --git a/demos/DemoRD1.m b/demos/DemoRD1.m deleted file mode 100644 index 5bb5f6b..0000000 --- a/demos/DemoRD1.m +++ /dev/null @@ -1,100 +0,0 @@ -% Demonstration of tomographic reconstruction from neutron tomography -% dataset (basalt sample) using Student t data fidelity -clear all -close all - -% adding paths -addpath('../data/'); -addpath('../main_func/'); addpath('../main_func/regularizers_CPU/'); -addpath('../supp/'); - -load('sino_basalt.mat') % load real neutron data - -size_det = size(sino_basalt, 1); % detector size -angSize = size(sino_basalt,2); % angles dim -recon_size = 650; % reconstruction size - -FBP = iradon(sino_basalt, rad2deg(angles),recon_size); -figure; imshow(FBP , [0, 0.45]); title ('FBP reconstruction'); -%% -% set projection/reconstruction geometry here -Z_slices = 1; -det_row_count = Z_slices; -proj_geom = astra_create_proj_geom('parallel3d', 1, 1, det_row_count, size_det, angles); -vol_geom = astra_create_vol_geom(recon_size,recon_size,Z_slices); -%% -fprintf('%s\n', 'Reconstruction using FISTA-LS without regularization...'); -clear params -params.proj_geom = proj_geom; % pass geometry to the function -params.vol_geom = vol_geom; -params.sino = sino_basalt; -params.iterFISTA = 50; -params.show = 0; -params.maxvalplot = 0.6; params.slice = 1; - -tic; [X_fista] = FISTA_REC(params); toc; -figure; imshow(X_fista , [0, 0.45]); title ('FISTA-LS reconstruction'); -%% -fprintf('%s\n', 'Reconstruction using FISTA-LS-TV...'); -clear params -params.proj_geom = proj_geom; % pass geometry to the function -params.vol_geom = vol_geom; -params.sino = sino_basalt; -params.iterFISTA = 60; -params.Regul_LambdaTV = 0.0003; % TV regularization parameter -params.show = 0; -params.maxvalplot = 0.6; params.slice = 1; - -tic; [X_fista_TV] = FISTA_REC(params); toc; -figure; imshow(X_fista_TV , [0, 0.45]); title ('FISTA-LS-TV reconstruction'); -%% -%% -fprintf('%s\n', 'Reconstruction using FISTA-GH-TV...'); -clear params -params.proj_geom = proj_geom; % pass geometry to the function -params.vol_geom = vol_geom; -params.sino = sino_basalt; -params.iterFISTA = 60; -params.Regul_LambdaTV = 0.0003; % TV regularization parameter -params.Ring_LambdaR_L1 = 0.001; % Soft-Thresh L1 ring variable parameter -params.Ring_Alpha = 20; % acceleration for ring variable -params.show = 0; -params.maxvalplot = 0.6; params.slice = 1; - -tic; [X_fista_GH_TV] = FISTA_REC(params); toc; -figure; imshow(X_fista_GH_TV , [0, 0.45]); title ('FISTA-GH-TV reconstruction'); -%% -%% -fprintf('%s\n', 'Reconstruction using FISTA-Student-TV...'); -clear params -params.proj_geom = proj_geom; % pass geometry to the function -params.vol_geom = vol_geom; -params.sino = sino_basalt; -params.iterFISTA = 50; -params.L_const = 3500; % Lipshitz constant -params.Regul_LambdaTV = 0.0003; % TV regularization parameter -params.fidelity = 'student'; % choosing Student t penalty -params.show = 1; -params.initialize = 1; % warm start with SIRT -params.maxvalplot = 0.6; params.slice = 1; - -tic; [X_fistaStudentTV] = FISTA_REC(params); toc; -figure; imshow(X_fistaStudentTV , [0, 0.45]); title ('FISTA-Student-TV reconstruction'); -%% - -fprintf('%s\n', 'Segmentation using OTSU method ...'); -level = graythresh(X_fista); -Segm_FISTA = im2bw(X_fista,level); -figure; imshow(Segm_FISTA, []); title ('Segmented FISTA-LS reconstruction'); - -level = graythresh(X_fista_TV); -Segm_FISTA_TV = im2bw(X_fista_TV,level); -figure; imshow(Segm_FISTA_TV, []); title ('Segmented FISTA-LS-TV reconstruction'); - -level = graythresh(X_fista_GH_TV); -BW_FISTA_GH_TV = im2bw(X_fista_GH_TV,level); -figure; imshow(BW_FISTA_GH_TV, []); title ('Segmented FISTA-GH-TV reconstruction'); - -level = graythresh(X_fistaStudentTV); -BW_FISTA_Student_TV = im2bw(X_fistaStudentTV,level); -figure; imshow(BW_FISTA_Student_TV, []); title ('Segmented FISTA-Student-LS reconstruction'); \ No newline at end of file diff --git a/demos/DemoRD2.m b/demos/DemoRD2.m deleted file mode 100644 index 717a55d..0000000 --- a/demos/DemoRD2.m +++ /dev/null @@ -1,188 +0,0 @@ -% Demonstration of tomographic 3D reconstruction from X-ray synchrotron -% dataset (dendrites) using various data fidelities -% ! 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/regularizers_GPU/NL_Regul/'); addpath('../main_func/regularizers_GPU/Diffus_HO/'); -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 -%% -% set projection/reconstruction geometry here -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); -%% -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-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 = 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, 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-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 = 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, 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-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 = 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, 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-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 = 12; -params.Regul_Lambda_FGPTV = 0.005; % TV regularization parameter for FGP-TV -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, 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 -% params.proj_geom = proj_geom; % pass geometry to the function -% params.vol_geom = vol_geom; -% params.sino = Sino3D(:,:,10); -% params.iterFISTA = 50; -% params.L_const = 0.01; % Lipshitz constant -% params.Regul_LambdaTV = 0.008; % TV regularization parameter for FISTA-TV -% params.fidelity = 'student'; % choosing Student t penalty -% params.weights = Weights3D(:,:,10); -% params.show = 0; -% params.initialize = 1; -% params.maxvalplot = 2.5; params.slice = 1; -% -% tic; [X_fistaStudentTV] = FISTA_REC(params); toc; -% figure; imshow(X_fistaStudentTV(:,:,1), [0, 2.5]); title ('FISTA-Student-TV reconstruction'); -%% diff --git a/demos/Demo_Phantom3D_Cone.m b/demos/Demo_Phantom3D_Cone.m index 3a9178b..a8f2c92 100644 --- a/demos/Demo_Phantom3D_Cone.m +++ b/demos/Demo_Phantom3D_Cone.m @@ -8,7 +8,6 @@ addpath('../data/'); addpath('../main_func/'); addpath('../main_func/regularizers_CPU/'); addpath('../main_func/regularizers_GPU/NL_Regul/'); addpath('../main_func/regularizers_GPU/Diffus_HO/'); addpath('../supp/'); - %% % build 3D phantom using TomoPhantom modelNo = 3; % see Phantom3DLibrary.dat file in TomoPhantom @@ -16,9 +15,11 @@ N = 256; % x-y-z size (cubic image) angles = 0:1.5:360; % angles vector in degrees angles_rad = angles*(pi/180); % conversion to radians det_size = round(sqrt(2)*N); % detector size -% in order to run functions you have to go to the directory: + +%---------TomoPhantom routines---------% pathTP = '/home/algol/Documents/MATLAB/TomoPhantom/functions/models/Phantom3DLibrary.dat'; % path to TomoPhantom parameters file TomoPhantom = buildPhantom3D(modelNo,N,pathTP); % generate 3D phantom +%--------------------------------------% %% % using ASTRA-toolbox to set the projection geometry (cone beam) % eg: astra.create_proj_geom('cone', 1.0 (resol), 1.0 (resol), detectorRowCount, detectorColCount, angles, originToSource, originToDetector) diff --git a/demos/Demo_Phantom3D_Parallel.m b/demos/Demo_Phantom3D_Parallel.m index 6a54450..9ef1f3c 100644 --- a/demos/Demo_Phantom3D_Parallel.m +++ b/demos/Demo_Phantom3D_Parallel.m @@ -9,18 +9,27 @@ addpath('../main_func/'); addpath('../main_func/regularizers_CPU/'); addpath('.. addpath('../supp/'); %% -% build 3D phantom using TomoPhantom and generate projection data +% Main reconstruction/data generation parameters modelNo = 2; % see Phantom3DLibrary.dat file in TomoPhantom N = 256; % x-y-z size (cubic image) angles = 1:0.5:180; % angles vector in degrees angles_rad = angles*(pi/180); % conversion to radians det_size = round(sqrt(2)*N); % detector size -% in order to run functions you have to go to the directory: + +%---------TomoPhantom routines---------% pathTP = '/home/algol/Documents/MATLAB/TomoPhantom/functions/models/Phantom3DLibrary.dat'; % path to TomoPhantom parameters file TomoPhantom = buildPhantom3D(modelNo,N,pathTP); % generate 3D phantom sino_tomophan3D = buildSino3D(modelNo, N, det_size, single(angles),pathTP); % generate ideal data +%--------------------------------------% % Adding noise and distortions if required -sino_artifacts = sino_add_artifacts(sino_tomophan3D,'rings'); +sino_tomophan3D = sino_add_artifacts(sino_tomophan3D,'rings'); +% adding Poisson noise +dose = 3e9; % photon flux (controls noise level) +multifactor = max(sino_tomophan3D(:)); +dataExp = dose.*exp(-sino_tomophan3D/multifactor); % noiseless raw data +dataRaw = astra_add_noise_to_sino(dataExp, dose); % pre-log noisy raw data (weights) +sino3D_log = log(dose./max(dataRaw,1))*multifactor; %log corrected data -> sinogram +clear dataExp sino_tomophan3D % %% % using ASTRA-toolbox to set the projection geometry (parallel beam) @@ -28,29 +37,34 @@ proj_geom = astra_create_proj_geom('parallel', 1, det_size, angles_rad); vol_geom = astra_create_vol_geom(N,N); %% fprintf('%s\n', 'Reconstructing with FBP using ASTRA-toolbox ...'); -for i = 1:k +reconASTRA_3D = zeros(size(TomoPhantom),'single'); +for k = 1:N vol_id = astra_mex_data2d('create', '-vol', vol_geom, 0); -proj_id = astra_mex_data2d('create', '-proj3d', proj_geom, sino_artifacts(:,:,k)); +proj_id = astra_mex_data2d('create', '-sino', proj_geom, sino3D_log(:,:,k)'); cfg = astra_struct('FBP_CUDA'); cfg.ProjectionDataId = proj_id; cfg.ReconstructionDataId = vol_id; cfg.option.MinConstraint = 0; alg_id = astra_mex_algorithm('create', cfg); -astra_mex_algorithm('iterate', alg_id, 15); -reconASTRA_3D = astra_mex_data2d('get', vol_id); +astra_mex_algorithm('iterate', alg_id, 1); +rec = astra_mex_data2d('get', vol_id); +reconASTRA_3D(:,:,k) = single(rec); end +figure; imshow(reconASTRA_3D(:,:,128), [0 1.3]); +%% %% -fprintf('%s\n', 'Reconstruction using FISTA-LS without regularization...'); +fprintf('%s\n', 'Reconstruction using OS-FISTA-PWLS without regularization...'); clear params % define parameters params.proj_geom = proj_geom; % pass geometry to the function params.vol_geom = vol_geom; -params.sino = single(sino_tomophan3D); % sinogram -params.iterFISTA = 5; %max number of outer iterations +params.sino = single(sino3D_log); % sinogram +params.iterFISTA = 12; %max number of outer iterations params.X_ideal = TomoPhantom; % ideal phantom +params.weights = dataRaw./max(dataRaw(:)); % statistical weight for PWLS +params.subsets = 12; % the number of subsets params.show = 1; % visualize reconstruction on each iteration -params.subsets = 12; -params.slice = round(N/2); params.maxvalplot = 1; +params.slice = 1; params.maxvalplot = 1.3; tic; [X_FISTA, output] = FISTA_REC(params); toc; error_FISTA = output.Resid_error; obj_FISTA = output.objective; @@ -63,4 +77,33 @@ subplot(1,2,2); imshow(Resid3D(:,:,params.slice),[0 0.1]); title('residual'); c figure(3); subplot(1,2,1); plot(error_FISTA); title('RMSE plot'); subplot(1,2,2); plot(obj_FISTA); title('Objective plot'); +%% +%% +fprintf('%s\n', 'Reconstruction using OS-FISTA-GH without FGP-TV regularization...'); +clear params +% define parameters +params.proj_geom = proj_geom; % pass geometry to the function +params.vol_geom = vol_geom; +params.sino = single(sino3D_log); % sinogram +params.iterFISTA = 15; %max number of outer iterations +params.X_ideal = TomoPhantom; % ideal phantom +params.weights = dataRaw./max(dataRaw(:)); % statistical weight for PWLS +params.subsets = 8; % the number of subsets +params.Regul_Lambda_FGPTV = 0.003; % TV regularization parameter for FGP-TV +params.Ring_LambdaR_L1 = 0.02; % Soft-Thresh L1 ring variable parameter +params.Ring_Alpha = 21; % to boost ring removal procedure +params.show = 1; % visualize reconstruction on each iteration +params.slice = 1; params.maxvalplot = 1.3; +tic; [X_FISTA_GH_TV, output] = FISTA_REC(params); toc; + +error_FISTA_GH_TV = output.Resid_error; obj_FISTA_GH_TV = output.objective; +fprintf('%s %.4f\n', 'Min RMSE for FISTA-PWLS reconstruction is:', min(error_FISTA_GH_TV(:))); + +Resid3D = (TomoPhantom - X_FISTA_GH_TV).^2; +figure(2); +subplot(1,2,1); imshow(X_FISTA_GH_TV(:,:,params.slice),[0 params.maxvalplot]); title('FISTA-LS reconstruction'); colorbar; +subplot(1,2,2); imshow(Resid3D(:,:,params.slice),[0 0.1]); title('residual'); colorbar; +figure(3); +subplot(1,2,1); plot(error_FISTA_GH_TV); title('RMSE plot'); +subplot(1,2,2); plot(obj_FISTA_GH_TV); title('Objective plot'); %% \ No newline at end of file diff --git a/demos/Demo_RealData3D_Parallel.m b/demos/Demo_RealData3D_Parallel.m new file mode 100644 index 0000000..e4c9eb0 --- /dev/null +++ b/demos/Demo_RealData3D_Parallel.m @@ -0,0 +1,186 @@ +% Demonstration of tomographic 3D reconstruction from X-ray synchrotron +% dataset (dendrites) using various data fidelities +% ! 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/regularizers_GPU/NL_Regul/'); addpath('../main_func/regularizers_GPU/Diffus_HO/'); +addpath('../supp/'); + +load('DendrRawData.mat') % load raw data of 3D dendritic set +angles_rad = angles*(pi/180); % conversion to radians +det_size = 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(det_size, 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 +%% +% set projection/reconstruction geometry here +proj_geom = astra_create_proj_geom('parallel', 1, det_size, angles_rad); +vol_geom = astra_create_vol_geom(recon_size,recon_size); +%% +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-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 = 12; +params.weights = Weights3D; +params.subsets = 16; % the number of ordered subsets +params.show = 1; +params.maxvalplot = 2.5; params.slice = 1; + +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-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 = 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, 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-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 = 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 = 1; + +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-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 = 12; +params.Regul_Lambda_FGPTV = 0.005; % TV regularization parameter for FGP-TV +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, 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 +% params.proj_geom = proj_geom; % pass geometry to the function +% params.vol_geom = vol_geom; +% params.sino = Sino3D(:,:,10); +% params.iterFISTA = 50; +% params.L_const = 0.01; % Lipshitz constant +% params.Regul_LambdaTV = 0.008; % TV regularization parameter for FISTA-TV +% params.fidelity = 'student'; % choosing Student t penalty +% params.weights = Weights3D(:,:,10); +% params.show = 0; +% params.initialize = 1; +% params.maxvalplot = 2.5; params.slice = 1; +% +% tic; [X_fistaStudentTV] = FISTA_REC(params); toc; +% figure; imshow(X_fistaStudentTV(:,:,1), [0, 2.5]); title ('FISTA-Student-TV reconstruction'); +%% diff --git a/main_func/FISTA_REC.m b/main_func/FISTA_REC.m index bdaeb18..fa98360 100644 --- a/main_func/FISTA_REC.m +++ b/main_func/FISTA_REC.m @@ -669,6 +669,10 @@ else counterInd = counterInd + numProjSub; end + if (i == 1) + r_old = r; + end + % working with a 'ring vector' if (lambdaR_L1 > 0) r = max(abs(r)-lambdaR_L1, 0).*sign(r); % soft-thresholding operator for ring vector diff --git a/main_func/regularizers_CPU/FGP_TV.c b/main_func/regularizers_CPU/FGP_TV.c index cfe5b9e..e0dae57 100644 --- a/main_func/regularizers_CPU/FGP_TV.c +++ b/main_func/regularizers_CPU/FGP_TV.c @@ -65,7 +65,7 @@ void mexFunction( A = (float *) mxGetData(prhs[0]); /*noisy image (2D/3D) */ lambda = (float) mxGetScalar(prhs[1]); /* regularization parameter */ iter = 50; /* default iterations number */ - epsil = 0.001; /* default tolerance constant */ + epsil = 0.0001; /* default tolerance constant */ methTV = 0; /* default isotropic TV penalty */ if ((nrhs == 3) || (nrhs == 4) || (nrhs == 5)) iter = (int) mxGetScalar(prhs[2]); /* iterations number */ @@ -88,7 +88,7 @@ void mexFunction( tk = 1.0f; tkp1=1.0f; - count = 1; + count = 0; re_old = 0.0f; if (number_of_dims == 2) { @@ -127,7 +127,7 @@ void mexFunction( } re = sqrt(re)/sqrt(re1); if (re < epsil) count++; - if (count > 3) { + if (count > 4) { Obj_func_CALC2D(A, D, funcvalA, lambda, dimX, dimY); break; } diff --git a/supp/sino_add_artifacts.m b/supp/sino_add_artifacts.m new file mode 100644 index 0000000..f601914 --- /dev/null +++ b/supp/sino_add_artifacts.m @@ -0,0 +1,33 @@ +function sino_artifacts = sino_add_artifacts(sino,artifact_type) +% function to add various distortions to the sinogram space, current +% version includes: random rings and zingers (streaks) +% Input: +% 1. sinogram +% 2. artifact type: 'rings' or 'zingers' (streaks) + + +[Detectors, anglesNumb, SlicesZ] = size(sino); +fprintf('%s %i %s %i %s %i %s \n', 'Sinogram has a dimension of', Detectors, 'detectors;', anglesNumb, 'projections;', SlicesZ, 'vertical slices.'); + +sino_artifacts = sino; + +if (strcmp(artifact_type,'rings')) + fprintf('%s \n', 'Adding rings...'); + NumRings = round(Detectors/20); % Number of rings relatively to the size of Detectors + IntenOff = linspace(0.05,0.5,NumRings); % the intensity of rings in the selected range + + for k = 1:SlicesZ + % generate random indices to propagate rings + RandInd = randperm(Detectors,Detectors); + for jj = 1:NumRings + ind_c = RandInd(jj); + sino_artifacts(ind_c,1:end,k) = sino_artifacts(ind_c,1:end,k) + IntenOff(jj).*sino_artifacts(ind_c,1:end,k); % generate a constant offset + end + + end +elseif (strcmp(artifact_type,'zingers')) + fprintf('%s \n', 'Adding zingers...'); +else + fprintf('%s \n', 'Nothing selected, the same sinogram returned...'); +end +end \ No newline at end of file -- cgit v1.2.3 From c91436873e48d531b9313f9c10fa5f89bcb90ab6 Mon Sep 17 00:00:00 2001 From: algol Date: Fri, 20 Oct 2017 16:16:32 +0100 Subject: demo_parallel update --- demos/Demo_Phantom3D_Parallel.m | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/demos/Demo_Phantom3D_Parallel.m b/demos/Demo_Phantom3D_Parallel.m index 9ef1f3c..402bdd2 100644 --- a/demos/Demo_Phantom3D_Parallel.m +++ b/demos/Demo_Phantom3D_Parallel.m @@ -32,6 +32,18 @@ sino3D_log = log(dose./max(dataRaw,1))*multifactor; %log corrected data -> sinog clear dataExp sino_tomophan3D % %% +%-------------Astra toolbox------------% +% one can generate data using ASTRA toolbox +proj_geom = astra_create_proj_geom('parallel', 1, det_size, angles_rad); +vol_geom = astra_create_vol_geom(N,N); +sino_ASTRA3D = zeros(det_size, length(angles), N, 'single'); +for i = 1:N +[sino_id, sinoT] = astra_create_sino_cuda(TomoPhantom(:,:,i), proj_geom, vol_geom); +sino_ASTRA3D(:,:,i) = sinoT'; +astra_mex_data2d('delete', sino_id); +end +%--------------------------------------% +%% % using ASTRA-toolbox to set the projection geometry (parallel beam) proj_geom = astra_create_proj_geom('parallel', 1, det_size, angles_rad); vol_geom = astra_create_vol_geom(N,N); -- cgit v1.2.3