From 1fbcb871342f491068056742d1ddf29ce2a70d45 Mon Sep 17 00:00:00 2001
From: Nicola Vigano <nicola.vigano@esrf.fr>
Date: Mon, 9 Jan 2017 17:34:17 +0100
Subject: [PATCH] Improved 3DTV absorption volume reconstruction, with KL data
divergence term, and l12 norm in the TV term
Signed-off-by: Nicola Vigano <nicola.vigano@esrf.fr>
---
5_reconstruction/gtAstra3DTV.m | 39 ++++++++++++++-----
5_reconstruction/gtAstraAbsorptionVolume.m | 34 ++++++++++------
.../gtAstraPrepareAbsorptionStack.m | 2 +-
.../gtRecAbsorptionDefaultParameters.m | 24 ++++++++++++
5_reconstruction/gtRecDefaultParameters.m | 4 +-
5 files changed, 78 insertions(+), 25 deletions(-)
create mode 100644 5_reconstruction/gtRecAbsorptionDefaultParameters.m
diff --git a/5_reconstruction/gtAstra3DTV.m b/5_reconstruction/gtAstra3DTV.m
index 6d92b6ed..19916ddf 100644
--- a/5_reconstruction/gtAstra3DTV.m
+++ b/5_reconstruction/gtAstra3DTV.m
@@ -33,7 +33,12 @@ function volume = gtAstra3DTV(proj, varargin)
% Version 002 26-09-2012 by LNervo
% Formatting, cleaning, commenting
- conf = struct('verbose', false, 'is_cone', false, 'epsilon', 1e-4);
+ conf = struct(...
+ 'verbose', false, ...
+ 'is_cone', false, ...
+ 'detector_norm', 'l2', ...
+ 'lambda', 1e-2, ...
+ 'epsilon', 1e-4);
conf = parse_pv_pairs(conf, varargin);
@@ -120,8 +125,13 @@ function volume = gtAstra3DTV(proj, varargin)
stats.add_task('primal_update', 'Update of Primal')
stats.add_task('grad_comp', 'Computing the Gradient')
stats.add_task('div_comp', 'Computing the Divergence')
+
res_initial = gtMathsNorm_l2(proj.stack);
residuals = [1; zeros(proj.num_iter, 1)];
+ if (strcmpi(conf.detector_norm, 'kl'))
+ data_term = 4 .* sigma .* proj.stack;
+ end
+
fprintf('\b\b: ')
for ii = 1:proj.num_iter
char_num = fprintf('%02d/%02d (res: %f)', ii-1, proj.num_iter, residuals(ii));
@@ -129,14 +139,24 @@ function volume = gtAstra3DTV(proj, varargin)
stats.tic('fp')
tomo_fp_tmp = fp(volume_enh);
stats.toc('fp')
- residuals(ii+1) = sqrt(gtMathsNorm_l2(tomo_fp_tmp - proj.stack) / res_initial);
+ residuals(ii+1) = gtMathsNorm_l2(tomo_fp_tmp - proj.stack) / res_initial;
stats.tic('dual_det_update')
- if (conf.epsilon >= 0)
- temp_dual_det = dual_det + sigma .* (tomo_fp_tmp - proj.stack);
- temp_dual_denom = temp_dual_det + (temp_dual_det == 0);
- dual_det = temp_dual_det .* max(1 - sigma .* conf.epsilon ./ temp_dual_denom, 0);
- else
- dual_det = (dual_det + sigma .* (tomo_fp_tmp - proj.stack)) .* sigma_1;
+ switch (lower(conf.detector_norm))
+ case 'l2'
+ if (conf.epsilon >= 0)
+ temp_dual_det = dual_det + sigma .* (tomo_fp_tmp - proj.stack);
+ temp_dual_denom = temp_dual_det + (temp_dual_det == 0);
+ dual_det = temp_dual_det .* max(1 - sigma .* conf.epsilon ./ temp_dual_denom, 0);
+ else
+ dual_det = (dual_det + sigma .* (tomo_fp_tmp - proj.stack)) .* sigma_1;
+ end
+ case 'kl'
+ temp_dual = dual_det + sigma .* tomo_fp_tmp;
+ dual_det = (1 + temp_dual - sqrt((temp_dual - 1) .^ 2 + data_term)) ./ 2;
+ otherwise
+ error([mfilename ':wrong_argument'], ...
+ 'Detector divergence term: "%s" not recognized', ...
+ conf.detector_norm)
end
stats.toc('dual_det_update')
@@ -145,7 +165,8 @@ function volume = gtAstra3DTV(proj, varargin)
stats.toc('grad_comp')
stats.tic('dual_grad_update')
temp_dual_tv = dual_tv + d;
- dual_tv = temp_dual_tv ./ max(1, abs(temp_dual_tv));
+ norm_dual_tv = max(1, sqrt(sum(temp_dual_tv .^ 2, 4)));
+ dual_tv = temp_dual_tv ./ norm_dual_tv(:, :, :, [1 1 1]);
stats.toc('dual_grad_update')
stats.tic('bp')
diff --git a/5_reconstruction/gtAstraAbsorptionVolume.m b/5_reconstruction/gtAstraAbsorptionVolume.m
index 07b10e8c..d61bc633 100644
--- a/5_reconstruction/gtAstraAbsorptionVolume.m
+++ b/5_reconstruction/gtAstraAbsorptionVolume.m
@@ -1,4 +1,4 @@
-function abs = gtAstraAbsorptionVolume(parameters)
+function [abs, proj] = gtAstraAbsorptionVolume(parameters, proj)
% GTASTRAABSORPTIONVOLUME reconstruct absorption volume using 3D SIRT or 2D FBP
% abs = gtAstraAbsorptionVolume(parameters)
% -----------------------------------------
@@ -13,24 +13,29 @@ function abs = gtAstraAbsorptionVolume(parameters)
if ~exist('parameters','var') || isempty(parameters)
- load('parameters.mat');
+ parameters = gtLoadParameters();
end
projname = fullfile(parameters.acq.dir, 'absorption_stack.mat');
- if (exist(projname, 'file'))
- check = inputwdefault('Absorption stack already exists - do you want to (re)-create it? [y/n]', 'n');
+ if (~exist('proj', 'var') || isempty(proj))
+ if (exist(projname, 'file'))
+ check = inputwdefault('Absorption stack already exists - do you want to (re)-create it? [y/n]', 'n');
- if strcmpi(check, 'n')
- fprintf('Loading absorption stack..')
- c = tic();
- proj = load(projname);
- fprintf('\b\b: Done. (%f seconds)\n', toc(c))
+ if strcmpi(check, 'n')
+ fprintf('Loading absorption stack..')
+ c = tic();
+ proj = load(projname);
+ fprintf('\b\b: Done. (%f seconds)\n', toc(c))
+ else
+ proj = gtAstraPrepareAbsorptionStack(parameters);
+ end
else
+ fprintf('Absorption stack doesn''t exist, yet. Building it..')
+ c = tic();
proj = gtAstraPrepareAbsorptionStack(parameters);
+ fprintf('\b\b: Done in %g seconds.\n', toc(c));
end
- else
- proj = gtAstraPrepareAbsorptionStack(parameters);
end
if (isfield(parameters.rec, 'absorption'))
@@ -50,7 +55,12 @@ function abs = gtAstraAbsorptionVolume(parameters)
case {'3DART', 'SIRT'}
abs = gtAstra3D(proj, 'is_cone', is_cone);
case {'3DTV'}
- abs = gtAstra3DTV(proj, 'is_cone', is_cone);
+ if (~isfield(rec_abs, 'options') || isempty(rec_abs.options))
+ rec_abs_options = {};
+ else
+ rec_abs_options = [fieldnames(rec_abs.options), struct2cell(rec_abs.options)]';
+ end
+ abs = gtAstra3DTV(proj, 'is_cone', is_cone, rec_abs_options{:});
case '2DFBP'
abs = gtAstra_FBP2D(parameters, proj);
otherwise
diff --git a/5_reconstruction/gtAstraPrepareAbsorptionStack.m b/5_reconstruction/gtAstraPrepareAbsorptionStack.m
index be9a1003..102bbfa4 100644
--- a/5_reconstruction/gtAstraPrepareAbsorptionStack.m
+++ b/5_reconstruction/gtAstraPrepareAbsorptionStack.m
@@ -62,7 +62,7 @@ imgs_nums = 0 : rec_abs.interval : totproj-1;
Omega_rad = angle_rad(imgs_nums + 1);
Omega_deg = Omega_rad * 180 / pi;
-bbpos = repmat(acq.bb, numel(Omega_deg), []);
+bbpos = repmat(acq.bb, [numel(Omega_deg), 1]);
geom = gtGeoProjForReconstruction([], Omega_deg', [], bbpos, [], ...
detgeo, labgeo, samgeo, recgeo, 'ASTRA_absorption');
geom(:, 10:12) = round(geom(:, 10:12));
diff --git a/5_reconstruction/gtRecAbsorptionDefaultParameters.m b/5_reconstruction/gtRecAbsorptionDefaultParameters.m
new file mode 100644
index 00000000..3ab191b1
--- /dev/null
+++ b/5_reconstruction/gtRecAbsorptionDefaultParameters.m
@@ -0,0 +1,24 @@
+function par_rec = gtRecAbsorptionDefaultParameters(algo)
+% FUNCTION par_rec = gtRecAbsorptionDefaultParameters(algo)
+% Arguments can be (default in {}):
+% - algo : {'SIRT'} | '2DFBP' | '3DTV'
+%
+
+ switch (upper(algo))
+ case {'SIRT', '2DFBP'}
+ par_rec = struct(...
+ 'algorithm', upper(algo), 'num_iter', 30, ...
+ 'interval', 10, 'padding', 5, 'options', []);
+
+ case '3DTV'
+ par_3DTV_rec_opts = struct( ...
+ 'verbose', false, ...
+ 'detector_norm', 'l2', ...
+ 'lambda', 1e-2, ...
+ 'epsilon', 1e-4 ...
+ );
+ par_rec = struct(...
+ 'algorithm', upper(algo), 'num_iter', 30, 'interval', 10, ...
+ 'padding', 5, 'options', par_3DTV_rec_opts);
+ end
+end
diff --git a/5_reconstruction/gtRecDefaultParameters.m b/5_reconstruction/gtRecDefaultParameters.m
index 17eecbe6..433429fb 100644
--- a/5_reconstruction/gtRecDefaultParameters.m
+++ b/5_reconstruction/gtRecDefaultParameters.m
@@ -11,9 +11,7 @@ function par_rec = gtRecDefaultParameters(varargin)
conf = parse_pv_pairs(conf, varargin);
% Absorption volume options
- par_rec.absorption = struct( ...
- 'algorithm', conf.absorption_algo, 'num_iter', 30, ...
- 'interval', 10, 'padding', 5, 'options', [] );
+ par_rec.absorption = gtRecAbsorptionDefaultParameters(conf.absorption_algo);
% Grain volumes options
par_rec.grains = gtRecGrainsDefaultParameters(conf.grains_algo);
--
GitLab