function image = sos_combine(images,noisecov,smap_norm)
%
% function image = sos_combine(images,noisecov,smap_norm)
%
% Simple sum-of-squares coil combination using optional noise covariance.
%
% Inputs:
%   images - 4-d matrix of images (4th dimension is coil channels)
%   noisecov - (optional) noise covariance matrix; if empty, uses identity
%   smap_norm - 3-d matrix for normalizing by sum-of-squares of
%       sensitivities; used to mitigate central darkening
%
% Outputs:
%   image - reconstructed 3-d volume
%
% Copyright 2012-10-09, Daniel Weller, University of Michigan

[Nx,Ny,Nz,NCha] = size(images);

if ~exist('noisecov','var') || isempty(noisecov), noisecov = eye(NCha); end

% eigenvalue-based inverse of noise covariance
[V,D] = eig(noisecov);
d = abs(diag(D));
dmin = eps(max(d));
dinvsqrt = 1./sqrt(d);
dinvsqrt(d <= dmin) = 0; % truncate eigenvalues
noisecovinvsqrt = V*diag(dinvsqrt)*V';

images = noisecovinvsqrt*reshape(images,[],NCha).';
image = reshape(sqrt(sum(abs(images).^2,1)),[Nx,Ny,Nz]);

% normalize
if exist('smap_norm','var') && ~isempty(smap_norm), image = image./smap_norm; end

end
