How do I create orthogonal basis based on two "almost" perpendicular vectors?
13 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
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.
0 comentarios
Respuesta aceptada
Cedric
el 1 de Feb. de 2013
Editada: Cedric
el 1 de Feb. de 2013
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 ;
0 comentarios
Más respuestas (2)
Jan
el 1 de Feb. de 2013
Editada: Jan
el 1 de Feb. de 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.
2 comentarios
Cedric
el 1 de Feb. de 2013
Editada: Cedric
el 1 de Feb. de 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.
Shashank Prasanna
el 1 de Feb. de 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)
0 comentarios
Ver también
Categorías
Más información sobre Software Development Tools en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!