How do I create orthogonal basis based on two "almost" perpendicular vectors?

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:

Z=[-1.02,1.53,-1.63]; Y=[2.39,-1.39,-2.8];

that are almost perpendicular, since their inner product is equal to 5e-4.

Then I find their cross product to create my 3rd basis:

X=cross(Y,Z);

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.

Answer 1

[Cedric]

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 ;

Answer 2

[Jan]

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:

1. Decide for one of the vectors to be kept exactly, called "1st" vector now
2. Create the third vector by the crossproduct
3. Rotate the second vector around the 3rd until it is orthogonal to the 1st one. This can be implemented as a cross product also.

N1 = Z ./ norm(Z);
N3 = NCross(Y, N1);
N2 = NCross(Z, N2);

with:

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.

[Doctor61]

Hi Jan, thanks. Did you mean

N2 = NCross(Z, N3);

?

[Cedric]

It is one of the following, as your statement seems to imply that you have to build X (or N3):

 N2 = NCross(Z, N3) ;
 N2 = NCross(N1, N3) ;

and the first operation ( N1 = Z ./ norm(Z)) is not necessary if you just need an orthogonal basis. It becomes necessary if you need an orthonormal basis though.

Answer 3

[Shashank Prasanna]

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.

pc(:,1)'*pc(:,2)