Skip to content
Snippets Groups Projects
Commit e20fdd24 authored by Laura Nervo's avatar Laura Nervo Committed by Nicola Vigano
Browse files

Commenting, formatting, cleaning... 


Signed-off-by: default avatarLaura Nervo <laura.nervo@esrf.fr>

git-svn-id: https://svn.code.sf.net/p/dct/code/trunk@809 4c865b51-4357-4376-afb4-474e03ccb993
parent cad6232d
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
5_reconstruction/gtAstra3D.m
+ 43
− 30
View file @ e20fdd24
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
5_reconstruction/gtAstraAbsorptionVolume.m
+ 17
− 16
View file @ e20fdd24
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
......
5_reconstruction/gtAstraPrepareAbsorptionStack.m
+ 25
− 17
View file @ e20fdd24
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
5_reconstruction/gtAstra_FBP2D.m
+ 27
− 18
View file @ e20fdd24
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
5_reconstruction/gtSetupReconstruction.m
+ 14
− 19
View file @ e20fdd24
......@@ -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
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment