Math: Combining a Translation, Rotation and Scale Matrix (edit: question solved by math magician Flavien :) )

I am currently continuing working my own double precision sandbox implementation (own double precision implementation of typical engine components) meant to run within Unity or later standalone (references to the Procedural Terrain Rendering How-To thread where we talk about procedural terrain/universes).
While the basic stuff (Object, GameObject, Component, Transform, Scene, Vector2d, Vector3d, Vector4d, Mathd, SystemTime etc etc) was more or less straightforward to implement, Quaterniond and the matrice classes Matrix2x2, Matrix3x3 and Matrix4x4 got a little more tricky. I think I am done with Quaterniond, Matrix2x2 and Matrix3x3. One remaining function within Matrix4x4 prevents me moving on.

As there are quite a few math geniouses around here by no doubt, or others who implemented an own engine as well, any help would be very much appreciated!

I’d like to add a “TRS” function that combines a Vector3d pos, Quaterniond q and Vector3d s to create a combined translation rotation and scale matrix. I have calculated the separate matrices within that function, but I am not sure how to merge these into the resulting TRS matrix. The translation components (m03,m13,m23,m33) are separated so I think at least for those components a simple matrix addition would work. However what about the components m00 to m22, where there is rotation and scale (m00,m11,m22) handled? Question to you math geniuses out there (which are around here by no doubt), what do you think would be the right way to merge that part?

See below the code (I normaly dont like posting code for such questions, but it helps to show the problem in that case).

        /// Creates a translation, rotation and scaling matrix.
        /// The returned matrix is such that places things at position pos, oriented in rotation q and scaled by s.
        /// </remarks>
        public static Matrix4x4Custom TRS(Vector3d pos, Quaterniond q, Vector3d s)
        {
            // Translation
            Matrix4x4Custom translation = new Matrix4x4Custom();
            translation.m03 = (float)pos.x;
            translation.m13 = (float)pos.y;
            translation.m23 = (float)pos.z;
            translation.m33 = 1.0f;

            // Rotation
            // Quaternion multiplication can be used to represent rotation. 
            // If a quaternion is represented by qw + i qx + j qy + k qz , then the equivalent matrix for rotation is:
            // Note: Small optimization would be to calcluate multiplications only once, similar to sqw, sqx, sqy, sqz.
            double sqw = q.w * q.w;
            double sqx = q.x * q.x;
            double sqy = q.y * q.y;
            double sqz = q.z * q.z;
            rotation.m00 = (float)(1 - 2 * sqy - 2 * sqz);
            rotation.m01 = (float)(2 * q.x * q.y - 2 * q.z * q.w);
            rotation.m02 = (float)(2 * q.x * q.z + 2 * q.y * q.w);
            rotation.m10 = (float)(2 * q.x * q.y + 2 * q.z * q.w);
            rotation.m11 = (float)(1 - 2 * sqx - 2 * sqz);
            rotation.m12 = (float)(2 * q.y * q.z - 2 * q.x * q.w);
            rotation.m20 = (float)(2 * q.x * q.z - 2 * q.y * q.w);
            rotation.m21 = (float)(2 * q.y * q.z + 2 * q.x * q.w);
            rotation.m22 = (float)(1 - 2 * sqx - 2 * sqy);
            rotation.m33 = 1.0f;

            // Scale
            Matrix4x4Custom scale = new Matrix4x4Custom();
            scale.m00 = (float)s.x;
            scale.m11 = (float)s.y;
            scale.m22 = (float)s.z;
            scale.m33 = 1.0f;

            // TODO How to merge these?
            Matrix4x4Custom result = new Matrix4x4Custom();
            result = ???
            return result;
        }

Note: There are a few float casts within, thats because I initially have implemented the matrices in float and not yet in double precision. I have not yet determined if it makes sense and if double precision would be needed (I think it does), and if there’d be an performance impact.

3 Likes

You should multiply all those matrices together, not add them. So the result would be: rotation * scale * translation, or the opposite ( translation * scale * rotation ) depending on the order of operations.

It’d be more performant to build the final matrix directly though. For that you’d simply multiply m00 m11 and m22 between the rotation and scale components; in m01 m02 m10 m12 m20 m21 you keep the current rotation factors, and in m03 m13 m23 you put the translation components.

9 Likes

Thanks a lot Flavien, very much appreciated! I changed the code accordingly, and also took your optimization hint into account. You are right besides three unnecessary matrix initializations I can indeed avoid the more costly two full matrix multiplication and just multiply the elements where necessary.

Code has shortened quite nicely. There would be still some minor room for improvement in calculating some of the q.? * q.? operations only once, but that’d be a small enhancement, and I think I’ll stick with more readable code at that point.

Thanks again for the hint Flavien!

PS: In case someone stumbles across the same topic here is the codeupdate.

   public static Matrix4x4Custom TRS(Vector3d pos, Quaterniond q, Vector3d s)
   {
       Matrix4x4Custom result = new Matrix4x4Custom();
       // Rotation and Scale
       // Quaternion multiplication can be used to represent rotation. 
       // If a quaternion is represented by qw + i qx + j qy + k qz , then the equivalent matrix for rotation is (including scale):
       // Remarks: https://forums.inovaestudios.com/t/math-combining-a-translation-rotation-and-scale-matrix-question-to-you-math-magicians/5194/2
       double sqw = q.w * q.w;
       double sqx = q.x * q.x;
       double sqy = q.y * q.y;
       double sqz = q.z * q.z;
       result.m00 = (float)(1 - 2 * sqy - 2 * sqz) * (float)s.x;
       result.m01 = (float)(2 * q.x * q.y - 2 * q.z * q.w);
       result.m02 = (float)(2 * q.x * q.z + 2 * q.y * q.w);
       result.m10 = (float)(2 * q.x * q.y + 2 * q.z * q.w);
       result.m11 = (float)(1 - 2 * sqx - 2 * sqz) * (float)s.y;
       result.m12 = (float)(2 * q.y * q.z - 2 * q.x * q.w);
       result.m20 = (float)(2 * q.x * q.z - 2 * q.y * q.w);
       result.m21 = (float)(2 * q.y * q.z + 2 * q.x * q.w);
       result.m22 = (float)(1 - 2 * sqx - 2 * sqy) * (float)s.z;
       // Translation
       result.m03 = (float)pos.x;
       result.m13 = (float)pos.y;
       result.m23 = (float)pos.z;
       result.m33 = 1.0f;
       // Return result
       return result;
   }
8 Likes

I love this community!

5 Likes