 function [psf, fwhm, tbeta] = res2beta(G, C, psf_fhandle, reg_beta, psf_offset, mask, target_fwhm)
%function [psf, fwhm, [tbeta]] = res2beta(G, C, psf_fhandle, reg_beta, psf_offset, mask, [target_fwhm])
% Builds and makes the system matrix and penalty function
% for resolution analysis.  This is implemented both
% using CG recon and according to the 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
%       target_fwhm             [Nbd,1]                 target fwhm [voxels]
% Out:
%       psf                     [Nbd,Nbetas,Noffsets]   point spread function evaluated at reg_beta
%       fwhm                    [Nbd,Nbetas,Noffsets]   full width half max evaluated at reg_beta
%       tbeta                   [Nbd,1]                 beta value for the target_fwhm
%
%  Copyright 2005-6-20, Valur Olafsson, The University of Michigan

if (nargin < 6) | (nargin > 7)
    help mfilename
    error 'Wrong number of input variables'
end
if (nargin == 7) & (nargout ~= 3)  % target_fwhm is an input varible
    help mfilename
    error 'Need to have 3 output arguments when using target_fwhm'
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
            psf(:,ii,kk) = psf_tmp(:);
            fwhm(1,ii,kk) = fwhm_tmp;

            % 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
            % 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
            psf(:,ii,kk) = psf_tmp(:);
            fwhm(1,ii,kk) = fwhm_tmp;

            % Clear variables
            clear psf_tmp var_tmp fwhm_tmp
        case 'qpwls_psf_st'
            % read psf into a temporary variable
            psfR_tmp = psf_tmp(1:end/2,:);
            psfI_tmp = psf_tmp(end/2+1:end,:);
            
            % Relocating the impulse response to correct coordinate
            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
            psf(:,ii,kk) = [psfR_tmp(:);psfI_tmp(:)];
            fwhm_tmp, fwhm(:,ii,kk) = fwhm_tmp(:);

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

if exist('target_fwhm','var')
    switch func2str(psf_fhandle)
    case 'qpwls_psf_st'
        if all((max(mean(fwhm,3),[],2) > target_fwhm) & (min(mean(fwhm,3),[],2) < target_fwhm))
            % x,y coordinate is fwhm averaged over the offsets
            betaR = griddata(mean(fwhm(1,:,:),3).',mean(fwhm(2,:,:),3).',reg_beta(:,1),target_fwhm(1),target_fwhm(2),'cubic');
            betaI = griddata(mean(fwhm(1,:,:),3).',mean(fwhm(2,:,:),3).',reg_beta(:,2),target_fwhm(1),target_fwhm(2),'cubic');
            
            % set output variable
            tbeta = [betaR;betaI];
        else
            warning 'target fwhm is not within the range of fwhm values calculated reg_beta, returning tbeta = [NaN;NaN]'
            tbeta = [nan;nan];
        end
    case 'qpwls_psf'
        if (max(mean(fwhm,3),[],2) > target_fwhm) & (min(mean(fwhm,3),[],2) < target_fwhm)
            % x,y coordinate is fwhm averaged over the offsets
            tbeta = interp1(mean(fwhm,3),reg_beta,target_fwhm(1),'cubic');
        else
            warning 'target fwhm is not within the range of fwhm values calculated reg_beta, returning tbeta = NaN'
            tbeta = nan;
        end
    case 'qpwls_psf_gram'
        if (max(mean(fwhm,3),[],2) > target_fwhm) & (min(mean(fwhm,3),[],2) < target_fwhm)
            % x,y coordinate is fwhm averaged over the offsets
            tbeta = interp1(mean(fwhm,3),reg_beta,target_fwhm(1),'cubic');
        else
            warning 'target fwhm is not within the range of fwhm values calculated reg_beta, returning tbeta = NaN'
            tbeta = nan;
        end
    end
end
