From e20fdd2424c2494fd5a854551edcb57b5634015a Mon Sep 17 00:00:00 2001
From: Laura Nervo <lnervo@esrf.fr>
Date: Wed, 26 Sep 2012 11:53:46 +0000
Subject: [PATCH] Commenting, formatting, cleaning...
Signed-off-by: Laura Nervo <laura.nervo@esrf.fr>
git-svn-id: https://svn.code.sf.net/p/dct/code/trunk@809 4c865b51-4357-4376-afb4-474e03ccb993
---
5_reconstruction/gtAstra3D.m | 73 +++++++++++--------
5_reconstruction/gtAstraAbsorptionVolume.m | 33 +++++----
.../gtAstraPrepareAbsorptionStack.m | 42 ++++++-----
5_reconstruction/gtAstra_FBP2D.m | 45 +++++++-----
5_reconstruction/gtSetupReconstruction.m | 33 ++++-----
5 files changed, 126 insertions(+), 100 deletions(-)
diff --git a/5_reconstruction/gtAstra3D.m b/5_reconstruction/gtAstra3D.m
index 3df7d41b..c6b3d5a7 100644
--- a/5_reconstruction/gtAstra3D.m
+++ b/5_reconstruction/gtAstra3D.m
@@ -1,28 +1,30 @@
-function volume = gtAstra3D(~, ~ ,~, proj)
+function volume = gtAstra3D(~, ~, ~, proj)
% GTASTRA3D 3D SIRT reconstruction
-% function volume = gtAstra3D(~, ~ ,~, proj);
-%
+% volume = gtAstra3D(~, ~, ~, proj)
+% ---------------------------------
%
-
-% INPUTS: proj structure containing fields:
-% .stack stack of projections
-% .geom projection geometry for Astra
-% geom is a N x 12 matrix, witn N = #projections
-% Each row is:
-% ( rayX, rayY, rayZ, dX, dY, dZ, uX, uY, uZ, vX, vY, vZ )
-% ray: the ray direction
-% d : the center of the detector
-% u : the vector from detector pixel (0,0) to (0,1)
-% v : the vector from detector pixel (0,0) to (1,0)
-% .num_rows vertical size of projection images
-% .num_cols horizontal size of projection images
-% .num_iter number of iterations
-% .vol_size_x x size of output volume
-% .vol_size_y y size of output volume
-% .vol_size_z z size of output volume
+% INPUT: proj structure containing fields:
+% .stack stack of projections
+% .geom projection geometry for Astra
+% geom is a N x 12 matrix, witn N = #projections
+% Each row is:
+% ( rayX, rayY, rayZ, dX, dY, dZ, uX, uY, uZ, vX, vY, vZ )
+% ray: the ray direction
+% d : the center of the detector
+% u : the vector from detector pixel (0,0) to (0,1)
+% v : the vector from detector pixel (0,0) to (1,0)
+% .num_rows vertical size of projection images
+% .num_cols horizontal size of projection images
+% .num_iter number of iterations
+% .vol_size_x x size of output volume
+% .vol_size_y y size of output volume
+% .vol_size_z z size of output volume
%
-% OUTPUT: volume <double vol_size_x, vol_size_y, vol_size_z>)
-
+% OUTPUT: volume <double vol_size_x, vol_size_y, vol_size_z>)
+%
+%
+% Version 002 26-09-2012 by LNervo
+% Formatting, cleaning, commenting
OK = gtCheckGpu();
@@ -31,36 +33,47 @@ if ~OK
mexc = MException('ASTRA:no_GPU','gtAstra3D: Submit this job on a GPU machine ! Exiting...');
throw(mexc)
else
-
- %% projection geometry (type, #rows, #columns, vectors)
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ % projection geometry (type, #rows, #columns, vectors)
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
proj_geom = astra_create_proj_geom('parallel3d_vec', proj.num_rows, proj.num_cols, proj.geom);
- %% create volume and projection data objects (x, y, z)
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ % create volume and projection data objects (x, y, z)
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% proj_data should be a 3d "matrix" of size [#columns,#projections,#rows]
vol_geom = astra_create_vol_geom(proj.vol_size_x, proj.vol_size_y, proj.vol_size_z);
proj_id = astra_mex_data3d('create', '-proj3d', proj_geom, proj.stack);
- %% allocate and zero memory for reconstruction
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ % allocate and zero memory for reconstruction
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
vol_id = astra_mex_data3d('create', '-vol', vol_geom, 0);
- %% prepare reconstruction algorithm
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ % prepare reconstruction algorithm
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
cfg = astra_struct('SIRT3D_CUDA');
cfg.ProjectionDataId = proj_id;
cfg.ReconstructionDataId = vol_id;
cfg.option.MinConstraint = 0;
alg_id = astra_mex_algorithm('create', cfg);
- %% reconstruct
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ % reconstruct
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%disp('reconstructing...');
astra_mex_algorithm('iterate', alg_id, proj.num_iter);
volume = astra_mex_data3d('get', vol_id);
% vol is now a 3d "matrix" of size [x,y,z]
- %% clean up
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ % clean up
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
astra_mex_algorithm('delete', alg_id);
astra_mex_data3d('delete', vol_id);
astra_mex_data3d('delete', proj_id);
end
-end
\ No newline at end of file
+end % end of function
diff --git a/5_reconstruction/gtAstraAbsorptionVolume.m b/5_reconstruction/gtAstraAbsorptionVolume.m
index 666497fb..6ece5f68 100644
--- a/5_reconstruction/gtAstraAbsorptionVolume.m
+++ b/5_reconstruction/gtAstraAbsorptionVolume.m
@@ -1,24 +1,22 @@
function abs = gtAstraAbsorptionVolume(parameters)
-% GTASTRAABSORPTIONVOLUME reconstruct absorption volume using 3D SIRT or
-% 2D FBP
+% GTASTRAABSORPTIONVOLUME reconstruct absorption volume using 3D SIRT or 2D FBP
+% abs = gtAstraAbsorptionVolume(parameters)
+% -----------------------------------------
+% Checks if absorption_stack.mat exists (will create it otherwise)
+% Checks if absorption_mask.mat exists (will create it otherwise)
+% Outputs reconstructed volume into file absorption.mat
%
-% function abs = gtAstraAbsorptionVolume(parameters)
+% parameters.rec.method : '3DART' / '2DFBP'
%
-% Checks if absorption_stack.mat exists (will create it otherwise)
-% Checks if absorption_mask.mat exists (will create it otherwise)
-% Outputs reconstructed volume into file absorption.mat
-
-if ~exist('parameters','var')
- load parameters
-end
-
-
-
-acq = parameters.acq;
+% Version 002 26-09-2012 by LNervo
+% Formatting, cleaning, commenting
+if ~exist('parameters','var') || isempty(parameters)
+ load('parameters.mat');
+end
-projname = fullfile(acq.dir, 'absorption_stack.mat');
+projname = fullfile(parameters.acq.dir, 'absorption_stack.mat');
proj = [];
if exist(projname, 'file')
@@ -39,9 +37,12 @@ switch parameters.rec.method
abs = gtAstra3D([],[],[],proj);
case '2DFBP'
abs = gtAstra_FBP2D(parameters,proj);
+ otherwise
+ disp('Method for reconstruction not known. Please check parameters file')
+ return;
end
-end
+end % end of function
diff --git a/5_reconstruction/gtAstraPrepareAbsorptionStack.m b/5_reconstruction/gtAstraPrepareAbsorptionStack.m
index 830e01b9..ce17e4a2 100644
--- a/5_reconstruction/gtAstraPrepareAbsorptionStack.m
+++ b/5_reconstruction/gtAstraPrepareAbsorptionStack.m
@@ -1,15 +1,21 @@
function proj = gtAstraPrepareAbsorptionStack(parameters)
-% GTASTRAPREPAREABSORPTIONSTACK assembles the 3D projection stack of
-% absorption images and saves them to file absorption_stack.mat
+% GTASTRAPREPAREABSORPTIONSTACK Assembles the 3D projection stack of absorption
+% images and saves them to file absorption_stack.mat
%
-% INPUTS: parameters <structure>
+% proj = gtAstraPrepareAbsorptionStack(parameters)
+% ------------------------------------------------
+% INPUT:
+% parameters <structure>
%
-% OUTPUTS: proj <structure>
-% fields:
-% - stack <double X,X,Z> 3D stack of absorption images
-% - geom <double Z,12 > Vector describing Astra preojection geometry
-% - Omega <Double Z> Rotation angles in degrees
-
+% OUTPUT:
+% proj <structure> with fields:
+% - stack <double X,X,Z> 3D stack of absorption images
+% - geom <double Z,12 > Vector describing Astra preojection geometry
+% - Omega <Double Z> Rotation angles in degrees
+%
+%
+% Version 002 26-09-2012 by LNervo
+% Formatting, cleaning, commenting
acq = parameters.acq;
rec = parameters.rec;
@@ -18,13 +24,13 @@ geo = parameters.labgeo;
if strcmp(acq.type, '360degree')
totproj = 2 * acq.nproj;
increment = 2 * pi / totproj;
- angle_rad = [0 : increment : 2*pi-increment];
+ angle_rad = 0 : increment : 2*pi-increment;
elseif strcmp(acq.type, '180degree')
totproj = acq.nproj;
increment = pi / totproj;
- angle_rad = [pi/2 : increment: 3/2*pi-increment];
+ angle_rad = pi/2 : increment: 3/2*pi-increment;
else
- error('scan type not recognised!')
+ error('Scan type not recognised!')
end
% projection geometry... Here we assume that the absorption images are
@@ -32,7 +38,7 @@ end
% The origin of the sample coordinate system is on the rotation axis - in
% the middle of the illuminated area
-detpos0 = [0, 0, 0];
+detpos0 = [0, 0, 0];
rotcomp = gtMathsRotationMatrixComp(geo.rotdir,'row');
@@ -43,7 +49,7 @@ geom = zeros(totproj/rec.absinterval, 12);
jj = 1;
-disp('building sinogram...');
+fprintf('Building sinogram...');
for ii = 0 : rec.absinterval : totproj-1
@@ -54,7 +60,7 @@ for ii = 0 : rec.absinterval : totproj-1
end
Omega(jj) = angle_rad(ii + 1);
- %%
+
stack(:, jj, :) = -log(edf_read(absname))';
Rottens = gtMathsRotationTensor(-Omega(jj)*180/pi, rotcomp);
detdiru = geo.detdiru * Rottens;
@@ -65,7 +71,9 @@ for ii = 0 : rec.absinterval : totproj-1
geom(jj, :) = [r(1), r(2), r(3), detpos(1), detpos(2), detpos(3), detdiru(1), detdiru(2), detdiru(3), detdirv(1), detdirv(2), detdirv(3)];
jj = jj + 1;
+
end
+fprintf('...done!\n')
proj.stack = stack;
proj.Omega = Omega;
@@ -77,8 +85,8 @@ proj.vol_size_x = acq.bb(3);
proj.vol_size_y = acq.bb(3);
proj.vol_size_z = acq.bb(4);
-disp('writing sinogram');
name = fullfile(acq.dir, 'absorption_stack.mat');
+disp(['Writing sinogram in ' name]);
save(name, '-struct', 'proj');
-end
\ No newline at end of file
+end % end of function
diff --git a/5_reconstruction/gtAstra_FBP2D.m b/5_reconstruction/gtAstra_FBP2D.m
index e6a09d51..cf52c380 100644
--- a/5_reconstruction/gtAstra_FBP2D.m
+++ b/5_reconstruction/gtAstra_FBP2D.m
@@ -1,28 +1,32 @@
-function abs=gtAstra_FBP2D(parameters, proj)
-
-acq=parameters.acq;
-rec=parameters.rec;
-geo=parameters.labgeo;
-
-if ~exist('proj', 'var')
- name=sprintf('%s/absorption_stack.mat',acq.dir);
+function abs = gtAstra_FBP2D(parameters, proj)
+% GTASTRA_FBP2D
+% abs = gtAstra_FBP2D(parameters, proj)
+% -------------------------------------
+%
+%
+% Version 002 26-09-2012 by LNervo
+% Formatting, cleaning, commenting
+
+
+if ~exist('proj','var')
+ name = fullfile(parameters.acq.dir,'absorption_stack.mat');
if exist(name,'file')
check = inputwdefault('Absorption stack already exists - do you want to (re)-create it ? [y/n]', 'n');
if strcmpi(check,'n')
disp('reading projections...');
- proj=load(name);
+ proj = load(name);
else
- proj=gtAstraPrepareAbsorptionStack(parameters);
+ proj = gtAstraPrepareAbsorptionStack(parameters);
end
else
-
- proj=gtAstraPrepareAbsorptionStack(parameters);
+ proj = gtAstraPrepareAbsorptionStack(parameters);
end
end
-
-%% now start FBP reconstruction of absorption volume
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% now start FBP reconstruction of absorption volume
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%proj_data=permute(proj.stack,[3 2 1]); % transforms stack of projection images to stack of sinograms
proj_data = proj.stack;
@@ -42,8 +46,9 @@ yvol=hsize;
zvol=vsize;
-
-%% projection geometry (type, #rows, #columns, vectors)
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% projection geometry (type, #rows, #columns, vectors)
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
proj_geom = astra_create_proj_geom('parallel', 1.0, hsize_zp, proj.Omega);
vol_geom = astra_create_vol_geom(hsize, hsize);
proj_id = astra_create_projector('linear', proj_geom, vol_geom);
@@ -52,7 +57,9 @@ sinogram_id = astra_mex_data2d('create','-sino', proj_geom, 0);
% allocate space to store reconstruction image
recon_id = astra_mex_data2d('create', '-vol', vol_geom, 0);
-% create FBP algorithm configuration
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% create FBP algorithm configuration
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
cfg = astra_struct('FBP_CUDA');
cfg.ProjectionDataId = sinogram_id;
cfg.ReconstructionDataId = recon_id;
@@ -83,11 +90,13 @@ for i=1:vsize
astra_mex_algorithm('delete', alg_id);
end
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% clean up
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
astra_mex_data2d('delete', sinogram_id);
astra_mex_data2d('delete', recon_id);
astra_mex_projector('delete', proj_id);
-
+end % end of function
diff --git a/5_reconstruction/gtSetupReconstruction.m b/5_reconstruction/gtSetupReconstruction.m
index 90fb7659..dda3c3fb 100644
--- a/5_reconstruction/gtSetupReconstruction.m
+++ b/5_reconstruction/gtSetupReconstruction.m
@@ -6,11 +6,12 @@ function gtSetupReconstruction(phaseID)
% ------------------------------
%
% INPUT:
-% phaseID = number of the phase <double> {default: 1}
+% phaseID = number of the phase <double> {1}
%
% Version 003 29-05-2012 by LNervo
-% Add comments in header about functions used;
-% Add possibility to recompile 'gtChangeThreshold3D
+% Add comments in header about functions used
+% Add possibility to recompile 'gtChangeThreshold3D'
+%
% Version 002 25-05-2012 by LNervo
% Commenting and formatting: fprintf -> disp; y/n -> [y/n];
%
@@ -29,17 +30,17 @@ if ~GPU_OK
disp('from rnice: oarsub -I -l "walltime=hh:mm:ss" -p "gpu = ''YES''"')
return;
end
+
+if ~exist('phaseID','var') || isempty(phaseID)
+ phaseID = 1;
+end
-check=inputwdefault('Do you want to recompile the function ''gtChangeThreshold3D''? [y/n]','n');
+check = inputwdefault('Do you want to recompile the function ''gtChangeThreshold3D''? [y/n]','n');
if strcmpi(check,'y')
gt_mcc('gtChangeThreshold3D');
end
-
-if ~exist('phaseID','var')
- phaseID = 1;
-end
-parameters=[];
+parameters = [];
load('parameters.mat');
parameters_savefile = fullfile(parameters.acq.dir, 'parameters.mat');
@@ -48,13 +49,12 @@ nof_phases = length(parameters.cryst);
if nof_phases > 1
phaseID = inputwdefaultnumeric('Enter phaseID', num2str(1));
end
-
phaseID_str = sprintf('%02d',phaseID);
-disp(['Setting up Reconstruction for phase ' num2str(phaseID) ' (' parameters.cryst(phaseID).name ')']);
+disp(['Setting up Reconstruction for phase ' phaseID_Str ' (' parameters.cryst(phaseID).name ')']);
-filetable = sprintf('%s_filetable', parameters.acq.name);
+filetable = [parameters.acq.name '_filetable'];
samplefile = fullfile(parameters.acq.dir, '4_grains', 'sample.mat');
if exist(samplefile, 'file')
sample = GtSample.loadFromFile(samplefile);
@@ -81,16 +81,11 @@ absorption_mask = fullfile(parameters.acq.dir, '5_reconstruction', 'volume_mask.
if ~exist(absorption_mask,'file')
check = inputwdefault('Can not find absorption mask - would you like to segment absorption volume ? [y/n]', 'y');
if strcmpi(check, 'y')
- disp('Please create file volume_mask.mat from absorption volume (right click in the volume to select a seed point)');
+ disp('Please create file volume_mask.mat from absorption volume (right click in the volume to select a seed point and then press ''Segment'')');
GtGuiThresholdVolume(abs);
end
end
-
-
-
-
-
% get reconstruction parameters
check = inputwdefault('Use default parameters for grain reconstruction ? [y/n]', 'y');
if ~strcmpi(check, 'y')
@@ -126,7 +121,7 @@ if strcmpi(check, 'y')
disp(OAR_parameters)
disp(['Using ' num2str(oar.njobs) ' OAR jobs']);
pause(0.5);
- OAR_make('gtChangeThreshold3D',first, last, oar.njobs, [parameters.acq.dir, sprintf(' %d', phaseID)], true, 'walltime', OAR_parameters.walltime);
+ OAR_make('gtChangeThreshold3D',first, last, oar.njobs, [parameters.acq.dir ' ' phaseID_str], true, 'walltime', OAR_parameters.walltime);
else
gtChangeThreshold3D(first, last, parameters.acq.dir, phaseID);
end
--
GitLab