function [voitc, weightv, weightr] = vtc_VOITimeCourse(hfile, voi, weight, fliplr)
% VTC::VOITimeCourse  - extract VOI time course data
%
% FORMAT:       voitc [, uvec, uvecr] = vtc.VOITimeCourse(voi [, weight, fliplr])
%
% Input fields:
%
%       voi         VOI file or coordinates (e.g. from VOI::BVCoords)
%       weight      (cell array of) Nx1 vector(s) with voxel weights
%                   give scalar 0 for unique, scalar 1 for no weighting,
%                   a scalar -1 for SVD after z-transform, or
%                   a scalar [Inf] to get a cell array of TxV arrays
%       fliplr      flip left/right (Z axes) for radiological convention
%
% Output fields:
%
%       voitc       TxV time course of voi(s)
%       uvec        unique VTC voxel indices within VOI, so that
%                   voi.VOI(i).Voxels(uvec{i}) leads to those coordinates
%       uvecr       reverse indexing (to find out which anatomical voxels
%                   fall into which functional voxels)

% Version:  v0.7c
% Build:    7100212
% Date:     Oct-02 2007, 12:43 PM CEST
% Author:   Jochen Weber, Brain Innovation, B.V., Maastricht, NL
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 2 || ...
    numel(hfile) ~= 1 || ...
   ~isBVQXfile(hfile, 'vtc') || ...
   (~all(isBVQXfile(voi(:), 'voi')) && ~isa(voi, 'double')) || ...
    isempty(voi)
    error( ...
        'BVQXfile:BadArguments', ...
        'Invalid call to %s.', ...
        mfilename ...
    );
end
bc = bvqxfile_getcont(hfile.L);
if numel(voi) == 1 && ...
    isBVQXfile(voi, 'voi')
    voic = bvqxfile_getcont(voi.L);
    if isempty(voic.VOI)
        error( ...
            'BVQXfile:BadArgument', ...
            'Invalid VOI object given.' ...
        );
    end
    numvois = numel(voic.VOI);
    vois = cell(1, numvois);
    for vc = 1:numvois
        vois{vc} = voi_BVCoords(voi, vc);
    end
else
    if ~any(size(voi)) == 3 || ...
        any(isinf(voi(:)) | isnan(voi(:)) | voi(:) < 0 | voi(:) > 255)
        error( ...
            'BVQXfile:BadArgument', ...
            'Invalid VOI coordinates given.' ...
        );
    end
    if size(voi, 1) == 3 && ...
        size(voi, 2) ~= 3
        voi = voi';
    end
    vois = {round(voi)};
end

% get VTC info
vres = bc.Resolution;
xstr = bc.XStart;
xend = bc.XEnd;
ystr = bc.YStart;
yend = bc.YEnd;
zstr = bc.ZStart;
zend = bc.ZEnd;
vtcsz = size(bc.VTCData);
numtp = vtcsz(1);
if vres * vtcsz(2) ~= (xend - xstr) || ...
    vres * vtcsz(3) ~= (yend - ystr) || ...
    vres * vtcsz(4) ~= (zend - zstr)
    error( ...
        'BVQXfile:BadObject', ...
        'Invalid Resolution/Start/End/Size combination.' ...
    );
end
voff = [xstr, ystr, zstr];

% check vois / weights
numvois = numel(vois);
if nargin < 3 || ...
    isempty(weight) || ...
   (~iscell(weight) && ~isa(weight, 'double'))
    weight = 1;
end
if ~iscell(weight)
    weight = {weight(:)};
end
if numel(weight) ~= numvois && ...
    numel(weight) == 1
    weight = weight(ones(1, numvois));
end
if numel(weight) ~= numvois
    error( ...
        'BVQXfile:BadArgument', ...
        'Number of weighting vectors mismatches VOIs.' ...
    );
end

% flipping
if nargin > 3 && ...
   (islogical(fliplr) || isa(fliplr, 'double')) && ...
   ~isempty(fliplr) && ...
    fliplr(1)
    fliplr = true;
else
    fliplr = false;
end

% make dimension check
cellout = false;
for vc = 1:numvois
    if numel(weight{vc}) == 1
        weight{vc} = repmat(weight{vc}, [size(vois{vc}, 1), 1]);
    end
    if any(isinf(weight{vc}))
        cellout = true;
    end
    if size(vois{vc}, 2) ~= 3 || ...
        size(weight{vc}, 2) ~= 1 || ...
        size(vois{vc}, 1) ~= size(weight{vc}, 1)
        error( ...
            'BVQXfile:InvalidArgument', ...
            'Invalid VOI/weight combination (dim mismatch).' ...
        );
    end

    % flip
    if fliplr
        vois{vc}(:, 3) = 256 - vois{vc}(:, 3);
    end

    % prepare coord lists
    vois{vc} = ...
        round(1 + (vois{vc} - repmat(voff, [size(vois{vc}, 1), 1])) / vres);

    % remove bad entries
    be = find( ...
        vois{vc}(:, 1) < 1 | vois{vc}(:, 1) > vtcsz(2) | ...
        vois{vc}(:, 2) < 1 | vois{vc}(:, 2) > vtcsz(3) | ...
        vois{vc}(:, 3) < 1 | vois{vc}(:, 3) > vtcsz(4));
    vois{vc}(be, :) = [];
    weight{vc}(be) = [];

end

% initialize voitc
weightv = {};
weightr = {};
if cellout
    voitc = {zeros(numtp, 1)};
    if numvois > 1
        voitc(2:numvois) = voitc(1);
    end
    weightv = cell(1, numvois);
    weightr = cell(1, numvois);
else
    voitc = zeros(numtp, numvois);
end

% iterate over vois
for vc = 1:numvois
    
    % extract voxel time courses
    voxs = vois{vc};
    dosvd = false;
    if all(weight{vc} == weight{vc}(1))
        if weight{vc} == 0
            doweight = false;
            voxs = unique(voxs, 'rows');
        elseif weight{vc}(1) == -1
            doweight = false;
            dosvd = true;
        elseif isinf(weight{vc}(1))
            doweight = false;
            dosvd = false;
        else
            doweight = true;
        end
    else
        doweight = true;
    end
    numvox = size(voxs, 1);
    vvtc = zeros(numtp, numvox);
    for xc = 1:numvox
        vox = voxs(xc, :);
        vvtc(:, xc) = bc.VTCData(:, vox(1), vox(2), vox(3));
    end
    
    % which method (weighting or unique)
    if doweight
        vvtc = vvtc .* repmat(weight{vc}', [numtp, 1]);
        voitc(:, vc) = sum(vvtc, 2) / sum(weight{vc});
    elseif dosvd
        [u{1:3}] = svd(ztrans(vvtc));
        u = u{1}(:, 1);
        u = u ./ std(u);
        cv = cov([mean(vvtc, 2), u]);
        if cv(1, 2) < 0
            u = -u;
        end
        voitc(:, vc) = u;
    else
        if cellout
            voitc{vc} = vvtc;
            [uvec{1:3}] = unique(voxs, 'rows');
            weightv{vc} = uvec{2};
            weightr{vc} = uvec{3};
        else
            voitc(:, vc) = mean(vvtc, 2);
        end
    end
end