[BUG] [bugfix-2.0.x] G02 mathematical bug in calculations

[kpter]

Description

Steps to Reproduce

1. Latest bugfix-2.0.x
2. Enable in configs: CNC and G2 G3 and SD card support (configs attached)
3. Build with Platform IO to Arduino 2560 r3 and RAMPS 1.4
4. Real example G-code that behavior explain the problem

; - 1 - BUG, Mistake! (Seen linear fast motions parallel to the X axis, not circles)
G00  Z0 
G00  X61.075 Y10 
G02  X64.325 Y10 R1.625 F500 
G02  X57.825 Y10 Z-0.125 R3.25 

; - 2 - GOOD (circles)
G00  Z0
G00  X21.075 Y10 
G02  X24.325 Y10 R1.625 F500 
G02  X17.825 Y10 Z-0.125 R3.25 

; - 3 - GOOD (circles)
G00  Z0 
G00  X61.075 Y10 
G02  X64.325 Y10 R1.62501 F500 
G02  X57.825 Y10 Z-0.125 R3.25001 

Additional Information

Analyzing the differences, I came to the conclusion that here the arithmetic error of calculations with roundings in the AVR platform leads to the fact that the semicircle radius is smaller than the calculated diameter of the semicircle along the coordinates.

I do not program in C - therefore I cannot implement a check myself that would correct the situation.

All my configs.zip

[kpter]

I made a temporary workaround by adding + 0.00001f

$ diff G2_G3.cpp Marlin/src/gcode/motion/G2_G3.cpp
281c281
<       const float r = parser.value_linear_units(),
---
>       const float r = parser.value_linear_units() + 0.00001f,

It is resolve my example.
But is it good?
Maybe there is a stronger solution.

[shitcreek]

@thinkyhead, what do you think of this issue?

[kpter]

I found another error in the implementation of G2/G3.
If enabled (as in my config) #define CNC_WORKSPACE_PLANES ,
the code G2 working unexpected (as fast line):

G0 X8 Y10 Z0
G18  ; change plan to ZX
G2 X15 Z-7 I7 F600
G1 X30
G1 Y18

This is a real good sample, I check it on Camotics Camo.zip .
But, when I replace in G2 X15 Z-7 I7 F600 only one letter I by J - this works good as in simulation attached ( Camo.zip).

[DerAndere1]

Regarding the last comment: According to DIN 66025 / ISO 6983 and beckhoff , G2 should accept parameters I (for X-offset) and K (for Z-offset) after G18 ; change to ZX plane

According to sample 3 of the beckhoff reference, this should happen:

G01 G18 X100 Y100 Z100 F6000
G02 I0 K50 X150 Z150 ; quarter circle in ZX plane 

So it does not make sense to me that Marlin works with parameter G02 J after G18.

[kpter]

In source code "parser.cpp" there are no implementation of processing of 'K' codes

    #if ENABLED(GCODE_MOTION_MODES)
      #if ENABLED(ARC_SUPPORT)
        case 'I': case 'J': case 'R':
          if (motion_mode_codenum != 2 && motion_mode_codenum != 3) return;
      #endif
      case 'P': case 'Q':
        if (motion_mode_codenum != 5) return;
      case 'X': case 'Y': case 'Z': case 'E': case 'F':
        if (motion_mode_codenum < 0) return;
        command_letter = 'G';

[edwilliams16]

The most obvious candidate for loss of precision is line 289 in G2_G3.cpp:

h = SQRT(sq(r) - sq(d * 0.5f)), // Distance to the arc pivot-point

to preserve precision this should be

h = SQRT((r - d * 0.5f) * (r + d * 0.5f)), // Distance to the arc pivot-point

[thinkyhead]

The PR linked above should correct the workspace plane parameters issue. We'll check G5 next to make sure it's also compliant.

[kpter]

The most obvious candidate for loss of precision is line 289 in G2_G3.cpp:

h = SQRT(sq(r) - sq(d * 0.5f)), // Distance to the arc pivot-point

to preserve precision this should be

h = SQRT((r - d * 0.5f) * (r + d * 0.5f)), // Distance to the arc pivot-point

No, I checked it out yesterday on my CNC with the latest firmware, as well as on the old firmware (I checked this too) - it does not change the situation.
Moreover, with the latest firmware it has become even worse. Now does not work all my three examples that I gave in the first post. Also, it not helped solve the problem by applying my old workaround (add +0.00001) in the new firmware.

The PR linked above should correct the workspace plane parameters issue. We'll check G5 next to make sure it's also compliant.

Yes, I checked it yesterday on my CNC with last firmware - now the CNC workspace plane works fine.

[kpter]

New changes in version
#define STRING_DISTRIBUTION_DATE "2019-10-08"
from previous version that I tested
#define STRING_DISTRIBUTION_DATE "2019-08-21"
led to deterioration.
So we must look for that has changed.
I have already verified that the form of expression is not affected. Therefore, we can return to the previous form (h = SQRT(sq(r) - sq(d * 0.5f))) , if it is readable better in source code.

[edwilliams16]

I've pointed out above two errors in the commits. Things were complicated by the fact that I had a week-old copy of the code when I made my comment and the code in the meantime had been refactored using struct.

(r - x)*(r + x) and r*r - x*x are obviously equivalent algebraically, but the former is less susceptible to rounding error, when (r -x) is small.

@thinkyhead

shouldn't this be
h = SQRT((r - len * 0.5f) * (r + len * 0.5f));

I'm not sure what it does using d It looks illegal.

Also, I believe
const xy_pos_t s = { d.x, -d.y }; // Inverse slope of the perpendicular bisector
should be
const xy_pos_t s = { -d.y, d.x }; // move rotated 90 degrees

The above error was introduced in the refactoring #15204```

 

[edwilliams16]

Here's a suggested rewrite of this code section:

    if (parser.seenval('R')) {
      const float r = parser.value_linear_units();
      if (r) {
        const xy_pos_t p1 = current_position, p2 = destination;
        if (p1 != p2) {
          const xy_pos_t d2 = (p2 - p1) * 0.5f;                 // XY vector to midpoint of move from current
          const float e = clockwise ^ (r < 0) ? -1 : 1,         // clockwise -1/1, counterclockwise 1/-1
                      len = d2.magnitude(),                     // Distance to mid-point of move from current
                      h = SQRT((r - len) * (r + len));          // Distance to the arc pivot-point from midpoint
          const xy_pos_t s = { -d2.y/len, d2.x/len };           // Unit vector along perpendicular bisector
          arc_offset = d2 + s * e * h;                          // The calculated offset
        }
      }
    }

I see one remaining problem. The code will crash with a floating point error if given an invalid G-code command with the radius less than half the distance to the destination point (geometrically impossible). How should that be handled? Ignore it? Do a G1 instead?

ab_float_t arc_offset = { 0, 0 };
    if (parser.seenval('R')) {
      const float r = parser.value_linear_units();
      if (r) {
        const xy_pos_t p1 = current_position, p2 = destination;
        if (p1 != p2) {
          const xy_pos_t d2 = (p2 - p1) * 0.5f;                 // XY vector to midpoint of move from current
          const float e = clockwise ^ (r < 0) ? -1 : 1,         // clockwise -1/1, counterclockwise 1/-1
                      len = d2.magnitude(),                     // Distance to mid-point of move from current
                      h2 = (r - len) * (r + len),
                      h = (h2  >= 0) ? SQRT(h2) : 0.0f;          // Distance to the arc pivot-point from midpoint
          const xy_pos_t s = { -d2.y/len, d2.x/len };           // Unit vector along perpendicular bisector
          arc_offset = d2 + s * e * h;                          // The calculated offset
        }
      }
    }

would fix it by replacing the impossible arc with a straight line, assuming the arc implementation can handle an infinite radius circle.

I can work up a pull request if you prefer, as this is no longer a one-liner.

I'm fighting git, but here's the diff:

-          const xy_pos_t d = p2 - p1, m = (p1 + p2) * 0.5f;   // XY distance and midpoint
-          const float e = clockwise ^ (r < 0) ? -1 : 1,       // clockwise -1/1, counterclockwise 1/-1
-                      len = d.magnitude(),                    // Total move length
-                      h = SQRT((r - d * 0.5f) * (r + d * 0.5f)); // Distance to the arc pivot-point
-          const xy_pos_t s = { d.x, -d.y };                   // Inverse Slope of the perpendicular bisector
-          arc_offset = m + s * RECIPROCAL(len) * e * h - p1;  // The calculated offset
+          const xy_pos_t d2 = (p2 - p1) * 0.5f;          // XY vector to midpoint of move from current
+          const float e = clockwise ^ (r < 0) ? -1 : 1,  // clockwise -1/1, counterclockwise 1/-1
+                      len = d2.magnitude(),              // Distance to mid-point of move from current
+                      h2 = (r - len) * (r + len),
+                      h = (h2  >= 0) ? SQRT(h2) : 0.0f;  // Distance to the arc pivot-point from midpoin
+          const xy_pos_t s = { -d2.y/len, d2.x/len};     // Unit vector along perpendicular bisector
+          arc_offset = d2 + s * e * h;                   // The calculated offset (mid-point if |r| < len)

[kpter]

Yes you are right!
I made these changes to the firmware 2019-10-08.
And it was all the work as I described in the first report. Then I added another amendment +0.00001f and it worked all my examples again.
Eventually, I adjusted the following lines

$ diff _Orig/G2_G3.cpp Marlin/src/gcode/motion/G2_G3.cpp
284c284
<       const float r = parser.value_linear_units();
---
>       const float r = parser.value_linear_units() +0.00001f;
291,292c291,292
<                       h = SQRT((r - d * 0.5f) * (r + d * 0.5f)); // Distance to the arc pivot-point
<           const xy_pos_t s = { d.x, -d.y };                   // Inverse Slope of the perpendicular bisector
---
>                       h = SQRT((r - len * 0.5f) * (r + len * 0.5f)); // Distance to the arc pivot-point
>           const xy_pos_t s = { -d.y, d.x };                   // Inverse Slope of the perpendicular bisector

Maybe insert check:
If r>d/2 than r=d/2 - it is better then my +0.00001f ?

[kpter]

Just now I saw that you have added the solution into your message.
I tested it on the CNC. It Works! Thank you!

I made these changes to the firmware 2019-10-08:

288,293c288,294
<           const xy_pos_t d = p2 - p1, m = (p1 + p2) * 0.5f;   // XY distance and midpoint
<           const float e = clockwise ^ (r < 0) ? -1 : 1,       // clockwise -1/1, counterclockwise 1/-1
<                       len = d.magnitude(),                    // Total move length
<                       h = SQRT((r - d * 0.5f) * (r + d * 0.5f)); // Distance to the arc pivot-point
<           const xy_pos_t s = { d.x, -d.y };                   // Inverse Slope of the perpendicular bisector
<           arc_offset = m + s * RECIPROCAL(len) * e * h - p1;  // The calculated offset
---
>           const xy_pos_t d2 = (p2 - p1) * 0.5f;          // XY vector to midpoint of move from current
>           const float e = clockwise ^ (r < 0) ? -1 : 1,  // clockwise -1/1, counterclockwise 1/-1
>                       len = d2.magnitude(),              // Distance to mid-point of move from current
>                       h2 = (r - len) * (r + len),
>                       h = (h2  >= 0) ? SQRT(h2) : 0.0f;  // Distance to the arc pivot-point from midpoin
>           const xy_pos_t s = { -d2.y/len, d2.x/len};     // Unit vector along perpendicular bisector
>           arc_offset = d2 + s * e * h;                   // The calculated offset (mid-point if |r| < len)

[Roxy-3D]

I see one remaining problem. The code will crash with a floating point error if given an invalid G-code command with the radius less than half the distance to the destination point (geometrically impossible). How should that be handled? Ignore it? Do a G1 instead?
...
would fix it by replacing the impossible arc with a straight line, assuming the arc implementation can handle an infinite radius circle.

I think I agree.... Replacing the impossible arc with a straight line probably makes sense. The extruder needs to be in the correct location at the end of the GCode command.

But that raises another question... The Extruder Axis needs to be in the 'correct' location also. So what do we do? My vote would be to extrude nothing along the straight line just because we don't know that there is room to do that. Or worse... The Slicer might be moving the extruder back to that location to fill in material and if we extrude along the straight line (instead of the arc) we could be causing a problem.

My suggestion is we artificially set the E-Axis to the 'correct' location at the end of the straight line move.

[edwilliams16]

@Roxy-3D @kpter

I mis-stated what my impossible arc fix above does. If the diameter of the arc is less than the distance to the destination point, it makes a minimum-radius arc (i.e. a semi-circle.) between the points.
Other than outright bugs in the slicer, this covers the most likely cause of the problem, namely that rounding of the slicer G-code output makes the radius of the semicircle slightly too small to fit between the points. If we assume this was what was intended, we don't need to play with the extruder.

[kpter]

I completely agree.
If this error occurs, it is negligible.
This error should not affect the correctness of the extrusion because it is clear that the error is only smaller than 0,00001 mm.
I vote for this implementation of check and correction h = (h2 >= 0) ? SQRT(h2) : 0.0f;.
And I have already checked these fixes on the real device. The problem is completely solved. You can add to the main code the corrections and close successfully my issue.

[edwilliams16]

OK. I created a pull request above.

[github-actions]

This issue has been automatically locked since there has not been any recent activity after it was closed. Please open a new issue for related bugs.