1 | function U = orthogonal(V, W); |
---|
2 | % ORTHOGONAL(V) constructs a unit pure quaternion perpendicular to V, and |
---|
3 | % perpendicular to W, if given. V (and W if given) must be pure quaternion |
---|
4 | % scalar(s). W need not be perpendicular to V, but it must not be parallel |
---|
5 | % (that is the scalar product of W and V must not have unit modulus). |
---|
6 | |
---|
7 | % Copyright © 2005 Stephen J. Sangwine and Nicolas Le Bihan. |
---|
8 | % See the file : Copyright.m for further details. |
---|
9 | |
---|
10 | % Modified 27 September 2005 by SJS to include comments on complexified |
---|
11 | % case. |
---|
12 | |
---|
13 | error(nargchk(1, 2, nargin)), error(nargoutchk(0, 1, nargout)) |
---|
14 | |
---|
15 | if size(V) ~= [1, 1] |
---|
16 | error('First argument must not be a vector or matrix.'); |
---|
17 | end |
---|
18 | |
---|
19 | if ~isa(V, 'quaternion') | ~ispure(V) |
---|
20 | error('First argument must be a pure quaternion.') |
---|
21 | end |
---|
22 | |
---|
23 | V = unit(V); % Ensure that V is a unit pure quaternion. |
---|
24 | |
---|
25 | if nargin == 2 |
---|
26 | if size(W) ~= [1, 1] |
---|
27 | error('Second argument must not be a vector or matrix.'); |
---|
28 | end |
---|
29 | |
---|
30 | if ~isa(W, 'quaternion') | ~ispure(W) |
---|
31 | error('Second argument must be a pure quaternion.') |
---|
32 | end |
---|
33 | |
---|
34 | % W must be distinct in direction from V, that is it must not be |
---|
35 | % parallel to V. We verify this by computing the dot (scalar) product |
---|
36 | % between V and unit(W), which will be +1 or -1 in the case where they |
---|
37 | % are parallel. |
---|
38 | % |
---|
39 | % Notice that, if the dot product is complex, the abs function ensures |
---|
40 | % a real result before 1 is subtracted. Thus 'parallel' in the case of |
---|
41 | % complex quaternions means having a dot product with a modulus of 1. |
---|
42 | % It is possible this restriction could be relaxed when it is better |
---|
43 | % understood, but since W is somewhat arbitrary, it is not difficult to |
---|
44 | % make a choice that satisfies the check below. |
---|
45 | |
---|
46 | if abs(abs(dot(V, unit(W))) - 1) <= eps |
---|
47 | error('Arguments must not be parallel to each other.'); |
---|
48 | end |
---|
49 | end |
---|
50 | |
---|
51 | % Now compute the required vector. The method used was published in: |
---|
52 | % |
---|
53 | % Ell, T. A. and Sangwine, S. J., "Decomposition of 2D Hypercomplex Fourier |
---|
54 | % Transforms into Pairs of Complex Fourier Transforms", in: Gabbouj, M. and |
---|
55 | % Kuosmanen (eds), "Signal Processing X, Theories and Applications", |
---|
56 | % Proceedings of EUSIPCO 2000, Tenth European Signal Processing Conference, |
---|
57 | % Tampere, Finland, 5-8 September 2000, Vol. II, 1061-1064, section 6. |
---|
58 | % |
---|
59 | % In the paper, a specific case was shown. Here we have a general method |
---|
60 | % for choosing one of the three basis vectors qi, qj, qk, making the choice |
---|
61 | % in such a way that if V is one of these vectors we will choose a |
---|
62 | % different one. |
---|
63 | |
---|
64 | % Note: This has not been thoroughly generalized to the complex case. If V |
---|
65 | % is a pure multivector, the result will be a pure multivector (that is a |
---|
66 | % complex pure quaternion), but if V is purely imaginary, the result will |
---|
67 | % be purely real. Since this code is used only if no value is given for W, |
---|
68 | % this is not a serious limitation, since the user is free to specify W and |
---|
69 | % thus obtain a complex result if needed. However, this issue merits |
---|
70 | % further study. SJS 27 September 2005. |
---|
71 | |
---|
72 | % Note added 3 October 2005. Perhaps we should test here not for equality, |
---|
73 | % but for parallelism, using the dot product. The code would read something |
---|
74 | % like: if dot(V, qi) == 1 | dot(V, qi) == i ? The problem is that if V is |
---|
75 | % complex, the dot product may be imaginary. What do we mean by parallel? |
---|
76 | % Whatever we mean, we want to avoid it and make a different choice. |
---|
77 | |
---|
78 | if nargin == 1 |
---|
79 | |
---|
80 | % W was not given, so we have to choose an arbitrary vector not |
---|
81 | % parallel to V. We make this choice in such a way that if V is |
---|
82 | % in the set {qi, qj, qk}, the result of the function is the |
---|
83 | % next element of the set (cyclically). This is not essential, |
---|
84 | % but since the result is arbitary, it is neater to make a |
---|
85 | % sensible choice, rather than an arbitrary one. |
---|
86 | |
---|
87 | if V == qi |
---|
88 | W = -qk; |
---|
89 | elseif V == qj |
---|
90 | W = -qi; |
---|
91 | else |
---|
92 | W = -qj; |
---|
93 | end; |
---|
94 | end |
---|
95 | |
---|
96 | % Now compute the result. We have to normalise this, because W may not be |
---|
97 | % perpendicular to V, and therefore the result may not have unit modulus. |
---|
98 | |
---|
99 | U = unit(cross(V, W)); |
---|