% function SPATIALBASESGEN constructs spatial basis vectors from the
% collection of dipole forward solutions. The spatial bases span approximately
% the same space as the dipole forward solutions.
%
% inputs: 1) ForwardSoln is the collection of forward solution for K dipoles,
%         which is either
%           a) a 3-D matrix with dimension [KxNx3] or [KxNx2], where each
%         of the K [Nx3] or [Nx2] matrix is the forward solution for a single
%         dipole with fixed but unknown moment. N is the number of sensor
%         channels
%         OR
%           b) a 2-D matrix with dimension [NxK], where each Nx1 column is
%         the forward solution for a single dipole with known moment.
%         2) PatchID is a [Kx1] vector with only 0 or 1 entries. The
%         locations of 1's in the vector corresponds to
%         the locations of dipoles that belong to that patch.
%         3) P - a positive scalar value which can have 2 different meanings
%            0<P<1, real value, corresponds to the fraction of energy that the
%            constructed spatial bases preserve from the space originally spanned by
%            the collection of dipole forward solutions. The number of
%            spatial bases is determined automatically.
%            P>=1, integer only, corresponds to the total number of bases.
%            We found that P = 3 is usually sufficient for a patch with area less
%            than 1000 mm2.
%
% output: An [NxM] matrix U, where the M columns of U are the spatial basis vectors
%         (M = P if P is an integer >= 1)
%
% Author: Tulaya Limpiti
% last updated: 5/25/2009

function U = SpatialBasesGen(ForwardSoln,PatchID,P)

% find dipoles that belong to the patch
dipoleLoc = find(PatchID);

% concatenate all relevant dipole forward solutions

if (length(size(ForwardSoln)) == 2)
    % [Nx1] forward solution with known moment
    Hr = ForwardSoln(:,dipoleLoc);

elseif (length(size(ForwardSoln)) == 3)
    % [Nx2] or [Nx3] forward solution with fixed but unknown moment
    Htmp = ForwardSoln(dipoleLoc,:,:);
    Hr = squeeze(Htmp(1,:,:));
    if (length(dipoleLoc) >1)
        for i = 2:length(dipoleLoc)
            Hr = [Hr, squeeze(Htmp(i,:,:))];
        end
    end
    clear Htmp
end

% construct spatial bases
[Uf,U,Uorth,Sr,Vr,Pr] = rnkreduce(Hr,P);
clear Sr Vr Uorth Uf