Asked by Doctor61
on 1 Feb 2013

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.

*No products are associated with this question.*

Answer by Cedric Wannaz
on 1 Feb 2013

Edited by Cedric Wannaz
on 1 Feb 2013

Accepted answer

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 by Jan Simon
on 1 Feb 2013

Edited by Jan Simon
on 1 Feb 2013

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:

- Decide for one of the vectors to be kept exactly, called "1st" vector now
- Create the third vector by the crossproduct
- 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.

Cedric Wannaz
on 1 Feb 2013

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 by Shashank Prasanna
on 1 Feb 2013

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)

Opportunities for recent engineering grads.

## 0 Comments