Shear transformation is incorrect when done on more than one axis

[ABotond]

Currently, the 2D shearing transformation matrix is defined as: $\begin{pmatrix} 1 & S_x \\ S_y & 1 \end{pmatrix}$.
However, as shearing on the x-axis is defined as $\begin{pmatrix} 1 & S_x \\ 0 & 1 \end{pmatrix}$, on the y-axis as $\begin{pmatrix} 1 & 0 \\ S_y & 1 \end{pmatrix}$, their composition is either $\begin{pmatrix} 1 + S_xS_y & S_x \\ S_y & 1 \end{pmatrix}$ or $\begin{pmatrix} 1 & S_x \\ S_y & 1 + S_xS_y \end{pmatrix}$.

Firstly, this means that if a shearing is done on both axis at the same time, the resulting operation is not a real shearing, as it is not area-preserving. Secondly, if two separate shearings are done after each other, the end result is not the same as if the transformation was done in one step.

Reproduction
A very minimal example:

import cv2
import monai.transforms
import numpy as np

if __name__ == "__main__":
    dd = {"image": np.expand_dims(cv2.imread("input.png", flags=cv2.IMREAD_GRAYSCALE), 0)}

    xy_transform= monai.transforms.Compose([
        monai.transforms.Affined(["image"], shear_params=(0.3, 0.0), padding_mode="zeros"),
        monai.transforms.Affined(["image"], shear_params=(0.0, 0.3), padding_mode="zeros"),
    ])

    yx_transform = monai.transforms.Compose([
        monai.transforms.Affined(["image"], shear_params=(0.0, 0.3), padding_mode="zeros"),
        monai.transforms.Affined(["image"], shear_params=(0.3, 0.0), padding_mode="zeros"),
    ])

    joined_transform = monai.transforms.Affined(["image"], shear_params=(0.3, 0.3), padding_mode="zeros")
    xy_image = xy_transform(dd)["image"]
    yx_image = yx_transform(dd)["image"]
    joined_image = joined_transform(dd)["image"]

    cv2.imwrite("variant-1.png", xy_image[0, :, :].astype(np.uint8))
    cv2.imwrite("variant-2.png", yx_image[0, :, :].astype(np.uint8))
    cv2.imwrite("variant-3.png", joined_image[0, :, :].astype(np.uint8))

Original	x-shear - y-shear	y-shear - x-shear	MONAI

Expected behavior
I expect the output of the composed shearing to be the same when I first shear only one axis, and then in a succinct step on another.

Proposed solution
In monai/transforms/utils.py in the _create_shear function line 988 should be changed either to
out[0, 1], out[1, 0], out[1, 1] = coefs[0], coefs[1], 1 + coefs[0] * coefs[1]
or to
out[0, 0], out[0, 1], out[1, 0] = 1 + coefs[0] * coefs[1], coefs[0], coefs[1]
depening on the desired order of the shearing operations.

The 3D version is also wrong, however, I did not derive the correct matrix right now.

Environment
I tested it on MONAI 0.9.0 and on MONAI 1.4.0

Edit: changed the images to make the differences easier to see.