decode_qform0

 Decode qform info from NIFTI-1 headers.
 _______________________________________________________________________
 Copyright (C) 2005 Wellcome Department of Imaging Neuroscience

0001 function M = decode_qform0(hdr)
0002 % Decode qform info from NIFTI-1 headers.
0003 % _______________________________________________________________________
0004 % Copyright (C) 2005 Wellcome Department of Imaging Neuroscience
0005 
0006 %
0007 % $Id: decode_qform0.m 253 2005-10-13 15:31:34Z guillaume $
0008 
0009 
0010 dim    = double(hdr.dim);
0011 pixdim = double(hdr.pixdim);
0012 if ~isfield(hdr,'magic') || hdr.qform_code <= 0,
0013     flp = nic_spm_flip_analyze_images;
0014     %disp('------------------------------------------------------');
0015     %disp('The images are in a form whereby it is not possible to');
0016     %disp('tell the left and right sides of the brain apart.');
0017     %if flp,
0018     %    disp('They are assumed to be stored right-handed.');
0019     %else
0020     %    disp('They are assumed to be stored left-handed.');
0021     %end;
0022     %disp('------------------------------------------------------');
0023 
0024     %R     = eye(4);
0025     n      = min(dim(1),3);
0026     vox    = [pixdim(2:(n+1)) ones(1,3-n)];
0027 
0028     if ~isfield(hdr,'origin') || ~any(hdr.origin(1:3)),
0029        origin = (dim(2:4)+1)/2;
0030     else
0031         origin = double(hdr.origin(1:3));
0032     end;
0033     off     = -vox.*origin;
0034     M       = [vox(1) 0 0 off(1) ; 0 vox(2) 0 off(2) ; 0 0 vox(3) off(3) ; 0 0 0 1];
0035 
0036     % Stuff for default orientations
0037     if flp, M = diag([-1 1 1 1])*M; end;
0038 else
0039 
0040     % Rotations from quaternions
0041     R = Q2M(double([hdr.quatern_b hdr.quatern_c hdr.quatern_d]));
0042 
0043     % Translations
0044     T = [eye(4,3) double([hdr.qoffset_x hdr.qoffset_y hdr.qoffset_z 1]')];
0045 
0046     % Zooms.  Note that flips are derived from the first
0047     % element of pixdim, which is normally unused.
0048     n = min(dim(1),3);
0049     Z = [pixdim(2:(n+1)) ones(1,4-n)];
0050     Z(Z<0) = 1;
0051     if pixdim(1)<0, Z(3) = -Z(3); end;
0052     Z = diag(Z);
0053 
0054     M = T*R*Z;
0055 
0056     % Convert from first voxel at [1,1,1]
0057     % to first voxel at [0,0,0]
0058     M = M * [eye(4,3) [-1 -1 -1 1]'];
0059 end;
0060 return;
0061