MATLAB Central Newsreader - Directional Cosine Matrix

Directional Cosine Matrix

Chia Khai Low — Thu, 19 Jun 2008 13:49:01 -0400

I have two vectors in different excel files. I have
successfully imported them and matched them by date. These
3d vectors are in cartesian coordinates. I want to find the
directional cosine matrix to convert the first vector into
the second vector.

Anyone with the method please post it here or show me links
on the math behind it. Once I have this DCMs(for every
instance) I can make Euler angles, visualizations and stuff.

Thanks and bless you.

ps. reply quickly :)

Re: Directional Cosine Matrix

Roger Stafford — Fri, 20 Jun 2008 06:54:02 -0400

"Chia Khai Low" <kaiserwulf@gmail.com> wrote in message <g3do4d$667
$1@fred.mathworks.com>...
> I have two vectors in different excel files. I have ........

  It is not entirely evident from your description what you asking, Chia. For
example, the statement, "These 3d vectors are in cartesian coordinates" is not
completely clear to me. I am guessing what you mean is that you have two n
by 3 arrays, A and B, in which the corresponding rows in them are x,y,z
(Cartesian) coordinates in two different coordinate systems, one presumably
rotated with respect to the other about a common origin. It is this 3 x 3
rotation matrix that you are apparently seeking as your "directional cosine
matrix". Have I guessed correctly?

  You didn't mention a possible translation in aligning the two sets of points.
If that were necessary, you could subtract the mean values of each so that the
sum of the coordinates in A and B would then be zero. The necessary
translation would then be the vector difference between these two mean value
vectors. I have therefore assumed below that the sums of each coordinate
column in A and B are zero: sum(A,1) = sum(B,1) = [0,0,0].

  I do not regard the method I am about to give as entirely satisfactory. It has
two fundamental weaknesses. The first weakness is that if A and B are not
accurate rotations of one another, the method may not give a perfect least
squares optimum rotation matrix for them. In fact, if they vary too far from
being rotations of one another, the method could conceivably fail altogether.

  The second weakness has to do with the eigenvalues (singular values in this
case) that are computed below. If it were to happen that two or three of them
were essentially equal, the method again falls apart.

  If we (temporarily) ignore these difficulties, do this:

[U1,S1,V1] = svd(A'*A,0); % Find eigenvalues and eigenvectors
[U2,S2,V2] = svd(B'*B,0); % for A'*A and B'*B
mn = inf;
II = [-1,-1,-1;-1,-1,1;-1,1,-1;-1,1,1;1,-1,-1;1,-1,1;1,1,-1;1,1,1];
for k = 1:8
  Rt = V1*diag(II(k,:))*V2'; % One of these eight
  X = B-A*Rt; % possibilities should be the right one
  tt = norm(X(:)); % Choose the closest fit to find it
  if tt < mn, mn = tt; R = Rt; end
end

Then R is the desired rotation matrix from the n points in A to those in B. We
should have essential equality, B = A*R, provided the points in B are actually
rotated versions of those in A.

  The reasoning here is that the eigenvalues of A'*A and B'*B should
(hopefully) be the same, and then the matrix V1*V2', with one of eight
possible combinations of sign changes in the columns of V1, should be the
required transformation from A to B. We have only to search through the
eight possibilities to find it.

  As stated above, if two or three of the eigenvalues in S1 or S2 should
happen to be equal, the logic in the previous paragraph fails. There would
have to be a search through infinitely many possible matrices between V1 and
V2', which of course is impossible. The method would break down in this
case.

  I should mention here that if some kind of reflection has occurred between
the points of A and those of B, this method will come up with an R matrix
whose determinant is -1, that is, it is a reflection, not a rotation, matrix.

  I am sure there exists a superior method of solving your problem that would
avoid the two aforementioned weaknesses, but I have not been able, up to
this time, to come up with a practical one. If I manage to think of a better
method, I will inform you, or perhaps some other kind soul will present one
here.

Roger Stafford

Re: Directional Cosine Matrix

Roger Stafford — Fri, 20 Jun 2008 15:55:05 -0400

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in
message <g3fk6a$co7$1@fred.mathworks.com>...
> It is not entirely evident from your description what you asking, Chia.
> ........

  A good night's sleep has furnished a much superior method of solving your
problem, Chia. This method has none of the disadvantages of the one I sent
you yesterday and is much simpler. There is no need to run through the
previous eight possibilities.

  Assume, as before, that A and B are n x 3 arrays in which all column sums
are zero. Then do this:

[U,S,V] = svd(A'*B);
R = U*V';

This R will be your desired unitary (direction cosine) matrix that, if applied to
A, will bring you to B: B = A*R.

  If B differs from such a transformation of A, then this solution yields an R
which I have reason to believe, (but can't prove,) is the best one in a least
squares sense.

  In the case that two or all three of the singular values in S are equal, there
will be an inherent indeterminacy, but this time it does not interfere with
finding a valid solution. In such cases there would be an infinite continuum
of equally good solutions.

  Again I remind you that the above solution, R, can in some cases have a
determinant of -1 instead of +1 and therefore be a reflection, rather than
rotation, matrix. In these cases the reflection matrix would provide a closer
fit than a true rotation matrix. It would nevertheless be unitary, that is, a
matrix of direction cosines.

Roger Stafford

Re: Directional Cosine Matrix

Chia Khai Low — Fri, 20 Jun 2008 21:47:02 -0400

Dear Roger,

Thank you so much for your comprehensive and yet quick
reply. And the follow-up too.

Have a great day.

Thanks again

Chia

Re: Directional Cosine Matrix

Roger Stafford — Mon, 23 Jun 2008 01:31:02 -0400

"Chia Khai Low" <kaiserwulf@gmail.com> wrote in message <g3h8gm$jd4
$1@fred.mathworks.com>...
> Thank you so much for your comprehensive and yet quick .........
> Chia

Hello Chia,

This is to let you know that I am now able to prove rigorously that the solution
I gave you in the second of my articles in this thread does in fact provide the
best orthogonal (rotation/reflection) matrix R in the least squares sense between
one set of points and another, which means that you can now use this method
with confidence. The proof involves the Schwartz-Cauchy inequality. As you
may recall, I was unable to come up with a proof earlier.

I hope for my peace of mind you will completely abandon the method I
proposed in the first article involving a search through eight possibilities. That
is not the right way to do the problem.

Roger Stafford