vmr_writeanalyzevol.m

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

function rvalue = vmr_WriteAnalyzeVol(hfile, filename)
% VMR::WriteAnalyzeVol  - write an Analyze image from the VMR volume
%
% FORMAT:       [success = ] vmr.WriteAnalyzeVol(filename)
%
% Input fields:
%
%       filename    volume filename
%
% Output fields:
%
%       success     true if successful

% Version:  v0.7b
% Build:    7090216
% Date:     Sep-02 2007, 4:58 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, 'vmr') || ...
   ~ischar(filename) || ...
    isempty(filename) || ...
    numel(filename) ~= length(filename) || ...
    length(filename) < 5 || ...
  (~strcmpi(filename(end-3:end), '.hdr') && ...
   ~strcmpi(filename(end-3:end), '.img'))
    error( ...
        'BVQXfile:BadArgument', ...
        'Invalid call to %s.', ...
        mfilename ...
    );
end
sc = bvqxfile_getcont(hfile.L);
bc = sc.C;

% default rvalue
rvalue = false;

% create struct version
[opath{1:2}] = fileparts(sc.F);
ofile = opath{2};

% filename
filename = filename(:)';
filename = [filename(1:end-3) 'img'];
hdrfname = [filename(1:end-3) 'hdr'];
matfname = [filename(1:end-3) 'mat'];

% dimension and datatype
xpix = bc.DimX;
ypix = bc.DimY;
zpix = bc.DimZ;
isiz = [xpix, ypix, zpix];

% X/Y/Z resolution
try
    xres = bc.FoVRows / bc.NRows;
    yres = bc.FoVCols / bc.NCols;
    zres = bc.SliceThickness + bc.GapThickness;
catch
    xres = 1;
    yres = 1;
    zres = 1;
end

% positional information
try
    xvec = xres .* [bc.RowDirX; bc.RowDirY; bc.RowDirZ];
    yvec = yres .* [bc.ColDirX; bc.ColDirY; bc.ColDirZ];
    ovec = [bc.Slice1CenterX; bc.Slice1CenterY; bc.Slice1CenterZ];
    lvec = [bc.SliceNCenterX; bc.SliceNCenterY; bc.SliceNCenterZ];
    zspan = bc.DimZ - 1;
    zvec = (1 / zspan) .* (lvec - ovec);
    if numel(bc.Trf) > 0
        znrm = 1 / sqrt(zvec(:)' * zvec(:));
        ovec(3) = ovec(3) * znrm;
        zvec = zvec .* znrm;
    end
    svec = ovec - ((xpix + 1) / 2) * xvec - ((ypix + 1) / 2) * yvec;
catch
    xvec = [0; 1; 0];
    yvec = [0; 0;-1];
    zvec = [1; 0; 0];
    svec = [-128; -128; 128];
end
% build mat
tmat = [[xvec, yvec, zvec], svec;  0, 0, 0, 1];

% older transformations, then discard info?
if isfield(bc, 'Trf')
    for tc = 1:length(bc.Trf)
        tt = bc.Trf(tc);
        if tt.NrOfSpatialTransformationValues && ...
           all(size(tt.TransformationValues) == 4)
            tmat = [0, 0, 1, -128; 1, 0, 0, -128; 0, -1, 0, 128; 0, 0, 0,1];
            break;
        end
    end
end

% radiological convention ?
if (isfield(bc, 'Convention') && ...
    bc.Convention(1)) || ...
   (isfield(bc, 'RadiologicalConvention') && ...
    strcmpi(bc.RadiologicalConvention(:)', 'yes'))
    tfmat = eye(4);
    tfmat(2, 2) = -1;
    tmat = tfmat * tmat;
else
    if ~isfield(bc, 'Convention') && ...
       ~isfield(bc, 'RadiologicalConvention')
        tfmat = eye(4);
        tfmat(2, 2) = -1;
        tmat = tfmat * tmat;
    end
end

% try volume creation
hdrc = cell(1, 1);
try
    hdr = BVQXfile('new:hdr');
    hdrc{1} = hdr;
    hdrbc = bvqxfile_getcont(hdr.L);
catch
    clearbvqxobjects(hdrc);
    error( ...
        'BVQXfile:InternalError', ...
        'Error creating Analyze header object: %s.', ...
        lasterr ...
    );
end

% set dims and data
hdrbc.ImgDim.Dim(1:4) = [3, isiz];
hdrbc.ImgDim.PixSpacing(2:4) = [xres, yres, zres];
if bc.VMR8bit && ...
    isempty(bc.VMRData16)
    hdrbc.ImgDim.DataType = 2;
    hdrbc.ImgDim.BitsPerPixel = 8;
    hdrbc.ImgDim.CalMaxDisplay = 255;
    vxd = bc.VMRData(:, :, :);
    otype = 'VMR';
else
    hdrbc.ImgDim.DataType = 132;
    hdrbc.ImgDim.BitsPerPixel = 16;
    hdrbc.ImgDim.CalMaxDisplay = 32767;
    if bc.VMR8bit
        vxd = bc.VMRData16(:, :, :);
    else
        vxd = bc.VMRData(:, :, :);
    end
    otype = 'V16';
end
hdrbc.ImgDim.CalMinDisplay = 0;
hdrbc.ImgDim.GLMax = hdrbc.ImgDim.CalMaxDisplay;
hdrbc.ImgDim.GLMin = 0;
hdrbc.DataHist.Description = sprintf('%s Volume of %s', otype, ofile);
hdrbc.VoxelData = vxd;
bvqxfile_setcont(hdr.L, hdrbc);

% save hdr/img
try
    aft_SaveAs(hdr, hdrfname);
    rvalue = hdr_SaveVoxelData(hdr);
catch
    warning( ...
        'BVQXfile:InternalError', ...
        'Error writing Analyze header/image %s.', ...
        hdrfname ...
    );
    return;
end
if ~rvalue
    return;
end

% save mat
eval('M=tmat;mat=tmat;save(matfname,''M'',''mat'',''-v6'');', '')
if exist(matfname, 'file') ~= 2
    rvalue = false;
end