% Multiplication of (generalized) inverse with data matrix
% If dim~=1, dim rows of the inverse are considered to belong to one location in the source space
% and for the solution only the absolute value for every location is output
% More measurement points (electrodes, sensors) than inverse coefficients are
% allowed (only the first measurement points are used), but not vice versa
%   This software is protected by german copyright and international treaties.             
%   Copyright 2004 Markus Junghöfer & Peter Peyk. All Rights Reserved.                     
%                                                                                          
%   THIS SOFTWARE AND DOCUMENTATION IS PROVIDED "AS IS," AND COPYRIGHT HOLDERS MAKE        
%   NO REPRESENTATIONS OR WARRANTIES, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO,    
%   WARRANTIES OF MERCHANTABILITY OR FITNESS FOR ANY PARTICULAR PURPOSE OR THAT THE USE    
%   OF THE SOFTWARE OR DOCUMENTATION WILL NOT INFRINGE ANY THIRD PARTY PATENTS, COPYRIGHTS,
%   TRADEMARKS OR OTHER RIGHTS. COPYRIGHT HOLDERS WILL NOT BE LIABLE FOR ANY DIRECT,       
%   INDIRECT, SPECIAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF ANY USE OF THE SOFTWARE      
%   OR DOCUMENTATION.                                                                      
% Different rows of rec correpond to different locations in the source space
% Different columns of rec correspond to different data points (time points)
% OH 29.10.97

function rec = inv_recon(inverse, data, dim);

[mi ni] = size(inverse); %mi=NDip;  ni=NChan
[md nd] = size(data);    %md=NChan; nd=NPoints

if (nargin==2)
	dim = 1;
end

if ni>md,
	disp('Fewer measurement points than inverse coefficients!!! (inv_recon) ');
	return;
end

if mod(mi,dim)~=0,
	disp('Dimension of sourcespace incompatible with number of solution points!!! (inv_recon) ');
	return;
end

if (nargin==2) | (dim==1),
	rec=inverse*data;
else
	cnt = 1;
	for i=1:dim:mi,	% Solution points (vectors)
		comp=inverse(i:i+dim-1,:)*data(1:ni,:);
		for j=1:nd	% Data points
			rec(cnt,j)=norm(comp(:,j));			
		end
		cnt=cnt+1;
	end
end
