We may start with orthogonal matrices, but when we begin to do operations on them small rounding errors can build up, causing the matrices to be slightly de-orthogonalised. Therefore, after many such operations, we may need a method to reorthonalise them.
where we doing a lot of calculations and errors can build up, for
An example of where reorthogonalising a matrix may be important is when we are recalculating at every frame in an animation, as described here.
Unfortunately I have not found an easy formula for doing this, here are the methods I am aware of:
- Recalculating Basis Vectors
- Deriving Correction Matrix
- Using SVD
None of these methods are ideal in that they either seem to have inaccuracies or they are very complicated, I would welcome any suggestions for better methods or improvements.
Recalculating Basis Vectors
The idea here is to remove the redundancy involved in a rotation matrix and then reconstruct the missing bits so that the matrix is perfectly orthogonal.
This is the simplest approach, but it is not ideal in that the redundant information might have helped us cancel out the rounding errors, instead this approach may degrade the accuracy of the result but the matrix will be orthogonal.
If we consider a rotation matrix as a set of basis vectors (as discussed on this page) then we need to make sure these basis vectors are unit length, divide by √(x² + y² + z²), and also mutually perpendicular. We could do that by deleting the third basis vector and recalculating it from the cross product of the other two. So far we have not made sure that the first two basis vectors are perpendicular with each other, so we the need to delete on of them and then recalculate from the other two.
Another approach could be to convert the matrix to a quaternion and then convert back to a matrix.
Either way its messy and we could be basing our results on the least accurate elements of the matrix.
Deriving Correction Matrix
The following derivation is evolved from this discussion Colin Plumb on newsgroup: sci.math
We want a correction matrix [C] which when multiplied with the original matrix [M] to give an orthogonal matrix [O] such that:
[O] = [C]*[M]
and since [O] is orthogonal then,
[O]-1 = [O]T
combining these equations gives:
[C]*[M] * ([C]*[M])T = [I]
where [I] = identity matrix
as explained here, ([A] * [B])T = [B]T * [A]T so substituting for [C] and [M] gives,
[C]*[M] * [M]T*[C]T = [I]
post multiplying both sides by ([C]T)-1 gives:
[C]*[M] * [M]T = ([C]T)-1
pre multiplying both sides by [C]-1 gives:
[M] * [M]T = [C]-1 ([C]T)-1
If we constrain [C] to be symmetric then [C]T = [C]
[M] * [M]T = [C]-2
[C] = ([M] * [M]T)-(1/2)
So using the taylor series here and using the first two terms we get:
x^(-1/2) = (3-x)/2
So substituting in above gives:
[C] = ([M] * [M]T)-(1/2)
= (3-([M] * [M]T))/2
since we can't subtract a matrix from a scalar I think it should really be written:
[C] = (3[I]-([M] * [M]T))/2
At least this seems to give appropriate results in the examples below:
There are some things about this derivation that make me uneasy, for instance constraining [C] to be symmetric seems to come out of nowhere. Vassilis has raised some issues with this method here, can anyone help me resolve them?
SVD Method
Hauke kindly suggested this method.
It is based on factoring the matrix into 3 parts:
[M] = [U][D][V]t
where:
- [M] = a matrix to be factored
- [U] = a matrix made up of a set of orthogonal 'output' vectors.
- [V] = a matrix made up of a set of orthogonal 'input' vectors.
- [D] = a diagonal matrix (non-diagonal terms are zero)
More information about SVD on this page.
This has separated the matrix our into a scaling part, the diagonal matrix [D] and the rotational part [U] and [V].
We want no (unity) scaling so it is tempting to set all the terms in the diagonal matrix to one:
[D] = [I] = identity matrix
which would give:
[M] = [U][V]t
However, because of the way that the SVD algorithm calculates [D], [U] and [V] it is possible that the terms of [D] might be -1 as well as +1. We therefore need to set the terms of [D] to +1 if the term was positive and -1 if the term was negative.
Alternatively we could recalculate the diagonal terms to give a +1 determinant.
D=diag([1 1 det(U*V^T)]).
Which means that our orthogonal matrix is:
[M] = U*diag([1 1 det(U*V^T])*V^T.
Examples
Orthogonal matrix
First test that this will preserve the orthogonal matrix, applying:
= (3[I]-([M] * [M]T))/2
Since the matrix is orthogonal then [M] * [M]T = [I] so,
| [C] = (3[I]-( |
|
))/2 |
| [C] = |
|
| [C] = |
|
so as required this wont alter the identity matrix as required:
Non-orthogonal matrix
now try this on a matrix that is not already orthogonal because it has a term 'error' added to m01, applying:
[C]= (3-([M] * [M]T))/2
| [C]= (3-( |
|
|
))/2 |
| [C]= |
|
| [C]= |
|
so if we multiply this be the original matrix we get:
| [O]= |
|
|
| [O]= |
|
| [O]= |
|
Is this more orthogonal than the original matrix, I can't really tell, it looks as though it might be?
