createvmr.m

toolbox of BVQX, This is the access between BV and matlab. It will help you to analysis data from BV

function vmr = createvmr(imafiles, flags)
% createvmr  - create a VMR from a set of IMA/DCM files
%
% FORMAT:       vmr = createvmr(imafiles [, flags])
%
% Input fields:
%
%       imafiles    FxC char or Fx1/1xF cell array with IMA filenames
%                   multiple-slice files must be sorted correctly
%       flags       optional settings
%        .flipped   override assume flipped flag of BVQXfile
%        .interp    method (if required) {'linear'}, 'cubic', 'spline'
%        .ismni     for Analyze, assumes the image to be in MNI space
%                   resulting in a full filled cube
%
% Output fields:
%
%       vmr         VMR object

% Version:  v0.7b
% Build:    7090508
% Date:     Sep-05 2007, 8:06 AM CEST
% Author:   Jochen Weber, Brain Innovation, B.V., Maastricht, NL
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 1 || ...
   (~ischar(imafiles) && ...
    ~iscell(imafiles))
    error( ...
        'BVQXtools:BadArguments', ...
        'Bad or missing argument.' ...
    );
end
if ischar(imafiles)
    imafiles = cellstr(imafiles);
end
if numel(imafiles{1}) < 5
    error( ...
        'BVQXtools:BadArguments', ...
        'Bad or missing argument.' ...
    );
end
if nargin < 2 || ...
    numel(flags) ~= 1 || ...
   ~isstruct(flags)
    flags = struct;
end
if ~isfield(flags, 'flipped') || ...
    isempty(flags.flipped) || ...
   (~isnumeric(flags.flipped) && ...
    ~islogical(flags.flipped))
    flags.flipped = BVQXfile(0, 'config', 'hdr', 'assumeflipped');
else
    flags.flipped = flags.flipped(1);
end
if ~isfield(flags, 'interp') || ...
   ~ischar(flags.interp) || ...
   ~any(strcmpi(flags.interp, {'cubic', 'linear', 'nearest', 'spline'}))
    flags.interp = 'linear';
else
    flags.interp = lower(flags.interp(:)');
end
if ~isfield(flags, 'ismni') || ...
    isempty(flags.ismni) || ...
   (~isnumeric(flags.ismni) && ...
    ~islogical(flags.ismni))
    flags.ismni = false;
else
    flags.ismni = false || flags.ismni(1);
end

% make sure files exist
for sc = 1:numel(imafiles)
    if ~ischar(imafiles{sc}) || ...
        exist(imafiles{sc}(:)', 'file') ~= 2
        error( ...
            'BVQXtools:FileNotFound', ...
            'File for slice %d not found.', ...
            sc ...
        );
    end
    imafiles{sc} = imafiles{sc}(:)';
end

% prepare array and set un-used file ID
intal = false;
vdt = uint16(0);
fid = 0;

% the next big scope decides on the type of file, reads in the
% file(s) and determines a few settings.
% after the scope, the following variables must be defined and
% set correctly to procede:
%
% DimX, DimY, DimZ   : dimensions of dataset
% gaps :             : gap size
% iop : double(4, 3) : [rowdir, coldir, slcdir, slcdircomp]
% pxs : double(1, 2) : pixel spacing within slice
% pxz :              : pixel spacing between slices
% sl1c / slnc        : slice 1 and N center coordinates (in world space)
% vdt                : data array with read voxel data

% what kind of files
cll = cell(1, 3);
if numel(imafiles) == 1
    
    % check extension
    exten = imafiles{1}(end-3:end);
    switch (lower(exten))
        
        % Analyze
        case {'.hdr', '.img', '.nii'}
            
            % for old analyze check HDR whether HDR exists
            gaps = 0;
            
            % read file
            afile = imafiles{1};
            if any(exten == upper(exten))
                hfile = [afile(1:end-4) '.HDR'];
            else
                hfile = [afile(1:end-4) '.hdr'];
            end
            if strcmpi(exten, '.img')
                afile = hfile;
            end
            if exist(afile, 'file') ~= 2
                error( ...
                    'BVQXtools:FileNotFound', ...
                    'Required Analyze file not found: %s.', ...
                    afile ...
                );
            end
            try
                oassfl = BVQXfile(0, 'config', 'hdr', 'assumeflipped');
                BVQXfile(0, 'config', 'hdr', 'assumeflipped', flags.flipped);
                afileo = BVQXfile(afile);
                cll{1} = afileo;
                aframe = afileo.CoordinateFrame;
                vdt = afileo.VoxelData(:, :, :);
                
                % get dims, resolution, and settings
                DimX = aframe.DimX;
                DimY = aframe.DimY;
                DimZ = aframe.DimZ;
                pxs = [aframe.ResX, aframe.ResY];
                pxz = aframe.ResZ;
                iop = [aframe.RowDir; aframe.ColDir; aframe.SlcDir];
                iop = iop ./ repmat(sqrt(sum(iop .* iop, 2)), [1, 3]);
                iop(4, :) = -cross(iop(1, :), iop(2, :));
                sl1c = aframe.Slice1Center;
                slnc = aframe.SliceNCenter;

            catch
                clearbvqxobjects(cll);
                BVQXfile(0, 'config', 'hdr', 'assumeflipped', oassfl);
                error( ...
                    'BVQXtools:ErrorReadingFile', ...
                    'Error reading specified Analyze file.' ...
                );
            end
            
            % exception : resampling is requested
            if flags.ismni
                
                % but check directions first (must be "clean")
                if abs(aframe.RowDir(1)) < 0.99 || ...
                    abs(aframe.ColDir(2)) < 0.99 || ...
                    abs(aframe.SlcDir(3)) < 0.99
                    flags.ismni = false;
                    warning( ...
                        'BVQXtools:BadOption', ...
                        'The given Analyze image seems not to be MNI.' ...
                    );
                end
            end
            
            % still requested
            if flags.ismni
                
                % get sampling matrix and sample
                smx = single(-127.5:127.5);
                try
                    [smy, smz, smx] = ndgrid(-smx, -smx, -smx);
                    smx = smx(:);
                    smy = smy(:);
                    smz = smz(:);
                    afileo = BVQXfile(afile);
                    afileo.VoxelData = vdt;
                    vds = single([]);
                    vds(256, 256, 256) = 0;
                    for vc = 1:1048576:numel(vds)
                        vcc = vc + 1048575;
                        vds(vc:vcc) = afileo.SampleCoords( ...
                            [smx(vc:vcc), smy(vc:vcc), smz(vc:vcc)], 1, ...
                            flags.interp, 'tal');
                    end
                    vds = reshape(vds, [256, 256, 256]);
                    vdt = vds;
                    vds = [];
                    clear vds;
                    intal = true;
                    BVQXfile(0, 'config', 'hdr', 'assumeflipped', oassfl);
                catch
                    clearbvqxobjects(cll);
                    BVQXfile(0, 'config', 'hdr', 'assumeflipped', oassfl);
                    error( ...
                        'BVQXtools:InternalError', ...
                        'Out of memory or sampling error: ''%s''.', ...
                        lasterr ...
                    );
                end
                
                % data limitting
                vdt = vdt - min(vdt(:));
                vdm = max(vdt(:));
                if vdm > 32767
                    vdt = uint16((32000 / vdm) .* vdt);
                else
                    vdt = uint16(vdt);
                end
                
                % re-setting values
                DimX = 256;
                DimY = 256;
                DimZ = 256;
                iop = [0, 1, 0; 0, 0, -1; 1, 0, 0; 1, 0, 0];
                pxs = [1, 1];
                pxz = 1;
                sl1c = [0, 0, -127.5];
                slnc = [0, 0,  127.5];
            end
            
        % PAR/REC
        case {'.par', '.rec'}
            try
                
                % read file
                p = readpar(imafiles{1});
                
                % get resolution
                DimX = p.SliceResolution(1);
                DimY = p.SliceResolution(2);
                DimZ = size(p.MatrixValues, 1);
                
                % get some settings
                try
                    angi = find(strcmp(p.MatrixHeaders, 'image_angulation_1'));
                    rotV = p.MatrixValues(1, [angi + 2, angi, angi + 1]);
                    rotX = [ ...
                        1, 0, 0; ...
                        0,  cos(rotV(1) * pi / 180), -sin(rotV(1) * pi / 180); ...
                        0,  sin(rotV(1) * pi / 180),  cos(rotV(1) * pi / 180)];
                    rotY = [ ...
                         cos(rotV(2) * pi / 180), 0,  sin(rotV(2) * pi / 180); ...
                        0, 1, 0;  ...
                        -sin(rotV(2) * pi / 180), 0,  cos(rotV(2) * pi / 180)]; ...
                    rotZ = [ ...
                         cos(rotV(3) * pi / 180), -sin(rotV(3) * pi / 180), 0; ...
                         sin(rotV(3) * pi / 180),  cos(rotV(3) * pi / 180), 0; ...
                         0, 0, 1]; ...
                    iop = (rotX * rotY * rotZ)';
                    iop(4, :) = -cross(iop(1, :), iop(2, :));
                    iop = iop ./ repmat(sqrt(sum(iop .* iop, 2)), [1, 3]);
                catch
                    iop = [0, 1, 0; 0, 0, -1; 1 0 0; 1 0 0];
                end
                
                % and some more settings
                resi = find(strcmp(p.MatrixHeaders, 'pixel_spacing_1'));
                if isempty(resi)
                    error('UNSUPPORTED_PARREC');
                end
                pxs  = p.MatrixValues(1, [resi, resi + 1]);
                slci = find(strcmp(p.MatrixHeaders, 'image_offcentre_1'));
                if isempty(slci)
                    error('UNSUPPORTED_PARREC');
                end
                sl1c = p.MatrixValues( 1 , [slci + 2, slci, slci + 1]);
                slnc = p.MatrixValues(end, [slci + 2, slci, slci + 1]);