ultima modifica
07/08/2013 15:44
function [x, iter, err, discr, times, varargout] = sgp_deblurring_boundary(psf, gn, varargin)
% sgp_deblurring - SGP algorithm for non-regularized deblurring with boundary effects correction
% This function solves an image deblurring problem by applying the SGP
% algorithm to the minimization of the generalized Kullback-Leibler
% divergence with no regularization [1], with boundary effects correction [2]:
%
% min KL(A*x + bg, gn)
% x in OMEGA
%
% where KL(u,v) is the generalized Kullback-Leibler divergence between
% vectors u and v, bg is the background, gn are the observed data and
% the feasible set OMEGA is either
% 'nonneg': x(i) >= 0;
% 'nonneg_eq': same as 'nonneg' plus a liner constraint of total flux
% conservation.
%
% Note: the PSF must be normalized. All columns of A must sum-up to 1.
%
% [1] S. Bonettini, R. Zanella, L. Zanni,
% "A scaled gradient projection method for constrained image deblurring",
% Inverse Problems 25(1), 2009, January, 015002.
% [2] M. Prato, R. Cavicchioli, L. Zanni, P. Boccacci, M. Bertero,
% "Efficient deconvolution methods for astronomical imaging: algorithms
% and IDL-GPU codes",
% Astronomy and Astrophysics 539, 2012, January, A133, 1-11.
%
% SYNOPSIS
% [x, iter, err, discr, time] = sgp_deblurring_boundary(psf, gn[, opts])
%
% MANDATORY INPUT
% psf (double array) - the Point Spread Function (PSF) of the blurring
% operator A. The routine internally defines how to
% compute A*x and transpose(A)*x.
% N.B. The psf array must sum-up to 1.
% gn (double array) - measured image
%
% OPTIONAL INPUT
% The following options must be provided as keyword/value pairs.
% 'OBJ' - Exact solution, for error calculation (double array)
% 'BG' - Background value (double)
% DEFAULT = 0
% 'INITIALIZATION' - Choice for starting point:
% 0 - all zero starting point
% 1 - random starting point
% 2 - initialization with gn
% 3 - initialization with
% ones(size(gn))*sum(gn(:) - bg) / numel(gn)
% x0 - user-provided starting point (double array)
% DEFAULT = 0
% 'MAXIT' - Maximum number of iterations (integer)
% DEFAULT = 1000
% 'VERBOSE' - Verbosity level (integer)
% 0 - silent
% 1 - print configuration parameters at startup
% 2 - print configuration parameters at startup and
% some information at each iteration
% DEFAULT = 0
% 'STOPCRITERION' - Choice for stopping rule (integer)
% 1 -> iter > MAXIT
% 2 -> ||x_k - x_(k-1)|| <= tol*||x_k|| OR iter > MAXIT
% 3 -> |KL_k - KL_(k-1)| <= tol*|KL_k| OR iter > MAXIT
% 4 -> (2/N)*KL_k <= tol OR iter > MAXIT
% DEFAULT = 1;
% 'TOL' - Tolerance used in the stopping criterion
% DEFAULT = 1e-4 if STOPCRITERION = 2 or 3
% DEFAULT = 1+1/mean(gn) if STOPCRITERION = 4
% 'TOLPOS' - Threshold against zero to detect non-positive values
% DEFAULT = 1e-12;
% 'M' - Nonmonotone lineasearch memory (positive integer)
% If M = 1 the algorithm is monotone
% DEFAULT = 1 (monotone)
% 'GAMMA' - Linesearch sufficient-decrease parameter (double)
% DEFAULT = 1e-4
% 'BETA' - Linesearch backtracking parameter (double)
% DEFAULT = 0.4
% 'ALPHA_MIN' - Lower bound for Barzilai-Borwein' steplength (double>0)
% DEFAULT = 1e-5
% 'ALPHA_MAX' - upper bound for Barzilai-Borwein' steplength (double>0)
% DEFAULT = 1e5
% 'MALPHA' - Memory length for alphaBB2 (positive integer)
% DEFAULT = 3
% 'TAUALPHA' - Alternating parameter for Barzilai-Borwein' steplength
% (double)
% DEFAULT = 0.5
% 'INITALPHA' - Initial value for Barzilai-Borwein' steplength (double)
% DEFAULT = 1.3
% 'BOUNDARYSIZE' - Enlargement parameter for boundary effects correction
% (nonnegative integer/vector)
% If it is a nonnegative integer, then the correction
% is the same in all the dimensions. If the parameter is
% a vector, it must have as many elements as the object's
% dimensions: in this case different correction can be
% done in different directions.
% For each coordinate direction, the nonnegative value
% represents of the corresponding PSF size that is used
% to enlarge the reconstruzion domanin.
% A value of 0 (zero) in a vector means NO boundary
% effect correction will be taken in that coordinate
% direction. Giving the scalar zero only completely
% disable boundary effects correction.
% EXPERIMENTAL, NOT IMPLEMENTED YET:
% If BOUNDARYSIZE = 'PowerOfTwo', then the
% domain will be enlarged in each direction up to the
% first power of 2 containing the sum of the object's and
% the PSF's domains.
% DEFAULT = 1
% 'OMEGA' - Constraints type (string):
% 'nonneg' : non negativity
% 'nonneg_eq': non negativity and flux conservation
% ( estimated flux is sum(gn(:) - bg) )
% DEFAULT: 'nonneg'
% 'SAVE' - Directory name + optionally problem name (string) for
% saving iterations data as controlled by the 'SAVEFREQ'
% option. If the option is not given, then the saving of
% intermediate results is disabled. If its value is the
% empty string, then the current directory is selected.
% Saved data include:
% - the reconstructed image x_k
% - the Puetter residual: (x_k - gn) ./ sqrt(x_k)
% DEFAULT: '.' (current directory)
% 'SAVEFREQ' - Saving frequency (integer or integer array): if it is a
% scalar, then the reconstructed image x_k and the Putter
% residual (x_k - gn)./ sqrt(x_k) are saved in the
% directory named in 'SAVE' option. If it is a positive
% integer vector of iteration counts, then the image
% and the residual are saved at those iterations
% only.
% If the value/s is/are not integer, then the integer
% part is taken. If they are not positive or the argument
% is present but empty, then a warning appears and the
% saving is disabled.
% DEFAULT = 1 (save at each iteration)
% 'RESTART' - Data filename (string) to continue a previously
% suspended computation. If the option is not given, then
% the deblurring process start from scratch, otherwise
% the state saved in the named file is restored and the
% computation continues from that point. If the option is
% given, but the value is the empty string, then the
% default filename is used to load the computation state.
% DEFAULT: 'SGPstate'
%
% WARNING: since this is an optional input, the user MUST
% provide the first two mandatory function inputs as
% empty matrices. Usually, the 'RESTART' option is given
% as THE ONLY optional input. However, the user have the
% opportunity to also give a subset of the other optional
% inputs, so that to overright the value of some of the
% parameters during computation. In this case, the user
% MUST BE VERY CAREFUL and know exactly what he/she is
% doing, because the resulting algorithm could no longer
% be convergent. Overrightable parameters are:
% 'MAXIT', 'M', 'GAMMA', 'BETA', 'ALPHA_MIN', 'ALPHA_MAX',
% 'MALPHA', 'TAUALPHA', 'VERBOSE', 'TOL', 'SAVE',
% 'SAVEFREQ', 'SUSPEND', 'SUSPENDFILE'
% 'SUSPEND' - Iteration count for computation suspension (positive
% integer). If it is given, then the deblurring process
% is suspended at the given iteration and the process
% state is saved to a MAT file as a whole for, later
% reloading. The value of this option can only be a
% non-negative integer and cannot be empty. A zero value
% disable the suspension even if the option in given.
% The state filename can optionally be given as the value
% of the 'SUSPENDFILE' option.
% WARNING: the iterative deblurring process is stopped
% at the given iteration, but the function is not halted.
% The current values of the output parameters are
% returned to the caller. If verb > 0 then the user is
% alerted of occurence of the suspension by a command
% window message.
% DEFAULT: 0
% 'SUSPENDFILE' - MAT-filename for saving process state (string). It can
% optionally include directory name + problem name. It is
% meaningful only if suspension is enabled, otherwise it
% remains unused. The current state of the deblurring
% process, as it is at the end of the iteration selected
% in the 'SUSPEND' option, is saved as a whole to the
% file for later restart.
% WARNING: suspension can only occur at the end of an
% iteration, including the first one. At the moment the
% whole workspace is saved, so the file can become very
% large.
% DEFAULT: 'SGPstate' (same as of the 'RESTART' option)
%
% OUTPUT
% x - Reconstructed data
% iter - Total number of iterations executed
% err - Error value at each iteration. If OBJ was not given,
% then err is the empty matrix.
% discr - Discrepancy value after each iteration:
% D = 2/numel(x_k) * KL( Ax_k + bg, gn)
% time - CPU time after each iteration
%
% ------------------------------------------------------------------------------
%
% This software is developed within the research project
%
% PRISMA - Optimization methods and software for inverse problems
% http://www.unife.it/prisma
%
% funded by the Italian Ministry for University and Research (MIUR), under
% the PRIN2008 initiative, grant n. 2008T5KA4L, 2010-2012. This software is
% part of the package "IRMA - Image Reconstruction in Microscopy and Astronomy"
% currently under development within the PRISMA project.
%
% Version: 1.0.3
% Date: May. 2012
% Authors:
% Riccardo Zanella, Gaetano Zanghirati
% Dept. of Mathematics, University of Ferrara, Italy
% riccardo.zanella@unife.it, g.zanghirati@unife.it
% Roberto Cavicchioli, Luca Zanni
% Dept. of Pure Appl. Math., Univ. of Modena and Reggio Emilia, Italy
% roberto.cavicchioli@unimore.it, luca.zanni@unimore.it
%
% Software homepage: http://www.unife.it/prisma/software
%
% Copyright (C) 2011 by R. Cavicchioli, R. Zanella, G. Zanghirati, L. Zanni.
% ------------------------------------------------------------------------------
% COPYRIGHT NOTIFICATION
%
% Permission to copy and modify this software and its documentation for
% internal research use is granted, provided that this notice is retained
% thereon and on all copies or modifications. The authors and their
% respective Universities makes no representations as to the suitability
% and operability of this software for any purpose. It is provided "as is"
% without express or implied warranty. Use of this software for commercial
% purposes is expressly prohibited without contacting the authors.
%
% This program is free software; you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by the
% Free Software Foundation; either version 3 of the License, or (at your
% option) any later version.
%
% This program is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License along
% with this program; if not, either visite http://www.gnu.org/licenses/
% or write to
% Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
% ==============================================================================
% start the clock
%t0 = cputime; % <<< MOVED RIGHT BEFORE THE "ARRAY SIZE FOR BOUNDARY CORRECTION" SECTION!!!
% Enable full verbose flag to see the statistics table, useful when tracking the behavior of
% the algorithm among different architectures, leave it false for time-aware reconstructions
fullVerboseFlag = false;
% use gradient = AT( M_s - gn ./( A(x) + bg ) )
% instead of gradient = ATMimg - AT( gn ./( A(x) + bg ) )
% no need for ATMimage, it will be cleared just after creating R_MASK
modifiedGradient = true;
% Modified weights:
% Weights_i = 1.0 if i \notin R
% Weights_i = 1.0 / (AT(M_s)) if i \in R
% You can then evaluate D^{-1} by "inverting" D
% Useless, except for tracking, leave it false
modifiedWeights = false;
restart = false; % flag for restarting: false -> start from scratch
suspendIter = 0; % iteration count for process suspension: 0 -> disable suspension
restartFile = 'SGPstate'; % default filename to load SGP state data for restart
suspendFile = 'SGPstate'; % default filename to save SGP state data for later restart
if (nargout > 5), RIflag = 1; else, RIflag = 0; end % Return the "Improvement Factor" information
% The 'RESTART' option is treated differently from the others
lengthVarArgIn = length(varargin);
lastOption = lengthVarArgIn;
i = 1;
while( i <= lastOption )
switch upper(varargin{i})
case 'RESTART'
restart = true;
if ( ~isempty(varargin{i+1}) )
if ( ischar(varargin{i+1}) )
restartFile = varargin{i+1};
else
error('Wrong data type for RESTART option: it must be a string');
end
end
varargin(i:i+1) = []; % removing this options from input list
lastOption = lastOption - 2;
otherwise
i = i + 2;
end
end
if ( restart )
% Continue a previously interrupted computation
if (strcmp(restartFile(end-4:end),'.mat')), ext = ''; else, ext = '.mat'; end
fprintf('\n >> Restarting computation from file %s%s',restartFile,ext);
varargin_new = varargin; % TO DO: remove varargin from saved state!!!
% >> TO DO: ADJUST THE RESTARTING OF THE COMPUTING TIME!!
% >> NOTE: check for t0 right before the "Array size for boundary correction" section
load( restartFile );
varargin = varargin_new;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Read the optional parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if (rem(length(varargin),2)==1)
error('Optional parameters should always go by pairs');
else
for i=1:2:(length(varargin)-1)
switch upper(varargin{i})
case 'MAXIT'
MAXIT = varargin{i+1};
case 'M'
M = varargin{i+1};
case 'GAMMA'
gamma = varargin{i+1};
case 'BETA'
beta = varargin{i+1};
case 'ALPHA_MIN'
alpha_min = varargin{i+1};
case 'ALPHA_MAX'
alpha_max = varargin{i+1};
case 'MALPHA'
Malpha = varargin{i+1};
case 'TAUALPHA'
tau = varargin{i+1};
case 'VERBOSE'
verb = varargin{i+1};
case 'TOL'
tol = varargin{i+1};
case 'TOLPOS'
tolpos = varargin{i+1};
case 'SAVE'
mysave = true;
if (~isempty(varargin{i+1}))
if (ischar(varargin{i+1}))
dest_dir = varargin{i+1};
else
error('Wrong data type for directory name of intermediate results');
end
else
dest_dir = '.';
end
% possible need of output dir handling
[success, mess] = mkdir(dest_dir);
if not(success)
error('%s: %s',dest_dir,mess);
end
if (verb > 0)
if (~isempty(mess))
fprintf('\n%s\n',mess);
end
end
case 'SAVEFREQ'
if (~isempty(varargin{i+1}))
savefreq = fix(varargin{i+1}(:));
if ( any(savefreq <= 0) )
warning('Non-positive saving frequency: saving disabled.');
mysave = false;
end
else
warning('Empty saving frequency: saving disabled.');
mysave = false;
end
case 'SUSPEND'
if ( isempty(varargin{i+1}) )
error('SUSPEND option requires an iteration number');
elseif (~isnumeric(varargin{i+1}) || varargin{i+1}~=abs(fix(varargin{i+1})))
error('Wrong data type for SUSPEND option: it must be a positive integer');
else
suspendIter = varargin{i+1};
end
case 'SUSPENDFILE'
if ( ~isempty(varargin{i+1}) )
if ( ischar(varargin{i+1}) )
suspendFile = varargin{i+1};
else
error('Wrong data type for SUSPENDFILE option: it must be a string');
end
end
otherwise
error(['Unrecognized or not changeable (after restart) option: ''' varargin{i} '''']);
end;
end;
end
else
% test for number of required parametres
if (nargin-lengthVarArgIn) ~= 2
error('Wrong number of required parameters');
end
%%%%%%%%%%%%%%%%%%%%%%%%
% SGP default parameters
%%%%%%%%%%%%%%%%%%%%%%%%
% SGPpar = par_init();
MAXIT = 1000; % maximum number of iterations
gamma = 1e-4; % for sufficient decrease
beta = 0.4; % backtracking parameter
M = 1; % memory in obj. function value (if M = 1 monotone)
alpha_min = 1e-5; % alpha lower bound
alpha_max = 1e5; % alpha upper bound
Malpha = 3; % alfaBB1 memory
tau = 0.5; % alternating parameter
initalpha = 1.3; % initial alpha
initflag = 0; % 0 -> initial x all zeros
errflag = false; % 0 -> no error calculation
err = []; % errors w.r.t. the object, if it is given
verb = 0; % verbosity: 0 -> silent
stopcrit = 1; % stopping criterion: 1 -> number of iterations
boundarycorr = 1; % flag for boundary effects correction
boundarysize = 1; % size parameter for boundary effects correction
bg = 0; % background value
omega = 'nonneg'; % non negativity constraints
AT = []; % function handle for A'*x
dest_dir = '.'; % where to store intermediate and final results
mysave = false; % true = save, false = don't save intermediate resutls
savefreq = 1; % iteration interval for storing intermediate results
savedLastIter = 0; % flag: SET BY THE PRORAM
sigma_tol = 1.0e-2; % zero-threshold for small elements in the input arrays
tolpos = 1.0e-12; % threshold against zero to detect non-positive values
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Read the optional parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if (rem(length(varargin),2)==1)
error('Optional parameters should always go by pairs');
else
for i=1:2:(length(varargin)-1)
switch upper(varargin{i})
% case 'AT'
% AT = varargin{i+1};
% warning('This version only works with the PSF (point spread function) matrix: it does not need AT.');
case 'OBJ'
obj = varargin{i+1};
errflag = true;
case 'INITIALIZATION'
if numel(varargin{i+1}) > 1 % initial x provided by user
initflag = 999;
x = varargin{i+1};
else
initflag = varargin{i+1};
end
case 'MAXIT'
MAXIT = varargin{i+1};
case 'M'
M = varargin{i+1};
case 'GAMMA'
gamma = varargin{i+1};
case 'BETA'
beta = varargin{i+1};
case 'ALPHA_MIN'
alpha_min = varargin{i+1};
case 'ALPHA_MAX'
alpha_max = varargin{i+1};
case 'MALPHA'
Malpha = varargin{i+1};
case 'TAUALPHA'
tau = varargin{i+1};
case 'INITALPHA'
initalpha = varargin{i+1};
case 'BOUNDARYSIZE'
boundarysize = varargin{i+1};
case 'VERBOSE'
verb = varargin{i+1};
case 'TOL'
tol = varargin{i+1};
case 'TOLPOS'
tolpos = varargin{i+1};
case 'STOPCRITERION'
stopcrit = varargin{i+1};
case 'BG'
bg = varargin{i+1};
case 'OMEGA'
omega = varargin{i+1};
case 'SAVE'
mysave = true;
if (~isempty(varargin{i+1}))
if (ischar(varargin{i+1}))
dest_dir = varargin{i+1};
else
error('Wrong directory/file name for intermediate results');
end
end
case 'SAVEFREQ'
if (~isempty(varargin{i+1}))
savefreq = fix( varargin{i+1}(:) );
if ( any(savefreq <= 0) )
warning('Non-positive saving frequency: saving disabled.');
mysave = false;
end
else
warning('Empty saving frequency: saving disabled.');
mysave = false;
end
case 'SUSPEND'
if ( isempty(varargin{i+1}) )
error('SUSPEND option requires an iteration number');
elseif (~isnumeric(varargin{i+1}) || varargin{i+1}~=abs(fix(varargin{i+1})))
error('Wrong data type for SUSPEND option: it must be a positive integer');
else
suspendIter = varargin{i+1};
end
case 'SUSPENDFILE'
if ( ~isempty(varargin{i+1}) )
if ( ischar(varargin{i+1}) )
suspendFile = varargin{i+1};
else
error('Wrong data type for SUSPENDFILE option: it must be a string');
end
end
otherwise
error(['Unrecognized option: ''' varargin{i} '''']);
end;
end;
end
%%%%%%%%%%%%%%%%%%
% function handles
%%%%%%%%%%%%%%%%%%
if ( isa(psf,'function_handle') )
error('This version only works with the PSF (point spread function) matrix: function handle is not allowed.');
% if ( isempty(AT) )
% error('Missing parameter: AT');
% end
% if( ~isa(AT,'function_handle') )
% error('AT is not a function handle');
% end
else
% Check column-normalization condition on A
% sumColsA = sum(A)';
% tolCheckA = 1.0e4*eps;
% checkA = find(abs(sumColsA-1) > tolCheckA);
% if (~isempty(checkA))
% errmsg = sprintf('\n\t%d %s\n\t%s %d:\n\t%s%d%s%e, %s%e',...
% length(checkA),...
% 'not-normalized columns found in blurring matrix A.',...
% 'The first one is column',checkA(1),...
% '|sum(A(:,',checkA(1),')) - 1| = ',...
% abs(sumColsA(checkA(1))-1),'tolerance = ',tolCheckA);
% error('Not-normalized blurring matrix A: %s%s',...
% 'provide a normalized A (see documentation).',errmsg);
% end
%AT = @(x) A'*x;
%A = @(x) A*x;
% Check normalization of the blurring operator A (PSF)
sumPSF = sum(psf(:));
checkPSF = abs(sumPSF-1); tolCheckPSF = 1.0e4*eps;
if (checkPSF > tolCheckPSF)
errmsg = sprintf('\n\t|sum(psf(:)) - 1| = %e, tolerance = %e\n',checkPSF,tolCheckPSF);
error('Not-normalized PSF: provide a normalized PSF (see documentation).%s',errmsg);
end
end
%%%%%%%%%%%%%%%%
% starting point >>> ATTENZIONE: sovrascritto piu' avanti!!! (vedi IDL)
%%%%%%%%%%%%%%%%
switch initflag
case 0 % all zeros
x = zeros(size(gn));
case 1 % random
x = randn(size(gn));
case 2 % gn
x = gn;
case 3 % same flux as gn - bg
x = sum(gn(:) - bg)/numel(gn)*ones(size(gn));
case 999 % x is explicitly given, check dimension
if not( size(x) == size(gn) )
error('Invalid size of the initial point.');
end
otherwise
error('Unknown initialization option.');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% every image is treated as a vector
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
obj_size = size(x);
img_size = size(gn);
psf_size = size(psf);
x = x(:);
gn = gn(:);
%keyboard
% start the clock
t0 = cputime;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% determine array size for boundary correction
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
max_size = max([obj_size; psf_size; img_size]);
% work_size = 2.^(floor(log2(max_size))+1);
work_size = img_size + psf_size;
obj_margin = floor( (work_size - obj_size)/2 ) + 1; % matlab's array indices start from 1
psf_margin = floor( (work_size - psf_size)/2 ) + 1; % matlab's array indices start from 1
img_margin = floor( (work_size - img_size)/2 ) + 1; % matlab's array indices start from 1
% WARNING: at this level, if A and gn differ in dimensions, then
% the input A MUST be sized work_size by the user
% if ( ~isa(A,'function_handle') )
psf_work = zeros(work_size);
psf_upperind = psf_margin + psf_size - 1;
psf_ranges = [psf_margin; psf_upperind]';
psf_work_ind = findIndices( psf_ranges, work_size );
psf_work( psf_work_ind ) = psf(:);
% if ( isequal(psf_size,obj_size) )
% obj_upperind = psf_upperind;
% obj_ranges = psf_ranges;
% obj_work_ind = psf_work_ind;
% else
% obj_upperind = obj_margin + obj_size - 1;
% obj_ranges = [obj_margin; obj_upperind]';
% obj_work_ind = findIndices( obj_ranges, work_size );
% end
% if ( isequal(img_size,obj_size) )
% img_upperind = obj_upperind;
% img_ranges = obj_ranges;
% img_work_ind = obj_work_ind;
% else
img_upperind = img_margin + img_size - 1;
img_ranges = [img_margin; img_upperind]';
img_work_ind = findIndices( img_ranges, work_size );
% end
% else
% warning('The blurring operator A is a function handle: you must ensure compliant size of the resulting vector');
% end
TF = fftn(fftshift(psf_work));
CTF = conj(TF);
A = @(x) Afunction(x,TF,work_size);
AT = @(x) Afunction(x,CTF,work_size);
%%%%%%%%%%%%%%%%
% stop criterion
%%%%%%%%%%%%%%%%
if not( ismember(stopcrit, [1 2 3 4]) )
error('Unknown stopping criterion: ''%d''',num2str(stopcrit));
end
switch stopcrit
case 1
tol=[];
case { 2 , 3 }
if not(exist('tol','var'))
tol=1e-4;
end
case 4
if not(exist('tol','var'))
tol=1+1/mean(gn);
end
end
if (verb > 0)
par_print();
end
%%%%%%%%%%%%%%
% data scaling
%%%%%%%%%%%%%%
scaling = max(gn(:));
sqrtscaling = sqrt(scaling);
gn = gn/scaling;
bg = bg/scaling;
x = x/scaling;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% change the null pixels of the observed image
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
vmin = min( gn(gn > 0) );
gn(gn<=0) = vmin*eps*eps;
%gn(gn < tolpos) = max(tolpos, vmin*eps*eps);
% embed gn in a larger array
gn_work = zeros(work_size);
gn_work = gn_work(:);
gn_work( img_work_ind ) = gn;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% some computations needed only once
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N = numel(gn); % pixels in the image
flux = sum(gn) - N * bg; % exact flux
iter = 1; % iteration counter
alpha = initalpha; % initial alpha
Valpha = alpha_max * ones(Malpha,1); % memory buffer for alpha
Fold = -1e30 * ones(M, 1); % memory buffer for obj. func.
Discr_coeff = 2/N*scaling; % discrepancy coefficient
ONE = ones(work_size);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% arrays for boundary effects removal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
img_mask = zeros(work_size);
img_mask( img_work_ind ) = 1; % Mask_S
img_mask = img_mask(:);
ATMimg = AT(img_mask); % conv(PSF,M_S)
ATMimg = ATMimg(:);
gtsigma_ind = find(ATMimg >= sigma_tol); % indices of conv(PSF,M_S) >= sigma, that is the domain R
if ( modifiedWeights )
Weight = ones(size(ATMimg));
else
Weight = zeros(size(ATMimg));
end
Weight(gtsigma_ind) = 1.0 ./ ATMimg(gtsigma_ind); % M_R / alpha
R_mask = zeros( size(img_mask) );
R_mask( gtsigma_ind ) = 1;
% save('img_mask.mat', 'img_mask');
% save('ATMimg.mat', 'ATMimg');
% save('R_mask.mat', 'R_mask');
if ( modifiedGradient )
clear ATMimg;
else
ATMimg = ATMimg .* R_mask;
end
AT = @(x)( R_mask .* Afunction(x,CTF,work_size) );
A = @(x)( img_mask .* Afunction(x,TF,work_size) );
%%%%%%%%%%%%%%%%
% starting point >>> ATTENZIONE: sovrascrive i precedenti (vedi IDL)
%%%%%%%%%%%%%%%%
x_work = zeros(work_size);
x_work = x_work(:);
x_work(gtsigma_ind) = flux / numel(gtsigma_ind);
%%%%%%%%%%%%%%%%%
% projection type
%%%%%%%%%%%%%%%%%
switch (omega)
case 'nonneg'
pflag = 0;
case 'nonneg_eq'
pflag = 1;
otherwise
error('projection %s is not implemented',omega);
end
%%%%%%%%%%%%%%%%%%%%%
% output dir handling
%%%%%%%%%%%%%%%%%%%%%
if mysave
[success, mess] = mkdir(dest_dir);
if not(success)
error('%s: %s',dest_dir,mess);
end
if not(isempty(mess))
fprintf('%s\n\n',mess);
end
end
%%%%%%%%%%%%%%%%%%%
% vector allocation
%%%%%%%%%%%%%%%%%%%
if errflag
err = zeros(MAXIT+1,1);
obj = obj(:);
obj = obj/scaling;
obj_sum = sum(obj.*obj);
obj_work = zeros( work_size );
obj_work = obj_work(:);
obj_work_ind = img_work_ind;
obj_work( obj_work_ind ) = obj;
if (RIflag)
KLobjxk = zeros(MAXIT+1,1);
objnzind = find(obj_work);
objnz = obj_work( objnzind );
% objsum = sum(obj);
objsum = sum(objnz);
end
end
discr = zeros(MAXIT+1,1);
times = zeros(MAXIT+1,1);
times(1) = 0;
%%%%%%%%%%%%%%
% start of SGP
%%%%%%%%%%%%%%
% projection of the initial point
switch (pflag)
case 0 % non negativity
x_work( x_work < 0 ) = 0;
case 1 % non negativity and flux conservation
% we have no diagonal scaling matrix yet, so
% we project using euclidean norm
[x_work, biter, siter, r] = projectDF(flux, x_work, ones(size(x_work)));
end
% error
if errflag
e = x_work( obj_work_ind ) - obj;
err(1) = sqrt(sum(e.*e)/obj_sum);
% Relative Improvement
if (RIflag) % Compute KL(obj,gn) just once
% KLobjgn = sum( obj.*log( obj ./ gn_work(obj_work_ind) ) + sum(gn_work(obj_work_ind)) - objsum; % = KL(obj,gn)
% KLobjxk(1) = sum( obj.* log( obj ./ x_work(obj_work_ind) ) ) + sum(x_work(obj_work_ind)) - objsum; % = KL(obj,xk)
KLobjgn = sum( objnz.*log( objnz ./ gn_work(objnzind) ) ) + sum(gn_work(objnzind)) - objsum; % = KL(obj,gn)
KLobjxk(1) = sum( objnz.* log( objnz ./ x_work(objnzind) ) ) + sum(x_work(objnzind)) - objsum; % = KL(obj,xk)
end
end
%% work_ind must be set to the characteristic function of the
%% reconstruction domain R, which must be larger than the object domain S.
%work_ind = gtsigma_ind;
% objective gradient and function value
x_tf_work = A(x_work);
den = x_tf_work + bg;
temp = gn_work./den;
% g = ONE - AT(temp); substituted by the following
if ( modifiedGradient )
g_work = AT( img_mask - temp );
else
g_work = ATMimg - AT(temp);
end
fv = sum( gn_work(img_work_ind) .* log(temp(img_work_ind)) ) + sum( x_tf_work(img_work_ind) ) - flux;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% bounds for the scaling matrices
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
y_work = (flux/(flux+N*bg)).*AT(gn_work);
y = y_work( gtsigma_ind ); % y_work( work_ind ); % DEVE VALERE SU R
X_low_bound = min(y(y>0)); % Lower bound for the scaling matrix
X_upp_bound = max(y); % Upper bound for the scaling matrix
if X_upp_bound/X_low_bound < 50
X_low_bound = X_low_bound/10;
X_upp_bound = X_upp_bound*10;
end
fprintf('low: %e upp: %e\n', X_low_bound, X_upp_bound);
% discrepancy
discr(1) = Discr_coeff * fv;
% scaling matrix
if initflag == 0
X = ones(size(x_work));
else
X = x_work;
% bounds
X( X < X_low_bound ) = X_low_bound;
X( X > X_upp_bound ) = X_upp_bound;
end
X = X .* Weight;
if pflag == 1
D = 1./X;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% tolerance for stop criterion
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
switch stopcrit
case 2
if verb > 1
fprintf('it %04d || x_k - x_(k-1) ||^2 / || x_k ||^2 %e \n',iter-1,0);
end
% instead of using sqrt each iteration
tol = tol*tol;
case 3
if verb > 1
fprintf('it %04d | f_k - f_(k-1) | / | f_k | %e \n',iter-1,0);
end
case 4
if verb > 1
fprintf('it %04d D_k %e \n',iter-1,discr(1));
end
end
if ( suspendIter == 1 )
if (verb > 0)
fprintf('\n >> Suspending computation to file %s.mat\n',suspendFile);
end
save( suspendFile );
loop = false;
else
loop = true;
end
end
flag = 0; lam = 0; fr = 0; gd = 0; alpha1 = 0; alpha2 = 0;
if (fullVerboseFlag)
fprintf('IT | fv | fr | gd | alpha | lambda | alpha1 | alpha2 | F | tauAlpha \n');
fprintf('% 4d %+e %+e %+e %+e %+e %+e %+e % 2d %+e\n', iter-1, fv, fr, gd, alpha, lam, alpha1, alpha2, flag, tau);
end
%%%%%%%%%%%
% main loop
%%%%%%%%%%%
while loop
% store alpha and objective function values
Valpha(1:Malpha-1) = Valpha(2:Malpha);
if (Malpha > 1), Fold(1:M-1) = Fold(2:M); end
Fold(M) = fv;
% compute descent direction
y_work = x_work - alpha*X.*g_work;
switch (pflag) % projection onto the feasible set
case 0 % non negativity
y_work( y_work < 0 ) = 0;
case 1 % non negativity and flux conservation
[y_work, biter, siter, r] = projectDF(flux, y_work.*D, D);
end
d_work = y_work - x_work;
% backtracking loop for linesearch
gd = dot(d_work,g_work);
lam = 1;
fcontinue = 1;
d_tf_work = A(d_work); % exploiting linearity
fr = max(Fold);
while fcontinue
xplus_work = x_work + lam*d_work;
x_tf_try = x_tf_work + lam*d_tf_work;
den = x_tf_try + bg;
temp = gn_work ./ den;
fv = sum( gn_work(img_work_ind) .* log(temp(img_work_ind)) ) + sum( x_tf_try(img_work_ind) ) - flux;
if ( fv <= fr + gamma * lam * gd || lam < 1e-12)
x_work = xplus_work; clear xplus_work;
sk = lam*d_work;
x_tf_work = x_tf_try; clear x_tf_try;
if ( modifiedGradient )
gtemp = AT( img_mask - temp );
else
gtemp = ATMimg - AT(temp);
end
yk = gtemp - g_work;
g_work = gtemp; clear gtemp;
fcontinue = 0;
else
lam = lam * beta;
end
end
if (fv >= fr) && (verb > 0)
disp(' Worning: fv >= fr');
end
% update the scaling matrix and the steplength
X = x_work;
X( X < X_low_bound ) = X_low_bound;
X( X > X_upp_bound ) = X_upp_bound;
X = X .* Weight;
if ( modifiedWeights )
D = 1./X;
else
D = zeros(size(X));
D( gtsigma_ind ) = 1./X( gtsigma_ind );
end
sk2 = sk.*D; yk2 = yk.*X;
%fprintf('internal yk''yk: %e\n', dot(yk,yk));
bk = dot(sk2,yk); ck = dot(yk2,sk);
if (bk <= 0)
alpha1 = min(10*alpha,alpha_max);
else
alpha1BB = sum(dot(sk2,sk2))/bk;
alpha1 = min(alpha_max, max(alpha_min, alpha1BB));
end
if (ck <= 0)
alpha2 = min(10*alpha,alpha_max);
else
alpha2BB = ck/sum(dot(yk2,yk2));
alpha2 = min(alpha_max, max(alpha_min, alpha2BB));
end
Valpha(Malpha) = alpha2;
if (iter <= 20)
flag = 2;
alpha = min(Valpha);
elseif (alpha2/alpha1 < tau)
flag = 2;
alpha = min(Valpha);
tau = tau*0.9;
else
flag = 1;
alpha = alpha1;
tau = tau*1.1;
end
alpha = double(single(alpha));
iter = iter + 1;
times(iter) = cputime - t0;
if errflag
e = x_work(img_work_ind) - obj;
err(iter) = sqrt(sum(e.*e)/obj_sum);
% Relative Improvement
if (RIflag)
% KLobjxk(iter) = sum( obj.* log( obj ./ x_work(obj_work_ind) ) ) + sum(x_work(obj_work_ind)) - objsum; % = KL(obj,xk)
KLobjxk(iter) = sum( objnz.* log( objnz ./ x_work(objnzind) ) ) + sum(x_work(objnzind)) - objsum; % = KL(obj,xk)
end
end
discr(iter) = Discr_coeff * fv;
%%%%%%%%%%%%%%%
% stop criteria
%%%%%%%%%%%%%%%
switch stopcrit
case 1
if verb > 1
fprintf('\nit %04d of %04d',iter-1,MAXIT);
end
case 2
normstep = dot(sk(img_work_ind),sk(img_work_ind)) / dot(x_work(img_work_ind),x_work(img_work_ind));
loop = (normstep > tol);
if verb > 1
fprintf('\nit %04d || x_k - x_(k-1) ||^2 / || x_k ||^2 %e tol %e',iter-1,normstep,tol);
end
case 3
reldecrease = abs(fv - Fold(M)) / abs(fv);
loop = (reldecrease > tol);
if verb > 1
fprintf('\nit %04d | f_k - f_(k-1) | / | f_k | %e tol %e',iter-1,reldecrease,tol);
end
case 4
loop = (discr(iter) > tol);
if verb > 1
fprintf('\nit %04d D_k %e tol %e',iter-1,discr(iter),tol);
end
end
if iter > MAXIT
loop = false;
end
%%%%%%%%
% images
%%%%%%%%
if mysave
if ( length(savefreq) > 1 && any(iter-1 == savefreq) ) ...
|| ( length(savefreq) == 1 && ~rem(iter-1, savefreq) )
if (verb > 1)
fprintf(' saving intermediate results...');
end
% reconstructed image
% filename=sprintf('%s%crec_%04d.fits',dest_dir,filesep,iter-1);
% fitswrite((reshape(x,obj_size))',filename);
filename = sprintf('%s%crec_%04d.mat',dest_dir,filesep,iter-1);
tmp = reshape( x_work(img_work_ind), img_size ) * scaling;
save(filename,'tmp');
% residual image
res = sqrtscaling*(den(img_work_ind) - gn)./sqrt(den(img_work_ind));
% filename=sprintf('%s%cputter_res_%04d.fits',dest_dir,filesep,iter-1);
% fitswrite((reshape(res,obj_size))',filename);
filename = sprintf('%s%cputter_res_%04d.mat',dest_dir,filesep,iter-1);
tmp = reshape(res,img_size);
save(filename,'tmp');
savedLastIter = 1;
else
savedLastIter = 0;
end
end
if ( iter == suspendIter )
if (verb > 0)
fprintf('\n >> Suspending computation to file %s.mat\n',suspendFile);
end
save( suspendFile );
break;
end
if (fullVerboseFlag)
fprintf('% 4d %+e %+e %+e %+e %+e %+e %+e % 2d %+e\n', iter-1, fv, fr, gd, alpha, lam, alpha1, alpha2, flag, tau);
end
end
x = reshape( x_work(img_work_ind), img_size ) * scaling;
%%%%%%%%%%%%%%%%%%%
% save final images
%%%%%%%%%%%%%%%%%%%
if mysave
if ( ~savedLastIter )
% reconstructed image
% filename=sprintf('%s%crec_%04d.fits',dest_dir,filesep,iter-1);
% fitswrite((reshape(x,obj_size))',filename);
if (verb > 1)
fprintf(' saving last iteration...\n');
end
filename = sprintf('%s%crec_%04d.mat',dest_dir,filesep,iter-1);
tmp = x; % it is already reshaped
save(filename,'tmp');
% residual image
res = sqrtscaling*(den(img_work_ind) - gn)./sqrt(den(img_work_ind));
% filename=sprintf('%s%cputter_res_%04d.fits',dest_dir,filesep,iter-1);
% fitswrite((reshape(res,obj_size))',filename);
filename = sprintf('%s%cputter_res_%04d.mat',dest_dir,filesep,iter-1);
tmp = reshape(res,obj_size);
save(filename,'tmp');
end
end
if errflag
err = err(1:iter);
% Relative Improvement
if (RIflag)
KLobjxk = KLobjxk(1:iter);
varargout{1} = KLobjxk / KLobjgn;
end
end
discr = discr(1:iter);
times = times(1:iter);
iter = iter - 1;
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function par_print()
% par_print - Utility inner function to print-out parameters' values
initflag = evalin('caller','initflag');
switch initflag
case 0
temp = 'zeros';
case 1
temp = 'randn';
case 2
temp = 'gn';
case 3
temp = 'const';
case 999
temp = 'input';
otherwise
error(['Unrecognized initflag ''' num2str(initflag) '''']);
end
fprintf('Initialization %s\n',temp);
fprintf('MAXIT %d\n',evalin('caller','MAXIT'));
M = evalin('caller','M');
if M > 1
temp = 'nonmonotone';
else
temp = 'monotone';
end
fprintf('M %d (%s)\n',M,temp);
fprintf('gamma %e\n',evalin('caller','gamma'));
fprintf('beta %g\n',evalin('caller','beta'));
fprintf('alpha_min %e\n',evalin('caller','alpha_min'));
fprintf('alpha_max %e\n',evalin('caller','alpha_max'));
fprintf('Malpha %d\n',evalin('caller','Malpha'));
fprintf('tau %g\n',evalin('caller','tau'));
fprintf('initalpha %g\n',evalin('caller','initalpha'));
fprintf('background %g\n',evalin('caller','bg'));
errflag = evalin('caller','errflag');
switch errflag
case 0
temp = 'false';
case 1
temp = 'true';
otherwise
error(['unrecognized errflag ''' num2str(errflag) '''']);
end
fprintf('err calc %s\n',temp);
fprintf('verbosity %d\n',evalin('caller','verb'));
stopcrit = evalin('caller','stopcrit');
switch stopcrit
case 1
temp = 'iterations';
case 2
temp = 'relstepx or iterations';
case 3
temp = 'relstepf or iterations';
case 4
temp = 'discr. val. or iterations';
end
fprintf('stop criterion %s\n',temp);
fprintf('tolerance %e\n',evalin('caller','tol'));
fprintf('constraints %s\n',evalin('caller','omega'))
mysave = evalin('caller','mysave');
if mysave
dest_dir = sprintf('%s%s%s',pwd,filesep,evalin('caller','dest_dir'));
else
dest_dir = 'no';
end
fprintf('save %s\n',dest_dir);
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [SGPpar] = par_init()
% par_init - Utility inner function to initialize parameters' values
%
%%%%%%%%%%%%%%%%%%%%%%%%
% SGP default parameters
%%%%%%%%%%%%%%%%%%%%%%%%
SGPpar = struct(...
'MAXIT', 1000, ... % maximum number of iterations
'gamma', 1e-4, ... % for sufficient decrease
'beta', 0.4, ... % backtracking parameter
'M', 1, ... % memory in obj. function value (if M = 1 monotone)
'alpha_min', 1e-5, ... % alpha lower bound
'alpha_max', 1e5, ... % alpha upper bound
'Malpha', 3, ... % alfaBB1 memory
'tau', 0.5, ... % alternating parameter
'initalpha', 1.3, ... % initial alpha
'initflag', 0, ... % 0 -> initial x all zeros
'errflag', false, ... % 0 -> no error calculation
'err', [], ... % errors w.r.t. the object, if it is given
'verb', 0, ... % 0 -> silent
'stopcrit', 1, ... % 1 -> number of iterations
'bg', 0, ... % background value
'omega', 'nonneg', ... % non negativity constraints
'AT', [], ... % function handle for A'*x
'dest_dir', '.', ... % where to store intermediate and final results
'mysave', false, ... % true = save, false = don't save intermediate resutls
'savefreq', 1, ... % iteration interval for storing intermediate results
'savedLastIter', 0); % flag: SET BY THE PRORAM
% ==============================================================================
% End of sgp_deblurring.m file - IRMA package
% ==============================================================================
