function [dipout] = minimumnormestimate(dip, grad, vol, dat, varargin);

% LORETA reconstructs the current density in a distributed source model
% using the LORETA algorithm.
%
% Use as
%   [dipout] = loreta(dip, grad, vol, dat, ...)
%
% The algorithm is described in detail in
% Pascual-Marqui RD, Michel CM, and Lehmann D (1994) Low resolution
% electromagnetic tomography: a new method for localizing electrical
% activity in the brain. Int. J. Psychophysiol. 18:49-65.

% Copyright (C) 2005, Robert Oostenveld
%
% $Log: loreta.m,v $
% Revision 1.3  2007/01/08 09:07:35  roboos
% multiple changes, no idea what exactly
%
% Revision 1.2  2005/10/05 06:32:59  roboos
% many changes, still work in progress
%
% Revision 1.1  2005/09/29 08:48:30  roboos
% new implementation, created at the EEGLAB workshop in Porto
%

% get the values of the optional arguments
lambda        = keyval('lambda',        varargin); if isempty(lambda  ),      lambda = 0;                   end
keepleadfield = keyval('keepleadfield', varargin); if isempty(keepleadfield), keepleadfield = 'no';         end
chaninv       = keyval('chaninv', varargin);       if isempty(chaninv),       chaninv = 'yes';              end
% convert the yes/no arguments to the corresponding logical values
keepleadfield = strcmp(keepleadfield, 'yes');
chaninv       = strcmp(chaninv, 'yes');

% ensure that these are row-vectors
dip.inside = dip.inside(:)';
dip.outside = dip.outside(:)';

Nchan = size(dat,1);
Ndip  = prod(dip.dim);

if Ndip~=(length(dip.inside)+length(dip.outside))
  error('ambiguity in the number of dipoles in the grid');
end

% concatenate the leadfield components of all sources into one large matrix
% the sources outside the brain will have a leadfield of zero
A = zeros(Nchan, 3*Ndip);
V = spalloc(3*Ndip, 3*Ndip, 3*Ndip);
sel = zeros(1, 3*Ndip);
for i=dip.inside;
  if ~isfield(dip, 'leadfield')
    lf = compute_leadfield(dip.pos(i,:), grad, vol);
  else
    lf = dip.leadfield{i};
  end
  lx = (i-1)*3+1;
  ly = (i-1)*3+2;
  lz = (i-1)*3+3;
  A(:,[lx ly lz]) = lf;
  % compute the norm of each leadfield column
  V(lx,lx) = norm(lf(:,1));
  V(ly,ly) = norm(lf(:,2));
  V(lz,lz) = norm(lf(:,3));
  if keepleadfield
    dipout.leadfield{i} = lf;
  end
  sel(lx) = 1;
  sel(ly) = 1;
  sel(lz) = 1;
  if keepleadfield
    dip.leadfield{i} = lf;
  end
end
sel = find(sel);

fprintf('computing discrete spatial 3-D laplacian\n');
B = loreta_lap3d(dip);

tic

% assign the a-priori source covariance R
Rinv = V' * B' * B * V;
if ~chaninv
  R = pinv(full(Rinv));
end

% assume an identity noise matrix, scaled with lambda
if lambda>0
  C    = sparse(eye(Nchan,Nchan)) * lambda;
  Cinv = sparse(eye(Nchan,Nchan)) / lambda;
else
  C    = sparse(zeros(Nchan,Nchan));
  Cinv = sparse(eye(Nchan,Nchan));
end

% compute the weights
if chaninv
  W = inv(A' * Cinv * A + Rinv) * A' * Cinv;
else
  W = R * A' * inv(A * R * A' + C);
end

toc

fprintf('computing LORETA source reconstruction\n');

% for each of the timebins, estimate the source strength
mom = W * dat;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% reassign the estimated source strength over the inside and outside dipoles
n = 1;
for i=dip.inside
  lx = (i-1)*3+1;
  ly = (i-1)*3+2;
  lz = (i-1)*3+3;
  dipout.mom{i} = mom([lx ly lz], :);
end
dipout.mom(dip.outside) = {nan};

% for convenience also compute power at each location
for i=dip.inside
  dipout.pow(i,:) = sum(dipout.mom{i}.^2, 1);
end
dipout.pow(dip.outside,:) = nan;

if keepleadfield
  dipout.leadfield(dip.inside) = dip.leadfield(dip.inside);
  dipout.leadfield(dip.outside) = {[]};
end

% add other descriptive information to the output source model
dipout.pos     = dip.pos;
dipout.inside  = dip.inside;
dipout.outside = dip.outside;