encode_qform0

 Encode an affine transform into qform
 _______________________________________________________________________
 Copyright (C) 2005 Wellcome Department of Imaging Neuroscience

0001 function hdr = encode_qform0(M,hdr)
0002 % Encode an affine transform into qform
0003 % _______________________________________________________________________
0004 % Copyright (C) 2005 Wellcome Department of Imaging Neuroscience
0005 
0006 %
0007 % $Id: encode_qform0.m 253 2005-10-13 15:31:34Z guillaume $
0008 
0009 
0010 % Convert from first voxel at [1,1,1] to first voxel at [0,0,0]
0011 M = M * [eye(4,3) [1 1 1 1]'];
0012 
0013 % Translations
0014 hdr.qoffset_x = M(1,4);
0015 hdr.qoffset_y = M(2,4);
0016 hdr.qoffset_z = M(3,4);
0017 
0018 % Rotations and zooms
0019 R         = M(1:3,1:3);
0020 vx        = sqrt(sum(M(1:3,1:3).^2));
0021 vx(vx==0) = 1;
0022 R         = R * diag(1./vx);
0023 
0024 % Ensure that R is O(3)
0025 [U,S,V] = svd(R);
0026 R       = U*V';
0027 if any(abs(diag(S)-1)>1e-3), warning('QFORM0 representation has been rounded.'); end;
0028 
0029 % Ensure that R is SO(3)
0030 if det(R)>0
0031     hdr.pixdim(1:4) = [ 1 vx];
0032 else
0033     R               = R*diag([1 1 -1]);
0034     hdr.pixdim(1:4) = [-1 vx];
0035 end;
0036 
0037 % Convert to quaternions
0038 Q             = M2Q(R);
0039 hdr.quatern_b = Q(1);
0040 hdr.quatern_c = Q(2);
0041 hdr.quatern_d = Q(3);
0042 
0043 if hdr.qform_code == 0, hdr.qform_code = 2; end;
0044 return;
0045