function ft_realtime_fmriviewer(cfg)

% FT_REALTIME_FMRIVIEWER allows for realtime visualization of the fMRI data stream
%
% Use as
%   ft_realtime_fmriviewer(cfg)
% with the following configuration options
%   cfg.bufferdata = whether to start on the 'first or 'last' data that is available (default = 'last')
%
% The source of the data is configured as
%   cfg.dataset       = string
% or alternatively to obtain more low-level control as
%   cfg.datafile      = string
%   cfg.headerfile    = string
%   cfg.eventfile     = string
%   cfg.dataformat    = string, default is determined automatic
%   cfg.headerformat  = string, default is determined automatic
%   cfg.eventformat   = string, default is determined automatic

% Copyright (C) 2012, Stefan Klanke
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

if ~isfield(cfg, 'headerfile'),    	cfg.headerfile = 'buffer://localhost:1972'; end
if ~isfield(cfg, 'datafile'),      	cfg.datafile = cfg.headerfile; end
if ~isfield(cfg, 'eventfile'),    	cfg.eventfile = cfg.headerfile; end
if ~isfield(cfg, 'headerformat'),   cfg.headerformat = []; end
if ~isfield(cfg, 'dataformat'),   	cfg.dataformat = []; end
if ~isfield(cfg, 'bufferdata'),     cfg.bufferdata = 'last';  end % first or last

% translate dataset into datafile+headerfile
cfg = ft_checkconfig(cfg, 'dataset2files', 'yes');
cfg = ft_checkconfig(cfg, 'required', {'datafile' 'headerfile'});

% ensure that the persistent variables related to caching are cleared
clear ft_read_header
% start by reading the header from the realtime buffer
hdr = ft_read_header(cfg.headerfile);

disp(hdr);

if isfield(hdr,'nifti_1')
	ddim = double(hdr.nifti_1.dim);
    width  = ddim(1);
    height = ddim(2);
    numSlices = ddim(3);
elseif isfield(hdr,'siemensap') && isstruct(hdr.siemensap)
    width = siemensap.sKSpace.lBaseResolution;
    phaseFOV = siemensap.sSliceArray.asSlice{1}.dPhaseFOV;
    readoutFOV = siemensap.sSliceArray.asSlice{1}.dReadoutFOV;
    height = width * phaseFOV / readoutFOV;
    numSlices = siemensap.sSliceArray.lSize;
else
  ft_warning('No protocol information found!')
  width = sqrt(hdr.nChans);
  height = width;
  numSlices = 1;
end
  
% open the figure and set the colormap
h = figure;
colormap(gray);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% this is the general BCI loop where realtime incoming data is handled
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

blocksize  = 1;
prevSample = 0;

% simple smoothing kernel
kern1d = [1 2 1];
kern = convn(kern1d'*kern1d,reshape(kern1d,1,1,3))
kern = kern/sum(kern(:));

while true
  
  % determine number of samples available in buffer
  hdr = ft_read_header(cfg.headerfile, 'headerformat', cfg.headerformat, 'cache', true);
  
  if hdr.nSamples > prevSample
    % determine the samples to process
    if strcmp(cfg.bufferdata, 'last')
      begsample  = hdr.nSamples*hdr.nTrials - blocksize + 1;
      endsample  = hdr.nSamples*hdr.nTrials;
    elseif strcmp(cfg.bufferdata, 'first')
      begsample  = prevSample+1;
      endsample  = prevSample+blocksize;
    else
      ft_error('unsupported value for cfg.bufferdata');
    end
    
    prevSample  = endsample;
    fprintf('\nprocessing scan %d\n', begsample);
    
    % read data from buffer (only the last scan)
    dat = ft_read_data(cfg.datafile, 'header', hdr, 'dataformat', cfg.dataformat, 'begsample', begsample, 'endsample', endsample);
    
    % read events from buffer
    ev = ft_read_event(cfg.eventfile, 'header', hdr);
    if ~isempty(ev)
      ev_sam = [ev.sample];
      ind = find(ev_sam==begsample-1);
      if ~isempty(ind)
        disp([ev(ind).type ev(ind).value])
      end
    end
    
    % go to the correct figure
    figure(h);
    
    % actually, you should reshape to 3D (remove the *)
    S = reshape(dat, [width  height*numSlices]);
    
    V = reshape(dat, [width  height numSlices]);
    tic;
    Vn = convn(V,kern);
    toc
    Vn = Vn(2:end-1,2:end-1,2:end-1);
    size(Vn)
    S = reshape(Vn, [width height*numSlices]);
    
    imagesc(S(:,:)',[0 2048]);
    
    % force Matlab to update the figure
    drawnow
  else
    pause(0.01);
  end % if enough new samples
end % while true
