Orthonormalize three vectors

There are two alternatives for rotational transformations: quaternions or matrices. I like matrices. But multiplying them many times produces numerical error which makes their vectors flow away from unit length and orthogonality. Therefore orthonormalization has to be performed eventually after n transformations.

The easy way is this:

mat3 RotMtx;
vec3 i = RotMtx[0];
vec3 j = RotMtx[1];
vec3 k = RotMtx[2];
vec3 Ox = normalize(i);
vec3 Oz = normalize(cross(i,j));
vec3 Oy = cross(Oz,Ox);
RotMtx[0] = Ox;
RotMtx[1] = Oy;
RotMtx[2] = Oz;

But this is not the ideal solution, because one axis’ direction is preserved (i), another changed (k->Oz), and the last one (Oy) does not really depend on it’s original direction at all - it may even form a right-hand vector triplet out of left-hand one.

The perfect solution, as I think, must be equal to all three input vectors: the direction of each one must change by the same angle to form an orthogonal triplet, and it shouldn’t be an iterative algorithm. But I just can’t find a solution!

So simply saying, the task is the following: we have three non-collinear vectors, a,b,c (let’s say, they are of unit length); we need to make them orthogonal by rotating each of them by the equal angle.

Here is what we have:
|a|=|b|=|c|=1;
a,b,c are not collinear
a -> A
b -> B
c -> C

Here is what we need:
dot(A,B)=0; (1)
dot(A,C)=0; (2)
dot(B,C)=0; (3)
|A|2=1; (4)
|B|2=1; (5)
|C|2=1; (6)
dot(a,A)=dot(b,B); (7)
dot(b,B)=dot(c,C); (8)

As you can see, I came up with 8 equations, but there are three unknown 3-dimensional vectors, which means 3*3=9 equations required. What do I miss?
And maybe there are some geometrical solution exists? Any tips, people?..:confused:

Gram–Schmidt orthonormalization.

Or just flush the whole thing and use quaternions or dual quaternions for transform composition.

But isn’t the Gram–Schmidt orthonormalization unfair to the change of input vectors’ directions just like in an “easy way” I show above?

Actually, I just developed some alternative solution. It is still rotates vectors by unequal angles to orthonormalize them, but the difference in those angles is not very significant. At least it gives the same triplet no matter which order the input vectors are given, so it’s kinda fair solution.

The idea is that the average direction of all three vectors

vec3 n=normalize(a+b+c));

should remain constant after a,b,c become orthonormal. If so, then the centroid of the triangle abc (which is based on those input vectors) will also become the circumcenter of the triangle ABC. And if A,B,C - are desired orthonormal vectors, then the ABC triangle will be equilateral with it’s circumcenter (point R) distanced from origin by “1/sqrt(3)”. So we find coordinates of this point first:

vec3 R = normalize(a+b+c)*0.57735026918962576450914878050196;

So we know the point laying at the desired triangle’ plane (the point is R) and we know that vector n is perpendicular to the triangle. Let’s scale the input vectors to make them touch the triangle’ plane:

a = a * ( dot(R,n) / dot(a,n) );
b = b * ( dot(R,n) / dot(b,n) );
c = c * ( dot(R,n) / dot(c,n) );

As the desired ABC triangle is based on orthonormal vectors of unit length, then we know that A,B,C should be at the equal distance from it’s circumcenter R, and this distance must be “sqrt(2/3)”. In other words, the projections of the desired orthonormal vectors onto the ABC plane must be equal to “sqrt(2/3)”. But the intersection points we just found are not necessarily at this distance right now, so we must bring them there:

a = R+normalize(axABC-R)*0.81649658092772603273242802490196;
b = R+normalize(bxABC-R)*0.81649658092772603273242802490196;
c = R+normalize(cxABC-R)*0.81649658092772603273242802490196;

At this point we’ve got three vectors of unit length (a,b,c), all making equal angle with the vector directed from origin to circumcenter of the triangle ABC (vector R). The last thing we have to correct is an angles between those vectors: to make them orthonormal we just need to ensure that angles between their projections onto ABC plane are 120 degrees - if so, than triangle abc will become equilateral (abc -> ABC). The rotation axis is obviously R - we can spin our a,b,c vectors around it and they will not change their length.
So now we will work with the projections of a,b,c vectors onto ABC plane (normalize them also):

vec3 pa = normalize(a-R);
vec3 pb = normalize(b-R);
vec3 pc = normalize(c-R);

Reminder: we must spin pa,pb,pc around R to ensure they spread around with 120 degree between each of them. Let’s check what are the angles between them now:

float papb = acos( dot(pa,pb) ) * sign(dot(cross(pa,pb),n));
float pbpc = acos( dot(pb,pc) ) * sign(dot(cross(pb,pc),n));
float pcpa = acos( dot(pc,pa) ) * sign(dot(cross(pc,pa),n));

The formulas are a bit complicated to distinguish between positive and negative angles as this is important.
Now the final task: we must spin pa, pb and pc by the corresponding error-correction angles to ensure those projections are spread 120 degrees apart from each other. I will omit the 3 simple equations I’ve got and show only the derived formulas for those angles:

float paRotAngle = (papb-pcpa)*0.33333333333333333333333333333333;
float pbRotAngle = (pbpc-papb)*0.33333333333333333333333333333333;
float pcRotAngle = (pcpa-pbpc)*0.33333333333333333333333333333333;

Let’s direct the pa,pb,pc wherever they should point by spinning them around normalized R vector:

vec3 RotAxis = normalize®;
pa = pa*cos(paRotAngle) + cross(RotAxis,pa)sin(paRotAngle);
pb = pb
cos(pbRotAngle) + cross(RotAxis,pb)sin(pbRotAngle);
pc = pc
cos(pcRotAngle) + cross(RotAxis,pc)*sin(pcRotAngle);

Now we are ready to calculate the resulting orthogonal vectors A,B,C:

vec3 A = R+pa0.81649658092772603273242802490196;
vec3 B = R+pb
0.81649658092772603273242802490196;
vec3 C = R+pc*0.81649658092772603273242802490196;

Done!

As the result, we got orthogonal vectors A,B,C of unit length out of non-collinear vectors of unit length a,b,c. The angles each of the a,b,c vectors were redirected are minimum possible to preserve the average direction of the vector triplet:

normalize(a+b+c) == normalize(A+B+C);

Gimme the cookie, I did good, didn’t I? :smiley:

The possibility that there may be multiple solutions?

If you replace equations 7 and 8 with

dot(a,A)=t (7)
dot(b,B)=t (8)
dot(c,C)=t (9)

Then for any given value of t, you have 9 equations in 9 unknowns. You’re (presumably) looking for the smallest value of t for which there exists a solution.

Geometric interpretation:

Equations 4 through 6 constrain A,B,C to lie on the unit sphere. Each of equations 7 through 9 constrain one of the vectors to lie on a cone whose axis is a, b or c and whose half-angle is acos(t). Between the two groups of equations, each vector is constrained to lie on a circle.

Equations 1 through 3 fix the distances between vectors (|A-B|, |B-C|, |C-A|).

Consider vectors A and B. Each is constrained to a circle by eqs 4+7 and 5+8 respectively, and the distance between them is fixed by eqs 1+4+5. This leaves one degree of freedom, resulting in a situation equivalent to a four-bar linkage. For any A, you typically have two valid values for B (intersection of a circle with a sphere).

The locus of C is now determined by eqs 2+3. The intersection of this locus with the circle defined by eqs 6+9 determines the solution. For small values of t, there will be no solutions (somewhere, a sqrt() will have a negative argument). For larger values there may be multiple solutions. At the limiting value of t, there would typically be one solution, where the locus defined by A and B is tangent to the circle defined by eqs 6+9.

[QUOTE=Yandersen;1263786]I just developed some alternative solution. …The idea is that the average direction of all three vectors

vec3 n=normalize(a+b+c));

should remain constant after a,b,c become orthonormal. If so, then the centroid of the triangle abc (which is based on those input vectors) will also become the circumcenter of the triangle ABC. And if A,B,C - are desired orthonormal vectors, then the ABC triangle will be equilateral with it’s circumcenter (point R) distanced from origin by “1/sqrt(3)”. So we find coordinates of this point first:

vec3 R = normalize(a+b+c)*0.57735026918962576450914878050196;

…[/QUOTE]

Fun exercise. I haven’t picked through all of your method, but here’s a few thoughts to consider on the first part. Suppose your axis vectors have drifted away from normal (extreme example to illustrate the point: if |a| = 2, |b| = 1, and |c| = 1). What does summing them give you? Not an average direction. Need to normalize them first. Then |a|+|b|+|c| gives you something.

Also, your 1/sqrt(3) presumes orthonormal already, which they aren’t necessarily (that’s the point of this exercise!). However, if you just assume the 3 normalized axes endpoints form a triangle, you can use barycentric coords to find the center:

vec3 R = 1/3 * ( normalize(a) + normalize(b) + normalize© )

But… Do you really need to deal with orthonormal bases drifting out of orthonormal and reorthonormalizing occasionally, …or could just use quaternions?

Right! Just an illustration: if we take an orthotriplet and rotate it by any angle around the central axis - it will match another orthotriplet, and all 3 vectors will be rotated by the same angle in result. So there are unlimited solutions if vectors are initially orthogonal. Howether, for non-orthogonal vectors I am not sure there could be more than one solution (or two, incluiding the mirror-variant from the other side of the origin). So it would really help to know what would be the minimum angle from which we can derive t. No idea how to pick the value. It just “feels” like the angle can be predicted based on angles between the original vectors and desired 90 degree angles, but no prove.

But it brings us back to the classic normalization approach where one axis normalized and preserved and the others dance around it. We can’t give preferences to any of the axises just to make a starting point and make further calculations based on that decision.

[QUOTE=Dark Photon;1263792]Suppose your axis vectors have drifted away from normal (extreme example to illustrate the point: if |a| = 2, |b| = 1, and |c| = 1). What does summing them give you?[/QUOTE]Sry, yeah, those a,b,c are of unit length, prenormalized as the very first step. Just forgot to mention that (said in first post only).

[QUOTE=Dark Photon;1263792]Also, your 1/sqrt(3) presumes orthonormal already, which they aren’t necessarily (that’s the point of this exercise!).[/QUOTE]You get it wrong. Once the vectors become orthonormal, the triangle will become equilateral. Then the circumcenter, centroid and all other magic points will match. So we just mark a point on a diagonal and claim it will be the circumcenter of the resulting triangle, and from there we build it then. As the result, we have R laying on the diagonal, and, later on, it will remain on the diagonal as well. That is the point.

[QUOTE=Dark Photon;1263792]But… Do you really need to deal with orthonormal bases drifting out of orthonormal and reorthonormalizing occasionally, …or could just use quaternions?[/QUOTE]I don’t like those! I like matrices, don’t you understand? I am a big fan of matrices, yes! :slight_smile:
Actually, that orthonormalization is not to be performed every time the matrix is multiplied. It can take a thousand cycles before error will become significant enough, even if single-precision floating point matrix is used. Here is the function code:

             #include "GLSL.hpp"   //GLSL-style stuff

using namespace GLSL;[FONT=monospace] //Expose sin, cos, sign and other standard math routines from GLSL.hpp

mat3 Orthonormalize(const mat3& m){
//Normalize axises to ensure they represent directions
dvec3 Ox = normalize(m[0]), Oy = normalize(m[1]), Oz = normalize(m[2]);
//Find direction to the centroid of the triangle based on those vectors
dvec3 n = normalize(Ox+Oy+Oz);
//Find the desired centroid location
dvec3 C = n*0.57735026918962576450914878050196;
//Find intersections of axises with the plane perpendicular to n
Ox = 0.57735026918962576450914878050196/dot(Ox,n);
Oy = 0.57735026918962576450914878050196/dot(Oy,n);
Oz = 0.57735026918962576450914878050196/dot(Oz,n);
//Project each axis onto plane perpendicular to n
dvec3 OxC = 0.81649658092772603273242802490196
normalize(Ox-C);
dvec3 OyC = 0.81649658092772603273242802490196
normalize(Oy-C);
dvec3 OzC = 0.81649658092772603273242802490196
normalize(Oz-C);
//Ensure each axis makes desired equal angle with n
Ox = C + OxC;
Oy = C + OyC;
Oz = C + OzC;
//Calculate angles between axises’ projections
double aOxOy = acos(dot(OxC,OyC)*1.5)*sign(dot(cross(OxC,OyC),n));
double aOyOz = acos(dot(OyC,OzC)*1.5)*sign(dot(cross(OyC,OzC),n));
double aOzOx = acos(dot(OzC,OxC)*1.5)*sign(dot(cross(OzC,OxC),n));
//Calculate rotation angles
double OxRotA = (aOxOy-aOzOx)*0.33333333333333333333333333333333;
double OyRotA = (aOyOz-aOxOy)*0.33333333333333333333333333333333;
double OzRotA = (aOzOx-aOyOz)0.33333333333333333333333333333333;
//Rotate each axis
Ox = C + OxC
cos(OxRotA) + cross(n,OxC)sin(OxRotA);
Oy = C + OyC
cos(OyRotA) + cross(n,OyC)sin(OyRotA);
Oz = C + OzC
cos(OzRotA) + cross(n,OzC)*sin(OzRotA);
//Return the resulting matrix
return mat3(Ox,Oy,Oz);
}[/FONT]

P.S.: GLSL.hpp is required.

This response suggests that you either didn’t read what I wrote, or didn’t understand it.

When solving a system of equations, the order in which you rearrange the equations doesn’t affect the set of solutions.

Oh, sry, GClements, I see, I get it wrong. So you suggest solving those 9 equations with solution containing the t parameter, then determine t[SUB]min[/SUB] from there and find the result for t[SUB]min[/SUB]. So the system you suggest, is this:

  1. |A-B|2 = 2
  2. |B-C|2 = 2
  3. |C-A|2 = 2
  4. |A|2 = 1
  5. |B|2 = 1
  6. |C|2 = 1
  7. dot(a,A)=t
  8. dot(b,B)=t
  9. dot(c,C)=t
    {10} t -> 0

And I have no idea how to solve it either mathematically or geometrically (thanks for sphere and circles model, btw, I happened not to think of it this way).
So in geometrical interpretation, the equations describe 3 circles laying on the unit sphere, and the equal radiuses of those circles are related to t: the smaller the t the wider the circles. So we go from t=1 to t=0 widening the circles until there will be 3 points on those 3 circles with sqrt(2) distance between each pair of those points.

But I see now why the whole thing is wrong. There is a valid solution may exist for non-ortho triplets, where the average of rotation axises for input vectors will be non-zero and as the result the system will turn, which is not desired. It’s like yes, all the vectors rotated by the same angle - just as we wanted - but at the end, all of them rotated primarily to one side, so the whole system get rotated to one side. In controversy, the alternative solution I gave above (with the code) ensures that the average direction stays the same, and i think this makes more sense.

[QUOTE=Yandersen;1263802]

You get it wrong. Once the vectors become orthonormal, the triangle will become equilateral. [/QUOTE]

Once they’re orthonormal, yes (orthogonal + normalized). But they’re only normalized at this point, not orthogonal. So they don’t form an equalateral triangle, and thus the center isn’t 1/sqrt(3) distance from the origin.
If they haven’t drifted too far out of orthogonal, then 1/sqrt(3) might be OK. But if they have, then…

Correct.

An analogous problem is: given three points on a plane, find an equilateral triangle of a given radius so that each vertex is the same distance to its associated point.

Solving such a system of implicit equations typically involves elimination: choosing an equation, re-arranging it with one variable on the left and all others on the right, substituting that into the next equation and so on. Much like solving a system of linear equations except that the re-arrangement typically involves finding roots of a polynomial rather than converting a*x+b=c to x=(c-b)/a. The fact that the polynomial must have a root constrains the solution; in this case, there will only be solutions for values of t which don’t result in a need to find a square root of a negative value.

Yes. Equal rotations doesn’t necessarily imply minimum “average” rotation.

Also, the formulation using dot-products is probably ill-conditioned. For an input which is almost orthonormal, dot(A,a) etc will be very close to one. This is likely to result in the equations being numerically unstable. If you want to assert co-linearity, it’s usually better to assert that the magnitude of the cross-product is zero, i.e. use |AxB|=0 rather than A.B=1.

[QUOTE=GClements;1263826]An analogous problem is: given three points on a plane, find an equilateral triangle of a given radius so that each vertex is the same distance to its associated point.[/QUOTE]The only thing is if the triangle lays on tips of the input vectors, then equal distances on the triangle plane will not imply equal rotation angles, actually. But I considered this approach and ended up with the function posted above. I choose plane to be the one perpendicular to the average direction, distanced by 1/sqrt(3) from the origin. Then I projected input vectors onto it and scaled those projections away from the intersection point to the sqrt(2/3) length. Then I simply spread the projections to make 120 degrees between each pair and get the result.

[QUOTE=GClements;1263826]Solving such a system of implicit equations typically involves elimination: choosing an equation, re-arranging it with one variable on the left and all others on the right, substituting that into the next equation and so on. Much like solving a system of linear equations except that the re-arrangement typically involves finding roots of a polynomial rather than converting a*x+b=c to x=(c-b)/a.[/QUOTE]It is still messy. Maybe better to use parameters to represent proportions for input vectors to build the required vectors out of those like if they were a basis vectors?
But, arggh, what is the point if the solution found already is not bad at all?

Dark Photon, If you think my solution is wrong, just copy-paste the function posted above and run it with few examples to see for yourself.

But it bothers me that the net sum of the rotation vectors is not zero - that means that the whole system is turned, actually. How is this possible if the average direction stays the same? Hm… Probably choosing the direction to centroid as the key vector is not the best choice. I also figured out that direction to circumcenter and orthocenter are not the best choices as those points may lay outside the triangle, so the system will definitely rotate to the side in such cases. Maybe the center of incircle will be the best choice?..

Orthogonal Sets of Vectors and the Gram-Schmidt Process