From 52ed0437bfbf5bb8a6fdafd65628c4460184fd2d Mon Sep 17 00:00:00 2001 From: Daniil Kazantsev Date: Tue, 3 Apr 2018 12:34:10 +0100 Subject: cleaning stage, leaving only FGP/ROF and related compilers, demos --- Wrappers/Matlab/demos/Demo_Phantom3D_Cone.m | 67 -------- Wrappers/Matlab/demos/Demo_Phantom3D_Parallel.m | 121 --------------- Wrappers/Matlab/demos/Demo_RealData3D_Parallel.m | 186 ----------------------- Wrappers/Matlab/demos/exportDemoRD2Data.m | 35 ----- 4 files changed, 409 deletions(-) delete mode 100644 Wrappers/Matlab/demos/Demo_Phantom3D_Cone.m delete mode 100644 Wrappers/Matlab/demos/Demo_Phantom3D_Parallel.m delete mode 100644 Wrappers/Matlab/demos/Demo_RealData3D_Parallel.m delete mode 100644 Wrappers/Matlab/demos/exportDemoRD2Data.m (limited to 'Wrappers/Matlab/demos') diff --git a/Wrappers/Matlab/demos/Demo_Phantom3D_Cone.m b/Wrappers/Matlab/demos/Demo_Phantom3D_Cone.m deleted file mode 100644 index a8f2c92..0000000 --- a/Wrappers/Matlab/demos/Demo_Phantom3D_Cone.m +++ /dev/null @@ -1,67 +0,0 @@ -% A demo script to reconstruct 3D synthetic data using FISTA method for -% CONE BEAM geometry -% requirements: ASTRA-toolbox and TomoPhantom toolbox - -close all;clc;clear 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/'); - -%% -% build 3D phantom using TomoPhantom -modelNo = 3; % see Phantom3DLibrary.dat file in TomoPhantom -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 - -%---------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) -vol_geom = astra_create_vol_geom(N,N,N); -proj_geom = astra_create_proj_geom('cone', 1.0, 1.0, N, det_size, angles_rad, 2000, 2160); -%% -% do forward projection using ASTRA -% inverse crime data generation -[sino_id, SinoCone3D] = astra_create_sino3d_cuda(TomoPhantom, proj_geom, vol_geom); -astra_mex_data3d('delete', sino_id); -%% -fprintf('%s\n', 'Reconstructing with CGLS using ASTRA-toolbox ...'); -vol_id = astra_mex_data3d('create', '-vol', vol_geom, 0); -proj_id = astra_mex_data3d('create', '-proj3d', proj_geom, SinoCone3D); -cfg = astra_struct('CGLS3D_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_data3d('get', vol_id); -%% -fprintf('%s\n', 'Reconstruction using FISTA-LS without regularization...'); -clear params -% define parameters -params.proj_geom = proj_geom; % pass geometry to the function -params.vol_geom = vol_geom; -params.sino = single(SinoCone3D); % sinogram -params.iterFISTA = 30; %max number of outer iterations -params.X_ideal = TomoPhantom; % ideal phantom -params.show = 1; % visualize reconstruction on each iteration -params.slice = round(N/2); params.maxvalplot = 1; -tic; [X_FISTA, output] = FISTA_REC(params); toc; - -error_FISTA = output.Resid_error; obj_FISTA = output.objective; -fprintf('%s %.4f\n', 'Min RMSE for FISTA-LS reconstruction is:', min(error_FISTA(:))); - -Resid3D = (TomoPhantom - X_FISTA).^2; -figure(2); -subplot(1,2,1); imshow(X_FISTA(:,:,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); title('RMSE plot'); colorbar; -subplot(1,2,2); plot(obj_FISTA); title('Objective plot'); colorbar; -%% \ No newline at end of file diff --git a/Wrappers/Matlab/demos/Demo_Phantom3D_Parallel.m b/Wrappers/Matlab/demos/Demo_Phantom3D_Parallel.m deleted file mode 100644 index 4219bd1..0000000 --- a/Wrappers/Matlab/demos/Demo_Phantom3D_Parallel.m +++ /dev/null @@ -1,121 +0,0 @@ -% A demo script to reconstruct 3D synthetic data using FISTA method for -% PARALLEL BEAM geometry -% requirements: ASTRA-toolbox and TomoPhantom toolbox - -close all;clc;clear; -% 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/'); - -%% -% 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 - -%---------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_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 -% -%% -%-------------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); -%% -fprintf('%s\n', 'Reconstructing with FBP using ASTRA-toolbox ...'); -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', '-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, 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 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(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 = 12; % the number of subsets -params.show = 1; % visualize reconstruction on each iteration -params.slice = 128; params.maxvalplot = 1.3; -tic; [X_FISTA, output] = FISTA_REC(params); toc; - -error_FISTA = output.Resid_error; obj_FISTA = output.objective; -fprintf('%s %.4f\n', 'Min RMSE for FISTA-PWLS reconstruction is:', min(error_FISTA(:))); - -Resid3D = (TomoPhantom - X_FISTA).^2; -figure(2); -subplot(1,2,1); imshow(X_FISTA(:,:,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); title('RMSE plot'); -subplot(1,2,2); plot(obj_FISTA); title('Objective plot'); -%% -%% -fprintf('%s\n', 'Reconstruction using OS-FISTA-GH with 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 weights for PWLS -params.subsets = 12; % the number of subsets -params.Regul_Lambda_FGPTV = 100; % 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 = 128; 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/Wrappers/Matlab/demos/Demo_RealData3D_Parallel.m b/Wrappers/Matlab/demos/Demo_RealData3D_Parallel.m deleted file mode 100644 index f82e0b0..0000000 --- a/Wrappers/Matlab/demos/Demo_RealData3D_Parallel.m +++ /dev/null @@ -1,186 +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 -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 = 18; -params.weights = Weights3D; -params.subsets = 8; % 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 = 18; -params.Regul_Lambda_FGPTV = 5.0000e+6; % TV regularization parameter for FGP-TV -params.weights = Weights3D; -params.subsets = 8; % the number of ordered subsets -params.show = 1; -params.maxvalplot = 2.5; params.slice = 10; - -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(:,:,10); -params.iterFISTA = 18; -params.Regul_Lambda_FGPTV = 5.0000e+6; % 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(:,:,10); -params.subsets = 8; % 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 = 5.0000e+6; % 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/Wrappers/Matlab/demos/exportDemoRD2Data.m b/Wrappers/Matlab/demos/exportDemoRD2Data.m deleted file mode 100644 index 028353b..0000000 --- a/Wrappers/Matlab/demos/exportDemoRD2Data.m +++ /dev/null @@ -1,35 +0,0 @@ -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