Skew calibration with calibration square? Maths help needed

[droftarts]

@oliof No worries, good meta gcode practice for me! I'm updating the main post for clarity, so might change in the next few minutes (but not the code). Let me know if you need it to do something different.

Ian

[oliof]

@droftarts I think just adding parameters so this can be put in as an M556.1 would be a great addition.

[droftarts]

@oliof I asked @dc42 that exact question, with that exact Gcode! However, I'm not even sure how many people are using skew compensation as to whether it's worth it, now that there is this workaround for diagonal lengths. Probably not a high priority. Feel free to add it to the firmware wishlist, though, see if it gets any traction.

It's also a question about whether this is a more accurate method of measuring skew than the existing method, though at least it doesn't need any extra vitamins. The accuracy of both methods depends on whether there are corner bulges and how uniform/good the print is. Any feedback on that would be useful.

Ian

[OwenD]

@droftarts
I'm not sure I agree with your drawing (as it applies to the calculations in Marlin).
If the X and Y axis are skewed you should end up with a rhombus.
i.e. all four sides are the same lengths, but the corner angles are no 90 degrees.

I came up with this macro
It requires RRF 3.5.0b1+ as it takes user input in M291
When you enter the measured length of the first diagonal it calculates what the other diagonal should be and puts that in as the default.

I have checked my calculations against this formula and they seem correct.

It also seems to work with the 100mm3 test piece I printed, although my printed Z dimensions do not match what they should be.
Someone has logged this as a bug, so I'm not sure if that's the cause

I think if we use M556 P1, the results can be directly applied.

I'd love if someone could review my methodology.

EDIT:
Code changed to correct mistake in calculating area and expected diagonals of rhombus

; set_skew_M556.g
; calculate the skew values based on diagonal measurements of a test object
;  https://www.thingiverse.com/thing:2972743/comments
; Mark the test point corners A,B,C,D as per the diagram below

; computation formulae from Marlin firmware
;  - Compute AB     : SQRT(2*AC*AC+2*BD*BD-4*AD*AD)/2
;  - XY_SKEW_FACTOR : TAN(PI/2-ACOS((AC*AC-AB*AB-AD*AD)/(2*AB*AD)))
; If desired, follow the same procedure for XZ and YZ.
;
; Use these diagrams for reference:
;
;    Y                     Z                     Z
;    ^     B-------C       ^     B-------C       ^     B-------C
;    |    /       /        |    /       /        |    /       /
;    |   /       /         |   /       /         |   /       /
;    |  A-------D          |  A-------D          |  A-------D
;    +-------------->X     +-------------->X     +-------------->Y
;     XY_SKEW_FACTOR        XZ_SKEW_FACTOR        YZ_SKEW_FACTOR

var xySkew = 0
var xzSkew =0
var yzSkew = 0

M291 S3 P"Measure X-Y Skew?" R"Measure X-Y"
G4 P200
; Maximum size limited to whichever axis is smallest
var max = min(move.axes[0].max,move.axes[1].max,move.axes[2].max)

M291 S6 P{"Enter side length of square (mm) max = " ^ var.max} R"Square" L10 H{var.max} F100.00 J1
var calculatedDiagonal = sqrt((input*input) + (input*input))
var A_D = input
var x2 = move.axes[0].max * move.axes[0].max
var y2 = move.axes[1].max * move.axes[1].max
; limit the diagonal to teh max possible iin the print volume
var maxDiagonal = sqrt(var.y2 + var.x2)

M291 S6 P"Enter A -> C measurement" R"A -> C"  R"X -> Y"  F{var.calculatedDiagonal} L10.0 H{var.maxDiagonal} J1
var A_C = input
; calculate area of rhombus given one side and diagonal
; 1/2 var.A_D sqrt(4*var.A_D*var.A_D)  - var.A_C * var.A_C
var area = 1/2* var.A_C * sqrt(4 * var.A_D * var.A_D  - var.A_C * var.A_C)
echo "Area = " ^ var.area
var difXY =  (var.area / var.A_C) * 2 ; calculates second diagonal of a rhombus given length of first
M291 S6 P"Enter B -> D measurement" R"B -> D" R"X -> Y" F{var.difXY} L10.0 H{var.maxDiagonal} J1
var B_D = input
var A_B = sqrt((2 * var.A_C * var.A_C) + (2 * var.B_D * var.B_D) - (4 * var.A_D * var.A_D))/2
set var.xySkew = tan(pi/2-acos((var.A_C*var.A_C-var.A_B*var.A_B-var.A_D*var.A_D)/(2*var.A_B*var.A_D)))
echo "Set M556 for X-Y Skew to M556 S1 X" ^ var.xySkew ^ " Ynn Znn"
G4 P200
M291 S4 P{"Apply skew (" ^ var.xySkew ^ ") to X-Y?"} R"Apply?" K{"Yes","No"} F0 J1 
if input = 0
	M556 S1 X{var.xySkew}
else
	set var.xySkew = 0 ; revert back to zero as not applied
G4 P200

; move to X Z 
M291 S3 P"Continue to X-Z Skew?" R"Measure X-Z"
G4 P200
M291 S6 P"Enter A -> C measurement" R"A -> C" R"X -> Z" F{var.calculatedDiagonal} L10.0 H{var.maxDiagonal} J1
set var.A_C = input
set var.area = 1/2* var.A_C * sqrt(4 * var.A_D * var.A_D  - var.A_C * var.A_C)
var difXZ =  (var.area / var.A_C) * 2 ; calculates second diagonal of a rhombus given length of first
echo var.difXZ
M291 S6 P"Enter B -> D measurement" R"B -> D" R"X -> Z" F{var.difXZ} L10.0 H{var.maxDiagonal} J1
set var.B_D = input
set var.A_B = sqrt((2 * var.A_C * var.A_C) + (2 * var.B_D * var.B_D) - (4 * var.A_D * var.A_D))/2
set var.xzSkew = tan(pi/2-acos((var.A_C*var.A_C-var.A_B*var.A_B-var.A_D*var.A_D)/(2*var.A_B*var.A_D)))
echo "Set M556 for X-Z skew to M556 S1 Xnn Ynn Z"^ var.xzSkew 
G4 P200
M291 S4 P{"Apply Skew (" ^var.xzSkew ^ ") to X-Z?"} R"Apply?" K{"Yes","No"} F0 J1 
if input = 0
	M556 S1 Z{var.xzSkew}
else
	set var.xzSkew = 0 ; revert to zero as not applied
G4 P200

; move to Y Z 
M291 S3 P"Continue to Y-Z Skew?" R"Measure Y-Z"
G4 P200
M291 S6 P"Enter A -> C measurement" R"A -> C" R"Y -> Z" F{var.calculatedDiagonal} L10.0 H{var.maxDiagonal} J1
set var.A_C = input
set var.area = 1/2* var.A_C * sqrt(4 * var.A_D * var.A_D  - var.A_C * var.A_C)
var difYZ =  var.area / var.A_C * 2
echo var.difYZ
M291 S6 P"Enter B -> D measurement" R"B -> D" R"Y -> Z" F{var.difYZ} L10.0 H{var.maxDiagonal} J1
set var.B_D = input
set var.A_B = sqrt((2 * var.A_C * var.A_C) + (2 * var.B_D * var.B_D) - (4 * var.A_D * var.A_D))/2
set var.yzSkew = tan(pi/2-acos((var.A_C*var.A_C-var.A_B*var.A_B-var.A_D*var.A_D)/(2*var.A_B*var.A_D)))
echo "Set M556 for Y-Z skew to M556 S1 Xnn Y" ^ var.yzSkew ^ "Znn" 
M291 S4 P{"Apply skew (" ^var.yzSkew ^ ") to Y-Z?"} R"Apply?" K{"Yes","No"} F0 J1 
if input = 0
	M556 S1 Y{var.yzSkew}
else
	set var.yzSkew = 0 ; revert back to zero as not applied

if (move.compensation.skew.tanXY != var.xySkew) || (move.compensation.skew.tanXZ != var.xzSkew) || (move.compensation.skew.tanYZ  != var.yzSkew)
	echo "Calculated skew settings should be:" 
	echo "M556 S1 X" ^ {var.xySkew}^ " Y" ^ {var.yzSkew} ^ " Z" ^ {var.xzSkew}
	echo "Actual setting"
	M556
	M291 S4 P{"One or more skew settings has not been applied - Apply Now?"} R"Settings not applied?" K{"Yes","No"} F0 J1 
	if input = 0
		M556 S1 X{var.xySkew} Y{var.yzSkew}	Z{var.xzSkew}
		M556
M291 S1 T10 P"Review results in console or send M556 to confirm" R"Done"

[OwenD]

@droftarts
I just ran your macro against mine and they essentially gave the same output when the same values were used, so I think I've misunderstood how you applied the measurements based on your drawing.

var AD = 180
 
var XY_AC = 244
var XY_BD = 265.57
var XY_skew_mm = 0
 
var XZ_AC = 254.5584
var XZ_BD = 254.5584
var XZ_skew_mm = 0
 
var YZ_AC = 254.5584
var YZ_BD = 254.5584
var YZ_skew_mm = 0
 

your output

My output

[OwenD]

Thinking on all this, I have to say that in many ways, the RRF method of using the measured skew at a given distance from a plane does probably make for an easier method.
It means you don't have to worry about over/under extrusion or other factors affecting the printed part size.
If I'm using diagonals then it only works if the printed part is exactly dimensionally accurate. Especially if it's based on a cube or square.
If I print something that's supposed to be 100x100x100 and it ends up 99.6 x 100 x 101.2 then my macro doesn't really work.
But if I print something approximately 110x110x110 and measure the skew at 100mm from the plane then the angles will be correct regardless of the actual dimension.
To me that probably gives you more chance of getting accurate skew settings without other possible inaccuracies in the system getting in the way.
Alternately we'd have to measure the actual dimensions on each face and calculate to diagonals of a parallelogram.
In truth it's probably taken longer to do the macro to use diagonals than it would have to make a measuring tool per the RRF docs.

As mentioned above, for some reason the diagonal measurement test piece I did was quite a bit out on the Z direction. (>1.5mm too high)
X & Y are near perfect (<0.2mm from dimension)
I have checked that Z steps are correct by measuring from frame to X gantry, carrying out a 100mm move and measuring again. Actual movement it exactly 100mm.
WTF?
Now I'm printing a plain 100mm cube to double check.
Down the rabbit hole I go

[droftarts]

@OwenD Thanks for chiming in! I was thinking about this earlier, and I think you are correct. I took the Marlin formulae as correct, and didn't think it through. The shape produced from misaligned axes would be a rhombus not a parallelogram. The two axes, assuming they have been configured correctly, should move the set distance in their plane, so the square would deform into a rhombus.

So, we shouldn't be calculating AB; AB = AD. Your calculations work in the same way as mine, taken from Marlin; AB is calculated where the point B is perpendicular to AD, at the distance AD. I haven't sat down and drawn out the two formula that Marlin uses, I took it that they were correct... but perhaps they are not! I did notice some small differences when calculating the skew factor from the Marlin formulae vs from the measured skew (the old way of doing it).

The interesting thing about the old way of doing this is that it is unambiguous. Take this:

If you measure up 150mm, you get a skew distance of 30mm. Skew factor is 30/150=0.2
If I run these diagonals and AD of 200 through yours I get 0.1999988, and through mine I also get 0.1999988, which is, arguably, very close. But it feels it should work out to 0.2!

However, as we're assuming the square deforms into a rhombus, the nice thing is that it makes an isosceles triangle with the diagonal measurements. So the trigonometry should be easy to do. But tomorrow, because it's a bit late here.

Ian

[OwenD]

@droftarts
I just realised I made a mistake in my calculations.
It only affected the calculation of the expected length of the diagonals and the area.
I was using the area of the square and not the rhombus.
If I use these dimensions with the RRF method and with my macro the results are essentially the same..

RRF Method

My macro

[oliof]

@droftarts I just meant to use the feature that allows you to give macros an unused gcode name and then calling that, without changes to the firmware (-;

[droftarts]

@oliof ah, yes, that could work too, once we decide on the correct calculations! I’ve been composing a longer reply for those for the last couple of days!

Ian

[droftarts]

I spent a bit of time looking at this, and have come to the conclusion that doing it the Marlin way is ... fine. A rhombus is a parallelogram with some extra constraints. If we knew that the shape created by misaligned axes was always going to be a rhombus, only the diagonals would need measuring, and everything else could be worked out from that. Using the Marlin formulae, it covers parallelograms AND rhombuses, though you need to provide the extra AD measurement. Realistically, I'd guess that diagonals can only be measured to an accuracy of perhaps 0.25mm (depending on accuracy of print/bulges/measuring device), so it's never going to be perfect, but certainly good enough.

And just because I wanted to know, I looked up what the two formulae Marlin uses are.

Compute AB : SQRT(2ACAC+2BDBD-4ADAD)/2

This is the parallelogram (and so also works for rhombi) formula to calculate the length of a side from the diagonals and the other side, see https://onlinemschool.com/math/formula/parallelogram/#h3
This can also be written as SQRT((AC*AC+BD*BD-AD*AD)/2)

XY_SKEW_FACTOR : TAN(PI/2-ACOS((ACAC-ABAB-ADAD)/(2AB*AD)))

ACOS((AC*AC-AB*AB-AD*AD)/(2*AB*AD)) This part is the 'Law of cosines' https://en.wikipedia.org/wiki/Law_of_cosines, and calculates the angle in radians at A given AB, AC and AD. See also https://www.calculatorsoup.com/calculators/geometry-plane/parallelogram.php
TAN(PI/2-A) gives the angle at A, ie between AE (the vertical point above A) and AB, and then using tan gives the skew factor.

@OwenD I had a look at your macro, and it works fine for me... in 3.5b1, because of the use of M291 S4 and S6! I also saw that you calculated some of the rhombus values, like area and the second diagonal from the first diagonal and AD, but then these weren't used in the calculations, I think?

@oliof Here's a generalised skew calculator macro. It also calculates the skew in mm, like the physical RRF method would measure, as a sense check. I'll do a M556.1 version when everyone agrees this is the best method, and I have some time!

; skew_calculator.g
;
; Bed Skew Compensation
;
; This feature corrects for misalignment in the XYZ axes.
;
; Take the following steps to get the bed skew in the XY plane:
; 1. Print a test square, e.g. https://www.thingiverse.com/thing:2563185 or https://www.thingiverse.com/thing:2972743
; 2. Mark the test point corners A,B,C,D as per the diagram below, if not already on the print
;
;    Y
;    ^  E--B-------C
;    |  | /       /
;    |  |/       /
;    |  A-------D
;    +-------------->X
;
; 3. Measure the diagonal A to C
; 4. Measure the diagonal B to D
; 5. Measure the edge A to D
; 6. Fill in the relevant measurements into the variables below, save and run the macro to show the result
; 7. If desired, follow the same procedure for XZ and YZ.
;
; Skew factors are computed automatically from these formulae, which may also be computed and set manually:
; Compute AB     : SQRT(2*AC*AC+2*BD*BD-4*AD*AD)/2
; XY_SKEW_FACTOR : TAN(PI/2-ACOS((AC*AC-AB*AB-AD*AD)/(2*AB*AD)))

var A_C = 320.15611
var B_D = 250
var A_D = 200

echo "Inputs: AC = " ^ var.A_C ^ ", BD = " ^ var.B_D ^ ", AD = " ^ var.A_D

var A_B = sqrt(2 * var.A_C * var.A_C + 2 * var.B_D * var.B_D - 4 * var.A_D * var.A_D)/2 ; if rhombus, A_B = A_D
var A_angle = acos((var.A_C * var.A_C - var.A_B * var.A_B - var.A_D * var.A_D)/(2 * var.A_B * var.A_D))
var xySkew = tan(pi/2-var.A_angle) ; skew
var A_E = var.A_B * sin(var.A_angle)
var E_B = sqrt(var.A_B * var.A_B - var.A_E * var.A_E)

echo "AB = " ^ var.A_B ^ ", XY skew = " ^ var.xySkew ^ ", Skew @ " ^ var.A_E ^ "mm = " ^ var.E_B ^ "mm"

Output:

M98 P"0:/macros/skew_calculator.g"
Inputs: AC = 320.1561, BD = 250, AD = 200
AB = 206.1552, XY skew = 0.2499996, Skew @ 200.0000mm = 49.99992mm

Ian

[dc42]

@droftarts maybe we should add a calculator for that to the existing calculators at reprapfirmware.org ?

[droftarts]

@dc42 I’ve added it to the wiki here: https://docs.duet3d.com/en/User_manual/Tuning/Orthogonal_axis_compensation#diagonal-measurement-method

Ian

[OwenD]

@droftarts said in Skew calibration with calibration square? Maths help needed:

I spent a bit of time looking at this, and have come to the conclusion that doing it the Marlin way is ... fine.
@OwenD I had a look at your macro, and it works fine for me... in 3.5b1, because of the use of M291 S4 and S6! I also saw that you calculated some of the rhombus values, like area and the second diagonal from the first diagonal and AD, but then these weren't used in the calculations, I think?

Thanks for your efforts

Yes, you're correct.
I used a formula to calculate what the second diagonal should be based on what's entered for the first one.
I only used the result to become the default displayed value in the M291 message box for the second diagonal.
This assumes that a skewed pair of axes should produce a rhombus.
My logic was that if my measurements for the second diagonal are very close to the displayed default then I'd use the default (or not) as I choose.
If it's significantly different then I should be asking why?
The actual calculation is done with whatever value is returned from the M291
However I felt, as you've indicated that if the piece being measured is not a rhombus then two A-D measurements are required.

I haven't yet tried "proving" the skew values on actual prints.
I did find that my YZ was significantly skewed (possibly knocked in a house move), but I couldn't bring myself to leave it that way, so I re-squared it.

[oliof]

@OwenD I completely agree that correcting skew mechanically is always the better option if at all possible.

[CalinFlorin86]

@droftarts
Maybe I am wrong , but the Marlin skew correction values are opposite than those for RRF. What I mean, as per old method of measuring the skew, if the angle was acute the values were negative ; if the angle is obtuse, the values were positive. As I see in the calculator, the acute angle has positive values and obtuse angle has negative values. With "angle" I mean the angle between AD and AB sides. Am I wrong?

[droftarts]

@CalinFlorin86 I think you are right! I'll have to do a couple of tests to check, but will update the documentation if it proves to be the case.

Ian