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 % ==============================================================================
ottobre 2014 »
ottobre
lumamegivesado
12345
6789101112
13141516171819
20212223242526
2728293031