 function kappa = calc_kappa(A, r2s_base, nhist)
%function kappa = calc_kappa(A, [r2s_base,] [nhist])
% This function finds the kappa values needed to fix
% spatially nonuniform penalties
%
% in
%       A                       system matrix
%       r2s_base                median baseline R2* used to find \beta (default 21.7 1/sec [3T-value])
%       nhist                   number of bins in histogram (default 100)
% out
%       kappa                   [[Nd]], uses same mask as A
%
% If A is of type 'fast_mr' then kappa is defined as:
%
%                  { \sum_{l=1}^{L} \Phi_l^2 \exp{-2 t_l R_n} }
%   kappa_n = \sqrt{ -----------------------------------------}
%                  {          \sum_{l=1}^{L} \Phi_l^2         }
%
% If A is of type 'diff_r2' then kappa is defined as:
%
%                        { \sum_{l=1}^{L} \Phi_l^2 t_l^2 \exp{-2 t_l R_n} }
%   kappa_n = |f_n| \sqrt{ -----------------------------------------------}
%                        {          \sum_{l=1}^{L} \Phi_l^2 t_l^2         }
%
% Copyright 2006-2-15, Valur Olafsson, The University of Michigan

% Check input arguments
if nargin > 4, error(sprintf('Incorrect number of input arguments for %s',mfilename)); end;

% Defaults
if ~isvar('r2s_base')
    r2s_base = mean([1/55.1e-3 1/39.6e-3]);  % From JMRI 9:531--38 (1999)
elseif mod(r2s_base,1) == 0 %CHECK IF CORRECT INPUT
    if r2s_base ~= 0
        warning 'r2s_base is an integer number...this is highly unlikely to occur'
        keyboard
    end
end
if ~isvar('nhist')
    nhist = 100;
elseif ~(mod(nhist,1) == 0)
    warning 'nhist should be an integer'
    keyboard
end


% Set chat = true for further computations
chat = false;

switch class(A)
case {'fast_mr','multi_echo_fastmr'}
    if isa(A,'multi_echo_fastmr')
        t = [];
        for ii = 1:length(A.a)
            t(:,ii) = A.a{ii}.tt;
        end
        t = t(:);
        A = A.a{1};
    else
        t = A.tt;
    end
    z = i*A.we;
    r2s = real(z);
    Nd = length(z);
    mask = A.G.arg.mask;
    Nt = length(t);
    if A.flgswt == 1
        swt = repmat(A.swt,Nt./length(A.swt),1);
    else
        swt = ones(Nt,1);
    end

    if sum(abs(r2s(:))) == 0,  % No R2* present
        kappa = mask;
    else
        r2sc = linspace(min(r2s),max(r2s),nhist);
        hedges = [r2sc-(max(r2s)-min(r2s))./(nhist-1)./2 r2sc(end)+(max(r2s)-min(r2s))./(nhist-1)./2];
        [hcount hindex] = histc(r2s,hedges);
        if chat % Show the histogram
            figure, hist(r2s,r2sc);
            prompt;
        end

        % Now to calculate numerator of kappa
        numer = sum(swt(:,ones(1,nhist)).^2.*exp(-2*t*r2sc),1);
        denom = sum(swt.^2.*exp(-2*t*r2s_base));
        dn = sqrt(numer./denom).';

        % Make kappa using results from the histogram
        kappa = dn(hindex);

        if chat % Compare histgram vs finer scaled plotting of r2s
            r2stmp = [1:.5:50];
            nr2s = length(r2stmp);
            numertmp = sum(swt(:,ones(1,nr2s)).^2.*exp(-2*t*r2stmp),1);
            % denomtmp = sum(swt.^2);
            % dntmp = sqrt(numertmp./denomtmp).';

            fh1 = figure; plot(r2stmp,numertmp,r2sc,numer); title('Numerator of kappa');
            xlabel('R_2^* values [1/s]')

            for ii = 1:Nd
                if mod(ii,100) == 0, disp(sprintf('ii = %i',ii)); end;
                en = zeros(Nd,1);
                en(ii) = 1;
                kappar(ii,1) = norm(A*en,2)/sqrt(sum(swt.^2.*exp(-2*t*r2s_base)));
            end
            
            fh2 = figure; im(221, embed(kappa,A.G.arg.mask),'Histogram Kappa'); colorbar;
            im(222, embed(kappar,A.G.arg.mask),'Real Kappa'); colorbar;
            im(223, embed(abs(kappa-kappar),A.G.arg.mask),'Kappa - Error'); colorbar;
            %im(224, embed(exp(-t(end)*r2s),mask),'kappa = e^{-t_{end}*R_2^*}'); colorbar;
            prompt
            close([fh1 fh2])
        end

        % Embedding kappa in mask
        kappa = embed(kappa,mask);
    end
case {'diff_r2','diff_r2v2','multi_echo_diffr2'}
    if isa(A,'multi_echo_diffr2')
        t = [];
        for ii = 1:length(A.a)
            t(:,ii) = A.a{ii}.A.tt;
        end
        t = t(:);
        A = A.a{1};
    else
        t = A.A.tt;
    end
    z = i*A.A.we;
    Nd = length(z);
    r2s = real(z);
    rho = A.rho;
    mask = A.A.G.arg.mask;
    %t = A.A.tt;
    Nt = length(t);
    if A.A.flgswt == 1
        swt = repmat(A.A.swt, Nt./length(A.A.swt),1);
    else
        swt = ones(Nt, 1);
    end

    if sum(abs(r2s)) == 0,
        error(sprintf('Missing R2* map in %s',mfilename));
    else
        r2sc = linspace(min(r2s),max(r2s),nhist);
        hedges = [r2sc-(max(r2s)-min(r2s))./(nhist-1)./2 r2sc(end)+(max(r2s)-min(r2s))./(nhist-1)./2];
        [hcount hindex] = histc(r2s,hedges);
        if chat % Show the histogram
            figure, hist(r2s,r2sc);
            prompt;
        end

        % Now to calculate numerator of kappa
        numer = sum(swt(:,ones(1,nhist)).^2.*t(:,ones(1,nhist)).^2.*exp(-2*t*r2sc),1);
        denom = sum(swt.^2.*t.^2.*exp(-2*t*r2s_base));
        dn = sqrt(numer./denom).';

        if 0
            % Extend rho outside the object
            % OLD METHOD
            %erode_suppmask = logical(imerode(embed(abs(rho)>.2*max(abs(rho)),mask),ones(3)));
            %rhos = phase_denoise_fun(embed(abs(rho),mask),erode_suppmask,-10,200); rhos = rhos(mask);
            % NEW METHOD *smarter in extending the kappa map*
            rhos = mri_phase_denoise(embed(abs(rho).*exp(i*abs(rho)./max(abs(rho(:)))),mask), ...
                     'l2b',-5,'order',1,'niter',200,'wi_ml',true,'pl',0)*max(abs(rho(:)));
            rhos = rhos(mask);
        else
            % NB! if this is enabled, get noisy edges in the reconstruction
            rhos = abs(rho)+eps;
        end

        % Make kappa using results from the histogram
        kappa = abs(rhos).*dn(hindex);

        if chat % Compare histgram vs finer scaled plotting of r2s
            r2stmp = [1:.5:50];
            nr2s = length(r2stmp);
            numertmp = sum(swt(:,ones(1,nr2s)).^2.*t(:,ones(1,nr2s)).^2.*exp(-2*t*r2stmp),1);
            % denomtmp = sum(swt.^2.*t.^2);
            % dntmp = sqrt(numertmp./denomtmp).';

            fh1 = figure; plot(r2stmp,numertmp,r2sc,numer); title('Numerator of kappa');
            xlabel('R_2^* values [1/s]')

            %fh3 = figure; plot(1:length(r2s),exp(-t(end)*r2s),'-x',1:length(r2s),dn(hindex),'-o');
            
            disp('Finding Non-histogramed kappa - SLOW!!!')
            Atmp = diff_r2(A.A,rhos);
            for ii = 1:Nd
                if mod(ii,100) == 0, disp(sprintf('ii = %i',ii)); end;
                en = zeros(Nd,1);
                en(ii) = 1;
                kappar(ii,1) = norm(Atmp*en,2)/sqrt(sum(swt.^2.*t.^2.*exp(-2*t*r2s_base)));
            end
            
            fh2 = figure; im(221, embed(kappa,mask),'Histogram Kappa'); colorbar;
            im(222, embed(kappar,mask),'Real Kappa'); colorbar;
            im(223, embed(abs(kappa-kappar),mask),'Kappa - Error'); colorbar;
            %im(224, embed(abs(rhos).*exp(-t(end)*r2s),mask),'kappa = |f|*e^{-t_{end}*R_2^*}'); colorbar;
            prompt
            close([fh1 fh2])
        end

        % Embedding kappa in mask
        kappa = embed(kappa,mask);
    end
otherwise
    error(sprintf('%s doesn''t work for %s',mfilename,class(A)));
end
