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 = pbcos(pbRotAngle) + cross(RotAxis,pb)*sin(pbRotAngle);*

pc = pccos(pcRotAngle) + cross(RotAxis,pc)*sin(pcRotAngle);

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

vec3 A = R+pa*0.81649658092772603273242802490196;*

vec3 B = R+pb0.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?