1 | function C = unadjoint(A, F) |
---|
2 | % UNADJOINT Given an adjoint matrix, constructs the original |
---|
3 | % quaternion matrix. |
---|
4 | % |
---|
5 | % unadjoint(A) or |
---|
6 | % unadjoint(A, 'complex') assumes A is a complex adjoint matrix. |
---|
7 | % unadjoint(A, 'real') assumes A is a real adjoint matrix. |
---|
8 | % |
---|
9 | % This function reverses the result of the ADJOINT function, so |
---|
10 | % that UNADJOINT(ADJOINT(A)) == A. |
---|
11 | |
---|
12 | % Copyright © 2005 Stephen J. Sangwine and Nicolas Le Bihan. |
---|
13 | % See the file : Copyright.m for further details. |
---|
14 | |
---|
15 | % References: see the code for ADJOINT. |
---|
16 | |
---|
17 | error(nargchk(1, 2, nargin)), error(nargoutchk(0, 1, nargout)) |
---|
18 | |
---|
19 | if nargin == 1 |
---|
20 | F = 'complex'; % Supply the default parameter value. |
---|
21 | end |
---|
22 | |
---|
23 | if ~strcmp(F, 'real') & ~strcmp(F, 'complex') |
---|
24 | error('Second parameter value must be ''real'' or ''complex''.') |
---|
25 | end |
---|
26 | |
---|
27 | % We have to be sure that A is a adjoint matrix, and not just some |
---|
28 | % arbitrary real or complex matrix. To do this we reconstruct the |
---|
29 | % original matrix, assuming A is an adjoint, then build the adjoint |
---|
30 | % from the reconstruction and compare it with A. If these are not |
---|
31 | % equal to within a tolerance, A was not an adjoint. We are also |
---|
32 | % able to eliminate some matrices as non-adjoints simply because of |
---|
33 | % their dimensions. When we build the adjoint C, we make an arbitrary |
---|
34 | % choice about which redundant parts of A to use. However, the unused |
---|
35 | % redundant parts are compared against the adjoint built from C and |
---|
36 | % if they are not close enough to the values used to build C an error |
---|
37 | % occurs. |
---|
38 | |
---|
39 | [r, c] = size(A); |
---|
40 | |
---|
41 | if strcmp(F, 'complex') |
---|
42 | |
---|
43 | r2 = r ./ 2; |
---|
44 | c2 = c ./ 2; |
---|
45 | |
---|
46 | if fix(r2) .* 2 ~= r | fix(c2) .* 2 ~= c |
---|
47 | error(['Matrix is incorrect size to be a ', F, ' adjoint']); |
---|
48 | end |
---|
49 | |
---|
50 | A1 = A(1 : r2, 1 : c2); |
---|
51 | A2 = A(1 : r2, c2 + 1 : c ); |
---|
52 | |
---|
53 | C = quaternion(real(A1), imag(A1), real(A2), imag(A2)); |
---|
54 | |
---|
55 | else % F must be 'real', since we checked it above. |
---|
56 | |
---|
57 | r4 = r ./ 4; |
---|
58 | c4 = c ./ 4; |
---|
59 | |
---|
60 | if fix(r4) .* 4 ~= r | fix(c4) .* 4 ~= c |
---|
61 | error(['Matrix is incorrect size to be a ', F, ' adjoint']); |
---|
62 | end |
---|
63 | |
---|
64 | |
---|
65 | W = A(1 : r4, 1 : c4); |
---|
66 | X = A(1 : r4, c4 + 1 : 2 .* c4); |
---|
67 | Y = A(1 : r4, 2 .* c4 + 1 : 3 .* c4); |
---|
68 | Z = A(1 : r4, 3 .* c4 + 1 : c ); |
---|
69 | |
---|
70 | C = quaternion(W, X, Y, Z); |
---|
71 | end |
---|
72 | |
---|
73 | T = adjoint(C, F); |
---|
74 | if any(size(T) ~= size(A)) |
---|
75 | error(['Matrix is incorrect size to be a ', F, ' adjoint']); |
---|
76 | end |
---|
77 | |
---|
78 | % Now we verify that A was indeed an adjoint matrix to within a tolerance. |
---|
79 | % We cannot test for exact equality, because A may be a product or the |
---|
80 | % result of some algorithm that theoretically yields an adjoint matrix, |
---|
81 | % but in practice yields an almost-adjoint because of rounding errors. |
---|
82 | |
---|
83 | RT = real(T); RA = real(A); |
---|
84 | IT = imag(T); IA = imag(A); % In the case of real quaternions, these will be 0. |
---|
85 | |
---|
86 | if any(any(abs(RT - RA) .* eps > abs(RA) | abs(IT - IA) .* eps > abs(IA))) |
---|
87 | error('Matrix is not (accurately) a valid adjoint matrix.'); |
---|
88 | end |
---|
89 | |
---|
90 | % Notes: the above test is not ideal. It is difficult to devise a simple |
---|
91 | % test that will work for complexified quaternions (because the modulus |
---|
92 | % can vanish). We can define a non-vanishing quasi-modulus and use that, |
---|
93 | % but it might not work correctly for all cases. For the moment the above |
---|
94 | % test will do. It will correctly raise an error in simple cases of bad |
---|
95 | % code. |
---|
96 | |
---|
97 | |
---|
98 | |
---|