@Paciente8159 — I wanted to double-check that this line is correct. Unlike the line in apply_leveling above, it doesn't include yz_skew_factor. Should it instead be…

const float temprx = raw[X_AXIS] + raw[Y_AXIS] * planner.xy_skew_factor
    + raw[Z_AXIS] * (planner.xz_skew_factor - (planner.xy_skew_factor * planner.yz_skew_factor));

…or something like it? Thanks for clearing this up. Probably no one has noticed this yet, if it is an issue.

Sorry for the late reply.
I checked my local repo code.
It's almost that (- sign and not a + sign when factoring the angles). It should be like this:

const float temprx = raw[X_AXIS] - (raw[Y_AXIS] * planner.xy_skew_factor) - (raw[Z_AXIS] * (planner.xz_skew_factor - (planner.xy_skew_factor * planner.yz_skew_factor)));

Edit
Unless you have changed the way the angles are calculated in the configuration file. It's a matter of perspective. :-)

@thinkyhead
I've checked the config file and after doing some math the signs must be negative.
For example if the internal angle (lower left at corner A) is smaller than 90º the skew factor will be positive. X will have to be decremented to compensate for the angle.

The code I'm pointing out is unskew — reversal of the skew.

What I currently have in the code…

FORCE_INLINE static void skew(float &cx, float &cy, const float &cz) {
  if (WITHIN(cx, X_MIN_POS + 1, X_MAX_POS) && WITHIN(cy, Y_MIN_POS + 1, Y_MAX_POS)) {
    const float sx = cx - (cy * xy_skew_factor) - (cz * (xz_skew_factor - (xy_skew_factor * yz_skew_factor))),
                sy = cy - (cz * yz_skew_factor);
    if (WITHIN(sx, X_MIN_POS, X_MAX_POS) && WITHIN(sy, Y_MIN_POS, Y_MAX_POS)) {
      cx = sx; cy = sy;
    }
  }
}

FORCE_INLINE static void unskew(float &cx, float &cy, const float &cz) {
  if (WITHIN(cx, X_MIN_POS, X_MAX_POS) && WITHIN(cy, Y_MIN_POS, Y_MAX_POS)) {
    const float ux = cx + cy * xy_skew_factor + cz * xz_skew_factor,
                uy = cy + cz * yz_skew_factor;
    if (WITHIN(ux, X_MIN_POS, X_MAX_POS) && WITHIN(uy, Y_MIN_POS, Y_MAX_POS)) {
      cx = ux; cy = uy;
    }
  }
}

The line that does the "unskew" for X uses your original code:

ux = cx + cy * XY + cz * XZ;

But I think maybe it should be this:

ux = cx + cy * XY + cz * (XZ - (XY * XZ));

Would you say it should actually be this?:

ux = cx + cy * XY + cz * (XZ + (XY * XZ));

Ahh OK.
The correct way is:

ux = cx + cy * XY + cz * XZ;
uy = cy + cz * YZ;

And the reverse operation of this is:

sx = cx - (cy * XY) - (cz * (XZ- (XY* YZ)));
sy = cy - (cz * YZ);

This is the result of applying matrix operations:
The unskew matrix - the matrix with the factors computed from the real machine distortion of what should be a perfect cube;
And the skew matrix - the math equation equivalent to the inverted unskew matrix.

You can see this in action in OpenSCAD by comenting the top multmatrix on this code snippet:

//Some random values for the skew factor
XY=0.2;
XZ=0.1;
YZ=-0.3;

multmatrix(m = [ [1, -XY, -(XZ-(XY*YZ)), 0],
                 [0, 1, -YZ, 0],
                 [0, 0, 1, 0],
                 [0, 0, 0,  1]
              ])
multmatrix(m = [ [1, XY, XZ, 0],
                 [0, 1, YZ, 0],
                 [0, 0, 1, 0],
                 [0, 0, 0,  1]
              ])
cube(size=[100,100,100],center=false);

This is why I say that bed leveling can affect skew and vice versa. The math is not that direct and changes on one axis affect all axis. Matrix operations are a more effective/correct way to perform 3D machine to world transformations although for small distortions this can be acceptable.

Makes sense, thanks!