function ft_omri_quality(cfg)

% FT_OMRI_QUALITY implements an online fMRI quality assurance stack

%

% Use as

% ft_omri_quality(cfg)

% where cfg is a structure with configuration settings.

%

% Configuration options are

% cfg.input = FieldTrip buffer containing raw scans (default='buffer://localhost:1972')

% cfg.numDummy = how many scans to ignore initially (default=0)

% cfg.showRawVariation = 1 to show variation in raw scans (default), 0 to show var. in processed scans

% cfg.clipVar = threshold to clip variation plot with as a fraction of signal magnitude (default=0.2)

% cfg.lambda = forgetting factor for the variaton plot (default=0.9)

% cfg.serial = serial port (default = /dev/ttyS0), set [] to disable motion reporting

% cfg.baudrate = serial port baudrate (default = 19200)

% cfg.maxAbs = threshold (mm) for absolute motion before 'A' is sent to serial port, default = Inf

% cfg.maxRel = threshold (mm) for relative motion before 'B' is sent to serial port, default = Inf

% Copyright (C) 2010, 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$

ft_defaults

ft_hastoolbox('spm8', 1);

if nargin<1

cfg = [];

end

if ~isfield(cfg, 'showRawVariation')

cfg.showRawVariation = 1;

end

if ~isfield(cfg, 'clipVar')

cfg.clipVar = 0.2;

end

if ~isfield(cfg, 'lambda')

cfg.lambda = 0.9;

end

if ~isfield(cfg, 'input')

cfg.input = 'buffer://localhost:1972';

end

if ~isfield(cfg, 'numDummy')

cfg.numDummy = 0;

end

if ~isfield(cfg, 'baudrate')

cfg.baudrate = 19200;

end

if ~isfield(cfg, 'serial')

cfg.serial = '/dev/ttyS0';

end

if ~isfield(cfg, 'maxAbs')

cfg.maxAbs = inf;

end

if ~isfield(cfg, 'maxRel')

cfg.maxRel = inf;

end

if ~isempty(cfg.serial)

try

serPort = serial(cfg.serial);

existPort = instrfind('Name',serPort.name, 'Status', 'open');

if ~isempty(existPort)

set(existPort, 'BaudRate', cfg.baudrate);

serPort = existPort;

else

set(serPort, 'BaudRate', cfg.baudrate);

fopen(serPort);

end

catch

% the "catch me" syntax is broken on MATLAB74, this fixes it

me = lasterror;

serPort = [];

disp(me.message);

end

else

serPort = [];

end

% Loop this forever (until user cancels)

while 1

clear ft_read_header

% start by reading the header from the realtime buffer

while 1

try

hdr = ft_read_header(cfg.input);

break;

catch

disp('Waiting for header');

pause(0.5);

end

end

% Ok, we got the header, try to make sense out of it

S = ft_omri_info_from_header(hdr);

if isempty(S)

ft_warning('No protocol information found!')

% restart loop

pause(0.5);

continue;

end

% reset motion estimates

motEst = [];

maxVal = 0;

maxDiff = 1e-6; % zero is not possible here (for imagesc range)

% Wait for numDummy scans (and drop them)

fprintf(1,'Waiting for %i dummy samples to come in...\n', cfg.numDummy);

while 1

threshold = struct('nsamples', cfg.numDummy);

newNum = ft_poll_buffer(cfg.input, threshold, 500);

if newNum.nsamples >= cfg.numDummy

break

end

pause(0.01);

end

fprintf(1,'Starting to process\n');

numTotal = cfg.numDummy;

numProper = 0;

discNum = 0;

% Loop this as long as the experiment runs with the same protocol (= data keeps coming in)

while 1

% determine number of samples available in buffer / wait for more than numTotal

threshold.nsamples = numTotal;

newNum = ft_poll_buffer(cfg.input, threshold, 500);

if newNum.nsamples < numTotal

% scanning seems to have stopped - re-read header to continue with next trial

break;

end

if newNum.nsamples == numTotal

% timeout -- go back to start of (inner) loop

drawnow;

continue;

end

% this is necessary for ft_read_data

hdr.nSamples = newNum.nsamples;

index = (cfg.numDummy + numProper) + 1;

fprintf('\nTrying to read %i. proper scan at sample index %d\n', numProper+1, index);

GrabSampleT = tic;

try

% read data from buffer (only the last scan)

dat = ft_read_data(cfg.input, 'header', hdr, 'begsample', index, 'endsample', index);

catch

ft_warning('Problems reading data - going back to poll operation...');

continue;

end

numProper = numProper + 1;

numTotal = numTotal + 1;

maxdat = max(dat);

if maxdat > maxVal

maxVal = maxdat;

maxDiff = cfg.clipVar * double(maxVal);

end

rawScan = single(reshape(dat, S.voxels));

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

% motion correction

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

if numProper == 1

% set up reference model for registration

flags = struct('mat', S.mat0);

fprintf(1,'Setting up first non-dummy scan as reference volume...\n');

RRM = ft_omri_align_init(rawScan, flags);

motEst = zeros(1,6);

procScan = single(rawScan);

lastPos = zeros(1,3);

else

fprintf('%-30s','Registration...');

tic;

[RRM, M, Mabs, procScan] = ft_omri_align_scan(RRM, rawScan);

toc

motPars = spm_imatrix(M);

motEst = [motEst; motPars(1:6).*[1 1 1 180/pi 180/pi 180/pi]];

if ~isempty(serPort)

curPos = motPars(1:3);

if any(abs(curPos) > cfg.maxAbs)

try

fprintf(serPort, 'A');

catch

% the "catch me" syntax is broken on MATLAB74, this fixes it

me = lasterror;

disp(me.message);

end

fprintf(1, 'A - too much absolute motion');

end

if any(abs(curPos - lastPos) > cfg.maxRel)

try

fprintf(serPort, 'B');

catch

% the "catch me" syntax is broken on MATLAB74, this fixes it

me = lasterror;

disp(me.message);

end

fprintf(1, 'B - too much relative motion');

end

lastPos = curPos;

end

end

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

% difference with previous image, after motion correction

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

if numProper == 1

diffScan = zeros(size(procScan));

if cfg.showRawVariation

lastScan = rawScan;

else

lastScan = procScan;

end

else

if cfg.showRawVariation

% Calculation based on raw scans

diffScan = max(cfg.lambda*diffScan, abs(rawScan - lastScan));

lastScan = rawScan;

else

% Calculation based on the processed scans

diffScan = max(cfg.lambda*diffScan, abs(procScan - lastScan));

lastScan = procScan;

end

discNum = cfg.lambda*discNum + 1;

end

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

% done with pre-processing, write output

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

fprintf('Done -- total time = %f\n', toc(GrabSampleT));

clippedDiff = max(diffScan, maxDiff);

% NOTE: the following line plots "rawScan" as the current volume

% If you'd like to see the processed scan (motion corrected)

% instead, exchange "rawScan" by "procScan"

ft_omri_quality_plot(motEst, rawScan, diffScan, maxVal, maxDiff);

% force Matlab to update the figure

drawnow

end % while true

end