[37] | 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 | |
---|