neeman@s5.csrd.uiuc.edu (Henry J. Neeman) (03/20/91)
Here's an interesting problem that I need to solve. Any ideas?
The basic idea: given an object in d-space defined by d points, construct a
rotation matrix that rotates the object to be perpendicular to one of the
coordinate axes. The method should be fast, but I don't care about minimizing
floating point error.
Suppose we have a (d-1)-dimensional structure S in d-space (e.g., line in 2D;
plane in 3D; hyperplane in 4D), described by d points; i.e., S={Pi|i=1..d}.
Note that S is "flat" (e.g., a line segment/triangle/tetrahedron/hyperpentatope
in 2/3/4/5-space) rather than some curvy thing like a spline; also that we
are interested in the unbounded structure (e.g., the line/plane/hyperplane)
rather than the bounded structure that the points Pi specifically define.
We want to rotate S such that the normal to S, N, is parallel to a particular
coordinate axis; for example, in 3D we might want to rotate the plane S such
that N is parallel to the z-axis and S is parallel to the xy-plane.
Now in 2D and 3D this is pretty straightforward and can be done efficiently:
by finding the slope of S in 2D or the normal from cross products in 3D. But
as we move into higher dimensions, it becomes a more interesting problem.
Let me tell you how I'm currently solving this: Let Pi=(pi1, pi2, ..., pid).
We can obtain the d-1 vectors V1, V2, ..., V(d-1) defining the unbounded
structure as
Vi = Pi - Pd = [vi1=pi1-pd1 vi2=pi2-pd2 ... vid=pid-pdd].
From the Vi's, we can easily obtain the normal to the d-1 structure S, N, from
the solution to the homogeneous equation A*N=0:
[ v11 v12 ... v1d ] * [n1] = [0]
[ v21 v22 ... v2d ] [n2] [0]
A=[ . . . ] N=[. ] [0]
[ : : : ] [: ] [0]
[v(d-1)1 v(d-1)2 ... v(d-1)d] [nd] [0]
That is, N is orthogonal to all vectors Vi, i=1..d-1, so the dot product
N.Vi is zero for all i=1..d-1. (Actually, N is a family of such vectors,
but we can choose its length arbitrarily for now; later we'll normalize.)
This system of equations has a solution, because the rank of A is d-1 and
the rank of the space is d.
Currently, I solve this system by Gaussian elimination with back substitution,
using 1.0 as my initial constant on the family of solutions (i.e., set Nd=1.0,
then back substitute). *This method is O(d^3), with a large constant of
proportionality.* My first question is: is there a faster way to solve this
system that's computationally simple (e.g., no iteration)?
Anyway: now that we have a set of d-1 vectors in our d-1 structure S,
plus the dth vector, which is S's normal N, we want to rotate S such that
it is parallel to d-1 of the coordinate axes and orthogonal to the dth;
i.e., we want N to be parallel to one of the coordinate axes.
Here's how I'm currently doing this: first, I construct an orthonormal basis
for S and N with the same orientation as S. I have N orthogonal to S, so any
vector which is orthogonal to N must lie in (unbounded) S. Denoting the rows
of A as Ai, we have:
Let U1=N.
for i=1..d-1 {
Replace Ai by Ui.
Solve A*U(i+1)=0 for U(i+1), which is orthogonal to N and thus
lies in S, and is also orthogonal to all Uj, j <= i.
}
This algorithm produces a set of mutually orthogonal vectors Ui, i=1..d, such
that U1=N is orthogonal to S and Ui, i=2..d, lie in (unbounded) S. Now we
normalize: Bi=Ui/|Ui|=[bi1 bi2 ... bid], i=1..d. Then the set B={Bi|i=1..d}
is an orthonormal basis for d-space with the same orientation as S. *Thus,
the total time complexity is O(d^4), with a large constant of proportionality.*
From B, we can construct a rotation matrix R which will rotate the orthonormal
basis vectors Bi into the orthonormal vectors of the coordinate axes, except
with respect to origins, as follows: Let C={Ci|i=1..d} be the set of coordinate
axis vectors (e.g., in 3-space, C1, C2 and C3 correspond to i, j and k
respectively). We construct R such that the rows of R are the vectors in the
orthonormal basis:
R = [B1.C1 B1.C2 ... B1.Cd] = [b11 b12 ... b1d]
[B2.C1 B2.C2 ... B2.Cd] [b21 b22 ... b1d]
[ : : : ] [ . . . ]
[ : : : ] [ : : : ]
[Bd.C1 Bd.C2 ... Bd.Cd] [bd1 bd2 ... bdd]
Then R is the rotation matrix from B to C.
My question is: is there a more efficient way to obtain a rotation of S's
normal into one of the coordinate axes? For example, is there some way to
obtain such a rotation just from S and its normal? It seems like there should
be a better way, but I don't know enough linear algebra or geometry to be able
to find it. I'm not overly concerned about reducing error, since my task is
to rotate an object that's "very diagonal" to the coordinate axes so that it
becomes "mostly parallel" to the coordinate axes.
Thanks for your help.
Henry Neeman
neeman@csrd.uiuc.edu