function [dataout] = lyngby_filter(datain,filtertype,kernelX,kernelY,kernelZ)

% lyngby_filter		- function for averaging volumes
%
%	function [dataout] = lyngby_filter(datain,filtertype,kernelX,kernelY,kernelZ) 
%
%	Input:	datain		a data matrix where the rows are the time frames
%                        	and the volume indices correspond to the colums.
%		filtertype	If set to 1, then use a Hann window. 
%                        	If set to 0, then use a uniform (square) window.
%            	kernelX		Size of averaging kernel in the x-direction is
%				2*kernelX+1. If set to zero do not average in this
%                        	direction.
%		kernelY		Size of averaging kernel in the y-direction is
%				2*kernelY+1. If set to zero do not average in this
%				direction.
%		kernelZ		Size of averaging kernel in the z-direction is
%				2*kernelZ+1. If set to zero do not average in this
%				direction.
%
%	Output: dataout		The filtered volume

%	Date		Author		Version
%	1996-09-17	Peter Toft	Initial Version
%	1996-09-25	Finn Aa Nielsen	Renamed: filtfMRI -> lyngby_filter
%	1996-11-06	fnielsen	Cosmetically changed

% cvs :  $Revision: 1.1.1.1 $

femary_ui_log('Filtering');

X=ROI_VOXELS(1,2)-ROI_VOXELS(1,1);
Y=ROI_VOXELS(2,2)-ROI_VOXELS(2,1);
Z=ROI_VOXELS(3,2)-ROI_VOXELS(3,1);

[NoFrames NoVoxels]=size(datain);

if NoVoxels~=X*Y*Z,
  error('Error in size of input matrix');
end

kernelX=round(kernelX);  %Just checking it is an integer
kernelY=round(kernelY);  %Just checking it is an integer
kernelZ=round(kernelZ);  %Just checking it is an integer

if kernelZ==0,  %Just filter in x and y if kernelX=0
  Filt3D=0;
else
  Filt3D=1;
end

if filtertype==0,
  filtX=ones(2*kernelX+1,1);
  filtY=ones(2*kernelY+1,1);
  filtZ=ones(2*kernelZ+1,1);
else 
  if filtertype==1;
    if kernelX>0,
      filtX=0.5*(1-cos(2*pi/(2*kernelX+2)*(1:2*kernelX+1)));
    end
    if kernelY>0,
      filtY=0.5*(1-cos(2*pi/(2*kernelY+2)*(1:2*kernelY+1)));
    end
    if kernelZ>0,
      filtZ=0.5*(1-cos(2*pi/(2*kernelZ+2)*(1:2*kernelZ+1)));
    end
  end
end

filtX=filtX/sum(filtX);
filtY=filtY/sum(filtY);
filtZ=filtZ/sum(filtZ);

for f=1:NoFrames,
% Filter in x-direction
  if kernelX>0
    for k=0:Z-1,
      for j=0:Y-1,
        ar=datain(f,(1:X)+X*(j+k*Y));
        arout=conv(ar,filtX);
        datain(f,(1:X)+X*(j+k*Y))=arout(kernelX+1:kernelX+X);
      end
    end
    % Filter in y-direction
    if kernelY>0,
      for k=0:Z-1      
        for i=1:X,
          ar=datain(f,i+X*((0:Y-1)+k*Y));
          arout=conv(ar,filtY);
          datain(f,i+X*((0:Y-1)+k*Y))=arout(kernelY+1:kernelY+Y);
        end      
      end
    end
  end
  % Filter in z-direction  
  if Filt3D==1 & kernelZ>0,
    for k=0:Z-1
      for i=1:X,
        for j=0:Y-1,
          ar=datain(f,i+X*(j+(0:Z-1)*Y));
          arout=conv(ar,filtX);
          datain(f,i+X*(j+(0:Z-1)*Y))=arout(kernelZ+1:kernelZ+Z);
        end
      end
    end
  end
end

dataout=datain;
