I am trying to create an orthogonal coordinate system based on two "almost" perpendicular vectors, which are deduced from medical images. I have two vectors, for example:
that are almost perpendicular, since their inner product is equal to 5e-4.
Then I find their cross product to create my 3rd basis:
Even his third vector is not completely orthogonal to Z and Y, as their inner products are in the order of -15 and -16, but I guess that is almost zero. In order to use this set of vectors as my orthogonal basis for a local coordinate system, I assume they should be almost completely perpendicular. I first thought I can do this by rounding my vectors to less decimal figures, but did not help. I guess I need to find a way to alter my initial vectors a little to make them more perpendicular, but I don't know how to do that.
I would appreciate any suggestions.
No products are associated with this question.
If you have no preference about which component to modify, you could perform something as simple as
>> Z(end) = -Z(1:2)*Y(1:2).'/Y(end) Z = -1.0200 1.5300 -1.6302 >> Z*Y.' ans = 0
but I would definitely implement a test and generate an error/warning if the new 3rd component of Z is too far from the original ..
EDIT: we would have to check that Y(end) is non-zero, so here is another version that doesn't have this flaw:
>> Yn = Y/norm(Y) ; >> Zcor = Z - (Z*Yn.')*Yn ;
Your input is:
Z = [-1.02, 1.53, -1.63]; Y = [2.39, -1.39, -2.8];
Now a typical procedure for creating an orthonormal tripod is:
N1 = Z ./ norm(Z); N3 = NCross(Y, N1); N2 = NCross(Z, N2);
function [NV, LenV, V] = NCross(L1, L2) % Normalized vector orthogonal to L1 and L2 % [NV, LenV, V] = NCross(L1, L2) % INPUT: % L1, L2: [1 x 3] or [T x 3] double arrays. The size of the 1st % dimension of L1 and L2 must match or one must have 1 row. % L1 and L2 needn't be normalized. % OUTPUT: % NV: [T x 3] array, normalized cross product between rows of the % input matrices: % LenV: [T x 1] array, length of the cross products. % V: [T x 3] array, cross vector before normalization. % % MATHEMATICAL DEFINITION: % V = cross(L1, L2) % LenV = norm(V, 2) % NV = V / LenV % % If input vectors are almost parallel, the normalization will be % weak, therefore all vectors with length smaller than SQRT(EPS) are set to Inf. % % Tested: Matlab 6.5, 7.7, 7.8, 7.13, WinXP/32, Win7/64 % Author: Jan Simon, Heidelberg, (C) 2009-2013 % $License: BSD $
if ~(isequal(2, ndims(L1), ndims(L2)) && ... (isequal(size(L1), size(L2)) || ... size(L1, 1) == 1 || size(L2, 1) == 1)) error('JSimon:NCross:BadInputSize', ... 'NCross needs matching dimensions for L1 and L2!'); end
%smallVal = 100.0 * eps; smallVal = 1.4901161193847656e-008; % SQRT(EPS) infVal = Inf;
V = [ ( L1(:, 2) .* L2(:, 3) - L1(:, 3) .* L2(:, 2) ), ... ( L1(:, 3) .* L2(:, 1) - L1(:, 1) .* L2(:, 3) ), ... ( L1(:, 1) .* L2(:, 2) - L1(:, 2) .* L2(:, 1) ) ];
LenV = sqrt(V(:, 1).*V(:, 1) + V(:, 2).*V(:, 2) + V(:, 3).*V(:, 3));
lowI = (LenV < smallVal); LenV(lowI) = 1; % Prevent division by zero NV = V ./ LenV(:, [1, 1, 1]); % Or BSXFUN of course LenV(lowI) = 0; badI = (isfinite(LenV) == 0); LenV(badI) = infVal;
% Set weak vectors to [Inf, Inf, Inf]: weakI = (lowI | badI); NV(weakI, :) = infVal; V(weakI, :) = infVal;
Sorry that the function NCross is so long, but catching the exceptions is prone to mistakes. Creating a normalized tripod can be important for numerical reasons, because strange effects can appear, if the lengths of the vectors differ by a factor of 1e13.
Upto 1e-16 is as close to a precision beyond which it becomes questionable due to finite precision arithmetic on machines. Again there are several ways to generate orthogonal basis, a popular method is Gram-Schmidt Process for which there are couple of File Exchange submissions.
However PCA or PRINCOMP (principal component analysis) always gives an orthogonal basis for your data, albeit a special basis which you can get as follows:
[pc,score,latent,tsquare] = princomp(X);
The variable pc has normalized vectors which are all necessarily orthogonal to each other.