I've been trying to implement transformations which are defined as always being performed around the current coordinate system (i.e., the coordinate system that's defined by applying Current Transformation Matrix).

I'm working in a row-major environment i.e., result = point * CTM.

In my particular definition of transform semantics, performing a -90º rotation and then a scale of Y by 2 gives the same result regardless of order as the operations are defined to operate on the coordinate system defined the current CTM.

My first approach was to implement an operation as follows: firstly, append the inverse of the current CTM which gets us the identity, then apply the transform in question and then re-apply the CTM. For example, if S = standard scale Y by 2, I will then append CTM^-1 * S * CTM, giving us result = point * CTM * CTM^-1 * S * CTM = point * S * CTM. This seemingly gives the correct result as S is performed in the identity coordinate system, thus it will scale along the Y axis.

But this strategy fails when applied again. For example, if R = rotate by 90º, I will then append (S * CTM)^-1 * R * (S * CTM), finally giving us result = point * (S * CTM) * (S * CTM)^-1 * R * (S * CTM) = point * R * S * CTM. Obviously this produces the wrong result, as the scale will be performed with the Y axis rotated. The root issue is that S was constructed assuming it was relative to the identity coordinate system but it's not anymore.

How would I go about constructing scale / rotations / translation that are defined on the coordinate system defined by the CTM? Do I have to derive matrices that perform those from scratch? For example, each transformation will be given the origin and unit vectors, in addition to the usual parameters for the particular transform.

Any help is greatly appreciated.

I've been trying to implement transformations which are defined as always being performed around the current coordinate system (i.e., the coordinate system that's defined by applying Current Transformation Matrix).

That's already a bit confusing. You are going to implement transformations in some affine space. So far, so good, but I am not sure what it means to performing them "around the current coordinate system". I'll ignore that for a minute.

I'm working in a row-major environment i.e., result = point * CTM.

So you are using a notation where vectors are rows and you apply transformations by multiplying a row by a matrix on the right. That's fine, but it's not what "row major" means. It refers to an order in which to store the elements of the matrix in memory.

In my particular definition of transform semantics, performing a -90º rotation and then a scale of Y by 2 gives the same result regardless of order as the operations are defined to operate on the coordinate system defined the current CTM.

Performing a -90º rotation and then scaling Y by 2 cannot possibly give you the same result than doing the operations in the reverse order, although I don't understand what you mean by the bit about the CTM.

The rest of the post is kind of confusing, but the confusion seems to stem from your understanding that rotating and scaling one axis can possibly commute. Perhaps you can try to explain things a bit more clearly. But if you could, chances are you wouldn't have a question. :)

CTM * ScaleMatrix will get one result, ScaleMatrix * CTM will get another. One will be applying the scale "inside" the coordinate system defined by CTM, the other will scale around the world axes. I'm not sure which is which offhand, but I'm pretty sure "inside" is what you want. Then a Y scale matrix will scale the local Y axis, so it will have the same effect whether it's done before or after a rotation.

Something that took me a long time to notice is that a row-major rotate/scale matrix is really just an array of axis vectors. In 3D, with a 3x3 matrix, the first 3 elements are the X axis (1,0,0 in an identity matrix), middle 3 are the Y axis (0,1,0 identity), and last 3 are the Z axis (0,0,1). After rotating, those axis vectors will be pointing in different directions, but they'll still be orthogonal to eachother. And the length of each axis vector is the scale for that axis. Instead of applying a Y scale matrix, you could cheat and simply multiply the middle 3 elements by the Y scale you want :)

In 2D, the first 2 elements are the X axis (1,0) and the second two are the Y axis (0,1), and the same rules apply. It's actually simpler, but a bit harder to see the beauty of it in 2D, at least for me.

I'm not sure which is which offhand, but I'm pretty sure "inside" is what you want. Then a Y scale matrix will scale the local Y axis, so it will have the same effect whether it's done before or after a rotation.

Correct: inside is what geometrically can be thought of as transforming the coordinate system. But it does not quite work when rotations get involved. The reason is that if you firstly scale, we get "scale * CTM" and then we rotate, the rotation will be prepended, i.e., "rotate * scale * CTM" and the scale no longer has the same meaning when it was applied (e.g., if scale was along Y but rotation was 90º, the scale would instead end up rotating the logical x-axis).

In any case, here are the functions that I've implemented which achieve the desired effect. Each of them takes the CTM to be able to construct the matrix according to the definition.

function scale(ctm, sx, sy) {
  scale_x = scale_along_vector(ctm, (1, 0), sx)
  scale_y = scale_along_vector(ctm, (0, 1), sy)
  return scale_x * scale_y
}

function scale_along_vector(ctm, p, k) {
  // k = scale
  mapped_origin = (0, 0) * ctm
  mapped_point = p * ctm
  v = normalise (mapped_point - mapped_origin)
  
  // 3x3 matrix in form of:
  // (a b 0)
  // (c d 0)
  // (0 0 1)
  scale_transform = {
    .a = 1 + (k - 1) * v.x^2
    .b = (k - 1) * v.x * v.y
    .c = (k - 1) * v.x * v.y
    .d = 1 + (k - 1) * v.y^2
  };

  translate_to_origin = translation(-mapped_origin.x, -mapped_origin.y)
  translate_back = translation(mapped_origin.x, mapped_origin.y)
  return translate_to_origin * scale_transform * translate_back
}

function rotate(ctm, radians) {
  // determine multiple (1.0 or -1.0) by mapping (1, 0) and
  // (0, 1) using the ctm and then determining whether the
  // rotation should be in opposite or normal direction
  multiplier = ... 
  
  rotation_transform = standard rotation matrix(radians * multiplier)
  translate_to_origin = translation(-mapped_origin.x, -mapped_origin.y)
  translate_back = translation(mapped_origin.x, mapped_origin.y)
  
  return translate_to_origin * rotation_transform * translate_back
}

function translate(ctm, tx, ty) {
  // compute unit vectors
  mapped_origin = (0, 0) * ctm
  xv = ((1, 0) * ctm) - mapped_origin
  yv = ((0, 1) * ctm) - mapped_origin
  
  // translate by tx, ty along unit vectors
  dx = xv.x * tx + yv.x * ty
  dy = xv.y * tx + yv.y * ty
  
  return translation matrix(dx, dy)
}