From 99cfddc7764a292c438868a02c9d512dabefbdc4 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Thu, 19 Oct 2017 12:47:02 +0100 Subject: moved code out of if-then-else fixed #5 --- main_func/FISTA_REC.m | 54 ++++++++++++++------------------------------------- 1 file changed, 15 insertions(+), 39 deletions(-) diff --git a/main_func/FISTA_REC.m b/main_func/FISTA_REC.m index bea1860..d177748 100644 --- a/main_func/FISTA_REC.m +++ b/main_func/FISTA_REC.m @@ -517,6 +517,18 @@ 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; @@ -528,18 +540,7 @@ else if (lambdaR_L1 > 0) % 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 - % 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 = zeros(Detectors, numProjSub, SlicesZ,'single'); % residual for a chosen subset for kkk = 1:numProjSub @@ -550,18 +551,7 @@ else 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 + % artifacts removal with Students t penalty residualSub = squeeze(weights(:,CurrSubIndeces,:)).*(sino_updt_Sub - squeeze(sino(:,CurrSubIndeces,:))); @@ -576,21 +566,7 @@ else objective(i) = ff; % for the objective function output else % 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 - [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 - residualSub = squeeze(weights(:,CurrSubIndeces,:)).*(sino_updt_Sub - squeeze(sino(:,CurrSubIndeces,:))); - objective(i) = 0.5*norm(residualSub(:)); % for the objective function output - end + 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 -- 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(-) 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