 function [fwhm, psf] = beta2res(G, C, psf_fhandle, reg_beta, psf_offset, mask)
%function [fwhm, [psf]] = beta2res(G, C, psf_fhandle, reg_beta, psf_offset, mask)
% Calculate resolution properties for a range of beta values.
% This is implemented using fast method from
% Res. Prop. paper by Fessler and Rogers
% in:
%       G                       [Nt,N]                  system matrix for psf eval
%       C                                               penalty matrix with beta=1(!!!!) for psf eval
%       psf_fhandle:                                    function handle to the psf eval routine
%         @qpwls_psf                                    psf eval using fft when using 'qpwls_pcg'
%         @qpwls_psf_gram                               psf eval using fft when using 'qpwls_pcg'
%         @qpwls_psf_st                                 psf eval using fft when using 'qpwls_pcg_st'
%       reg_beta                [Nbetas,Nbd]            betas used in psf evaluation
%       psf_offset              [Noffsets,Nbd]          impulse offsets
%       mask                    [Nx,Ny]                 logical array
% out:
%       fwhm                    [Nbd,Nbetas,Noffsets]   full width half max evaluated at reg_beta
%       psf                     [Nbd,Nbetas,Noffsets]   point spread function evaluated at reg_beta
%
%  Copyright 2007-5-06, Valur Olafsson, The University of Michigan

if ~(nargin == 6)
    help mfilename
    error 'Wrong number of input variables'
end

dims = size(mask);

%%%%%
% Loop through all betas and local impulses
for kk = 1:size(psf_offset,1)
    offset = psf_offset(kk,:);
    for ii = 1:size(reg_beta,1)
        Cbeta = reg_beta(ii,:);
        
        % Finding the impulse response
        [psf_tmp, var_tmp, fwhm_tmp] = psf_fhandle(G, C, Cbeta, mask, 1, offset);
        
        switch func2str(psf_fhandle)
        case 'qpwls_psf'
            % Relocating the impulse response to correct coordinate
            % NO NEED FOR THIS FOR THE CURRENT IMPLIMENTATION OF qpwls_psf
            %psf_tmp = circshift(psf_tmp,offset);

            % Only show the results from voxels within mask
            psf_tmp = embed(psf_tmp(mask),mask);
            
            % Setting the output variables
            fwhm(1,ii,kk) = fwhm_tmp;
            if nargout == 2
                psf(:,ii,kk) = psf_tmp(:);
            end

            % Clear variables
            clear psf_tmp var_tmp fwhm_tmp
        case 'qpwls_psf_gram'
            warning 'This has not been tested'

            % Relocating the impulse response to correct coordinate
            % ONLY NEEDED FOR OLD VERSION
            %psf_tmp = circshift(psf_tmp,offset);

            % Only show the results from voxels within mask
            psf_tmp = embed(psf_tmp(mask),mask);
            
            % Setting the output variables
            fwhm(1,ii,kk) = fwhm_tmp;
            if nargout == 2
                psf(:,ii,kk) = psf_tmp(:);
            end

            % Clear variables
            clear psf_tmp var_tmp fwhm_tmp
        case {'qpwls_psf_st'}
            % read psf into a temporary variable
            psfR_tmp = psf_tmp(:,:,1);
            psfI_tmp = psf_tmp(:,:,2);
            
            % Relocating the impulse response to correct coordinate
            % ONLY NEEDED FOR OLD VERSION
%            psfR_tmp = circshift(psfR_tmp,offset);
%            psfI_tmp = circshift(psfI_tmp,offset);

            % Only show the results from voxels within mask
            psfR_tmp = embed(psfR_tmp(mask),mask);
            psfI_tmp = embed(psfI_tmp(mask),mask);
            
            % Setting the output variables
            fwhm(:,ii,kk) = fwhm_tmp(:);
            if nargout == 2
                psf(:,ii,kk) = [psfR_tmp(:);psfI_tmp(:)];
            end

            % clearing variables
            clear psf_tmp psfR_tmp psfI_tmp var_tmp fwhm_tmp
        otherwise
            error 'unknown psf function handle'
        end
    end
end
