/
spm_imatrix.m
48 lines (41 loc) · 1.61 KB
/
spm_imatrix.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
function P = spm_imatrix(M)
% Return the parameters for creating an affine transformation matrix
% FORMAT P = spm_imatrix(M)
% M - Affine transformation matrix
% P - Parameters (see spm_matrix for definitions)
%__________________________________________________________________________
%
% See also: spm_matrix.m
%__________________________________________________________________________
% Copyright (C) 1996-2011 Wellcome Trust Centre for Neuroimaging
% John Ashburner & Stefan Kiebel
% $Id: spm_imatrix.m 4414 2011-08-01 17:51:40Z guillaume $
%-Translations and Zooms
%--------------------------------------------------------------------------
R = M(1:3,1:3);
C = chol(R'*R);
P = [M(1:3,4)' 0 0 0 diag(C)' 0 0 0];
if det(R)<0, P(7)=-P(7); end % Fix for -ve determinants
%-Shears
%--------------------------------------------------------------------------
C = diag(diag(C))\C;
P(10:12) = C([4 7 8]);
R0 = spm_matrix([0 0 0 0 0 0 P(7:12)]);
R0 = R0(1:3,1:3);
R1 = R/R0;
%-This just leaves rotations in matrix R1
%--------------------------------------------------------------------------
%[ c5*c6, c5*s6, s5]
%[-s4*s5*c6-c4*s6, -s4*s5*s6+c4*c6, s4*c5]
%[-c4*s5*c6+s4*s6, -c4*s5*s6-s4*c6, c4*c5]
% There may be slight rounding errors making x>1 or x<-1.
rang = @(x) min(max(x, -1), 1);
P(5) = asin(rang(R1(1,3)));
if (abs(P(5))-pi/2)^2 < 1e-9
P(4) = 0;
P(6) = atan2(-rang(R1(2,1)), rang(-R1(3,1)/R1(1,3)));
else
c = cos(P(5));
P(4) = atan2(rang(R1(2,3)/c), rang(R1(3,3)/c));
P(6) = atan2(rang(R1(1,2)/c), rang(R1(1,1)/c));
end