From 0847a315ce744e52be3dade398fb16c58323084e Mon Sep 17 00:00:00 2001 From: Daniil Kazantsev Date: Wed, 18 Oct 2017 13:44:57 +0100 Subject: linked demos to TomoPhantom, cone beam demo, some FISTA modificastions --- demos/Demo1.m | 174 ---------------------------------------- demos/Demo_Phantom3D_Cone.m | 66 +++++++++++++++ demos/Demo_Phantom3D_Parallel.m | 49 +++++++++++ 3 files changed, 115 insertions(+), 174 deletions(-) delete mode 100644 demos/Demo1.m create mode 100644 demos/Demo_Phantom3D_Cone.m create mode 100644 demos/Demo_Phantom3D_Parallel.m (limited to 'demos') diff --git a/demos/Demo1.m b/demos/Demo1.m deleted file mode 100644 index 15e2e5b..0000000 --- a/demos/Demo1.m +++ /dev/null @@ -1,174 +0,0 @@ -% Demonstration of tomographic reconstruction from noisy and corrupted by -% artifacts undersampled projection data using Students't penalty -% Optimisation problem is solved using FISTA algorithm (see Beck & Teboulle) - -% see Readme file for instructions -%% -% compile MEX-files ones -% cd .. -% cd main_func -% compile_mex -% cd .. -% cd demos -%% - -close all;clc;clear all; -% adding paths -addpath('../data/'); -addpath('../main_func/'); addpath('../main_func/regularizers_CPU/'); -addpath('../supp/'); - -load phantom_bone512.mat % load the phantom -load my_red_yellowMAP.mat % load the colormap -% load sino1.mat; % load noisy sinogram - -N = 512; % the size of the tomographic image NxN -theta = 1:1:180; % acquisition angles (in parallel beam from 0 to Pi) -theta_rad = theta*(pi/180); % conversion to radians -P = 2*ceil(N/sqrt(2))+1; % the size of the detector array -ROI = find(phantom > 0); - -% using ASTRA to set the projection geometry -% potentially parallel geometry can be replaced with a divergent one -Z_slices = 1; -det_row_count = Z_slices; -proj_geom = astra_create_proj_geom('parallel3d', 1, 1, det_row_count, P, theta_rad); -vol_geom = astra_create_vol_geom(N,N,Z_slices); - -zing_rings_add; % generating data, adding zingers and stripes -%% -fprintf('%s\n', 'Direct reconstruction using FBP...'); -FBP_1 = iradon(sino_zing_rings', theta, N); - -fprintf('%s %.4f\n', 'RMSE for FBP reconstruction:', RMSE(FBP_1(:), phantom(:))); - -figure(1); -subplot_tight(1,2,1, [0.05 0.05]); imshow(FBP_1,[0 0.6]); title('FBP reconstruction of noisy and corrupted by artifacts sinogram'); colorbar; -subplot_tight(1,2,2, [0.05 0.05]); imshow((phantom - FBP_1).^2,[0 0.1]); title('residual: (ideal phantom - FBP)^2'); colorbar; -colormap(cmapnew); - -%% -fprintf('%s\n', 'Reconstruction using 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 = sino_zing_rings; % sinogram -params.iterFISTA = 45; %max number of outer iterations -params.X_ideal = phantom; % ideal phantom -params.ROI = ROI; % phantom region-of-interest -params.show = 1; % visualize reconstruction on each iteration -params.slice = 1; params.maxvalplot = 0.6; -params.weights = Dweights; % statistical weighting -tic; [X_FISTA, output] = FISTA_REC(params); toc; - -fprintf('%s %.4f\n', 'Min RMSE for FISTA-PWLS reconstruction is:', min(error_FISTA(:))); -error_FISTA = output.Resid_error; obj_FISTA = output.objective; - -figure(2); clf -%set(gcf, 'Position', get(0,'Screensize')); -subplot(1,2,1, [0.05 0.05]); imshow(X_FISTA,[0 0.6]); title('FISTA-PWLS reconstruction'); colorbar; -subplot(1,2,2, [0.05 0.05]); imshow((phantom - X_FISTA).^2,[0 0.1]); title('residual'); colorbar; -colormap(cmapnew); -figure(3); clf -subplot(1,2,1, [0.05 0.05]); plot(error_FISTA); title('RMSE plot'); colorbar; -subplot(1,2,2, [0.05 0.05]); plot(obj_FISTA); title('Objective plot'); colorbar; -colormap(cmapnew); -%% -fprintf('%s\n', 'Reconstruction using FISTA-PWLS-TV...'); -clear params -% define parameters -params.proj_geom = proj_geom; % pass geometry to the function -params.vol_geom = vol_geom; -params.sino = sino_zing_rings; -params.iterFISTA = 45; % max number of outer iterations -params.Regul_LambdaTV = 0.0015; % regularization parameter for TV problem -params.X_ideal = phantom; % ideal phantom -params.ROI = ROI; % phantom region-of-interest -params.weights = Dweights; % statistical weighting -params.show = 1; % visualize reconstruction on each iteration -params.slice = 1; params.maxvalplot = 0.6; -tic; [X_FISTA_TV, output] = FISTA_REC(params); toc; - -fprintf('%s %.4f\n', 'Min RMSE for FISTA-PWLS-TV reconstruction is:', min(error_FISTA_TV(:))); -error_FISTA_TV = output.Resid_error; obj_FISTA_TV = output.objective; - -figure(4); clf -subplot(1,2,1, [0.05 0.05]); imshow(X_FISTA_TV,[0 0.6]); title('FISTA-PWLS-TV reconstruction'); colorbar; -subplot(1,2,2, [0.05 0.05]); imshow((phantom - X_FISTA_TV).^2,[0 0.1]); title('residual'); colorbar; -colormap(cmapnew); -figure(5); clf -subplot(1,2,1, [0.05 0.05]); plot(error_FISTA_TV); title('RMSE plot'); colorbar; -subplot(1,2,2, [0.05 0.05]); plot(obj_FISTA_TV); title('Objective plot'); colorbar; -colormap(cmapnew); -%% -fprintf('%s\n', 'Reconstruction using FISTA-GH-TV...'); -clear params -% define parameters -params.proj_geom = proj_geom; % pass geometry to the function -params.vol_geom = vol_geom; -params.sino = sino_zing_rings; -params.iterFISTA = 50; % max number of outer iterations -params.Regul_LambdaTV = 0.0015; % regularization parameter for TV problem -params.X_ideal = phantom; % ideal phantom -params.ROI = ROI; % phantom region-of-interest -params.weights = Dweights; % statistical weighting -params.Ring_LambdaR_L1 = 0.002; % parameter to sparsify the "rings vector" -params.Ring_Alpha = 20; % to accelerate ring-removal procedure -params.show = 0; % visualize reconstruction on each iteration -params.slice = 1; params.maxvalplot = 0.6; -tic; [X_FISTA_GH_TV, output] = FISTA_REC(params); toc; - -fprintf('%s %.4f\n', 'Min RMSE for FISTA-GH-TV reconstruction is:', min(error_FISTA_GH_TV(:))); -error_FISTA_GH_TV = output.Resid_error; obj_FISTA_GH_TV = output.objective; - -figure(6); clf -subplot(1,2,1, [0.05 0.05]); imshow(X_FISTA_GH_TV,[0 0.6]); title('FISTA-GH-TV reconstruction'); colorbar; -subplot(1,2,2, [0.05 0.05]);imshow((phantom - X_FISTA_GH_TV).^2,[0 0.1]); title('residual'); colorbar; -colormap(cmapnew); - -figure(7); clf -subplot(1,2,1, [0.05 0.05]); plot(error_FISTA_GH_TV); title('RMSE plot'); colorbar; -subplot(1,2,2, [0.05 0.05]); plot(obj_FISTA_GH_TV); title('Objective plot'); colorbar; -colormap(cmapnew); -%% -fprintf('%s\n', 'Reconstruction using FISTA-Student-TV...'); -clear params -% define parameters -params.proj_geom = proj_geom; % pass geometry to the function -params.vol_geom = vol_geom; -params.sino = sino_zing_rings; -params.iterFISTA = 55; % max number of outer iterations -params.L_const = 0.1; % Lipshitz constant (can be chosen manually to accelerate convergence) -params.Regul_LambdaTV = 0.00152; % regularization parameter for TV problem -params.X_ideal = phantom; % ideal phantom -params.ROI = ROI; % phantom region-of-interest -params.weights = Dweights; % statistical weighting -params.fidelity = 'student'; % selecting students t fidelity -params.show = 1; % visualize reconstruction on each iteration -params.slice = 1; params.maxvalplot = 0.6; -params.initilize = 1; % warm start with SIRT -tic; [X_FISTA_student_TV, output] = FISTA_REC(params); toc; - -fprintf('%s %.4f\n', 'Min RMSE for FISTA-Student-TV reconstruction is:', min(error_FISTA_student_TV(:))); -error_FISTA_student_TV = output.Resid_error; obj_FISTA_student_TV = output.objective; - -figure(8); -set(gcf, 'Position', get(0,'Screensize')); -subplot(1,2,1, [0.05 0.05]); imshow(X_FISTA_student_TV,[0 0.6]); title('FISTA-Student-TV reconstruction'); colorbar; -subplot(1,2,2, [0.05 0.05]); imshow((phantom - X_FISTA_student_TV).^2,[0 0.1]); title('residual'); colorbar; -colormap(cmapnew); - -figure(9); -subplot(1,2,1, [0.05 0.05]); plot(error_FISTA_student_TV); title('RMSE plot'); colorbar; -subplot(1,2,2, [0.05 0.05]); plot(obj_FISTA_student_TV); title('Objective plot'); colorbar; -colormap(cmapnew); -%% -% print all RMSE's -fprintf('%s\n', '--------------------------------------------'); -fprintf('%s %.4f\n', 'RMSE for FBP reconstruction:', RMSE(FBP_1(:), phantom(:))); -fprintf('%s %.4f\n', 'Min RMSE for FISTA-PWLS reconstruction:', min(error_FISTA(:))); -fprintf('%s %.4f\n', 'Min RMSE for FISTA-PWLS-TV reconstruction:', min(error_FISTA_TV(:))); -fprintf('%s %.4f\n', 'Min RMSE for FISTA-GH-TV reconstruction:', min(error_FISTA_GH_TV(:))); -fprintf('%s %.4f\n', 'Min RMSE for FISTA-Student-TV reconstruction:', min(error_FISTA_student_TV(:))); -% \ No newline at end of file diff --git a/demos/Demo_Phantom3D_Cone.m b/demos/Demo_Phantom3D_Cone.m new file mode 100644 index 0000000..6419386 --- /dev/null +++ b/demos/Demo_Phantom3D_Cone.m @@ -0,0 +1,66 @@ +% 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 +% in order to run functions you have to go to the directory: +cd /home/algol/Documents/MATLAB/TomoPhantom/functions/ +TomoPhantom = buildPhantom3D(modelNo,N); % 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/demos/Demo_Phantom3D_Parallel.m b/demos/Demo_Phantom3D_Parallel.m new file mode 100644 index 0000000..fd8096a --- /dev/null +++ b/demos/Demo_Phantom3D_Parallel.m @@ -0,0 +1,49 @@ +% 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 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 and generate projection data +modelNo = 3; % 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: +cd /home/algol/Documents/MATLAB/TomoPhantom/functions/ +TomoPhantom = buildPhantom3D(modelNo,N); % generate 3D phantom +sino_tomophan3D = buildSino3D(modelNo, N, det_size, single(angles)); % generate data +%% +% 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', '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(sino_tomophan3D); % sinogram +params.iterFISTA = 5; %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-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'); colorbar; +subplot(1,2,2); plot(obj_FISTA); title('Objective plot'); colorbar; +%% \ No newline at end of file -- cgit v1.2.3 From cb8ef11f00e897b6f5b3049126dd32baa3c50cf9 Mon Sep 17 00:00:00 2001 From: Daniil Kazantsev Date: Wed, 18 Oct 2017 21:46:29 +0100 Subject: ordered subsets fix for GH term --- demos/Demo_Phantom3D_Parallel.m | 5 +- main_func/FISTA_REC.m | 179 +++++++++++++++++++++++----------------- 2 files changed, 108 insertions(+), 76 deletions(-) (limited to 'demos') diff --git a/demos/Demo_Phantom3D_Parallel.m b/demos/Demo_Phantom3D_Parallel.m index fd8096a..ac9827c 100644 --- a/demos/Demo_Phantom3D_Parallel.m +++ b/demos/Demo_Phantom3D_Parallel.m @@ -33,6 +33,7 @@ params.sino = single(sino_tomophan3D); % sinogram params.iterFISTA = 5; %max number of outer iterations params.X_ideal = TomoPhantom; % ideal phantom params.show = 1; % visualize reconstruction on each iteration +params.subsets = 12; params.slice = round(N/2); params.maxvalplot = 1; tic; [X_FISTA, output] = FISTA_REC(params); toc; @@ -44,6 +45,6 @@ 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; +subplot(1,2,1); plot(error_FISTA); title('RMSE plot'); +subplot(1,2,2); plot(obj_FISTA); title('Objective plot'); %% \ No newline at end of file diff --git a/main_func/FISTA_REC.m b/main_func/FISTA_REC.m index dde0e73..bea1860 100644 --- a/main_func/FISTA_REC.m +++ b/main_func/FISTA_REC.m @@ -106,14 +106,14 @@ if (isfield(params,'L_const')) else % using Power method (PM) to establish L constant fprintf('%s %s %s \n', 'Calculating Lipshitz constant for',proj_geom.type, 'beam geometry...'); - if (strcmp(proj_geom.type,'parallel') || strcmp(proj_geom.type,'fanflat') || strcmp(proj_geom.type,'fanflat_vec')) - % for 2D geometry we can do just one selected slice + if (strcmp(proj_geom.type,'parallel') || strcmp(proj_geom.type,'fanflat') || strcmp(proj_geom.type,'fanflat_vec')) + % for 2D geometry we can do just one selected slice niter = 15; % number of iteration for the PM x1 = rand(N,N,1); sqweight = sqrt(weights(:,:,1)); [sino_id, y] = astra_create_sino_cuda(x1, proj_geom, vol_geom); y = sqweight.*y'; - astra_mex_data2d('delete', sino_id); + astra_mex_data2d('delete', sino_id); for i = 1:niter [x1] = astra_create_backprojection_cuda((sqweight.*y)', proj_geom, vol_geom); s = norm(x1(:)); @@ -121,9 +121,9 @@ else [sino_id, y] = astra_create_sino_cuda(x1, proj_geom, vol_geom); y = sqweight.*y'; astra_mex_data2d('delete', sino_id); - end + end elseif (strcmp(proj_geom.type,'cone') || strcmp(proj_geom.type,'parallel3d') || strcmp(proj_geom.type,'parallel3d_vec') || strcmp(proj_geom.type,'cone_vec')) - % 3D geometry + % 3D geometry niter = 8; % number of iteration for PM x1 = rand(N,N,SlicesZ); sqweight = sqrt(weights); @@ -268,7 +268,7 @@ if (isfield(params,'initialize')) X = params.initialize; if ((size(X,1) ~= N) || (size(X,2) ~= N) || (size(X,3) ~= SlicesZ)) error('%s \n', 'The initialized volume has different dimensions!'); - end + end else X = zeros(N,N,SlicesZ, 'single'); % storage for the solution end @@ -320,10 +320,11 @@ if (subsets == 0) t_old = t; r_old = r; - % if the geometry is 2D use slice-by-slice projection-backprojection routine - if (strcmp(proj_geom.type,'parallel') || strcmp(proj_geom.type,'fanflat') || strcmp(proj_geom.type,'fanflat_vec')) + + if (strcmp(proj_geom.type,'parallel') || strcmp(proj_geom.type,'fanflat') || strcmp(proj_geom.type,'fanflat_vec')) + % if geometry is 2D use slice-by-slice projection-backprojection routine sino_updt = zeros(size(sino),'single'); - for kkk = 1:SlicesZ + for kkk = 1:SlicesZ [sino_id, sinoT] = astra_create_sino_cuda(X_t(:,:,kkk), proj_geom, vol_geom); sino_updt(:,:,kkk) = sinoT'; astra_mex_data2d('delete', sino_id); @@ -359,12 +360,11 @@ if (subsets == 0) else % no ring removal (LS model) residual = weights.*(sino_updt - sino); - objective(i) = (0.5*sum(residual(:).^2)); % for the objective function output + objective(i) = 0.5*norm(residual(:)); % for the objective function output end - % if the geometry is 2D use slice-by-slice projection-backprojection routine - if (strcmp(proj_geom.type,'parallel') || strcmp(proj_geom.type,'fanflat') || strcmp(proj_geom.type,'fanflat_vec')) + if (strcmp(proj_geom.type,'parallel') || strcmp(proj_geom.type,'fanflat') || strcmp(proj_geom.type,'fanflat_vec')) x_temp = zeros(size(X),'single'); for kkk = 1:SlicesZ [x_temp(:,:,kkk)] = astra_create_backprojection_cuda(squeeze(residual(:,:,kkk))', proj_geom, vol_geom); @@ -373,9 +373,9 @@ if (subsets == 0) [id, x_temp] = astra_create_backprojection3d_cuda(residual, proj_geom, vol_geom); astra_mex_data3d('delete', id); end - X = X_t - (1/L_const).*x_temp; + X = X_t - (1/L_const).*x_temp; - % ----------------Regularization part------------------------ + % ----------------Regularization part------------------------% if (lambdaFGP_TV > 0) % FGP-TV regularization if ((strcmp('2D', Dimension) == 1)) @@ -484,94 +484,128 @@ if (subsets == 0) end end else - % Ordered Subsets (OS) FISTA reconstruction routine (normally one order of magnitude faster than classical) + % Ordered Subsets (OS) FISTA reconstruction routine (normally one order of magnitude faster than the classical version) t = 1; X_t = X; - proj_geomSUB = proj_geom; - + proj_geomSUB = proj_geom; r = zeros(Detectors,SlicesZ, 'single'); % 2D array (for 3D data) of sparse "ring" vectors r_x = r; % another ring variable residual2 = zeros(size(sino),'single'); + sino_updt_FULL = zeros(size(sino),'single'); % Outer FISTA iterations loop - for i = 1:iterFISTA + for i = 1:iterFISTA - % With OS approach it becomes trickier to correlate independent subsets, hence additional work is required - % one solution is to work with a full sinogram at times - if ((i >= 3) && (lambdaR_L1 > 0)) - [sino_id2, sino_updt2] = astra_create_sino3d_cuda(X, proj_geom, vol_geom); - astra_mex_data3d('delete', sino_id2); + if ((i > 1) && (lambdaR_L1 > 0)) + % in order to make Group-Huber fidelity work with ordered subsets + % we still need to work with full sinogram + + % the offset variable must be calculated for the whole + % updated sinogram - sino_updt_FULL + for kkk = 1:anglesNumb + residual2(:,kkk,:) = squeeze(weights(:,kkk,:)).*(squeeze(sino_updt_FULL(:,kkk,:)) - (squeeze(sino(:,kkk,:)) - alpha_ring.*r_x)); + end + + r_old = r; + vec = sum(residual2,2); + if (SlicesZ > 1) + vec = squeeze(vec(:,1,:)); + end + r = r_x - (1./L_const).*vec; % update ring variable end % subsets loop counterInd = 1; for ss = 1:subsets X_old = X; - t_old = t; - r_old = r; + t_old = t; numProjSub = binsDiscr(ss); % the number of projections per subset CurrSubIndeces = IndicesReorg(counterInd:(counterInd + numProjSub - 1)); % extract indeces attached to the subset proj_geomSUB.ProjectionAngles = angles(CurrSubIndeces); + sino_updt_Sub = zeros(Detectors, numProjSub, SlicesZ,'single'); if (lambdaR_L1 > 0) - - % the ring removal part (Group-Huber fidelity) - % first 2 iterations do additional work reconstructing whole dataset to ensure - % the stablility - if (i < 3) - [sino_id2, sino_updt2] = astra_create_sino3d_cuda(X_t, proj_geom, vol_geom); - astra_mex_data3d('delete', sino_id2); + % Group-Huber fidelity (ring removal) + if (strcmp(proj_geom.type,'parallel') || strcmp(proj_geom.type,'fanflat') || strcmp(proj_geom.type,'fanflat_vec')) + % if geometry is 2D use slice-by-slice projection-backprojection routine + for kkk = 1:SlicesZ + [sino_id, sinoT] = astra_create_sino_cuda(X_t(:,:,kkk), proj_geomSUB, vol_geom); + sino_updt_Sub(:,:,kkk) = sinoT'; + astra_mex_data2d('delete', sino_id); + end else - [sino_id, sino_updt] = astra_create_sino3d_cuda(X_t, proj_geomSUB, vol_geom); - end - - for kkk = 1:anglesNumb - residual2(:,kkk,:) = squeeze(weights(:,kkk,:)).*(squeeze(sino_updt2(:,kkk,:)) - (squeeze(sino(:,kkk,:)) - alpha_ring.*r_x)); + % for 3D geometry (watch the GPU memory overflow in earlier ASTRA versions < 1.8) + [sino_id, sino_updt_Sub] = astra_create_sino3d_cuda(X_t, proj_geomSUB, vol_geom); + astra_mex_data3d('delete', sino_id); end - residual = zeros(Detectors, numProjSub, SlicesZ,'single'); + residualSub = zeros(Detectors, numProjSub, SlicesZ,'single'); % residual for a chosen subset for kkk = 1:numProjSub indC = CurrSubIndeces(kkk); - if (i < 3) - residual(:,kkk,:) = squeeze(residual2(:,indC,:)); - else - residual(:,kkk,:) = squeeze(weights(:,indC,:)).*(squeeze(sino_updt(:,kkk,:)) - (squeeze(sino(:,indC,:)) - alpha_ring.*r_x)); + residualSub(:,kkk,:) = squeeze(weights(:,indC,:)).*(squeeze(sino_updt_Sub(:,kkk,:)) - (squeeze(sino(:,indC,:)) - alpha_ring.*r_x)); + sino_updt_FULL(:,indC,:) = squeeze(sino_updt_Sub(:,kkk,:)); % filling the full sinogram + end + + elseif (studentt > 0) + % student t data fidelity + if (strcmp(proj_geom.type,'parallel') || strcmp(proj_geom.type,'fanflat') || strcmp(proj_geom.type,'fanflat_vec')) + % if geometry is 2D use slice-by-slice projection-backprojection routine + for kkk = 1:SlicesZ + [sino_id, sinoT] = astra_create_sino_cuda(X_t(:,:,kkk), proj_geomSUB, vol_geom); + sino_updt_Sub(:,:,kkk) = sinoT'; + astra_mex_data2d('delete', sino_id); end + else + % for 3D geometry (watch the GPU memory overflow in earlier ASTRA versions < 1.8) + [sino_id, sino_updt_Sub] = astra_create_sino3d_cuda(X_t, proj_geomSUB, vol_geom); + astra_mex_data3d('delete', sino_id); end - vec = sum(residual2,2); - if (SlicesZ > 1) - vec = squeeze(vec(:,1,:)); + + % artifacts removal with Students t penalty + residualSub = squeeze(weights(:,CurrSubIndeces,:)).*(sino_updt_Sub - squeeze(sino(:,CurrSubIndeces,:))); + + for kkk = 1:SlicesZ + res_vec = reshape(residualSub(:,:,kkk), Detectors*numProjSub, 1); % 1D vectorized sinogram + %s = 100; + %gr = (2)*res_vec./(s*2 + conj(res_vec).*res_vec); + [ff, gr] = studentst(res_vec, 1); + residualSub(:,:,kkk) = reshape(gr, Detectors, numProjSub); end - r = r_x - (1./L_const).*vec; + objective(i) = ff; % for the objective function output else - [sino_id, sino_updt] = astra_create_sino3d_cuda(X_t, proj_geomSUB, vol_geom); - - if (studentt == 1) - % artifacts removal with Students t penalty - residual = squeeze(weights(:,CurrSubIndeces,:)).*(sino_updt - squeeze(sino(:,CurrSubIndeces,:))); - + % PWLS model + if (strcmp(proj_geom.type,'parallel') || strcmp(proj_geom.type,'fanflat') || strcmp(proj_geom.type,'fanflat_vec')) + % if geometry is 2D use slice-by-slice projection-backprojection routine for kkk = 1:SlicesZ - res_vec = reshape(residual(:,:,kkk), Detectors*numProjSub, 1); % 1D vectorized sinogram - %s = 100; - %gr = (2)*res_vec./(s*2 + conj(res_vec).*res_vec); - [ff, gr] = studentst(res_vec, 1); - residual(:,:,kkk) = reshape(gr, Detectors, numProjSub); + [sino_id, sinoT] = astra_create_sino_cuda(X_t(:,:,kkk), proj_geomSUB, vol_geom); + sino_updt_Sub(:,:,kkk) = sinoT'; + astra_mex_data2d('delete', sino_id); end - objective(i) = ff; % for the objective function output else - % no ring removal (LS model) - residual = squeeze(weights(:,CurrSubIndeces,:)).*(sino_updt - squeeze(sino(:,CurrSubIndeces,:))); + % for 3D geometry (watch the GPU memory overflow in earlier ASTRA versions < 1.8) + [sino_id, sino_updt_Sub] = astra_create_sino3d_cuda(X_t, proj_geomSUB, vol_geom); + astra_mex_data3d('delete', sino_id); end + residualSub = squeeze(weights(:,CurrSubIndeces,:)).*(sino_updt_Sub - squeeze(sino(:,CurrSubIndeces,:))); + objective(i) = 0.5*norm(residualSub(:)); % for the objective function output end - [id, x_temp] = astra_create_backprojection3d_cuda(residual, proj_geomSUB, vol_geom); - X = X_t - (1/L_const).*x_temp; - astra_mex_data3d('delete', sino_id); - astra_mex_data3d('delete', id); + if (strcmp(proj_geom.type,'parallel') || strcmp(proj_geom.type,'fanflat') || strcmp(proj_geom.type,'fanflat_vec')) + % if geometry is 2D use slice-by-slice projection-backprojection routine + x_temp = zeros(size(X),'single'); + for kkk = 1:SlicesZ + [x_temp(:,:,kkk)] = astra_create_backprojection_cuda(squeeze(residualSub(:,:,kkk))', proj_geomSUB, vol_geom); + end + else + [id, x_temp] = astra_create_backprojection3d_cuda(residualSub, proj_geomSUB, vol_geom); + astra_mex_data3d('delete', id); + end + + X = X_t - (1/L_const).*x_temp; - % regularization + % ----------------Regularization part------------------------% if (lambdaFGP_TV > 0) % FGP-TV regularization if ((strcmp('2D', Dimension) == 1)) @@ -653,20 +687,17 @@ else end end - if (lambdaR_L1 > 0) - r = max(abs(r)-lambdaR_L1, 0).*sign(r); % soft-thresholding operator for ring vector - end - t = (1 + sqrt(1 + 4*t^2))/2; % updating t X_t = X + ((t_old-1)/t).*(X - X_old); % updating X - - if (lambdaR_L1 > 0) - r_x = r + ((t_old-1)/t).*(r - r_old); % updating r - end - counterInd = counterInd + numProjSub; 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 + r_x = r + ((t_old-1)/t).*(r - r_old); % updating r + end + if (show == 1) figure(10); imshow(X(:,:,slice), [0 maxvalplot]); if (lambdaR_L1 > 0) -- cgit v1.2.3 From ff543b21e175d49b95ea6da9c03f21a1789f6833 Mon Sep 17 00:00:00 2001 From: Daniil Kazantsev Date: Thu, 19 Oct 2017 13:08:13 +0100 Subject: OS loop now fixed --- demos/Demo_Phantom3D_Cone.m | 4 ++-- demos/Demo_Phantom3D_Parallel.m | 24 +++++++++++++++++++---- main_func/FISTA_REC.m | 43 +++++++++++++++++++++-------------------- 3 files changed, 44 insertions(+), 27 deletions(-) (limited to 'demos') diff --git a/demos/Demo_Phantom3D_Cone.m b/demos/Demo_Phantom3D_Cone.m index 6419386..3a9178b 100644 --- a/demos/Demo_Phantom3D_Cone.m +++ b/demos/Demo_Phantom3D_Cone.m @@ -17,8 +17,8 @@ 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: -cd /home/algol/Documents/MATLAB/TomoPhantom/functions/ -TomoPhantom = buildPhantom3D(modelNo,N); % generate 3D phantom +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 ac9827c..6a54450 100644 --- a/demos/Demo_Phantom3D_Parallel.m +++ b/demos/Demo_Phantom3D_Parallel.m @@ -10,19 +10,35 @@ addpath('../supp/'); %% % build 3D phantom using TomoPhantom and generate projection data -modelNo = 3; % see Phantom3DLibrary.dat file in TomoPhantom +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: -cd /home/algol/Documents/MATLAB/TomoPhantom/functions/ -TomoPhantom = buildPhantom3D(modelNo,N); % generate 3D phantom -sino_tomophan3D = buildSino3D(modelNo, N, det_size, single(angles)); % generate data +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'); +% %% % 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 ...'); +for i = 1:k +vol_id = astra_mex_data2d('create', '-vol', vol_geom, 0); +proj_id = astra_mex_data2d('create', '-proj3d', proj_geom, sino_artifacts(:,:,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); +end %% fprintf('%s\n', 'Reconstruction using FISTA-LS without regularization...'); clear params diff --git a/main_func/FISTA_REC.m b/main_func/FISTA_REC.m index d177748..bdaeb18 100644 --- a/main_func/FISTA_REC.m +++ b/main_func/FISTA_REC.m @@ -487,7 +487,7 @@ else % Ordered Subsets (OS) FISTA reconstruction routine (normally one order of magnitude faster than the classical version) t = 1; X_t = X; - proj_geomSUB = proj_geom; + proj_geomSUB = proj_geom; r = zeros(Detectors,SlicesZ, 'single'); % 2D array (for 3D data) of sparse "ring" vectors r_x = r; % another ring variable @@ -495,7 +495,7 @@ else sino_updt_FULL = zeros(size(sino),'single'); % Outer FISTA iterations loop - for i = 1:iterFISTA + for i = 1:iterFISTA if ((i > 1) && (lambdaR_L1 > 0)) % in order to make Group-Huber fidelity work with ordered subsets @@ -517,31 +517,31 @@ else % subsets loop counterInd = 1; - if (strcmp(proj_geom.type,'parallel') || strcmp(proj_geom.type,'fanflat') || strcmp(proj_geom.type,'fanflat_vec')) - % if geometry is 2D use slice-by-slice projection-backprojection routine - for kkk = 1:SlicesZ - [sino_id, sinoT] = astra_create_sino_cuda(X_t(:,:,kkk), proj_geomSUB, vol_geom); - sino_updt_Sub(:,:,kkk) = sinoT'; - astra_mex_data2d('delete', sino_id); - end - else - % for 3D geometry (watch the GPU memory overflow in earlier ASTRA versions < 1.8) - [sino_id, sino_updt_Sub] = astra_create_sino3d_cuda(X_t, proj_geomSUB, vol_geom); - astra_mex_data3d('delete', sino_id); - end for ss = 1:subsets X_old = X; - t_old = t; + t_old = t; numProjSub = binsDiscr(ss); % the number of projections per subset + sino_updt_Sub = zeros(Detectors, numProjSub, SlicesZ,'single'); CurrSubIndeces = IndicesReorg(counterInd:(counterInd + numProjSub - 1)); % extract indeces attached to the subset proj_geomSUB.ProjectionAngles = angles(CurrSubIndeces); - sino_updt_Sub = zeros(Detectors, numProjSub, SlicesZ,'single'); + + if (strcmp(proj_geom.type,'parallel') || strcmp(proj_geom.type,'fanflat') || strcmp(proj_geom.type,'fanflat_vec')) + % if geometry is 2D use slice-by-slice projection-backprojection routine + for kkk = 1:SlicesZ + [sino_id, sinoT] = astra_create_sino_cuda(X_t(:,:,kkk), proj_geomSUB, vol_geom); + sino_updt_Sub(:,:,kkk) = sinoT'; + astra_mex_data2d('delete', sino_id); + end + else + % for 3D geometry (watch the GPU memory overflow in earlier ASTRA versions < 1.8) + [sino_id, sino_updt_Sub] = astra_create_sino3d_cuda(X_t, proj_geomSUB, vol_geom); + astra_mex_data3d('delete', sino_id); + end if (lambdaR_L1 > 0) % Group-Huber fidelity (ring removal) - residualSub = zeros(Detectors, numProjSub, SlicesZ,'single'); % residual for a chosen subset for kkk = 1:numProjSub indC = CurrSubIndeces(kkk); @@ -551,8 +551,6 @@ else elseif (studentt > 0) % student t data fidelity - - % artifacts removal with Students t penalty residualSub = squeeze(weights(:,CurrSubIndeces,:)).*(sino_updt_Sub - squeeze(sino(:,CurrSubIndeces,:))); @@ -566,8 +564,11 @@ else objective(i) = ff; % for the objective function output else % PWLS model - + residualSub = squeeze(weights(:,CurrSubIndeces,:)).*(sino_updt_Sub - squeeze(sino(:,CurrSubIndeces,:))); + objective(i) = 0.5*norm(residualSub(:)); % for the objective function output + end + % perform backprojection of a subset if (strcmp(proj_geom.type,'parallel') || strcmp(proj_geom.type,'fanflat') || strcmp(proj_geom.type,'fanflat_vec')) % if geometry is 2D use slice-by-slice projection-backprojection routine x_temp = zeros(size(X),'single'); @@ -579,7 +580,7 @@ else astra_mex_data3d('delete', id); end - X = X_t - (1/L_const).*x_temp; + X = X_t - (1/L_const).*x_temp; % ----------------Regularization part------------------------% if (lambdaFGP_TV > 0) -- cgit v1.2.3