I'm trying to build my own little engine to get more familiar with the concepts of 3D programming.

I have a transform class that on each frame it creates a Scaling Matrix (S), a Rotation Matrix from a Quaternion (R) and concatenates them together (S*R). Once i have SR, I insert the translation values into the bottom of the three columns. So i end up with a transformation matrix that looks like:

[SR SR SR 0]
[SR SR SR 0]
[SR SR SR 0]
[tx ty tz 1]

This works perfectly in all cases except when rotating an object that has a non-uniform scale.

For example a unit cube with ScaleX = 4, ScaleY = 2, ScaleZ = 1 will give me a rectangular box that is 4 times as wide as the depth and twice as high as the depth.

If i then translate this around, the box stays the same and looks normal.

The problem happens whenever I try to rotate this scaled box. The shape itself becomes distorted and it appears as though the Scale factors are affecting the object on the World X,Y,Z axis rather than the local X,Y,Z axis of the object.

I've done some pretty extensive research through a variety of textbooks (Eberly, Moller/Hoffman, Phar etc) and there isn't a ton there to go off of. Online, most of the answers say to avoid non-uniform scaling which I understand the desire to avoid it, but I'd still like to figure out how to support it.

The only thing I can think off is that when constructing a Scale Matrix:

[sx 0  0  0]
[0  sy 0  0]
[0  0  sz 0]
[0  0  0  1]

This is scaling along the World Axis instead of the object's local Direction, Up and Right vectors or it's local Z, Y, X axis.

Does anyone have any tips or ideas on how to handle construction a transformation matrix that allows for non-uniform scaling and rotation?

Thanks!