Initial commit

10 files changed, 3909 insertions, 0 deletions

diff --git a/Marlin/src/feature/bedlevel/abl/abl.cpp b/Marlin/src/feature/bedlevel/abl/abl.cpp
new file mode 100644
index 0000000..3fb0cfc
--- /dev/null
+++ b/Marlin/src/feature/bedlevel/abl/abl.cpp
@@ -0,0 +1,421 @@
+/**
+ * Marlin 3D Printer Firmware
+ * Copyright (c) 2020 MarlinFirmware [https://github.com/MarlinFirmware/Marlin]
+ *
+ * Based on Sprinter and grbl.
+ * Copyright (c) 2011 Camiel Gubbels / Erik van der Zalm
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <https://www.gnu.org/licenses/>.
+ *
+ */
+
+#include "../../../inc/MarlinConfig.h"
+
+#if ENABLED(AUTO_BED_LEVELING_BILINEAR)
+
+#include "../bedlevel.h"
+
+#include "../../../module/motion.h"
+
+#define DEBUG_OUT ENABLED(DEBUG_LEVELING_FEATURE)
+#include "../../../core/debug_out.h"
+
+#if ENABLED(EXTENSIBLE_UI)
+  #include "../../../lcd/extui/ui_api.h"
+#endif
+
+xy_pos_t bilinear_grid_spacing, bilinear_start;
+xy_float_t bilinear_grid_factor;
+bed_mesh_t z_values;
+
+/**
+ * Extrapolate a single point from its neighbors
+ */
+static void extrapolate_one_point(const uint8_t x, const uint8_t y, const int8_t xdir, const int8_t ydir) {
+  if (!isnan(z_values[x][y])) return;
+  if (DEBUGGING(LEVELING)) {
+    DEBUG_ECHOPGM("Extrapolate [");
+    if (x < 10) DEBUG_CHAR(' ');
+    DEBUG_ECHO((int)x);
+    DEBUG_CHAR(xdir ? (xdir > 0 ? '+' : '-') : ' ');
+    DEBUG_CHAR(' ');
+    if (y < 10) DEBUG_CHAR(' ');
+    DEBUG_ECHO((int)y);
+    DEBUG_CHAR(ydir ? (ydir > 0 ? '+' : '-') : ' ');
+    DEBUG_ECHOLNPGM("]");
+  }
+
+  // Get X neighbors, Y neighbors, and XY neighbors
+  const uint8_t x1 = x + xdir, y1 = y + ydir, x2 = x1 + xdir, y2 = y1 + ydir;
+  float a1 = z_values[x1][y ], a2 = z_values[x2][y ],
+        b1 = z_values[x ][y1], b2 = z_values[x ][y2],
+        c1 = z_values[x1][y1], c2 = z_values[x2][y2];
+
+  // Treat far unprobed points as zero, near as equal to far
+  if (isnan(a2)) a2 = 0.0;
+  if (isnan(a1)) a1 = a2;
+  if (isnan(b2)) b2 = 0.0;
+  if (isnan(b1)) b1 = b2;
+  if (isnan(c2)) c2 = 0.0;
+  if (isnan(c1)) c1 = c2;
+
+  const float a = 2 * a1 - a2, b = 2 * b1 - b2, c = 2 * c1 - c2;
+
+  // Take the average instead of the median
+  z_values[x][y] = (a + b + c) / 3.0;
+  TERN_(EXTENSIBLE_UI, ExtUI::onMeshUpdate(x, y, z_values[x][y]));
+
+  // Median is robust (ignores outliers).
+  // z_values[x][y] = (a < b) ? ((b < c) ? b : (c < a) ? a : c)
+  //                                : ((c < b) ? b : (a < c) ? a : c);
+}
+
+//Enable this if your SCARA uses 180° of total area
+//#define EXTRAPOLATE_FROM_EDGE
+
+#if ENABLED(EXTRAPOLATE_FROM_EDGE)
+  #if GRID_MAX_POINTS_X < GRID_MAX_POINTS_Y
+    #define HALF_IN_X
+  #elif GRID_MAX_POINTS_Y < GRID_MAX_POINTS_X
+    #define HALF_IN_Y
+  #endif
+#endif
+
+/**
+ * Fill in the unprobed points (corners of circular print surface)
+ * using linear extrapolation, away from the center.
+ */
+void extrapolate_unprobed_bed_level() {
+  #ifdef HALF_IN_X
+    constexpr uint8_t ctrx2 = 0, xlen = GRID_MAX_POINTS_X - 1;
+  #else
+    constexpr uint8_t ctrx1 = (GRID_MAX_POINTS_X - 1) / 2, // left-of-center
+                      ctrx2 = (GRID_MAX_POINTS_X) / 2,     // right-of-center
+                      xlen = ctrx1;
+  #endif
+
+  #ifdef HALF_IN_Y
+    constexpr uint8_t ctry2 = 0, ylen = GRID_MAX_POINTS_Y - 1;
+  #else
+    constexpr uint8_t ctry1 = (GRID_MAX_POINTS_Y - 1) / 2, // top-of-center
+                      ctry2 = (GRID_MAX_POINTS_Y) / 2,     // bottom-of-center
+                      ylen = ctry1;
+  #endif
+
+  LOOP_LE_N(xo, xlen)
+    LOOP_LE_N(yo, ylen) {
+      uint8_t x2 = ctrx2 + xo, y2 = ctry2 + yo;
+      #ifndef HALF_IN_X
+        const uint8_t x1 = ctrx1 - xo;
+      #endif
+      #ifndef HALF_IN_Y
+        const uint8_t y1 = ctry1 - yo;
+        #ifndef HALF_IN_X
+          extrapolate_one_point(x1, y1, +1, +1);   //  left-below + +
+        #endif
+        extrapolate_one_point(x2, y1, -1, +1);     // right-below - +
+      #endif
+      #ifndef HALF_IN_X
+        extrapolate_one_point(x1, y2, +1, -1);     //  left-above + -
+      #endif
+      extrapolate_one_point(x2, y2, -1, -1);       // right-above - -
+    }
+
+}
+
+void print_bilinear_leveling_grid() {
+  SERIAL_ECHOLNPGM("Bilinear Leveling Grid:");
+  print_2d_array(GRID_MAX_POINTS_X, GRID_MAX_POINTS_Y, 3,
+    [](const uint8_t ix, const uint8_t iy) { return z_values[ix][iy]; }
+  );
+}
+
+#if ENABLED(ABL_BILINEAR_SUBDIVISION)
+
+  #define ABL_GRID_POINTS_VIRT_X (GRID_MAX_POINTS_X - 1) * (BILINEAR_SUBDIVISIONS) + 1
+  #define ABL_GRID_POINTS_VIRT_Y (GRID_MAX_POINTS_Y - 1) * (BILINEAR_SUBDIVISIONS) + 1
+  #define ABL_TEMP_POINTS_X (GRID_MAX_POINTS_X + 2)
+  #define ABL_TEMP_POINTS_Y (GRID_MAX_POINTS_Y + 2)
+  float z_values_virt[ABL_GRID_POINTS_VIRT_X][ABL_GRID_POINTS_VIRT_Y];
+  xy_pos_t bilinear_grid_spacing_virt;
+  xy_float_t bilinear_grid_factor_virt;
+
+  void print_bilinear_leveling_grid_virt() {
+    SERIAL_ECHOLNPGM("Subdivided with CATMULL ROM Leveling Grid:");
+    print_2d_array(ABL_GRID_POINTS_VIRT_X, ABL_GRID_POINTS_VIRT_Y, 5,
+      [](const uint8_t ix, const uint8_t iy) { return z_values_virt[ix][iy]; }
+    );
+  }
+
+  #define LINEAR_EXTRAPOLATION(E, I) ((E) * 2 - (I))
+  float bed_level_virt_coord(const uint8_t x, const uint8_t y) {
+    uint8_t ep = 0, ip = 1;
+    if (x > GRID_MAX_POINTS_X + 1 || y > GRID_MAX_POINTS_Y + 1) {
+      // The requested point requires extrapolating two points beyond the mesh.
+      // These values are only requested for the edges of the mesh, which are always an actual mesh point,
+      // and do not require interpolation. When interpolation is not needed, this "Mesh + 2" point is
+      // cancelled out in bed_level_virt_cmr and does not impact the result. Return 0.0 rather than
+      // making this function more complex by extrapolating two points.
+      return 0.0;
+    }
+    if (!x || x == ABL_TEMP_POINTS_X - 1) {
+      if (x) {
+        ep = GRID_MAX_POINTS_X - 1;
+        ip = GRID_MAX_POINTS_X - 2;
+      }
+      if (WITHIN(y, 1, ABL_TEMP_POINTS_Y - 2))
+        return LINEAR_EXTRAPOLATION(
+          z_values[ep][y - 1],
+          z_values[ip][y - 1]
+        );
+      else
+        return LINEAR_EXTRAPOLATION(
+          bed_level_virt_coord(ep + 1, y),
+          bed_level_virt_coord(ip + 1, y)
+        );
+    }
+    if (!y || y == ABL_TEMP_POINTS_Y - 1) {
+      if (y) {
+        ep = GRID_MAX_POINTS_Y - 1;
+        ip = GRID_MAX_POINTS_Y - 2;
+      }
+      if (WITHIN(x, 1, ABL_TEMP_POINTS_X - 2))
+        return LINEAR_EXTRAPOLATION(
+          z_values[x - 1][ep],
+          z_values[x - 1][ip]
+        );
+      else
+        return LINEAR_EXTRAPOLATION(
+          bed_level_virt_coord(x, ep + 1),
+          bed_level_virt_coord(x, ip + 1)
+        );
+    }
+    return z_values[x - 1][y - 1];
+  }
+
+  static float bed_level_virt_cmr(const float p[4], const uint8_t i, const float t) {
+    return (
+        p[i-1] * -t * sq(1 - t)
+      + p[i]   * (2 - 5 * sq(t) + 3 * t * sq(t))
+      + p[i+1] * t * (1 + 4 * t - 3 * sq(t))
+      - p[i+2] * sq(t) * (1 - t)
+    ) * 0.5f;
+  }
+
+  static float bed_level_virt_2cmr(const uint8_t x, const uint8_t y, const float &tx, const float &ty) {
+    float row[4], column[4];
+    LOOP_L_N(i, 4) {
+      LOOP_L_N(j, 4) {
+        column[j] = bed_level_virt_coord(i + x - 1, j + y - 1);
+      }
+      row[i] = bed_level_virt_cmr(column, 1, ty);
+    }
+    return bed_level_virt_cmr(row, 1, tx);
+  }
+
+  void bed_level_virt_interpolate() {
+    bilinear_grid_spacing_virt = bilinear_grid_spacing / (BILINEAR_SUBDIVISIONS);
+    bilinear_grid_factor_virt = bilinear_grid_spacing_virt.reciprocal();
+    LOOP_L_N(y, GRID_MAX_POINTS_Y)
+      LOOP_L_N(x, GRID_MAX_POINTS_X)
+        LOOP_L_N(ty, BILINEAR_SUBDIVISIONS)
+          LOOP_L_N(tx, BILINEAR_SUBDIVISIONS) {
+            if ((ty && y == (GRID_MAX_POINTS_Y) - 1) || (tx && x == (GRID_MAX_POINTS_X) - 1))
+              continue;
+            z_values_virt[x * (BILINEAR_SUBDIVISIONS) + tx][y * (BILINEAR_SUBDIVISIONS) + ty] =
+              bed_level_virt_2cmr(
+                x + 1,
+                y + 1,
+                (float)tx / (BILINEAR_SUBDIVISIONS),
+                (float)ty / (BILINEAR_SUBDIVISIONS)
+              );
+          }
+  }
+#endif // ABL_BILINEAR_SUBDIVISION
+
+// Refresh after other values have been updated
+void refresh_bed_level() {
+  bilinear_grid_factor = bilinear_grid_spacing.reciprocal();
+  TERN_(ABL_BILINEAR_SUBDIVISION, bed_level_virt_interpolate());
+}
+
+#if ENABLED(ABL_BILINEAR_SUBDIVISION)
+  #define ABL_BG_SPACING(A) bilinear_grid_spacing_virt.A
+  #define ABL_BG_FACTOR(A)  bilinear_grid_factor_virt.A
+  #define ABL_BG_POINTS_X   ABL_GRID_POINTS_VIRT_X
+  #define ABL_BG_POINTS_Y   ABL_GRID_POINTS_VIRT_Y
+  #define ABL_BG_GRID(X,Y)  z_values_virt[X][Y]
+#else
+  #define ABL_BG_SPACING(A) bilinear_grid_spacing.A
+  #define ABL_BG_FACTOR(A)  bilinear_grid_factor.A
+  #define ABL_BG_POINTS_X   GRID_MAX_POINTS_X
+  #define ABL_BG_POINTS_Y   GRID_MAX_POINTS_Y
+  #define ABL_BG_GRID(X,Y)  z_values[X][Y]
+#endif
+
+// Get the Z adjustment for non-linear bed leveling
+float bilinear_z_offset(const xy_pos_t &raw) {
+
+  static float z1, d2, z3, d4, L, D;
+
+  static xy_pos_t prev { -999.999, -999.999 }, ratio;
+
+  // Whole units for the grid line indices. Constrained within bounds.
+  static xy_int8_t thisg, nextg, lastg { -99, -99 };
+
+  // XY relative to the probed area
+  xy_pos_t rel = raw - bilinear_start.asFloat();
+
+  #if ENABLED(EXTRAPOLATE_BEYOND_GRID)
+    #define FAR_EDGE_OR_BOX 2   // Keep using the last grid box
+  #else
+    #define FAR_EDGE_OR_BOX 1   // Just use the grid far edge
+  #endif
+
+  if (prev.x != rel.x) {
+    prev.x = rel.x;
+    ratio.x = rel.x * ABL_BG_FACTOR(x);
+    const float gx = constrain(FLOOR(ratio.x), 0, ABL_BG_POINTS_X - (FAR_EDGE_OR_BOX));
+    ratio.x -= gx;      // Subtract whole to get the ratio within the grid box
+
+    #if DISABLED(EXTRAPOLATE_BEYOND_GRID)
+      // Beyond the grid maintain height at grid edges
+      NOLESS(ratio.x, 0); // Never <0 (>1 is ok when nextg.x==thisg.x)
+    #endif
+
+    thisg.x = gx;
+    nextg.x = _MIN(thisg.x + 1, ABL_BG_POINTS_X - 1);
+  }
+
+  if (prev.y != rel.y || lastg.x != thisg.x) {
+
+    if (prev.y != rel.y) {
+      prev.y = rel.y;
+      ratio.y = rel.y * ABL_BG_FACTOR(y);
+      const float gy = constrain(FLOOR(ratio.y), 0, ABL_BG_POINTS_Y - (FAR_EDGE_OR_BOX));
+      ratio.y -= gy;
+
+      #if DISABLED(EXTRAPOLATE_BEYOND_GRID)
+        // Beyond the grid maintain height at grid edges
+        NOLESS(ratio.y, 0); // Never < 0.0. (> 1.0 is ok when nextg.y==thisg.y.)
+      #endif
+
+      thisg.y = gy;
+      nextg.y = _MIN(thisg.y + 1, ABL_BG_POINTS_Y - 1);
+    }
+
+    if (lastg != thisg) {
+      lastg = thisg;
+      // Z at the box corners
+      z1 = ABL_BG_GRID(thisg.x, thisg.y);       // left-front
+      d2 = ABL_BG_GRID(thisg.x, nextg.y) - z1;  // left-back (delta)
+      z3 = ABL_BG_GRID(nextg.x, thisg.y);       // right-front
+      d4 = ABL_BG_GRID(nextg.x, nextg.y) - z3;  // right-back (delta)
+    }
+
+    // Bilinear interpolate. Needed since rel.y or thisg.x has changed.
+                L = z1 + d2 * ratio.y;   // Linear interp. LF -> LB
+    const float R = z3 + d4 * ratio.y;   // Linear interp. RF -> RB
+
+    D = R - L;
+  }
+
+  const float offset = L + ratio.x * D;   // the offset almost always changes
+
+  /*
+  static float last_offset = 0;
+  if (ABS(last_offset - offset) > 0.2) {
+    SERIAL_ECHOLNPAIR("Sudden Shift at x=", rel.x, " / ", bilinear_grid_spacing.x, " -> thisg.x=", thisg.x);
+    SERIAL_ECHOLNPAIR(" y=", rel.y, " / ", bilinear_grid_spacing.y, " -> thisg.y=", thisg.y);
+    SERIAL_ECHOLNPAIR(" ratio.x=", ratio.x, " ratio.y=", ratio.y);
+    SERIAL_ECHOLNPAIR(" z1=", z1, " z2=", z2, " z3=", z3, " z4=", z4);
+    SERIAL_ECHOLNPAIR(" L=", L, " R=", R, " offset=", offset);
+  }
+  last_offset = offset;
+  //*/
+
+  return offset;
+}
+
+#if IS_CARTESIAN && DISABLED(SEGMENT_LEVELED_MOVES)
+
+  #define CELL_INDEX(A,V) ((V - bilinear_start.A) * ABL_BG_FACTOR(A))
+
+  /**
+   * Prepare a bilinear-leveled linear move on Cartesian,
+   * splitting the move where it crosses grid borders.
+   */
+  void bilinear_line_to_destination(const feedRate_t &scaled_fr_mm_s, uint16_t x_splits, uint16_t y_splits) {
+    // Get current and destination cells for this line
+    xy_int_t c1 { CELL_INDEX(x, current_position.x), CELL_INDEX(y, current_position.y) },
+             c2 { CELL_INDEX(x, destination.x), CELL_INDEX(y, destination.y) };
+    LIMIT(c1.x, 0, ABL_BG_POINTS_X - 2);
+    LIMIT(c1.y, 0, ABL_BG_POINTS_Y - 2);
+    LIMIT(c2.x, 0, ABL_BG_POINTS_X - 2);
+    LIMIT(c2.y, 0, ABL_BG_POINTS_Y - 2);
+
+    // Start and end in the same cell? No split needed.
+    if (c1 == c2) {
+      current_position = destination;
+      line_to_current_position(scaled_fr_mm_s);
+      return;
+    }
+
+    #define LINE_SEGMENT_END(A) (current_position.A + (destination.A - current_position.A) * normalized_dist)
+
+    float normalized_dist;
+    xyze_pos_t end;
+    const xy_int8_t gc { _MAX(c1.x, c2.x), _MAX(c1.y, c2.y) };
+
+    // Crosses on the X and not already split on this X?
+    // The x_splits flags are insurance against rounding errors.
+    if (c2.x != c1.x && TEST(x_splits, gc.x)) {
+      // Split on the X grid line
+      CBI(x_splits, gc.x);
+      end = destination;
+      destination.x = bilinear_start.x + ABL_BG_SPACING(x) * gc.x;
+      normalized_dist = (destination.x - current_position.x) / (end.x - current_position.x);
+      destination.y = LINE_SEGMENT_END(y);
+    }
+    // Crosses on the Y and not already split on this Y?
+    else if (c2.y != c1.y && TEST(y_splits, gc.y)) {
+      // Split on the Y grid line
+      CBI(y_splits, gc.y);
+      end = destination;
+      destination.y = bilinear_start.y + ABL_BG_SPACING(y) * gc.y;
+      normalized_dist = (destination.y - current_position.y) / (end.y - current_position.y);
+      destination.x = LINE_SEGMENT_END(x);
+    }
+    else {
+      // Must already have been split on these border(s)
+      // This should be a rare case.
+      current_position = destination;
+      line_to_current_position(scaled_fr_mm_s);
+      return;
+    }
+
+    destination.z = LINE_SEGMENT_END(z);
+    destination.e = LINE_SEGMENT_END(e);
+
+    // Do the split and look for more borders
+    bilinear_line_to_destination(scaled_fr_mm_s, x_splits, y_splits);
+
+    // Restore destination from stack
+    destination = end;
+    bilinear_line_to_destination(scaled_fr_mm_s, x_splits, y_splits);
+  }
+
+#endif // IS_CARTESIAN && !SEGMENT_LEVELED_MOVES
+
+#endif // AUTO_BED_LEVELING_BILINEAR
diff --git a/Marlin/src/feature/bedlevel/abl/abl.h b/Marlin/src/feature/bedlevel/abl/abl.h
new file mode 100644
index 0000000..bbe2411
--- /dev/null
+++ b/Marlin/src/feature/bedlevel/abl/abl.h
@@ -0,0 +1,45 @@
+/**
+ * Marlin 3D Printer Firmware
+ * Copyright (c) 2020 MarlinFirmware [https://github.com/MarlinFirmware/Marlin]
+ *
+ * Based on Sprinter and grbl.
+ * Copyright (c) 2011 Camiel Gubbels / Erik van der Zalm
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <https://www.gnu.org/licenses/>.
+ *
+ */
+#pragma once
+
+#include "../../../inc/MarlinConfigPre.h"
+
+extern xy_pos_t bilinear_grid_spacing, bilinear_start;
+extern xy_float_t bilinear_grid_factor;
+extern bed_mesh_t z_values;
+float bilinear_z_offset(const xy_pos_t &raw);
+
+void extrapolate_unprobed_bed_level();
+void print_bilinear_leveling_grid();
+void refresh_bed_level();
+#if ENABLED(ABL_BILINEAR_SUBDIVISION)
+  void print_bilinear_leveling_grid_virt();
+  void bed_level_virt_interpolate();
+#endif
+
+#if IS_CARTESIAN && DISABLED(SEGMENT_LEVELED_MOVES)
+  void bilinear_line_to_destination(const feedRate_t &scaled_fr_mm_s, uint16_t x_splits=0xFFFF, uint16_t y_splits=0xFFFF);
+#endif
+
+#define _GET_MESH_X(I) float(bilinear_start.x + (I) * bilinear_grid_spacing.x)
+#define _GET_MESH_Y(J) float(bilinear_start.y + (J) * bilinear_grid_spacing.y)
+#define Z_VALUES_ARR  z_values
diff --git a/Marlin/src/feature/bedlevel/bedlevel.cpp b/Marlin/src/feature/bedlevel/bedlevel.cpp
new file mode 100644
index 0000000..0e0b87e
--- /dev/null
+++ b/Marlin/src/feature/bedlevel/bedlevel.cpp
@@ -0,0 +1,239 @@
+/**
+ * Marlin 3D Printer Firmware
+ * Copyright (c) 2020 MarlinFirmware [https://github.com/MarlinFirmware/Marlin]
+ *
+ * Based on Sprinter and grbl.
+ * Copyright (c) 2011 Camiel Gubbels / Erik van der Zalm
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <https://www.gnu.org/licenses/>.
+ *
+ */
+
+#include "../../inc/MarlinConfig.h"
+
+#if HAS_LEVELING
+
+#include "bedlevel.h"
+#include "../../module/planner.h"
+
+#if EITHER(MESH_BED_LEVELING, PROBE_MANUALLY)
+  #include "../../module/motion.h"
+#endif
+
+#if ENABLED(PROBE_MANUALLY)
+  bool g29_in_progress = false;
+#endif
+
+#if ENABLED(LCD_BED_LEVELING)
+  #include "../../lcd/marlinui.h"
+#endif
+
+#define DEBUG_OUT ENABLED(DEBUG_LEVELING_FEATURE)
+#include "../../core/debug_out.h"
+
+#if ENABLED(EXTENSIBLE_UI)
+  #include "../../lcd/extui/ui_api.h"
+#endif
+
+bool leveling_is_valid() {
+  return TERN1(MESH_BED_LEVELING,          mbl.has_mesh())
+      && TERN1(AUTO_BED_LEVELING_BILINEAR, !!bilinear_grid_spacing.x)
+      && TERN1(AUTO_BED_LEVELING_UBL,      ubl.mesh_is_valid());
+}
+
+/**
+ * Turn bed leveling on or off, fixing the current
+ * position as-needed.
+ *
+ * Disable: Current position = physical position
+ *  Enable: Current position = "unleveled" physical position
+ */
+void set_bed_leveling_enabled(const bool enable/*=true*/) {
+
+  const bool can_change = TERN1(AUTO_BED_LEVELING_BILINEAR, !enable || leveling_is_valid());
+
+  if (can_change && enable != planner.leveling_active) {
+
+    planner.synchronize();
+
+    #if ENABLED(AUTO_BED_LEVELING_BILINEAR)
+      // Force bilinear_z_offset to re-calculate next time
+      const xyz_pos_t reset { -9999.999, -9999.999, 0 };
+      (void)bilinear_z_offset(reset);
+    #endif
+
+    if (planner.leveling_active) {      // leveling from on to off
+      if (DEBUGGING(LEVELING)) DEBUG_POS("Leveling ON", current_position);
+      // change unleveled current_position to physical current_position without moving steppers.
+      planner.apply_leveling(current_position);
+      planner.leveling_active = false;  // disable only AFTER calling apply_leveling
+      if (DEBUGGING(LEVELING)) DEBUG_POS("...Now OFF", current_position);
+    }
+    else {                              // leveling from off to on
+      if (DEBUGGING(LEVELING)) DEBUG_POS("Leveling OFF", current_position);
+      planner.leveling_active = true;   // enable BEFORE calling unapply_leveling, otherwise ignored
+      // change physical current_position to unleveled current_position without moving steppers.
+      planner.unapply_leveling(current_position);
+      if (DEBUGGING(LEVELING)) DEBUG_POS("...Now ON", current_position);
+    }
+
+    sync_plan_position();
+  }
+}
+
+TemporaryBedLevelingState::TemporaryBedLevelingState(const bool enable) : saved(planner.leveling_active) {
+  set_bed_leveling_enabled(enable);
+}
+
+#if ENABLED(ENABLE_LEVELING_FADE_HEIGHT)
+
+  void set_z_fade_height(const float zfh, const bool do_report/*=true*/) {
+
+    if (planner.z_fade_height == zfh) return;
+
+    const bool leveling_was_active = planner.leveling_active;
+    set_bed_leveling_enabled(false);
+
+    planner.set_z_fade_height(zfh);
+
+    if (leveling_was_active) {
+      const xyz_pos_t oldpos = current_position;
+      set_bed_leveling_enabled(true);
+      if (do_report && oldpos != current_position)
+        report_current_position();
+    }
+  }
+
+#endif // ENABLE_LEVELING_FADE_HEIGHT
+
+/**
+ * Reset calibration results to zero.
+ */
+void reset_bed_level() {
+  if (DEBUGGING(LEVELING)) DEBUG_ECHOLNPGM("reset_bed_level");
+  #if ENABLED(AUTO_BED_LEVELING_UBL)
+    ubl.reset();
+  #else
+    set_bed_leveling_enabled(false);
+    #if ENABLED(MESH_BED_LEVELING)
+      mbl.reset();
+    #elif ENABLED(AUTO_BED_LEVELING_BILINEAR)
+      bilinear_start.reset();
+      bilinear_grid_spacing.reset();
+      GRID_LOOP(x, y) {
+        z_values[x][y] = NAN;
+        TERN_(EXTENSIBLE_UI, ExtUI::onMeshUpdate(x, y, 0));
+      }
+    #elif ABL_PLANAR
+      planner.bed_level_matrix.set_to_identity();
+    #endif
+  #endif
+}
+
+#if EITHER(AUTO_BED_LEVELING_BILINEAR, MESH_BED_LEVELING)
+
+  /**
+   * Enable to produce output in JSON format suitable
+   * for SCAD or JavaScript mesh visualizers.
+   *
+   * Visualize meshes in OpenSCAD using the included script.
+   *
+   *   buildroot/shared/scripts/MarlinMesh.scad
+   */
+  //#define SCAD_MESH_OUTPUT
+
+  /**
+   * Print calibration results for plotting or manual frame adjustment.
+   */
+  void print_2d_array(const uint8_t sx, const uint8_t sy, const uint8_t precision, element_2d_fn fn) {
+    #ifndef SCAD_MESH_OUTPUT
+      LOOP_L_N(x, sx) {
+        serial_spaces(precision + (x < 10 ? 3 : 2));
+        SERIAL_ECHO(int(x));
+      }
+      SERIAL_EOL();
+    #endif
+    #ifdef SCAD_MESH_OUTPUT
+      SERIAL_ECHOLNPGM("measured_z = ["); // open 2D array
+    #endif
+    LOOP_L_N(y, sy) {
+      #ifdef SCAD_MESH_OUTPUT
+        SERIAL_ECHOPGM(" [");             // open sub-array
+      #else
+        if (y < 10) SERIAL_CHAR(' ');
+        SERIAL_ECHO(int(y));
+      #endif
+      LOOP_L_N(x, sx) {
+        SERIAL_CHAR(' ');
+        const float offset = fn(x, y);
+        if (!isnan(offset)) {
+          if (offset >= 0) SERIAL_CHAR('+');
+          SERIAL_ECHO_F(offset, int(precision));
+        }
+        else {
+          #ifdef SCAD_MESH_OUTPUT
+            for (uint8_t i = 3; i < precision + 3; i++)
+              SERIAL_CHAR(' ');
+            SERIAL_ECHOPGM("NAN");
+          #else
+            LOOP_L_N(i, precision + 3)
+              SERIAL_CHAR(i ? '=' : ' ');
+          #endif
+        }
+        #ifdef SCAD_MESH_OUTPUT
+          if (x < sx - 1) SERIAL_CHAR(',');
+        #endif
+      }
+      #ifdef SCAD_MESH_OUTPUT
+        SERIAL_CHAR(' ', ']');            // close sub-array
+        if (y < sy - 1) SERIAL_CHAR(',');
+      #endif
+      SERIAL_EOL();
+    }
+    #ifdef SCAD_MESH_OUTPUT
+      SERIAL_ECHOPGM("];");               // close 2D array
+    #endif
+    SERIAL_EOL();
+  }
+
+#endif // AUTO_BED_LEVELING_BILINEAR || MESH_BED_LEVELING
+
+#if EITHER(MESH_BED_LEVELING, PROBE_MANUALLY)
+
+  void _manual_goto_xy(const xy_pos_t &pos) {
+
+    #ifdef MANUAL_PROBE_START_Z
+      constexpr float startz = _MAX(0, MANUAL_PROBE_START_Z);
+      #if MANUAL_PROBE_HEIGHT > 0
+        do_blocking_move_to_xy_z(pos, MANUAL_PROBE_HEIGHT);
+        do_blocking_move_to_z(startz);
+      #else
+        do_blocking_move_to_xy_z(pos, startz);
+      #endif
+    #elif MANUAL_PROBE_HEIGHT > 0
+      const float prev_z = current_position.z;
+      do_blocking_move_to_xy_z(pos, MANUAL_PROBE_HEIGHT);
+      do_blocking_move_to_z(prev_z);
+    #else
+      do_blocking_move_to_xy(pos);
+    #endif
+
+    current_position = pos;
+
+    TERN_(LCD_BED_LEVELING, ui.wait_for_move = false);
+  }
+
+#endif
+
+#endif // HAS_LEVELING
diff --git a/Marlin/src/feature/bedlevel/bedlevel.h b/Marlin/src/feature/bedlevel/bedlevel.h
new file mode 100644
index 0000000..a33f08a
--- /dev/null
+++ b/Marlin/src/feature/bedlevel/bedlevel.h
@@ -0,0 +1,102 @@
+/**
+ * Marlin 3D Printer Firmware
+ * Copyright (c) 2020 MarlinFirmware [https://github.com/MarlinFirmware/Marlin]
+ *
+ * Based on Sprinter and grbl.
+ * Copyright (c) 2011 Camiel Gubbels / Erik van der Zalm
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <https://www.gnu.org/licenses/>.
+ *
+ */
+#pragma once
+
+#include "../../inc/MarlinConfigPre.h"
+
+#if EITHER(RESTORE_LEVELING_AFTER_G28, ENABLE_LEVELING_AFTER_G28)
+  #define G28_L0_ENSURES_LEVELING_OFF 1
+#endif
+
+#if ENABLED(PROBE_MANUALLY)
+  extern bool g29_in_progress;
+#else
+  constexpr bool g29_in_progress = false;
+#endif
+
+bool leveling_is_valid();
+void set_bed_leveling_enabled(const bool enable=true);
+void reset_bed_level();
+
+#if ENABLED(ENABLE_LEVELING_FADE_HEIGHT)
+  void set_z_fade_height(const float zfh, const bool do_report=true);
+#endif
+
+#if EITHER(MESH_BED_LEVELING, PROBE_MANUALLY)
+  void _manual_goto_xy(const xy_pos_t &pos);
+#endif
+
+/**
+ * A class to save and change the bed leveling state,
+ * then restore it when it goes out of scope.
+ */
+class TemporaryBedLevelingState {
+  bool saved;
+  public:
+    TemporaryBedLevelingState(const bool enable);
+    ~TemporaryBedLevelingState() { set_bed_leveling_enabled(saved); }
+};
+#define TEMPORARY_BED_LEVELING_STATE(enable) const TemporaryBedLevelingState tbls(enable)
+
+#if HAS_MESH
+
+  typedef float bed_mesh_t[GRID_MAX_POINTS_X][GRID_MAX_POINTS_Y];
+
+  #if ENABLED(AUTO_BED_LEVELING_BILINEAR)
+    #include "abl/abl.h"
+  #elif ENABLED(AUTO_BED_LEVELING_UBL)
+    #include "ubl/ubl.h"
+  #elif ENABLED(MESH_BED_LEVELING)
+    #include "mbl/mesh_bed_leveling.h"
+  #endif
+
+  #define Z_VALUES(X,Y) Z_VALUES_ARR[X][Y]
+  #define _GET_MESH_POS(M) { _GET_MESH_X(M.a), _GET_MESH_Y(M.b) }
+
+  #if EITHER(AUTO_BED_LEVELING_BILINEAR, MESH_BED_LEVELING)
+
+    #include <stdint.h>
+
+    typedef float (*element_2d_fn)(const uint8_t, const uint8_t);
+
+    /**
+     * Print calibration results for plotting or manual frame adjustment.
+     */
+    void print_2d_array(const uint8_t sx, const uint8_t sy, const uint8_t precision, element_2d_fn fn);
+
+  #endif
+
+  struct mesh_index_pair {
+    xy_int8_t pos;
+    float distance;   // When populated, the distance from the search location
+    void invalidate() { pos = -1; }
+    bool valid() const { return pos.x >= 0 && pos.y >= 0; }
+    #if ENABLED(AUTO_BED_LEVELING_UBL)
+      xy_pos_t meshpos() {
+        return { ubl.mesh_index_to_xpos(pos.x), ubl.mesh_index_to_ypos(pos.y) };
+      }
+    #endif
+    operator xy_int8_t&() { return pos; }
+    operator const xy_int8_t&() const { return pos; }
+  };
+
+#endif
diff --git a/Marlin/src/feature/bedlevel/mbl/mesh_bed_leveling.cpp b/Marlin/src/feature/bedlevel/mbl/mesh_bed_leveling.cpp
new file mode 100644
index 0000000..ec5b95c
--- /dev/null
+++ b/Marlin/src/feature/bedlevel/mbl/mesh_bed_leveling.cpp
@@ -0,0 +1,133 @@
+/**
+ * Marlin 3D Printer Firmware
+ * Copyright (c) 2020 MarlinFirmware [https://github.com/MarlinFirmware/Marlin]
+ *
+ * Based on Sprinter and grbl.
+ * Copyright (c) 2011 Camiel Gubbels / Erik van der Zalm
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <https://www.gnu.org/licenses/>.
+ *
+ */
+
+#include "../../../inc/MarlinConfig.h"
+
+#if ENABLED(MESH_BED_LEVELING)
+
+  #include "../bedlevel.h"
+
+  #include "../../../module/motion.h"
+
+  #if ENABLED(EXTENSIBLE_UI)
+    #include "../../../lcd/extui/ui_api.h"
+  #endif
+
+  mesh_bed_leveling mbl;
+
+  float mesh_bed_leveling::z_offset,
+        mesh_bed_leveling::z_values[GRID_MAX_POINTS_X][GRID_MAX_POINTS_Y],
+        mesh_bed_leveling::index_to_xpos[GRID_MAX_POINTS_X],
+        mesh_bed_leveling::index_to_ypos[GRID_MAX_POINTS_Y];
+
+  mesh_bed_leveling::mesh_bed_leveling() {
+    LOOP_L_N(i, GRID_MAX_POINTS_X)
+      index_to_xpos[i] = MESH_MIN_X + i * (MESH_X_DIST);
+    LOOP_L_N(i, GRID_MAX_POINTS_Y)
+      index_to_ypos[i] = MESH_MIN_Y + i * (MESH_Y_DIST);
+    reset();
+  }
+
+  void mesh_bed_leveling::reset() {
+    z_offset = 0;
+    ZERO(z_values);
+    #if ENABLED(EXTENSIBLE_UI)
+      GRID_LOOP(x, y) ExtUI::onMeshUpdate(x, y, 0);
+    #endif
+  }
+
+  #if IS_CARTESIAN && DISABLED(SEGMENT_LEVELED_MOVES)
+
+    /**
+     * Prepare a mesh-leveled linear move in a Cartesian setup,
+     * splitting the move where it crosses mesh borders.
+     */
+    void mesh_bed_leveling::line_to_destination(const feedRate_t &scaled_fr_mm_s, uint8_t x_splits, uint8_t y_splits) {
+      // Get current and destination cells for this line
+      xy_int8_t scel = cell_indexes(current_position), ecel = cell_indexes(destination);
+      NOMORE(scel.x, GRID_MAX_POINTS_X - 2);
+      NOMORE(scel.y, GRID_MAX_POINTS_Y - 2);
+      NOMORE(ecel.x, GRID_MAX_POINTS_X - 2);
+      NOMORE(ecel.y, GRID_MAX_POINTS_Y - 2);
+
+      // Start and end in the same cell? No split needed.
+      if (scel == ecel) {
+        current_position = destination;
+        line_to_current_position(scaled_fr_mm_s);
+        return;
+      }
+
+      #define MBL_SEGMENT_END(A) (current_position.A + (destination.A - current_position.A) * normalized_dist)
+
+      float normalized_dist;
+      xyze_pos_t dest;
+      const int8_t gcx = _MAX(scel.x, ecel.x), gcy = _MAX(scel.y, ecel.y);
+
+      // Crosses on the X and not already split on this X?
+      // The x_splits flags are insurance against rounding errors.
+      if (ecel.x != scel.x && TEST(x_splits, gcx)) {
+        // Split on the X grid line
+        CBI(x_splits, gcx);
+        dest = destination;
+        destination.x = index_to_xpos[gcx];
+        normalized_dist = (destination.x - current_position.x) / (dest.x - current_position.x);
+        destination.y = MBL_SEGMENT_END(y);
+      }
+      // Crosses on the Y and not already split on this Y?
+      else if (ecel.y != scel.y && TEST(y_splits, gcy)) {
+        // Split on the Y grid line
+        CBI(y_splits, gcy);
+        dest = destination;
+        destination.y = index_to_ypos[gcy];
+        normalized_dist = (destination.y - current_position.y) / (dest.y - current_position.y);
+        destination.x = MBL_SEGMENT_END(x);
+      }
+      else {
+        // Must already have been split on these border(s)
+        // This should be a rare case.
+        current_position = destination;
+        line_to_current_position(scaled_fr_mm_s);
+        return;
+      }
+
+      destination.z = MBL_SEGMENT_END(z);
+      destination.e = MBL_SEGMENT_END(e);
+
+      // Do the split and look for more borders
+      line_to_destination(scaled_fr_mm_s, x_splits, y_splits);
+
+      // Restore destination from stack
+      destination = dest;
+      line_to_destination(scaled_fr_mm_s, x_splits, y_splits);
+    }
+
+  #endif // IS_CARTESIAN && !SEGMENT_LEVELED_MOVES
+
+  void mesh_bed_leveling::report_mesh() {
+    SERIAL_ECHOPAIR_F(STRINGIFY(GRID_MAX_POINTS_X) "x" STRINGIFY(GRID_MAX_POINTS_Y) " mesh. Z offset: ", z_offset, 5);
+    SERIAL_ECHOLNPGM("\nMeasured points:");
+    print_2d_array(GRID_MAX_POINTS_X, GRID_MAX_POINTS_Y, 5,
+      [](const uint8_t ix, const uint8_t iy) { return z_values[ix][iy]; }
+    );
+  }
+
+#endif // MESH_BED_LEVELING
diff --git a/Marlin/src/feature/bedlevel/mbl/mesh_bed_leveling.h b/Marlin/src/feature/bedlevel/mbl/mesh_bed_leveling.h
new file mode 100644
index 0000000..ade7a93
--- /dev/null
+++ b/Marlin/src/feature/bedlevel/mbl/mesh_bed_leveling.h
@@ -0,0 +1,127 @@
+/**
+ * Marlin 3D Printer Firmware
+ * Copyright (c) 2020 MarlinFirmware [https://github.com/MarlinFirmware/Marlin]
+ *
+ * Based on Sprinter and grbl.
+ * Copyright (c) 2011 Camiel Gubbels / Erik van der Zalm
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <https://www.gnu.org/licenses/>.
+ *
+ */
+#pragma once
+
+#include "../../../inc/MarlinConfig.h"
+
+enum MeshLevelingState : char {
+  MeshReport,     // G29 S0
+  MeshStart,      // G29 S1
+  MeshNext,       // G29 S2
+  MeshSet,        // G29 S3
+  MeshSetZOffset, // G29 S4
+  MeshReset       // G29 S5
+};
+
+#define MESH_X_DIST (float(MESH_MAX_X - (MESH_MIN_X)) / float(GRID_MAX_POINTS_X - 1))
+#define MESH_Y_DIST (float(MESH_MAX_Y - (MESH_MIN_Y)) / float(GRID_MAX_POINTS_Y - 1))
+#define _GET_MESH_X(I) mbl.index_to_xpos[I]
+#define _GET_MESH_Y(J) mbl.index_to_ypos[J]
+#define Z_VALUES_ARR mbl.z_values
+
+class mesh_bed_leveling {
+public:
+  static float z_offset,
+               z_values[GRID_MAX_POINTS_X][GRID_MAX_POINTS_Y],
+               index_to_xpos[GRID_MAX_POINTS_X],
+               index_to_ypos[GRID_MAX_POINTS_Y];
+
+  mesh_bed_leveling();
+
+  static void report_mesh();
+
+  static void reset();
+
+  FORCE_INLINE static bool has_mesh() {
+    GRID_LOOP(x, y) if (z_values[x][y]) return true;
+    return false;
+  }
+
+  static void set_z(const int8_t px, const int8_t py, const float &z) { z_values[px][py] = z; }
+
+  static inline void zigzag(const int8_t index, int8_t &px, int8_t &py) {
+    px = index % (GRID_MAX_POINTS_X);
+    py = index / (GRID_MAX_POINTS_X);
+    if (py & 1) px = (GRID_MAX_POINTS_X - 1) - px; // Zig zag
+  }
+
+  static void set_zigzag_z(const int8_t index, const float &z) {
+    int8_t px, py;
+    zigzag(index, px, py);
+    set_z(px, py, z);
+  }
+
+  static int8_t cell_index_x(const float &x) {
+    int8_t cx = (x - (MESH_MIN_X)) * RECIPROCAL(MESH_X_DIST);
+    return constrain(cx, 0, (GRID_MAX_POINTS_X) - 2);
+  }
+  static int8_t cell_index_y(const float &y) {
+    int8_t cy = (y - (MESH_MIN_Y)) * RECIPROCAL(MESH_Y_DIST);
+    return constrain(cy, 0, (GRID_MAX_POINTS_Y) - 2);
+  }
+  static inline xy_int8_t cell_indexes(const float &x, const float &y) {
+    return { cell_index_x(x), cell_index_y(y) };
+  }
+  static inline xy_int8_t cell_indexes(const xy_pos_t &xy) { return cell_indexes(xy.x, xy.y); }
+
+  static int8_t probe_index_x(const float &x) {
+    int8_t px = (x - (MESH_MIN_X) + 0.5f * (MESH_X_DIST)) * RECIPROCAL(MESH_X_DIST);
+    return WITHIN(px, 0, GRID_MAX_POINTS_X - 1) ? px : -1;
+  }
+  static int8_t probe_index_y(const float &y) {
+    int8_t py = (y - (MESH_MIN_Y) + 0.5f * (MESH_Y_DIST)) * RECIPROCAL(MESH_Y_DIST);
+    return WITHIN(py, 0, GRID_MAX_POINTS_Y - 1) ? py : -1;
+  }
+  static inline xy_int8_t probe_indexes(const float &x, const float &y) {
+    return { probe_index_x(x), probe_index_y(y) };
+  }
+  static inline xy_int8_t probe_indexes(const xy_pos_t &xy) { return probe_indexes(xy.x, xy.y); }
+
+  static float calc_z0(const float &a0, const float &a1, const float &z1, const float &a2, const float &z2) {
+    const float delta_z = (z2 - z1) / (a2 - a1),
+                delta_a = a0 - a1;
+    return z1 + delta_a * delta_z;
+  }
+
+  static float get_z(const xy_pos_t &pos
+    #if ENABLED(ENABLE_LEVELING_FADE_HEIGHT)
+      , const float &factor=1.0f
+    #endif
+  ) {
+    #if DISABLED(ENABLE_LEVELING_FADE_HEIGHT)
+      constexpr float factor = 1.0f;
+    #endif
+    const xy_int8_t ind = cell_indexes(pos);
+    const float x1 = index_to_xpos[ind.x], x2 = index_to_xpos[ind.x+1],
+                y1 = index_to_xpos[ind.y], y2 = index_to_xpos[ind.y+1],
+                z1 = calc_z0(pos.x, x1, z_values[ind.x][ind.y  ], x2, z_values[ind.x+1][ind.y  ]),
+                z2 = calc_z0(pos.x, x1, z_values[ind.x][ind.y+1], x2, z_values[ind.x+1][ind.y+1]);
+
+    return z_offset + calc_z0(pos.y, y1, z1, y2, z2) * factor;
+  }
+
+  #if IS_CARTESIAN && DISABLED(SEGMENT_LEVELED_MOVES)
+    static void line_to_destination(const feedRate_t &scaled_fr_mm_s, uint8_t x_splits=0xFF, uint8_t y_splits=0xFF);
+  #endif
+};
+
+extern mesh_bed_leveling mbl;
diff --git a/Marlin/src/feature/bedlevel/ubl/ubl.cpp b/Marlin/src/feature/bedlevel/ubl/ubl.cpp
new file mode 100644
index 0000000..b0640e5
--- /dev/null
+++ b/Marlin/src/feature/bedlevel/ubl/ubl.cpp
@@ -0,0 +1,257 @@
+/**
+ * Marlin 3D Printer Firmware
+ * Copyright (c) 2020 MarlinFirmware [https://github.com/MarlinFirmware/Marlin]
+ *
+ * Based on Sprinter and grbl.
+ * Copyright (c) 2011 Camiel Gubbels / Erik van der Zalm
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <https://www.gnu.org/licenses/>.
+ *
+ */
+
+#include "../../../inc/MarlinConfig.h"
+
+#if ENABLED(AUTO_BED_LEVELING_UBL)
+
+  #include "../bedlevel.h"
+
+  unified_bed_leveling ubl;
+
+  #include "../../../MarlinCore.h"
+  #include "../../../gcode/gcode.h"
+
+  #include "../../../module/settings.h"
+  #include "../../../module/planner.h"
+  #include "../../../module/motion.h"
+  #include "../../../module/probe.h"
+
+  #if ENABLED(EXTENSIBLE_UI)
+    #include "../../../lcd/extui/ui_api.h"
+  #endif
+
+  #include "math.h"
+
+  void unified_bed_leveling::echo_name() { SERIAL_ECHOPGM("Unified Bed Leveling"); }
+
+  void unified_bed_leveling::report_current_mesh() {
+    if (!leveling_is_valid()) return;
+    SERIAL_ECHO_MSG("  G29 I999");
+    GRID_LOOP(x, y)
+      if (!isnan(z_values[x][y])) {
+        SERIAL_ECHO_START();
+        SERIAL_ECHOPAIR("  M421 I", int(x), " J", int(y));
+        SERIAL_ECHOLNPAIR_F_P(SP_Z_STR, z_values[x][y], 4);
+        serial_delay(75); // Prevent Printrun from exploding
+      }
+  }
+
+  void unified_bed_leveling::report_state() {
+    echo_name();
+    SERIAL_ECHO_TERNARY(planner.leveling_active, " System v" UBL_VERSION " ", "", "in", "active\n");
+    serial_delay(50);
+  }
+
+  int8_t unified_bed_leveling::storage_slot;
+
+  float unified_bed_leveling::z_values[GRID_MAX_POINTS_X][GRID_MAX_POINTS_Y];
+
+  #define _GRIDPOS(A,N) (MESH_MIN_##A + N * (MESH_##A##_DIST))
+
+  const float
+  unified_bed_leveling::_mesh_index_to_xpos[GRID_MAX_POINTS_X] PROGMEM = ARRAY_N(GRID_MAX_POINTS_X,
+    _GRIDPOS(X,  0), _GRIDPOS(X,  1), _GRIDPOS(X,  2), _GRIDPOS(X,  3),
+    _GRIDPOS(X,  4), _GRIDPOS(X,  5), _GRIDPOS(X,  6), _GRIDPOS(X,  7),
+    _GRIDPOS(X,  8), _GRIDPOS(X,  9), _GRIDPOS(X, 10), _GRIDPOS(X, 11),
+    _GRIDPOS(X, 12), _GRIDPOS(X, 13), _GRIDPOS(X, 14), _GRIDPOS(X, 15)
+  ),
+  unified_bed_leveling::_mesh_index_to_ypos[GRID_MAX_POINTS_Y] PROGMEM = ARRAY_N(GRID_MAX_POINTS_Y,
+    _GRIDPOS(Y,  0), _GRIDPOS(Y,  1), _GRIDPOS(Y,  2), _GRIDPOS(Y,  3),
+    _GRIDPOS(Y,  4), _GRIDPOS(Y,  5), _GRIDPOS(Y,  6), _GRIDPOS(Y,  7),
+    _GRIDPOS(Y,  8), _GRIDPOS(Y,  9), _GRIDPOS(Y, 10), _GRIDPOS(Y, 11),
+    _GRIDPOS(Y, 12), _GRIDPOS(Y, 13), _GRIDPOS(Y, 14), _GRIDPOS(Y, 15)
+  );
+
+  volatile int16_t unified_bed_leveling::encoder_diff;
+
+  unified_bed_leveling::unified_bed_leveling() { reset(); }
+
+  void unified_bed_leveling::reset() {
+    const bool was_enabled = planner.leveling_active;
+    set_bed_leveling_enabled(false);
+    storage_slot = -1;
+    ZERO(z_values);
+    #if ENABLED(EXTENSIBLE_UI)
+      GRID_LOOP(x, y) ExtUI::onMeshUpdate(x, y, 0);
+    #endif
+    if (was_enabled) report_current_position();
+  }
+
+  void unified_bed_leveling::invalidate() {
+    set_bed_leveling_enabled(false);
+    set_all_mesh_points_to_value(NAN);
+  }
+
+  void unified_bed_leveling::set_all_mesh_points_to_value(const float value) {
+    GRID_LOOP(x, y) {
+      z_values[x][y] = value;
+      TERN_(EXTENSIBLE_UI, ExtUI::onMeshUpdate(x, y, value));
+    }
+  }
+
+  #if ENABLED(OPTIMIZED_MESH_STORAGE)
+
+    constexpr float mesh_store_scaling = 1000;
+    constexpr int16_t Z_STEPS_NAN = INT16_MAX;
+
+    void unified_bed_leveling::set_store_from_mesh(const bed_mesh_t &in_values, mesh_store_t &stored_values) {
+      auto z_to_store = [](const float &z) {
+        if (isnan(z)) return Z_STEPS_NAN;
+        const int32_t z_scaled = TRUNC(z * mesh_store_scaling);
+        if (z_scaled == Z_STEPS_NAN || !WITHIN(z_scaled, INT16_MIN, INT16_MAX))
+          return Z_STEPS_NAN; // If Z is out of range, return our custom 'NaN'
+        return int16_t(z_scaled);
+      };
+      GRID_LOOP(x, y) stored_values[x][y] = z_to_store(in_values[x][y]);
+    }
+
+    void unified_bed_leveling::set_mesh_from_store(const mesh_store_t &stored_values, bed_mesh_t &out_values) {
+      auto store_to_z = [](const int16_t z_scaled) {
+        return z_scaled == Z_STEPS_NAN ? NAN : z_scaled / mesh_store_scaling;
+      };
+      GRID_LOOP(x, y) out_values[x][y] = store_to_z(stored_values[x][y]);
+    }
+
+  #endif // OPTIMIZED_MESH_STORAGE
+
+  static void serial_echo_xy(const uint8_t sp, const int16_t x, const int16_t y) {
+    SERIAL_ECHO_SP(sp);
+    SERIAL_CHAR('(');
+    if (x < 100) { SERIAL_CHAR(' '); if (x < 10) SERIAL_CHAR(' '); }
+    SERIAL_ECHO(x);
+    SERIAL_CHAR(',');
+    if (y < 100) { SERIAL_CHAR(' '); if (y < 10) SERIAL_CHAR(' '); }
+    SERIAL_ECHO(y);
+    SERIAL_CHAR(')');
+    serial_delay(5);
+  }
+
+  static void serial_echo_column_labels(const uint8_t sp) {
+    SERIAL_ECHO_SP(7);
+    LOOP_L_N(i, GRID_MAX_POINTS_X) {
+      if (i < 10) SERIAL_CHAR(' ');
+      SERIAL_ECHO((int)i);
+      SERIAL_ECHO_SP(sp);
+    }
+    serial_delay(10);
+  }
+
+  /**
+   * Produce one of these mesh maps:
+   *   0: Human-readable
+   *   1: CSV format for spreadsheet import
+   *   2: TODO: Display on Graphical LCD
+   *   4: Compact Human-Readable
+   */
+  void unified_bed_leveling::display_map(const int map_type) {
+    const bool was = gcode.set_autoreport_paused(true);
+
+    constexpr uint8_t eachsp = 1 + 6 + 1,                           // [-3.567]
+                      twixt = eachsp * (GRID_MAX_POINTS_X) - 9 * 2; // Leading 4sp, Coordinates 9sp each
+
+    const bool human = !(map_type & 0x3), csv = map_type == 1, lcd = map_type == 2, comp = map_type & 0x4;
+
+    SERIAL_ECHOPGM("\nBed Topography Report");
+    if (human) {
+      SERIAL_ECHOLNPGM(":\n");
+      serial_echo_xy(4, MESH_MIN_X, MESH_MAX_Y);
+      serial_echo_xy(twixt, MESH_MAX_X, MESH_MAX_Y);
+      SERIAL_EOL();
+      serial_echo_column_labels(eachsp - 2);
+    }
+    else {
+      SERIAL_ECHOPGM(" for ");
+      serialprintPGM(csv ? PSTR("CSV:\n") : PSTR("LCD:\n"));
+    }
+
+    // Add XY probe offset from extruder because probe.probe_at_point() subtracts them when
+    // moving to the XY position to be measured. This ensures better agreement between
+    // the current Z position after G28 and the mesh values.
+    const xy_int8_t curr = closest_indexes(xy_pos_t(current_position) + probe.offset_xy);
+
+    if (!lcd) SERIAL_EOL();
+    for (int8_t j = GRID_MAX_POINTS_Y - 1; j >= 0; j--) {
+
+      // Row Label (J index)
+      if (human) {
+        if (j < 10) SERIAL_CHAR(' ');
+        SERIAL_ECHO(j);
+        SERIAL_ECHOPGM(" |");
+      }
+
+      // Row Values (I indexes)
+      LOOP_L_N(i, GRID_MAX_POINTS_X) {
+
+        // Opening Brace or Space
+        const bool is_current = i == curr.x && j == curr.y;
+        if (human) SERIAL_CHAR(is_current ? '[' : ' ');
+
+        // Z Value at current I, J
+        const float f = z_values[i][j];
+        if (lcd) {
+          // TODO: Display on Graphical LCD
+        }
+        else if (isnan(f))
+          serialprintPGM(human ? PSTR("  .   ") : PSTR("NAN"));
+        else if (human || csv) {
+          if (human && f >= 0.0) SERIAL_CHAR(f > 0 ? '+' : ' ');  // Space for positive ('-' for negative)
+          SERIAL_ECHO_F(f, 3);                                    // Positive: 5 digits, Negative: 6 digits
+        }
+        if (csv && i < GRID_MAX_POINTS_X - 1) SERIAL_CHAR('\t');
+
+        // Closing Brace or Space
+        if (human) SERIAL_CHAR(is_current ? ']' : ' ');
+
+        SERIAL_FLUSHTX();
+        idle_no_sleep();
+      }
+      if (!lcd) SERIAL_EOL();
+
+      // A blank line between rows (unless compact)
+      if (j && human && !comp) SERIAL_ECHOLNPGM("   |");
+    }
+
+    if (human) {
+      serial_echo_column_labels(eachsp - 2);
+      SERIAL_EOL();
+      serial_echo_xy(4, MESH_MIN_X, MESH_MIN_Y);
+      serial_echo_xy(twixt, MESH_MAX_X, MESH_MIN_Y);
+      SERIAL_EOL();
+      SERIAL_EOL();
+    }
+
+    gcode.set_autoreport_paused(was);
+  }
+
+  bool unified_bed_leveling::sanity_check() {
+    uint8_t error_flag = 0;
+
+    if (settings.calc_num_meshes() < 1) {
+      SERIAL_ECHOLNPGM("?Mesh too big for EEPROM.");
+      error_flag++;
+    }
+
+    return !!error_flag;
+  }
+
+#endif // AUTO_BED_LEVELING_UBL
diff --git a/Marlin/src/feature/bedlevel/ubl/ubl.h b/Marlin/src/feature/bedlevel/ubl/ubl.h
new file mode 100644
index 0000000..876063c
--- /dev/null
+++ b/Marlin/src/feature/bedlevel/ubl/ubl.h
@@ -0,0 +1,328 @@
+/**
+ * Marlin 3D Printer Firmware
+ * Copyright (c) 2020 MarlinFirmware [https://github.com/MarlinFirmware/Marlin]
+ *
+ * Based on Sprinter and grbl.
+ * Copyright (c) 2011 Camiel Gubbels / Erik van der Zalm
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <https://www.gnu.org/licenses/>.
+ *
+ */
+#pragma once
+
+//#define UBL_DEVEL_DEBUGGING
+
+#include "../../../module/motion.h"
+
+#define DEBUG_OUT ENABLED(DEBUG_LEVELING_FEATURE)
+#include "../../../core/debug_out.h"
+
+#define UBL_VERSION "1.01"
+#define UBL_OK false
+#define UBL_ERR true
+
+enum MeshPointType : char { INVALID, REAL, SET_IN_BITMAP };
+
+// External references
+
+struct mesh_index_pair;
+
+#define MESH_X_DIST (float(MESH_MAX_X - (MESH_MIN_X)) / float(GRID_MAX_POINTS_X - 1))
+#define MESH_Y_DIST (float(MESH_MAX_Y - (MESH_MIN_Y)) / float(GRID_MAX_POINTS_Y - 1))
+
+#if ENABLED(OPTIMIZED_MESH_STORAGE)
+  typedef int16_t mesh_store_t[GRID_MAX_POINTS_X][GRID_MAX_POINTS_Y];
+#endif
+
+class unified_bed_leveling {
+  private:
+
+    static int    g29_verbose_level,
+                  g29_phase_value,
+                  g29_repetition_cnt,
+                  g29_storage_slot,
+                  g29_map_type;
+    static bool   g29_c_flag;
+    static float  g29_card_thickness,
+                  g29_constant;
+    static xy_pos_t g29_pos;
+    static xy_bool_t xy_seen;
+
+    #if HAS_BED_PROBE
+      static int  g29_grid_size;
+    #endif
+
+    #if IS_NEWPANEL
+      static void move_z_with_encoder(const float &multiplier);
+      static float measure_point_with_encoder();
+      static float measure_business_card_thickness();
+      static void manually_probe_remaining_mesh(const xy_pos_t&, const float&, const float&, const bool) _O0;
+      static void fine_tune_mesh(const xy_pos_t &pos, const bool do_ubl_mesh_map) _O0;
+    #endif
+
+    static bool g29_parameter_parsing() _O0;
+    static void shift_mesh_height();
+    static void probe_entire_mesh(const xy_pos_t &near, const bool do_ubl_mesh_map, const bool stow_probe, const bool do_furthest) _O0;
+    static void tilt_mesh_based_on_3pts(const float &z1, const float &z2, const float &z3);
+    static void tilt_mesh_based_on_probed_grid(const bool do_ubl_mesh_map);
+    static bool smart_fill_one(const uint8_t x, const uint8_t y, const int8_t xdir, const int8_t ydir);
+    static inline bool smart_fill_one(const xy_uint8_t &pos, const xy_uint8_t &dir) {
+      return smart_fill_one(pos.x, pos.y, dir.x, dir.y);
+    }
+    static void smart_fill_mesh();
+
+    #if ENABLED(UBL_DEVEL_DEBUGGING)
+      static void g29_what_command();
+      static void g29_eeprom_dump();
+      static void g29_compare_current_mesh_to_stored_mesh();
+    #endif
+
+  public:
+
+    static void echo_name();
+    static void report_current_mesh();
+    static void report_state();
+    static void save_ubl_active_state_and_disable();
+    static void restore_ubl_active_state_and_leave();
+    static void display_map(const int) _O0;
+    static mesh_index_pair find_closest_mesh_point_of_type(const MeshPointType, const xy_pos_t&, const bool=false, MeshFlags *done_flags=nullptr) _O0;
+    static mesh_index_pair find_furthest_invalid_mesh_point() _O0;
+    static void reset();
+    static void invalidate();
+    static void set_all_mesh_points_to_value(const float value);
+    static void adjust_mesh_to_mean(const bool cflag, const float value);
+    static bool sanity_check();
+
+    static void G29() _O0;                          // O0 for no optimization
+    static void smart_fill_wlsf(const float &) _O2; // O2 gives smaller code than Os on A2560
+
+    static int8_t storage_slot;
+
+    static bed_mesh_t z_values;
+    #if ENABLED(OPTIMIZED_MESH_STORAGE)
+      static void set_store_from_mesh(const bed_mesh_t &in_values, mesh_store_t &stored_values);
+      static void set_mesh_from_store(const mesh_store_t &stored_values, bed_mesh_t &out_values);
+    #endif
+    static const float _mesh_index_to_xpos[GRID_MAX_POINTS_X],
+                       _mesh_index_to_ypos[GRID_MAX_POINTS_Y];
+
+    #if HAS_LCD_MENU
+      static bool lcd_map_control;
+      static void steppers_were_disabled();
+    #else
+      static inline void steppers_were_disabled() {}
+    #endif
+
+    static volatile int16_t encoder_diff; // Volatile because buttons may changed it at interrupt time
+
+    unified_bed_leveling();
+
+    FORCE_INLINE static void set_z(const int8_t px, const int8_t py, const float &z) { z_values[px][py] = z; }
+
+    static int8_t cell_index_x_raw(const float &x) {
+      return FLOOR((x - (MESH_MIN_X)) * RECIPROCAL(MESH_X_DIST));
+    }
+
+    static int8_t cell_index_y_raw(const float &y) {
+      return FLOOR((y - (MESH_MIN_Y)) * RECIPROCAL(MESH_Y_DIST));
+    }
+
+    static int8_t cell_index_x_valid(const float &x) {
+      return WITHIN(cell_index_x_raw(x), 0, (GRID_MAX_POINTS_X - 2));
+    }
+
+    static int8_t cell_index_y_valid(const float &y) {
+      return WITHIN(cell_index_y_raw(y), 0, (GRID_MAX_POINTS_Y - 2));
+    }
+
+    static int8_t cell_index_x(const float &x) {
+      return constrain(cell_index_x_raw(x), 0, (GRID_MAX_POINTS_X) - 2);
+    }
+
+    static int8_t cell_index_y(const float &y) {
+      return constrain(cell_index_y_raw(y), 0, (GRID_MAX_POINTS_Y) - 2);
+    }
+
+    static inline xy_int8_t cell_indexes(const float &x, const float &y) {
+      return { cell_index_x(x), cell_index_y(y) };
+    }
+    static inline xy_int8_t cell_indexes(const xy_pos_t &xy) { return cell_indexes(xy.x, xy.y); }
+
+    static int8_t closest_x_index(const float &x) {
+      const int8_t px = (x - (MESH_MIN_X) + (MESH_X_DIST) * 0.5) * RECIPROCAL(MESH_X_DIST);
+      return WITHIN(px, 0, GRID_MAX_POINTS_X - 1) ? px : -1;
+    }
+    static int8_t closest_y_index(const float &y) {
+      const int8_t py = (y - (MESH_MIN_Y) + (MESH_Y_DIST) * 0.5) * RECIPROCAL(MESH_Y_DIST);
+      return WITHIN(py, 0, GRID_MAX_POINTS_Y - 1) ? py : -1;
+    }
+    static inline xy_int8_t closest_indexes(const xy_pos_t &xy) {
+      return { closest_x_index(xy.x), closest_y_index(xy.y) };
+    }
+
+    /**
+     *                           z2   --|
+     *                 z0        |      |
+     *                  |        |      + (z2-z1)
+     *   z1             |        |      |
+     * ---+-------------+--------+--  --|
+     *   a1            a0        a2
+     *    |<---delta_a---------->|
+     *
+     *  calc_z0 is the basis for all the Mesh Based correction. It is used to
+     *  find the expected Z Height at a position between two known Z-Height locations.
+     *
+     *  It is fairly expensive with its 4 floating point additions and 2 floating point
+     *  multiplications.
+     */
+    FORCE_INLINE static float calc_z0(const float &a0, const float &a1, const float &z1, const float &a2, const float &z2) {
+      return z1 + (z2 - z1) * (a0 - a1) / (a2 - a1);
+    }
+
+    #ifdef UBL_Z_RAISE_WHEN_OFF_MESH
+      #define _UBL_OUTER_Z_RAISE UBL_Z_RAISE_WHEN_OFF_MESH
+    #else
+      #define _UBL_OUTER_Z_RAISE NAN
+    #endif
+
+    /**
+     * z_correction_for_x_on_horizontal_mesh_line is an optimization for
+     * the case where the printer is making a vertical line that only crosses horizontal mesh lines.
+     */
+    static inline float z_correction_for_x_on_horizontal_mesh_line(const float &rx0, const int x1_i, const int yi) {
+      if (!WITHIN(x1_i, 0, GRID_MAX_POINTS_X - 1) || !WITHIN(yi, 0, GRID_MAX_POINTS_Y - 1)) {
+
+        if (DEBUGGING(LEVELING)) {
+          if (WITHIN(x1_i, 0, GRID_MAX_POINTS_X - 1)) DEBUG_ECHOPGM("yi"); else DEBUG_ECHOPGM("x1_i");
+          DEBUG_ECHOLNPAIR(" out of bounds in z_correction_for_x_on_horizontal_mesh_line(rx0=", rx0, ",x1_i=", x1_i, ",yi=", yi, ")");
+        }
+
+        // The requested location is off the mesh. Return UBL_Z_RAISE_WHEN_OFF_MESH or NAN.
+        return _UBL_OUTER_Z_RAISE;
+      }
+
+      const float xratio = (rx0 - mesh_index_to_xpos(x1_i)) * RECIPROCAL(MESH_X_DIST),
+                  z1 = z_values[x1_i][yi];
+
+      return z1 + xratio * (z_values[_MIN(x1_i, GRID_MAX_POINTS_X - 2) + 1][yi] - z1); // Don't allow x1_i+1 to be past the end of the array
+                                                                                      // If it is, it is clamped to the last element of the
+                                                                                      // z_values[][] array and no correction is applied.
+    }
+
+    //
+    // See comments above for z_correction_for_x_on_horizontal_mesh_line
+    //
+    static inline float z_correction_for_y_on_vertical_mesh_line(const float &ry0, const int xi, const int y1_i) {
+      if (!WITHIN(xi, 0, GRID_MAX_POINTS_X - 1) || !WITHIN(y1_i, 0, GRID_MAX_POINTS_Y - 1)) {
+
+        if (DEBUGGING(LEVELING)) {
+          if (WITHIN(xi, 0, GRID_MAX_POINTS_X - 1)) DEBUG_ECHOPGM("y1_i"); else DEBUG_ECHOPGM("xi");
+          DEBUG_ECHOLNPAIR(" out of bounds in z_correction_for_y_on_vertical_mesh_line(ry0=", ry0, ", xi=", xi, ", y1_i=", y1_i, ")");
+        }
+
+        // The requested location is off the mesh. Return UBL_Z_RAISE_WHEN_OFF_MESH or NAN.
+        return _UBL_OUTER_Z_RAISE;
+      }
+
+      const float yratio = (ry0 - mesh_index_to_ypos(y1_i)) * RECIPROCAL(MESH_Y_DIST),
+                  z1 = z_values[xi][y1_i];
+
+      return z1 + yratio * (z_values[xi][_MIN(y1_i, GRID_MAX_POINTS_Y - 2) + 1] - z1); // Don't allow y1_i+1 to be past the end of the array
+                                                                                      // If it is, it is clamped to the last element of the
+                                                                                      // z_values[][] array and no correction is applied.
+    }
+
+    /**
+     * This is the generic Z-Correction. It works anywhere within a Mesh Cell. It first
+     * does a linear interpolation along both of the bounding X-Mesh-Lines to find the
+     * Z-Height at both ends. Then it does a linear interpolation of these heights based
+     * on the Y position within the cell.
+     */
+    static float get_z_correction(const float &rx0, const float &ry0) {
+      const int8_t cx = cell_index_x(rx0), cy = cell_index_y(ry0); // return values are clamped
+
+      /**
+       * Check if the requested location is off the mesh.  If so, and
+       * UBL_Z_RAISE_WHEN_OFF_MESH is specified, that value is returned.
+       */
+      #ifdef UBL_Z_RAISE_WHEN_OFF_MESH
+        if (!WITHIN(rx0, MESH_MIN_X, MESH_MAX_X) || !WITHIN(ry0, MESH_MIN_Y, MESH_MAX_Y))
+          return UBL_Z_RAISE_WHEN_OFF_MESH;
+      #endif
+
+      const float z1 = calc_z0(rx0,
+                               mesh_index_to_xpos(cx), z_values[cx][cy],
+                               mesh_index_to_xpos(cx + 1), z_values[_MIN(cx, GRID_MAX_POINTS_X - 2) + 1][cy]);
+
+      const float z2 = calc_z0(rx0,
+                               mesh_index_to_xpos(cx), z_values[cx][_MIN(cy, GRID_MAX_POINTS_Y - 2) + 1],
+                               mesh_index_to_xpos(cx + 1), z_values[_MIN(cx, GRID_MAX_POINTS_X - 2) + 1][_MIN(cy, GRID_MAX_POINTS_Y - 2) + 1]);
+
+      float z0 = calc_z0(ry0,
+                         mesh_index_to_ypos(cy), z1,
+                         mesh_index_to_ypos(cy + 1), z2);
+
+      if (DEBUGGING(MESH_ADJUST)) {
+        DEBUG_ECHOPAIR(" raw get_z_correction(", rx0);
+        DEBUG_CHAR(','); DEBUG_ECHO(ry0);
+        DEBUG_ECHOPAIR_F(") = ", z0, 6);
+        DEBUG_ECHOLNPAIR_F(" >>>---> ", z0, 6);
+      }
+
+      if (isnan(z0)) { // if part of the Mesh is undefined, it will show up as NAN
+        z0 = 0.0;      // in ubl.z_values[][] and propagate through the
+                       // calculations. If our correction is NAN, we throw it out
+                       // because part of the Mesh is undefined and we don't have the
+                       // information we need to complete the height correction.
+
+        if (DEBUGGING(MESH_ADJUST)) {
+          DEBUG_ECHOPAIR("??? Yikes!  NAN in get_z_correction(", rx0);
+          DEBUG_CHAR(',');
+          DEBUG_ECHO(ry0);
+          DEBUG_CHAR(')');
+          DEBUG_EOL();
+        }
+      }
+      return z0;
+    }
+    static inline float get_z_correction(const xy_pos_t &pos) { return get_z_correction(pos.x, pos.y); }
+
+    static inline float mesh_index_to_xpos(const uint8_t i) {
+      return i < GRID_MAX_POINTS_X ? pgm_read_float(&_mesh_index_to_xpos[i]) : MESH_MIN_X + i * (MESH_X_DIST);
+    }
+    static inline float mesh_index_to_ypos(const uint8_t i) {
+      return i < GRID_MAX_POINTS_Y ? pgm_read_float(&_mesh_index_to_ypos[i]) : MESH_MIN_Y + i * (MESH_Y_DIST);
+    }
+
+    #if UBL_SEGMENTED
+      static bool line_to_destination_segmented(const feedRate_t &scaled_fr_mm_s);
+    #else
+      static void line_to_destination_cartesian(const feedRate_t &scaled_fr_mm_s, const uint8_t e);
+    #endif
+
+    static inline bool mesh_is_valid() {
+      GRID_LOOP(x, y) if (isnan(z_values[x][y])) return false;
+      return true;
+    }
+
+}; // class unified_bed_leveling
+
+extern unified_bed_leveling ubl;
+
+#define _GET_MESH_X(I) ubl.mesh_index_to_xpos(I)
+#define _GET_MESH_Y(J) ubl.mesh_index_to_ypos(J)
+#define Z_VALUES_ARR ubl.z_values
+
+// Prevent debugging propagating to other files
+#include "../../../core/debug_out.h"
diff --git a/Marlin/src/feature/bedlevel/ubl/ubl_G29.cpp b/Marlin/src/feature/bedlevel/ubl/ubl_G29.cpp
new file mode 100644
index 0000000..41d2a36
--- /dev/null
+++ b/Marlin/src/feature/bedlevel/ubl/ubl_G29.cpp
@@ -0,0 +1,1783 @@
+/**
+ * Marlin 3D Printer Firmware
+ * Copyright (c) 2020 MarlinFirmware [https://github.com/MarlinFirmware/Marlin]
+ *
+ * Based on Sprinter and grbl.
+ * Copyright (c) 2011 Camiel Gubbels / Erik van der Zalm
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <https://www.gnu.org/licenses/>.
+ *
+ */
+
+#include "../../../inc/MarlinConfig.h"
+
+#if ENABLED(AUTO_BED_LEVELING_UBL)
+
+  #include "../bedlevel.h"
+
+  #include "../../../MarlinCore.h"
+  #include "../../../HAL/shared/eeprom_api.h"
+  #include "../../../libs/hex_print.h"
+  #include "../../../module/settings.h"
+  #include "../../../lcd/marlinui.h"
+  #include "../../../module/stepper.h"
+  #include "../../../module/planner.h"
+  #include "../../../module/motion.h"
+  #include "../../../module/probe.h"
+  #include "../../../gcode/gcode.h"
+  #include "../../../libs/least_squares_fit.h"
+
+  #if HAS_MULTI_HOTEND
+    #include "../../../module/tool_change.h"
+  #endif
+
+  #define DEBUG_OUT ENABLED(DEBUG_LEVELING_FEATURE)
+  #include "../../../core/debug_out.h"
+
+  #if ENABLED(EXTENSIBLE_UI)
+    #include "../../../lcd/extui/ui_api.h"
+  #endif
+
+  #include <math.h>
+
+  #define UBL_G29_P31
+
+  #if HAS_LCD_MENU
+
+    bool unified_bed_leveling::lcd_map_control = false;
+
+    void unified_bed_leveling::steppers_were_disabled() {
+      if (lcd_map_control) {
+        lcd_map_control = false;
+        ui.defer_status_screen(false);
+      }
+    }
+
+    void ubl_map_screen();
+
+  #endif
+
+  #define SIZE_OF_LITTLE_RAISE 1
+  #define BIG_RAISE_NOT_NEEDED 0
+
+  int    unified_bed_leveling::g29_verbose_level,
+         unified_bed_leveling::g29_phase_value,
+         unified_bed_leveling::g29_repetition_cnt,
+         unified_bed_leveling::g29_storage_slot = 0,
+         unified_bed_leveling::g29_map_type;
+  bool   unified_bed_leveling::g29_c_flag;
+  float  unified_bed_leveling::g29_card_thickness = 0,
+         unified_bed_leveling::g29_constant = 0;
+  xy_bool_t unified_bed_leveling::xy_seen;
+  xy_pos_t unified_bed_leveling::g29_pos;
+
+  #if HAS_BED_PROBE
+    int  unified_bed_leveling::g29_grid_size;
+  #endif
+
+  /**
+   *   G29: Unified Bed Leveling by Roxy
+   *
+   *   Parameters understood by this leveling system:
+   *
+   *   A     Activate   Activate the Unified Bed Leveling system.
+   *
+   *   B #   Business   Use the 'Business Card' mode of the Manual Probe subsystem with P2.
+   *                    Note: A non-compressible Spark Gap feeler gauge is recommended over a business card.
+   *                    In this mode of G29 P2, a business or index card is used as a shim that the nozzle can
+   *                    grab onto as it is lowered. In principle, the nozzle-bed distance is the same when the
+   *                    same resistance is felt in the shim. You can omit the numerical value on first invocation
+   *                    of G29 P2 B to measure shim thickness. Subsequent use of 'B' will apply the previously-
+   *                    measured thickness by default.
+   *
+   *   C     Continue   G29 P1 C continues the generation of a partially-constructed Mesh without invalidating
+   *                    previous measurements.
+   *
+   *   C                G29 P2 C tells the Manual Probe subsystem to not use the current nozzle
+   *                    location in its search for the closest unmeasured Mesh Point. Instead, attempt to
+   *                    start at one end of the uprobed points and Continue sequentially.
+   *
+   *                    G29 P3 C specifies the Constant for the fill. Otherwise, uses a "reasonable" value.
+   *
+   *   C     Current    G29 Z C uses the Current location (instead of bed center or nearest edge).
+   *
+   *   D     Disable    Disable the Unified Bed Leveling system.
+   *
+   *   E     Stow_probe Stow the probe after each sampled point.
+   *
+   *   F #   Fade       Fade the amount of Mesh Based Compensation over a specified height. At the
+   *                    specified height, no correction is applied and natural printer kenimatics take over. If no
+   *                    number is specified for the command, 10mm is assumed to be reasonable.
+   *
+   *   H #   Height     With P2, 'H' specifies the Height to raise the nozzle after each manual probe of the bed.
+   *                    If omitted, the nozzle will raise by Z_CLEARANCE_BETWEEN_PROBES.
+   *
+   *   H #   Offset     With P4, 'H' specifies the Offset above the mesh height to place the nozzle.
+   *                    If omitted, Z_CLEARANCE_BETWEEN_PROBES will be used.
+   *
+   *   I #   Invalidate Invalidate the specified number of Mesh Points near the given 'X' 'Y'. If X or Y are omitted,
+   *                    the nozzle location is used. If no 'I' value is given, only the point nearest to the location
+   *                    is invalidated. Use 'T' to produce a map afterward. This command is useful to invalidate a
+   *                    portion of the Mesh so it can be adjusted using other UBL tools. When attempting to invalidate
+   *                    an isolated bad mesh point, the 'T' option shows the nozzle position in the Mesh with (#). You
+   *                    can move the nozzle around and use this feature to select the center of the area (or cell) to
+   *                    invalidate.
+   *
+   *   J #   Grid       Perform a Grid Based Leveling of the current Mesh using a grid with n points on a side.
+   *                    Not specifying a grid size will invoke the 3-Point leveling function.
+   *
+   *   L     Load       Load Mesh from the previously activated location in the EEPROM.
+   *
+   *   L #   Load       Load Mesh from the specified location in the EEPROM. Set this location as activated
+   *                    for subsequent Load and Store operations.
+   *
+   *   The P or Phase commands are used for the bulk of the work to setup a Mesh. In general, your Mesh will
+   *   start off being initialized with a G29 P0 or a G29 P1. Further refinement of the Mesh happens with
+   *   each additional Phase that processes it.
+   *
+   *   P0    Phase 0    Zero Mesh Data and turn off the Mesh Compensation System. This reverts the
+   *                    3D Printer to the same state it was in before the Unified Bed Leveling Compensation
+   *                    was turned on. Setting the entire Mesh to Zero is a special case that allows
+   *                    a subsequent G or T leveling operation for backward compatibility.
+   *
+   *   P1    Phase 1    Invalidate entire Mesh and continue with automatic generation of the Mesh data using
+   *                    the Z-Probe. Usually the probe can't reach all areas that the nozzle can reach. For delta
+   *                    printers only the areas where the probe and nozzle can both reach will be automatically probed.
+   *
+   *                    Unreachable points will be handled in Phase 2 and Phase 3.
+   *
+   *                    Use 'C' to leave the previous mesh intact and automatically probe needed points. This allows you
+   *                    to invalidate parts of the Mesh but still use Automatic Probing.
+   *
+   *                    The 'X' and 'Y' parameters prioritize where to try and measure points. If omitted, the current
+   *                    probe position is used.
+   *
+   *                    Use 'T' (Topology) to generate a report of mesh generation.
+   *
+   *                    P1 will suspend Mesh generation if the controller button is held down. Note that you may need
+   *                    to press and hold the switch for several seconds if moves are underway.
+   *
+   *   P2    Phase 2    Probe unreachable points.
+   *
+   *                    Use 'H' to set the height between Mesh points. If omitted, Z_CLEARANCE_BETWEEN_PROBES is used.
+   *                    Smaller values will be quicker. Move the nozzle down till it barely touches the bed. Make sure the
+   *                    nozzle is clean and unobstructed. Use caution and move slowly. This can damage your printer!
+   *                    (Uses SIZE_OF_LITTLE_RAISE mm if the nozzle is moving less than BIG_RAISE_NOT_NEEDED mm.)
+   *
+   *                    The 'H' value can be negative if the Mesh dips in a large area. Press and hold the
+   *                    controller button to terminate the current Phase 2 command. You can then re-issue "G29 P 2"
+   *                    with an 'H' parameter more suitable for the area you're manually probing. Note that the command
+   *                    tries to start in a corner of the bed where movement will be predictable. Override the distance
+   *                    calculation location with the X and Y parameters. You can print a Mesh Map (G29 T) to see where
+   *                    the mesh is invalidated and where the nozzle needs to move to complete the command. Use 'C' to
+   *                    indicate that the search should be based on the current position.
+   *
+   *                    The 'B' parameter for this command is described above. It places the manual probe subsystem into
+   *                    Business Card mode where the thickness of a business card is measured and then used to accurately
+   *                    set the nozzle height in all manual probing for the duration of the command. A Business card can
+   *                    be used, but you'll get better results with a flexible Shim that doesn't compress. This makes it
+   *                    easier to produce similar amounts of force and get more accurate measurements. Google if you're
+   *                    not sure how to use a shim.
+   *
+   *                    The 'T' (Map) parameter helps track Mesh building progress.
+   *
+   *                    NOTE: P2 requires an LCD controller!
+   *
+   *   P3    Phase 3    Fill the unpopulated regions of the Mesh with a fixed value. There are two different paths to
+   *                    go down:
+   *
+   *                    - If a 'C' constant is specified, the closest invalid mesh points to the nozzle will be filled,
+   *                      and a repeat count can then also be specified with 'R'.
+   *
+   *                    - Leaving out 'C' invokes Smart Fill, which scans the mesh from the edges inward looking for
+   *                      invalid mesh points. Adjacent points are used to determine the bed slope. If the bed is sloped
+   *                      upward from the invalid point, it takes the value of the nearest point. If sloped downward, it's
+   *                      replaced by a value that puts all three points in a line. This version of G29 P3 is a quick, easy
+   *                      and (usually) safe way to populate unprobed mesh regions before continuing to G26 Mesh Validation
+   *                      Pattern. Note that this populates the mesh with unverified values. Pay attention and use caution.
+   *
+   *   P4    Phase 4    Fine tune the Mesh. The Delta Mesh Compensation System assumes the existence of
+   *                    an LCD Panel. It is possible to fine tune the mesh without an LCD Panel using
+   *                    G42 and M421. See the UBL documentation for further details.
+   *
+   *                    Phase 4 is meant to be used with G26 Mesh Validation to fine tune the mesh by direct editing
+   *                    of Mesh Points. Raise and lower points to fine tune the mesh until it gives consistently reliable
+   *                    adhesion.
+   *
+   *                    P4 moves to the closest Mesh Point (and/or the given X Y), raises the nozzle above the mesh height
+   *                    by the given 'H' offset (or default 0), and waits while the controller is used to adjust the nozzle
+   *                    height. On click the displayed height is saved in the mesh.
+   *
+   *                    Start Phase 4 at a specific location with X and Y. Adjust a specific number of Mesh Points with
+   *                    the 'R' (Repeat) parameter. (If 'R' is left out, the whole matrix is assumed.) This command can be
+   *                    terminated early (e.g., after editing the area of interest) by pressing and holding the encoder button.
+   *
+   *                    The general form is G29 P4 [R points] [X position] [Y position]
+   *
+   *                    The H [offset] parameter is useful if a shim is used to fine-tune the mesh. For a 0.4mm shim the
+   *                    command would be G29 P4 H0.4. The nozzle is moved to the shim height, you adjust height to the shim,
+   *                    and on click the height minus the shim thickness will be saved in the mesh.
+   *
+   *                    !!Use with caution, as a very poor mesh could cause the nozzle to crash into the bed!!
+   *
+   *                    NOTE:  P4 is not available unless you have LCD support enabled!
+   *
+   *   P5    Phase 5    Find Mean Mesh Height and Standard Deviation. Typically, it is easier to use and
+   *                    work with the Mesh if it is Mean Adjusted. You can specify a C parameter to
+   *                    Correct the Mesh to a 0.00 Mean Height. Adding a C parameter will automatically
+   *                    execute a G29 P6 C <mean height>.
+   *
+   *   P6    Phase 6    Shift Mesh height. The entire Mesh's height is adjusted by the height specified
+   *                    with the C parameter. Being able to adjust the height of a Mesh is useful tool. It
+   *                    can be used to compensate for poorly calibrated Z-Probes and other errors. Ideally,
+   *                    you should have the Mesh adjusted for a Mean Height of 0.00 and the Z-Probe measuring
+   *                    0.000 at the Z Home location.
+   *
+   *   Q     Test       Load specified Test Pattern to assist in checking correct operation of system. This
+   *                    command is not anticipated to be of much value to the typical user. It is intended
+   *                    for developers to help them verify correct operation of the Unified Bed Leveling System.
+   *
+   *   R #   Repeat     Repeat this command the specified number of times. If no number is specified the
+   *                    command will be repeated GRID_MAX_POINTS_X * GRID_MAX_POINTS_Y times.
+   *
+   *   S     Store      Store the current Mesh in the Activated area of the EEPROM. It will also store the
+   *                    current state of the Unified Bed Leveling system in the EEPROM.
+   *
+   *   S #   Store      Store the current Mesh at the specified location in EEPROM. Activate this location
+   *                    for subsequent Load and Store operations. Valid storage slot numbers begin at 0 and
+   *                    extend to a limit related to the available EEPROM storage.
+   *
+   *   S -1  Store      Print the current Mesh as G-code that can be used to restore the mesh anytime.
+   *
+   *   T     Topology   Display the Mesh Map Topology.
+   *                    'T' can be used alone (e.g., G29 T) or in combination with most of the other commands.
+   *                    This option works with all Phase commands (e.g., G29 P4 R 5 T X 50 Y100 C -.1 O)
+   *                    This parameter can also specify a Map Type. T0 (the default) is user-readable. T1
+   *                    is suitable to paste into a spreadsheet for a 3D graph of the mesh.
+   *
+   *   U     Unlevel    Perform a probe of the outer perimeter to assist in physically leveling unlevel beds.
+   *                    Only used for G29 P1 T U. This speeds up the probing of the edge of the bed. Useful
+   *                    when the entire bed doesn't need to be probed because it will be adjusted.
+   *
+   *   V #   Verbosity  Set the verbosity level (0-4) for extra details. (Default 0)
+   *
+   *   X #              X Location for this command
+   *
+   *   Y #              Y Location for this command
+   *
+   * With UBL_DEVEL_DEBUGGING:
+   *
+   *   K #  Kompare     Kompare current Mesh with stored Mesh #, replacing current Mesh with the result.
+   *                    This command literally performs a diff between two Meshes.
+   *
+   *   Q-1  Dump EEPROM Dump the UBL contents stored in EEPROM as HEX format. Useful for developers to help
+   *                    verify correct operation of the UBL.
+   *
+   *   W    What?       Display valuable UBL data.
+   *
+   *
+   *   Release Notes:
+   *   You MUST do M502, M500 to initialize the storage. Failure to do this will cause all
+   *   kinds of problems. Enabling EEPROM Storage is required.
+   *
+   *   When you do a G28 and G29 P1 to automatically build your first mesh, you are going to notice that
+   *   UBL probes points increasingly further from the starting location. (The starting location defaults
+   *   to the center of the bed.) In contrast, ABL and MBL follow a zigzag pattern. The spiral pattern is
+   *   especially better for Delta printers, since it populates the center of the mesh first, allowing for
+   *   a quicker test print to verify settings. You don't need to populate the entire mesh to use it.
+   *   After all, you don't want to spend a lot of time generating a mesh only to realize the resolution
+   *   or probe offsets are incorrect. Mesh-generation gathers points starting closest to the nozzle unless
+   *   an (X,Y) coordinate pair is given.
+   *
+   *   Unified Bed Leveling uses a lot of EEPROM storage to hold its data, and it takes some effort to get
+   *   the mesh just right. To prevent this valuable data from being destroyed as the EEPROM structure
+   *   evolves, UBL stores all mesh data at the end of EEPROM.
+   *
+   *   UBL is founded on Edward Patel's Mesh Bed Leveling code. A big 'Thanks!' to him and the creators of
+   *   3-Point and Grid Based leveling. Combining their contributions we now have the functionality and
+   *   features of all three systems combined.
+   */
+
+  void unified_bed_leveling::G29() {
+
+    bool probe_deployed = false;
+    if (g29_parameter_parsing()) return; // Abort on parameter error
+
+    const int8_t p_val = parser.intval('P', -1);
+    const bool may_move = p_val == 1 || p_val == 2 || p_val == 4 || parser.seen('J');
+    TERN_(HAS_MULTI_HOTEND, const uint8_t old_tool_index = active_extruder);
+
+    // Check for commands that require the printer to be homed
+    if (may_move) {
+      planner.synchronize();
+      // Send 'N' to force homing before G29 (internal only)
+      if (axes_should_home() || parser.seen('N')) gcode.home_all_axes();
+      TERN_(HAS_MULTI_HOTEND, if (active_extruder) tool_change(0));
+    }
+
+    // Invalidate Mesh Points. This command is a little bit asymmetrical because
+    // it directly specifies the repetition count and does not use the 'R' parameter.
+    if (parser.seen('I')) {
+      uint8_t cnt = 0;
+      g29_repetition_cnt = parser.has_value() ? parser.value_int() : 1;
+      if (g29_repetition_cnt >= GRID_MAX_POINTS) {
+        set_all_mesh_points_to_value(NAN);
+      }
+      else {
+        while (g29_repetition_cnt--) {
+          if (cnt > 20) { cnt = 0; idle(); }
+          const mesh_index_pair closest = find_closest_mesh_point_of_type(REAL, g29_pos);
+          const xy_int8_t &cpos = closest.pos;
+          if (cpos.x < 0) {
+            // No more REAL mesh points to invalidate, so we ASSUME the user
+            // meant to invalidate the ENTIRE mesh, which cannot be done with
+            // find_closest_mesh_point loop which only returns REAL points.
+            set_all_mesh_points_to_value(NAN);
+            SERIAL_ECHOLNPGM("Entire Mesh invalidated.\n");
+            break;            // No more invalid Mesh Points to populate
+          }
+          z_values[cpos.x][cpos.y] = NAN;
+          TERN_(EXTENSIBLE_UI, ExtUI::onMeshUpdate(cpos, 0.0f));
+          cnt++;
+        }
+      }
+      SERIAL_ECHOLNPGM("Locations invalidated.\n");
+    }
+
+    if (parser.seen('Q')) {
+      const int test_pattern = parser.has_value() ? parser.value_int() : -99;
+      if (!WITHIN(test_pattern, -1, 2)) {
+        SERIAL_ECHOLNPGM("Invalid test_pattern value. (-1 to 2)\n");
+        return;
+      }
+      SERIAL_ECHOLNPGM("Loading test_pattern values.\n");
+      switch (test_pattern) {
+
+        #if ENABLED(UBL_DEVEL_DEBUGGING)
+          case -1:
+            g29_eeprom_dump();
+            break;
+        #endif
+
+        case 0:
+          GRID_LOOP(x, y) {                                     // Create a bowl shape similar to a poorly-calibrated Delta
+            const float p1 = 0.5f * (GRID_MAX_POINTS_X) - x,
+                        p2 = 0.5f * (GRID_MAX_POINTS_Y) - y;
+            z_values[x][y] += 2.0f * HYPOT(p1, p2);
+            TERN_(EXTENSIBLE_UI, ExtUI::onMeshUpdate(x, y, z_values[x][y]));
+          }
+          break;
+
+        case 1:
+          LOOP_L_N(x, GRID_MAX_POINTS_X) {                     // Create a diagonal line several Mesh cells thick that is raised
+            z_values[x][x] += 9.999f;
+            z_values[x][x + (x < (GRID_MAX_POINTS_Y) - 1) ? 1 : -1] += 9.999f; // We want the altered line several mesh points thick
+            #if ENABLED(EXTENSIBLE_UI)
+              ExtUI::onMeshUpdate(x, x, z_values[x][x]);
+              ExtUI::onMeshUpdate(x, (x + (x < (GRID_MAX_POINTS_Y) - 1) ? 1 : -1), z_values[x][x + (x < (GRID_MAX_POINTS_Y) - 1) ? 1 : -1]);
+            #endif
+
+          }
+          break;
+
+        case 2:
+          // Allow the user to specify the height because 10mm is a little extreme in some cases.
+          for (uint8_t x = (GRID_MAX_POINTS_X) / 3; x < 2 * (GRID_MAX_POINTS_X) / 3; x++)     // Create a rectangular raised area in
+            for (uint8_t y = (GRID_MAX_POINTS_Y) / 3; y < 2 * (GRID_MAX_POINTS_Y) / 3; y++) { // the center of the bed
+              z_values[x][y] += parser.seen('C') ? g29_constant : 9.99f;
+              TERN_(EXTENSIBLE_UI, ExtUI::onMeshUpdate(x, y, z_values[x][y]));
+            }
+          break;
+      }
+    }
+
+    #if HAS_BED_PROBE
+
+      if (parser.seen('J')) {
+        save_ubl_active_state_and_disable();
+        tilt_mesh_based_on_probed_grid(g29_grid_size == 0); // Zero size does 3-Point
+        restore_ubl_active_state_and_leave();
+        #if ENABLED(UBL_G29_J_RECENTER)
+          do_blocking_move_to_xy(0.5f * ((MESH_MIN_X) + (MESH_MAX_X)), 0.5f * ((MESH_MIN_Y) + (MESH_MAX_Y)));
+        #endif
+        report_current_position();
+        probe_deployed = true;
+      }
+
+    #endif // HAS_BED_PROBE
+
+    if (parser.seen('P')) {
+      if (WITHIN(g29_phase_value, 0, 1) && storage_slot == -1) {
+        storage_slot = 0;
+        SERIAL_ECHOLNPGM("Default storage slot 0 selected.");
+      }
+
+      switch (g29_phase_value) {
+        case 0:
+          //
+          // Zero Mesh Data
+          //
+          reset();
+          SERIAL_ECHOLNPGM("Mesh zeroed.");
+          break;
+
+        #if HAS_BED_PROBE
+
+          case 1: {
+            //
+            // Invalidate Entire Mesh and Automatically Probe Mesh in areas that can be reached by the probe
+            //
+            if (!parser.seen('C')) {
+              invalidate();
+              SERIAL_ECHOLNPGM("Mesh invalidated. Probing mesh.");
+            }
+            if (g29_verbose_level > 1) {
+              SERIAL_ECHOPAIR("Probing around (", g29_pos.x);
+              SERIAL_CHAR(',');
+              SERIAL_DECIMAL(g29_pos.y);
+              SERIAL_ECHOLNPGM(").\n");
+            }
+            const xy_pos_t near_probe_xy = g29_pos + probe.offset_xy;
+            probe_entire_mesh(near_probe_xy, parser.seen('T'), parser.seen('E'), parser.seen('U'));
+
+            report_current_position();
+            probe_deployed = true;
+          } break;
+
+        #endif // HAS_BED_PROBE
+
+        case 2: {
+          #if HAS_LCD_MENU
+            //
+            // Manually Probe Mesh in areas that can't be reached by the probe
+            //
+            SERIAL_ECHOLNPGM("Manually probing unreachable points.");
+            do_z_clearance(Z_CLEARANCE_BETWEEN_PROBES);
+
+            if (parser.seen('C') && !xy_seen) {
+
+              /**
+               * Use a good default location for the path.
+               * The flipped > and < operators in these comparisons is intentional.
+               * It should cause the probed points to follow a nice path on Cartesian printers.
+               * It may make sense to have Delta printers default to the center of the bed.
+               * Until that is decided, this can be forced with the X and Y parameters.
+               */
+              g29_pos.set(
+                #if IS_KINEMATIC
+                  X_HOME_POS, Y_HOME_POS
+                #else
+                  probe.offset_xy.x > 0 ? X_BED_SIZE : 0,
+                  probe.offset_xy.y < 0 ? Y_BED_SIZE : 0
+                #endif
+              );
+            }
+
+            if (parser.seen('B')) {
+              g29_card_thickness = parser.has_value() ? parser.value_float() : measure_business_card_thickness();
+              if (ABS(g29_card_thickness) > 1.5f) {
+                SERIAL_ECHOLNPGM("?Error in Business Card measurement.");
+                return;
+              }
+              probe_deployed = true;
+            }
+
+            if (!position_is_reachable(g29_pos)) {
+              SERIAL_ECHOLNPGM("XY outside printable radius.");
+              return;
+            }
+
+            const float height = parser.floatval('H', Z_CLEARANCE_BETWEEN_PROBES);
+            manually_probe_remaining_mesh(g29_pos, height, g29_card_thickness, parser.seen('T'));
+
+            SERIAL_ECHOLNPGM("G29 P2 finished.");
+
+            report_current_position();
+
+          #else
+
+            SERIAL_ECHOLNPGM("?P2 is only available when an LCD is present.");
+            return;
+
+          #endif
+        } break;
+
+        case 3: {
+          /**
+           * Populate invalid mesh areas. Proceed with caution.
+           * Two choices are available:
+           *   - Specify a constant with the 'C' parameter.
+           *   - Allow 'G29 P3' to choose a 'reasonable' constant.
+           */
+
+          if (g29_c_flag) {
+            if (g29_repetition_cnt >= GRID_MAX_POINTS) {
+              set_all_mesh_points_to_value(g29_constant);
+            }
+            else {
+              while (g29_repetition_cnt--) {  // this only populates reachable mesh points near
+                const mesh_index_pair closest = find_closest_mesh_point_of_type(INVALID, g29_pos);
+                const xy_int8_t &cpos = closest.pos;
+                if (cpos.x < 0) {
+                  // No more REAL INVALID mesh points to populate, so we ASSUME
+                  // user meant to populate ALL INVALID mesh points to value
+                  GRID_LOOP(x, y) if (isnan(z_values[x][y])) z_values[x][y] = g29_constant;
+                  break; // No more invalid Mesh Points to populate
+                }
+                else {
+                  z_values[cpos.x][cpos.y] = g29_constant;
+                  TERN_(EXTENSIBLE_UI, ExtUI::onMeshUpdate(cpos, g29_constant));
+                }
+              }
+            }
+          }
+          else {
+            const float cvf = parser.value_float();
+            switch ((int)TRUNC(cvf * 10.0f) - 30) {   // 3.1 -> 1
+              #if ENABLED(UBL_G29_P31)
+                case 1: {
+
+                  // P3.1  use least squares fit to fill missing mesh values
+                  // P3.10 zero weighting for distance, all grid points equal, best fit tilted plane
+                  // P3.11 10X weighting for nearest grid points versus farthest grid points
+                  // P3.12 100X distance weighting
+                  // P3.13 1000X distance weighting, approaches simple average of nearest points
+
+                  const float weight_power  = (cvf - 3.10f) * 100.0f,  // 3.12345 -> 2.345
+                              weight_factor = weight_power ? POW(10.0f, weight_power) : 0;
+                  smart_fill_wlsf(weight_factor);
+                }
+                break;
+              #endif
+              case 0:   // P3 or P3.0
+              default:  // and anything P3.x that's not P3.1
+                smart_fill_mesh();  // Do a 'Smart' fill using nearby known values
+                break;
+            }
+          }
+          break;
+        }
+
+        case 4: // Fine Tune (i.e., Edit) the Mesh
+          #if HAS_LCD_MENU
+            fine_tune_mesh(g29_pos, parser.seen('T'));
+          #else
+            SERIAL_ECHOLNPGM("?P4 is only available when an LCD is present.");
+            return;
+          #endif
+          break;
+
+        case 5: adjust_mesh_to_mean(g29_c_flag, g29_constant); break;
+
+        case 6: shift_mesh_height(); break;
+      }
+    }
+
+    #if ENABLED(UBL_DEVEL_DEBUGGING)
+
+      //
+      // Much of the 'What?' command can be eliminated. But until we are fully debugged, it is
+      // good to have the extra information. Soon... we prune this to just a few items
+      //
+      if (parser.seen('W')) g29_what_command();
+
+      //
+      // When we are fully debugged, this may go away. But there are some valid
+      // use cases for the users. So we can wait and see what to do with it.
+      //
+
+      if (parser.seen('K')) // Kompare Current Mesh Data to Specified Stored Mesh
+        g29_compare_current_mesh_to_stored_mesh();
+
+    #endif // UBL_DEVEL_DEBUGGING
+
+
+    //
+    // Load a Mesh from the EEPROM
+    //
+
+    if (parser.seen('L')) {     // Load Current Mesh Data
+      g29_storage_slot = parser.has_value() ? parser.value_int() : storage_slot;
+
+      int16_t a = settings.calc_num_meshes();
+
+      if (!a) {
+        SERIAL_ECHOLNPGM("?EEPROM storage not available.");
+        return;
+      }
+
+      if (!WITHIN(g29_storage_slot, 0, a - 1)) {
+        SERIAL_ECHOLNPAIR("?Invalid storage slot.\n?Use 0 to ", a - 1);
+        return;
+      }
+
+      settings.load_mesh(g29_storage_slot);
+      storage_slot = g29_storage_slot;
+
+      SERIAL_ECHOLNPGM("Done.");
+    }
+
+    //
+    // Store a Mesh in the EEPROM
+    //
+
+    if (parser.seen('S')) {     // Store (or Save) Current Mesh Data
+      g29_storage_slot = parser.has_value() ? parser.value_int() : storage_slot;
+
+      if (g29_storage_slot == -1)                     // Special case, the user wants to 'Export' the mesh to the
+        return report_current_mesh();                 // host program to be saved on the user's computer
+
+      int16_t a = settings.calc_num_meshes();
+
+      if (!a) {
+        SERIAL_ECHOLNPGM("?EEPROM storage not available.");
+        goto LEAVE;
+      }
+
+      if (!WITHIN(g29_storage_slot, 0, a - 1)) {
+        SERIAL_ECHOLNPAIR("?Invalid storage slot.\n?Use 0 to ", a - 1);
+        goto LEAVE;
+      }
+
+      settings.store_mesh(g29_storage_slot);
+      storage_slot = g29_storage_slot;
+
+      SERIAL_ECHOLNPGM("Done.");
+    }
+
+    if (parser.seen('T'))
+      display_map(g29_map_type);
+
+    LEAVE:
+
+    #if HAS_LCD_MENU
+      ui.reset_alert_level();
+      ui.quick_feedback();
+      ui.reset_status();
+      ui.release();
+    #endif
+
+    #ifdef Z_PROBE_END_SCRIPT
+      if (DEBUGGING(LEVELING)) DEBUG_ECHOLNPAIR("Z Probe End Script: ", Z_PROBE_END_SCRIPT);
+      if (probe_deployed) {
+        planner.synchronize();
+        gcode.process_subcommands_now_P(PSTR(Z_PROBE_END_SCRIPT));
+      }
+    #else
+      UNUSED(probe_deployed);
+    #endif
+
+    TERN_(HAS_MULTI_HOTEND, tool_change(old_tool_index));
+    return;
+  }
+
+  void unified_bed_leveling::adjust_mesh_to_mean(const bool cflag, const float value) {
+    float sum = 0;
+    int n = 0;
+    GRID_LOOP(x, y)
+      if (!isnan(z_values[x][y])) {
+        sum += z_values[x][y];
+        n++;
+      }
+
+    const float mean = sum / n;
+
+    //
+    // Sum the squares of difference from mean
+    //
+    float sum_of_diff_squared = 0;
+    GRID_LOOP(x, y)
+      if (!isnan(z_values[x][y]))
+        sum_of_diff_squared += sq(z_values[x][y] - mean);
+
+    SERIAL_ECHOLNPAIR("# of samples: ", n);
+    SERIAL_ECHOLNPAIR_F("Mean Mesh Height: ", mean, 6);
+
+    const float sigma = SQRT(sum_of_diff_squared / (n + 1));
+    SERIAL_ECHOLNPAIR_F("Standard Deviation: ", sigma, 6);
+
+    if (cflag)
+      GRID_LOOP(x, y)
+        if (!isnan(z_values[x][y])) {
+          z_values[x][y] -= mean + value;
+          TERN_(EXTENSIBLE_UI, ExtUI::onMeshUpdate(x, y, z_values[x][y]));
+        }
+  }
+
+  void unified_bed_leveling::shift_mesh_height() {
+    GRID_LOOP(x, y)
+      if (!isnan(z_values[x][y])) {
+        z_values[x][y] += g29_constant;
+        TERN_(EXTENSIBLE_UI, ExtUI::onMeshUpdate(x, y, z_values[x][y]));
+      }
+  }
+
+  #if HAS_BED_PROBE
+    /**
+     * Probe all invalidated locations of the mesh that can be reached by the probe.
+     * This attempts to fill in locations closest to the nozzle's start location first.
+     */
+    void unified_bed_leveling::probe_entire_mesh(const xy_pos_t &nearby, const bool do_ubl_mesh_map, const bool stow_probe, const bool do_furthest) {
+      probe.deploy(); // Deploy before ui.capture() to allow for PAUSE_BEFORE_DEPLOY_STOW
+
+      TERN_(HAS_LCD_MENU, ui.capture());
+
+      save_ubl_active_state_and_disable();  // No bed level correction so only raw data is obtained
+      uint8_t count = GRID_MAX_POINTS;
+
+      mesh_index_pair best;
+      TERN_(EXTENSIBLE_UI, ExtUI::onMeshUpdate(best.pos, ExtUI::MESH_START));
+      do {
+        if (do_ubl_mesh_map) display_map(g29_map_type);
+
+        const int point_num = (GRID_MAX_POINTS) - count + 1;
+        SERIAL_ECHOLNPAIR("Probing mesh point ", point_num, "/", int(GRID_MAX_POINTS), ".");
+        TERN_(HAS_DISPLAY, ui.status_printf_P(0, PSTR(S_FMT " %i/%i"), GET_TEXT(MSG_PROBING_MESH), point_num, int(GRID_MAX_POINTS)));
+
+        #if HAS_LCD_MENU
+          if (ui.button_pressed()) {
+            ui.quick_feedback(false); // Preserve button state for click-and-hold
+            SERIAL_ECHOLNPGM("\nMesh only partially populated.\n");
+            ui.wait_for_release();
+            ui.quick_feedback();
+            ui.release();
+            probe.stow(); // Release UI before stow to allow for PAUSE_BEFORE_DEPLOY_STOW
+            return restore_ubl_active_state_and_leave();
+          }
+        #endif
+
+        best = do_furthest
+          ? find_furthest_invalid_mesh_point()
+          : find_closest_mesh_point_of_type(INVALID, nearby, true);
+
+        if (best.pos.x >= 0) {    // mesh point found and is reachable by probe
+          TERN_(EXTENSIBLE_UI, ExtUI::onMeshUpdate(best.pos, ExtUI::PROBE_START));
+          const float measured_z = probe.probe_at_point(
+                        best.meshpos(),
+                        stow_probe ? PROBE_PT_STOW : PROBE_PT_RAISE, g29_verbose_level
+                      );
+          z_values[best.pos.x][best.pos.y] = measured_z;
+          #if ENABLED(EXTENSIBLE_UI)
+            ExtUI::onMeshUpdate(best.pos, ExtUI::PROBE_FINISH);
+            ExtUI::onMeshUpdate(best.pos, measured_z);
+          #endif
+        }
+        SERIAL_FLUSH(); // Prevent host M105 buffer overrun.
+
+      } while (best.pos.x >= 0 && --count);
+
+      TERN_(EXTENSIBLE_UI, ExtUI::onMeshUpdate(best.pos, ExtUI::MESH_FINISH));
+
+      // Release UI during stow to allow for PAUSE_BEFORE_DEPLOY_STOW
+      TERN_(HAS_LCD_MENU, ui.release());
+      probe.stow();
+      TERN_(HAS_LCD_MENU, ui.capture());
+
+      probe.move_z_after_probing();
+
+      restore_ubl_active_state_and_leave();
+
+      do_blocking_move_to_xy(
+        constrain(nearby.x - probe.offset_xy.x, MESH_MIN_X, MESH_MAX_X),
+        constrain(nearby.y - probe.offset_xy.y, MESH_MIN_Y, MESH_MAX_Y)
+      );
+    }
+
+  #endif // HAS_BED_PROBE
+
+  #if HAS_LCD_MENU
+
+    typedef void (*clickFunc_t)();
+
+    bool click_and_hold(const clickFunc_t func=nullptr) {
+      if (ui.button_pressed()) {
+        ui.quick_feedback(false);         // Preserve button state for click-and-hold
+        const millis_t nxt = millis() + 1500UL;
+        while (ui.button_pressed()) {     // Loop while the encoder is pressed. Uses hardware flag!
+          idle();                         // idle, of course
+          if (ELAPSED(millis(), nxt)) {   // After 1.5 seconds
+            ui.quick_feedback();
+            if (func) (*func)();
+            ui.wait_for_release();
+            return true;
+          }
+        }
+      }
+      serial_delay(15);
+      return false;
+    }
+
+    void unified_bed_leveling::move_z_with_encoder(const float &multiplier) {
+      ui.wait_for_release();
+      while (!ui.button_pressed()) {
+        idle();
+        gcode.reset_stepper_timeout(); // Keep steppers powered
+        if (encoder_diff) {
+          do_blocking_move_to_z(current_position.z + float(encoder_diff) * multiplier);
+          encoder_diff = 0;
+        }
+      }
+    }
+
+    float unified_bed_leveling::measure_point_with_encoder() {
+      KEEPALIVE_STATE(PAUSED_FOR_USER);
+      move_z_with_encoder(0.01f);
+      return current_position.z;
+    }
+
+    static void echo_and_take_a_measurement() { SERIAL_ECHOLNPGM(" and take a measurement."); }
+
+    float unified_bed_leveling::measure_business_card_thickness() {
+      ui.capture();
+      save_ubl_active_state_and_disable();   // Disable bed level correction for probing
+
+      do_blocking_move_to(0.5f * (MESH_MAX_X - (MESH_MIN_X)), 0.5f * (MESH_MAX_Y - (MESH_MIN_Y)), MANUAL_PROBE_START_Z);
+        //, _MIN(planner.settings.max_feedrate_mm_s[X_AXIS], planner.settings.max_feedrate_mm_s[Y_AXIS]) * 0.5f);
+      planner.synchronize();
+
+      SERIAL_ECHOPGM("Place shim under nozzle");
+      LCD_MESSAGEPGM(MSG_UBL_BC_INSERT);
+      ui.return_to_status();
+      echo_and_take_a_measurement();
+
+      const float z1 = measure_point_with_encoder();
+      do_blocking_move_to_z(current_position.z + SIZE_OF_LITTLE_RAISE);
+      planner.synchronize();
+
+      SERIAL_ECHOPGM("Remove shim");
+      LCD_MESSAGEPGM(MSG_UBL_BC_REMOVE);
+      echo_and_take_a_measurement();
+
+      const float z2 = measure_point_with_encoder();
+      do_blocking_move_to_z(current_position.z + Z_CLEARANCE_BETWEEN_PROBES);
+
+      const float thickness = ABS(z1 - z2);
+
+      if (g29_verbose_level > 1) {
+        SERIAL_ECHOPAIR_F("Business Card is ", thickness, 4);
+        SERIAL_ECHOLNPGM("mm thick.");
+      }
+
+      restore_ubl_active_state_and_leave();
+
+      return thickness;
+    }
+
+    void unified_bed_leveling::manually_probe_remaining_mesh(const xy_pos_t &pos, const float &z_clearance, const float &thick, const bool do_ubl_mesh_map) {
+      ui.capture();
+
+      save_ubl_active_state_and_disable();  // No bed level correction so only raw data is obtained
+      do_blocking_move_to_xy_z(current_position, z_clearance);
+
+      ui.return_to_status();
+
+      mesh_index_pair location;
+      const xy_int8_t &lpos = location.pos;
+      do {
+        location = find_closest_mesh_point_of_type(INVALID, pos);
+        // It doesn't matter if the probe can't reach the NAN location. This is a manual probe.
+        if (!location.valid()) continue;
+
+        const xyz_pos_t ppos = {
+          mesh_index_to_xpos(lpos.x),
+          mesh_index_to_ypos(lpos.y),
+          Z_CLEARANCE_BETWEEN_PROBES
+        };
+
+        if (!position_is_reachable(ppos)) break; // SHOULD NOT OCCUR (find_closest_mesh_point only returns reachable points)
+
+        LCD_MESSAGEPGM(MSG_UBL_MOVING_TO_NEXT);
+
+        do_blocking_move_to(ppos);
+        do_z_clearance(z_clearance);
+
+        KEEPALIVE_STATE(PAUSED_FOR_USER);
+        ui.capture();
+
+        if (do_ubl_mesh_map) display_map(g29_map_type);  // show user where we're probing
+
+        serialprintPGM(parser.seen('B') ? GET_TEXT(MSG_UBL_BC_INSERT) : GET_TEXT(MSG_UBL_BC_INSERT2));
+
+        const float z_step = 0.01f;                         // existing behavior: 0.01mm per click, occasionally step
+        //const float z_step = planner.steps_to_mm[Z_AXIS]; // approx one step each click
+
+        move_z_with_encoder(z_step);
+
+        if (click_and_hold()) {
+          SERIAL_ECHOLNPGM("\nMesh only partially populated.");
+          do_z_clearance(Z_CLEARANCE_DEPLOY_PROBE);
+          return restore_ubl_active_state_and_leave();
+        }
+
+        z_values[lpos.x][lpos.y] = current_position.z - thick;
+        TERN_(EXTENSIBLE_UI, ExtUI::onMeshUpdate(location, z_values[lpos.x][lpos.y]));
+
+        if (g29_verbose_level > 2)
+          SERIAL_ECHOLNPAIR_F("Mesh Point Measured at: ", z_values[lpos.x][lpos.y], 6);
+        SERIAL_FLUSH(); // Prevent host M105 buffer overrun.
+      } while (location.valid());
+
+      if (do_ubl_mesh_map) display_map(g29_map_type);  // show user where we're probing
+
+      restore_ubl_active_state_and_leave();
+      do_blocking_move_to_xy_z(pos, Z_CLEARANCE_DEPLOY_PROBE);
+    }
+
+    inline void set_message_with_feedback(PGM_P const msg_P) {
+      ui.set_status_P(msg_P);
+      ui.quick_feedback();
+    }
+
+    void abort_fine_tune() {
+      ui.return_to_status();
+      do_z_clearance(Z_CLEARANCE_BETWEEN_PROBES);
+      set_message_with_feedback(GET_TEXT(MSG_EDITING_STOPPED));
+    }
+
+    void unified_bed_leveling::fine_tune_mesh(const xy_pos_t &pos, const bool do_ubl_mesh_map) {
+      if (!parser.seen('R'))      // fine_tune_mesh() is special. If no repetition count flag is specified
+        g29_repetition_cnt = 1;   // do exactly one mesh location. Otherwise use what the parser decided.
+
+      #if ENABLED(UBL_MESH_EDIT_MOVES_Z)
+        const float h_offset = parser.seenval('H') ? parser.value_linear_units() : MANUAL_PROBE_START_Z;
+        if (!WITHIN(h_offset, 0, 10)) {
+          SERIAL_ECHOLNPGM("Offset out of bounds. (0 to 10mm)\n");
+          return;
+        }
+      #endif
+
+      mesh_index_pair location;
+
+      if (!position_is_reachable(pos)) {
+        SERIAL_ECHOLNPGM("(X,Y) outside printable radius.");
+        return;
+      }
+
+      save_ubl_active_state_and_disable();
+
+      LCD_MESSAGEPGM(MSG_UBL_FINE_TUNE_MESH);
+      ui.capture();                                               // Take over control of the LCD encoder
+
+      do_blocking_move_to_xy_z(pos, Z_CLEARANCE_BETWEEN_PROBES);  // Move to the given XY with probe clearance
+
+      MeshFlags done_flags{0};
+      const xy_int8_t &lpos = location.pos;
+
+      #if IS_TFTGLCD_PANEL
+        lcd_mesh_edit_setup(0);                             // Change current screen before calling ui.ubl_plot
+        safe_delay(50);
+      #endif
+
+      do {
+        location = find_closest_mesh_point_of_type(SET_IN_BITMAP, pos, false, &done_flags);
+
+        if (lpos.x < 0) break;                              // Stop when there are no more reachable points
+
+        done_flags.mark(lpos);                              // Mark this location as 'adjusted' so a new
+                                                            // location is used on the next loop
+        const xyz_pos_t raw = {
+          mesh_index_to_xpos(lpos.x),
+          mesh_index_to_ypos(lpos.y),
+          Z_CLEARANCE_BETWEEN_PROBES
+        };
+
+        if (!position_is_reachable(raw)) break;             // SHOULD NOT OCCUR (find_closest_mesh_point_of_type only returns reachable)
+
+        do_blocking_move_to(raw);                           // Move the nozzle to the edit point with probe clearance
+
+        TERN_(UBL_MESH_EDIT_MOVES_Z, do_blocking_move_to_z(h_offset)); // Move Z to the given 'H' offset before editing
+
+        KEEPALIVE_STATE(PAUSED_FOR_USER);
+
+        if (do_ubl_mesh_map) display_map(g29_map_type);     // Display the current point
+
+        #if IS_TFTGLCD_PANEL
+          ui.ubl_plot(lpos.x, lpos.y);   // update plot screen
+        #endif
+
+        ui.refresh();
+
+        float new_z = z_values[lpos.x][lpos.y];
+        if (isnan(new_z)) new_z = 0;                        // Invalid points begin at 0
+        new_z = FLOOR(new_z * 1000) * 0.001f;               // Chop off digits after the 1000ths place
+
+        lcd_mesh_edit_setup(new_z);
+
+        SET_SOFT_ENDSTOP_LOOSE(true);
+
+        do {
+          idle();
+          new_z = lcd_mesh_edit();
+          TERN_(UBL_MESH_EDIT_MOVES_Z, do_blocking_move_to_z(h_offset + new_z)); // Move the nozzle as the point is edited
+          SERIAL_FLUSH();                                   // Prevent host M105 buffer overrun.
+        } while (!ui.button_pressed());
+
+        SET_SOFT_ENDSTOP_LOOSE(false);
+
+        if (!lcd_map_control) ui.return_to_status();        // Just editing a single point? Return to status
+
+        if (click_and_hold(abort_fine_tune)) break;         // Button held down? Abort editing
+
+        z_values[lpos.x][lpos.y] = new_z;                   // Save the updated Z value
+        TERN_(EXTENSIBLE_UI, ExtUI::onMeshUpdate(location, new_z));
+
+        serial_delay(20);                                   // No switch noise
+        ui.refresh();
+
+      } while (lpos.x >= 0 && --g29_repetition_cnt > 0);
+
+      if (do_ubl_mesh_map) display_map(g29_map_type);
+      restore_ubl_active_state_and_leave();
+
+      do_blocking_move_to_xy_z(pos, Z_CLEARANCE_BETWEEN_PROBES);
+
+      LCD_MESSAGEPGM(MSG_UBL_DONE_EDITING_MESH);
+      SERIAL_ECHOLNPGM("Done Editing Mesh");
+
+      if (lcd_map_control)
+        ui.goto_screen(ubl_map_screen);
+      else
+        ui.return_to_status();
+    }
+
+  #endif // HAS_LCD_MENU
+
+  bool unified_bed_leveling::g29_parameter_parsing() {
+    bool err_flag = false;
+
+    TERN_(HAS_LCD_MENU, set_message_with_feedback(GET_TEXT(MSG_UBL_DOING_G29)));
+
+    g29_constant = 0;
+    g29_repetition_cnt = 0;
+
+    if (parser.seen('R')) {
+      g29_repetition_cnt = parser.has_value() ? parser.value_int() : GRID_MAX_POINTS;
+      NOMORE(g29_repetition_cnt, GRID_MAX_POINTS);
+      if (g29_repetition_cnt < 1) {
+        SERIAL_ECHOLNPGM("?(R)epetition count invalid (1+).\n");
+        return UBL_ERR;
+      }
+    }
+
+    g29_verbose_level = parser.seen('V') ? parser.value_int() : 0;
+    if (!WITHIN(g29_verbose_level, 0, 4)) {
+      SERIAL_ECHOLNPGM("?(V)erbose level implausible (0-4).\n");
+      err_flag = true;
+    }
+
+    if (parser.seen('P')) {
+      const int pv = parser.value_int();
+      #if !HAS_BED_PROBE
+        if (pv == 1) {
+          SERIAL_ECHOLNPGM("G29 P1 requires a probe.\n");
+          err_flag = true;
+        }
+        else
+      #endif
+        {
+          g29_phase_value = pv;
+          if (!WITHIN(g29_phase_value, 0, 6)) {
+            SERIAL_ECHOLNPGM("?(P)hase value invalid (0-6).\n");
+            err_flag = true;
+          }
+        }
+    }
+
+    if (parser.seen('J')) {
+      #if HAS_BED_PROBE
+        g29_grid_size = parser.has_value() ? parser.value_int() : 0;
+        if (g29_grid_size && !WITHIN(g29_grid_size, 2, 9)) {
+          SERIAL_ECHOLNPGM("?Invalid grid size (J) specified (2-9).\n");
+          err_flag = true;
+        }
+      #else
+        SERIAL_ECHOLNPGM("G29 J action requires a probe.\n");
+        err_flag = true;
+      #endif
+    }
+
+    xy_seen.x = parser.seenval('X');
+    float sx = xy_seen.x ? parser.value_float() : current_position.x;
+    xy_seen.y = parser.seenval('Y');
+    float sy = xy_seen.y ? parser.value_float() : current_position.y;
+
+    if (xy_seen.x != xy_seen.y) {
+      SERIAL_ECHOLNPGM("Both X & Y locations must be specified.\n");
+      err_flag = true;
+    }
+
+    // If X or Y are not valid, use center of the bed values
+    if (!WITHIN(sx, X_MIN_BED, X_MAX_BED)) sx = X_CENTER;
+    if (!WITHIN(sy, Y_MIN_BED, Y_MAX_BED)) sy = Y_CENTER;
+
+    if (err_flag) return UBL_ERR;
+
+    g29_pos.set(sx, sy);
+
+    /**
+     * Activate or deactivate UBL
+     * Note: UBL's G29 restores the state set here when done.
+     *       Leveling is being enabled here with old data, possibly
+     *       none. Error handling should disable for safety...
+     */
+    if (parser.seen('A')) {
+      if (parser.seen('D')) {
+        SERIAL_ECHOLNPGM("?Can't activate and deactivate at the same time.\n");
+        return UBL_ERR;
+      }
+      set_bed_leveling_enabled(true);
+      report_state();
+    }
+    else if (parser.seen('D')) {
+      set_bed_leveling_enabled(false);
+      report_state();
+    }
+
+    // Set global 'C' flag and its value
+    if ((g29_c_flag = parser.seen('C')))
+      g29_constant = parser.value_float();
+
+    #if ENABLED(ENABLE_LEVELING_FADE_HEIGHT)
+      if (parser.seenval('F')) {
+        const float fh = parser.value_float();
+        if (!WITHIN(fh, 0, 100)) {
+          SERIAL_ECHOLNPGM("?(F)ade height for Bed Level Correction not plausible.\n");
+          return UBL_ERR;
+        }
+        set_z_fade_height(fh);
+      }
+    #endif
+
+    g29_map_type = parser.intval('T');
+    if (!WITHIN(g29_map_type, 0, 2)) {
+      SERIAL_ECHOLNPGM("Invalid map type.\n");
+      return UBL_ERR;
+    }
+    return UBL_OK;
+  }
+
+  static uint8_t ubl_state_at_invocation = 0;
+
+  #if ENABLED(UBL_DEVEL_DEBUGGING)
+    static uint8_t ubl_state_recursion_chk = 0;
+  #endif
+
+  void unified_bed_leveling::save_ubl_active_state_and_disable() {
+    #if ENABLED(UBL_DEVEL_DEBUGGING)
+      ubl_state_recursion_chk++;
+      if (ubl_state_recursion_chk != 1) {
+        SERIAL_ECHOLNPGM("save_ubl_active_state_and_disabled() called multiple times in a row.");
+        TERN_(HAS_LCD_MENU, set_message_with_feedback(GET_TEXT(MSG_UBL_SAVE_ERROR)));
+        return;
+      }
+    #endif
+    ubl_state_at_invocation = planner.leveling_active;
+    set_bed_leveling_enabled(false);
+  }
+
+  void unified_bed_leveling::restore_ubl_active_state_and_leave() {
+    TERN_(HAS_LCD_MENU, ui.release());
+    #if ENABLED(UBL_DEVEL_DEBUGGING)
+      if (--ubl_state_recursion_chk) {
+        SERIAL_ECHOLNPGM("restore_ubl_active_state_and_leave() called too many times.");
+        TERN_(HAS_LCD_MENU, set_message_with_feedback(GET_TEXT(MSG_UBL_RESTORE_ERROR)));
+        return;
+      }
+    #endif
+    set_bed_leveling_enabled(ubl_state_at_invocation);
+  }
+
+  mesh_index_pair unified_bed_leveling::find_furthest_invalid_mesh_point() {
+
+    bool found_a_NAN = false, found_a_real = false;
+
+    mesh_index_pair farthest { -1, -1, -99999.99 };
+
+    GRID_LOOP(i, j) {
+      if (!isnan(z_values[i][j])) continue;  // Skip valid mesh points
+
+      // Skip unreachable points
+      if (!probe.can_reach(mesh_index_to_xpos(i), mesh_index_to_ypos(j)))
+        continue;
+
+      found_a_NAN = true;
+
+      xy_int8_t nearby { -1, -1 };
+      float d1, d2 = 99999.9f;
+      GRID_LOOP(k, l) {
+        if (isnan(z_values[k][l])) continue;
+
+        found_a_real = true;
+
+        // Add in a random weighting factor that scrambles the probing of the
+        // last half of the mesh (when every unprobed mesh point is one index
+        // from a probed location).
+
+        d1 = HYPOT(i - k, j - l) + (1.0f / ((millis() % 47) + 13));
+
+        if (d1 < d2) {    // Invalid mesh point (i,j) is closer to the defined point (k,l)
+          d2 = d1;
+          nearby.set(i, j);
+        }
+      }
+
+      //
+      // At this point d2 should have the near defined mesh point to invalid mesh point (i,j)
+      //
+
+      if (found_a_real && nearby.x >= 0 && d2 > farthest.distance) {
+        farthest.pos = nearby; // Found an invalid location farther from the defined mesh point
+        farthest.distance = d2;
+      }
+    } // GRID_LOOP
+
+    if (!found_a_real && found_a_NAN) {        // if the mesh is totally unpopulated, start the probing
+      farthest.pos.set((GRID_MAX_POINTS_X) / 2, (GRID_MAX_POINTS_Y) / 2);
+      farthest.distance = 1;
+    }
+    return farthest;
+  }
+
+  mesh_index_pair unified_bed_leveling::find_closest_mesh_point_of_type(const MeshPointType type, const xy_pos_t &pos, const bool probe_relative/*=false*/, MeshFlags *done_flags/*=nullptr*/) {
+    mesh_index_pair closest;
+    closest.invalidate();
+    closest.distance = -99999.9f;
+
+    // Get the reference position, either nozzle or probe
+    const xy_pos_t ref = probe_relative ? pos + probe.offset_xy : pos;
+
+    float best_so_far = 99999.99f;
+
+    GRID_LOOP(i, j) {
+      if ( (type == (isnan(z_values[i][j]) ? INVALID : REAL))
+        || (type == SET_IN_BITMAP && !done_flags->marked(i, j))
+      ) {
+        // Found a Mesh Point of the specified type!
+        const xy_pos_t mpos = { mesh_index_to_xpos(i), mesh_index_to_ypos(j) };
+
+        // If using the probe as the reference there are some unreachable locations.
+        // Also for round beds, there are grid points outside the bed the nozzle can't reach.
+        // Prune them from the list and ignore them till the next Phase (manual nozzle probing).
+
+        if (!(probe_relative ? probe.can_reach(mpos) : position_is_reachable(mpos)))
+          continue;
+
+        // Reachable. Check if it's the best_so_far location to the nozzle.
+
+        const xy_pos_t diff = current_position - mpos;
+        const float distance = (ref - mpos).magnitude() + diff.magnitude() * 0.1f;
+
+        // factor in the distance from the current location for the normal case
+        // so the nozzle isn't running all over the bed.
+        if (distance < best_so_far) {
+          best_so_far = distance;   // Found a closer location with the desired value type.
+          closest.pos.set(i, j);
+          closest.distance = best_so_far;
+        }
+      }
+    } // GRID_LOOP
+
+    return closest;
+  }
+
+  /**
+   * 'Smart Fill': Scan from the outward edges of the mesh towards the center.
+   * If an invalid location is found, use the next two points (if valid) to
+   * calculate a 'reasonable' value for the unprobed mesh point.
+   */
+
+  bool unified_bed_leveling::smart_fill_one(const uint8_t x, const uint8_t y, const int8_t xdir, const int8_t ydir) {
+    const float v = z_values[x][y];
+    if (isnan(v)) {                           // A NAN...
+      const int8_t dx = x + xdir, dy = y + ydir;
+      const float v1 = z_values[dx][dy];
+      if (!isnan(v1)) {                       // ...next to a pair of real values?
+        const float v2 = z_values[dx + xdir][dy + ydir];
+        if (!isnan(v2)) {
+          z_values[x][y] = v1 < v2 ? v1 : v1 + v1 - v2;
+          TERN_(EXTENSIBLE_UI, ExtUI::onMeshUpdate(x, y, z_values[x][y]));
+          return true;
+        }
+      }
+    }
+    return false;
+  }
+
+  typedef struct { uint8_t sx, ex, sy, ey; bool yfirst; } smart_fill_info;
+
+  void unified_bed_leveling::smart_fill_mesh() {
+    static const smart_fill_info
+      info0 PROGMEM = { 0, GRID_MAX_POINTS_X,      0, GRID_MAX_POINTS_Y - 2,  false },  // Bottom of the mesh looking up
+      info1 PROGMEM = { 0, GRID_MAX_POINTS_X,      GRID_MAX_POINTS_Y - 1, 0,  false },  // Top of the mesh looking down
+      info2 PROGMEM = { 0, GRID_MAX_POINTS_X - 2,  0, GRID_MAX_POINTS_Y,      true  },  // Left side of the mesh looking right
+      info3 PROGMEM = { GRID_MAX_POINTS_X - 1, 0,  0, GRID_MAX_POINTS_Y,      true  };  // Right side of the mesh looking left
+    static const smart_fill_info * const info[] PROGMEM = { &info0, &info1, &info2, &info3 };
+
+    LOOP_L_N(i, COUNT(info)) {
+      const smart_fill_info *f = (smart_fill_info*)pgm_read_ptr(&info[i]);
+      const int8_t sx = pgm_read_byte(&f->sx), sy = pgm_read_byte(&f->sy),
+                   ex = pgm_read_byte(&f->ex), ey = pgm_read_byte(&f->ey);
+      if (pgm_read_byte(&f->yfirst)) {
+        const int8_t dir = ex > sx ? 1 : -1;
+        for (uint8_t y = sy; y != ey; ++y)
+          for (uint8_t x = sx; x != ex; x += dir)
+            if (smart_fill_one(x, y, dir, 0)) break;
+      }
+      else {
+        const int8_t dir = ey > sy ? 1 : -1;
+         for (uint8_t x = sx; x != ex; ++x)
+          for (uint8_t y = sy; y != ey; y += dir)
+            if (smart_fill_one(x, y, 0, dir)) break;
+      }
+    }
+  }
+
+  #if HAS_BED_PROBE
+
+    //#define VALIDATE_MESH_TILT
+
+    #include "../../../libs/vector_3.h"
+
+    void unified_bed_leveling::tilt_mesh_based_on_probed_grid(const bool do_3_pt_leveling) {
+      const float x_min = probe.min_x(), x_max = probe.max_x(),
+                  y_min = probe.min_y(), y_max = probe.max_y(),
+                  dx = (x_max - x_min) / (g29_grid_size - 1),
+                  dy = (y_max - y_min) / (g29_grid_size - 1);
+
+      xy_float_t points[3];
+      probe.get_three_points(points);
+
+      float measured_z;
+      bool abort_flag = false;
+
+      #ifdef VALIDATE_MESH_TILT
+        float z1, z2, z3;  // Needed for algorithm validation below
+      #endif
+
+      struct linear_fit_data lsf_results;
+      incremental_LSF_reset(&lsf_results);
+
+      if (do_3_pt_leveling) {
+        SERIAL_ECHOLNPGM("Tilting mesh (1/3)");
+        TERN_(HAS_DISPLAY, ui.status_printf_P(0, PSTR(S_FMT " 1/3"), GET_TEXT(MSG_LCD_TILTING_MESH)));
+
+        measured_z = probe.probe_at_point(points[0], PROBE_PT_RAISE, g29_verbose_level);
+        if (isnan(measured_z))
+          abort_flag = true;
+        else {
+          measured_z -= get_z_correction(points[0]);
+          #ifdef VALIDATE_MESH_TILT
+            z1 = measured_z;
+          #endif
+          if (g29_verbose_level > 3) {
+            serial_spaces(16);
+            SERIAL_ECHOLNPAIR("Corrected_Z=", measured_z);
+          }
+          incremental_LSF(&lsf_results, points[0], measured_z);
+        }
+
+        if (!abort_flag) {
+          SERIAL_ECHOLNPGM("Tilting mesh (2/3)");
+          TERN_(HAS_DISPLAY, ui.status_printf_P(0, PSTR(S_FMT " 2/3"), GET_TEXT(MSG_LCD_TILTING_MESH)));
+
+          measured_z = probe.probe_at_point(points[1], PROBE_PT_RAISE, g29_verbose_level);
+          #ifdef VALIDATE_MESH_TILT
+            z2 = measured_z;
+          #endif
+          if (isnan(measured_z))
+            abort_flag = true;
+          else {
+            measured_z -= get_z_correction(points[1]);
+            if (g29_verbose_level > 3) {
+              serial_spaces(16);
+              SERIAL_ECHOLNPAIR("Corrected_Z=", measured_z);
+            }
+            incremental_LSF(&lsf_results, points[1], measured_z);
+          }
+        }
+
+        if (!abort_flag) {
+          SERIAL_ECHOLNPGM("Tilting mesh (3/3)");
+          TERN_(HAS_DISPLAY, ui.status_printf_P(0, PSTR(S_FMT " 3/3"), GET_TEXT(MSG_LCD_TILTING_MESH)));
+
+          measured_z = probe.probe_at_point(points[2], PROBE_PT_STOW, g29_verbose_level);
+          #ifdef VALIDATE_MESH_TILT
+            z3 = measured_z;
+          #endif
+          if (isnan(measured_z))
+            abort_flag = true;
+          else {
+            measured_z -= get_z_correction(points[2]);
+            if (g29_verbose_level > 3) {
+              serial_spaces(16);
+              SERIAL_ECHOLNPAIR("Corrected_Z=", measured_z);
+            }
+            incremental_LSF(&lsf_results, points[2], measured_z);
+          }
+        }
+
+        probe.stow();
+        probe.move_z_after_probing();
+
+        if (abort_flag) {
+          SERIAL_ECHOLNPGM("?Error probing point. Aborting operation.");
+          return;
+        }
+      }
+      else { // !do_3_pt_leveling
+
+        bool zig_zag = false;
+
+        const uint16_t total_points = sq(g29_grid_size);
+        uint16_t point_num = 1;
+
+        xy_pos_t rpos;
+        LOOP_L_N(ix, g29_grid_size) {
+          rpos.x = x_min + ix * dx;
+          LOOP_L_N(iy, g29_grid_size) {
+            rpos.y = y_min + dy * (zig_zag ? g29_grid_size - 1 - iy : iy);
+
+            if (!abort_flag) {
+              SERIAL_ECHOLNPAIR("Tilting mesh point ", point_num, "/", total_points, "\n");
+              TERN_(HAS_DISPLAY, ui.status_printf_P(0, PSTR(S_FMT " %i/%i"), GET_TEXT(MSG_LCD_TILTING_MESH), point_num, total_points));
+
+              measured_z = probe.probe_at_point(rpos, parser.seen('E') ? PROBE_PT_STOW : PROBE_PT_RAISE, g29_verbose_level); // TODO: Needs error handling
+
+              abort_flag = isnan(measured_z);
+
+              #if ENABLED(DEBUG_LEVELING_FEATURE)
+                if (DEBUGGING(LEVELING)) {
+                  const xy_pos_t lpos = rpos.asLogical();
+                  DEBUG_CHAR('(');
+                  DEBUG_ECHO_F(rpos.x, 7);
+                  DEBUG_CHAR(',');
+                  DEBUG_ECHO_F(rpos.y, 7);
+                  DEBUG_ECHOPAIR_F(")   logical: (", lpos.x, 7);
+                  DEBUG_CHAR(',');
+                  DEBUG_ECHO_F(lpos.y, 7);
+                  DEBUG_ECHOPAIR_F(")   measured: ", measured_z, 7);
+                  DEBUG_ECHOPAIR_F("   correction: ", get_z_correction(rpos), 7);
+                }
+              #endif
+
+              measured_z -= get_z_correction(rpos) /* + probe.offset.z */ ;
+
+              if (DEBUGGING(LEVELING)) DEBUG_ECHOLNPAIR_F("   final >>>---> ", measured_z, 7);
+
+              if (g29_verbose_level > 3) {
+                serial_spaces(16);
+                SERIAL_ECHOLNPAIR("Corrected_Z=", measured_z);
+              }
+              incremental_LSF(&lsf_results, rpos, measured_z);
+            }
+
+            point_num++;
+          }
+
+          zig_zag ^= true;
+        }
+      }
+      probe.stow();
+      probe.move_z_after_probing();
+
+      if (abort_flag || finish_incremental_LSF(&lsf_results)) {
+        SERIAL_ECHOPGM("Could not complete LSF!");
+        return;
+      }
+
+      vector_3 normal = vector_3(lsf_results.A, lsf_results.B, 1).get_normal();
+
+      if (g29_verbose_level > 2) {
+        SERIAL_ECHOPAIR_F("bed plane normal = [", normal.x, 7);
+        SERIAL_CHAR(',');
+        SERIAL_ECHO_F(normal.y, 7);
+        SERIAL_CHAR(',');
+        SERIAL_ECHO_F(normal.z, 7);
+        SERIAL_ECHOLNPGM("]");
+      }
+
+      matrix_3x3 rotation = matrix_3x3::create_look_at(vector_3(lsf_results.A, lsf_results.B, 1));
+
+      GRID_LOOP(i, j) {
+        float mx = mesh_index_to_xpos(i),
+              my = mesh_index_to_ypos(j),
+              mz = z_values[i][j];
+
+        if (DEBUGGING(LEVELING)) {
+          DEBUG_ECHOPAIR_F("before rotation = [", mx, 7);
+          DEBUG_CHAR(',');
+          DEBUG_ECHO_F(my, 7);
+          DEBUG_CHAR(',');
+          DEBUG_ECHO_F(mz, 7);
+          DEBUG_ECHOPGM("]   ---> ");
+          DEBUG_DELAY(20);
+        }
+
+        apply_rotation_xyz(rotation, mx, my, mz);
+
+        if (DEBUGGING(LEVELING)) {
+          DEBUG_ECHOPAIR_F("after rotation = [", mx, 7);
+          DEBUG_CHAR(',');
+          DEBUG_ECHO_F(my, 7);
+          DEBUG_CHAR(',');
+          DEBUG_ECHO_F(mz, 7);
+          DEBUG_ECHOLNPGM("]");
+          DEBUG_DELAY(20);
+        }
+
+        z_values[i][j] = mz - lsf_results.D;
+        TERN_(EXTENSIBLE_UI, ExtUI::onMeshUpdate(i, j, z_values[i][j]));
+      }
+
+      if (DEBUGGING(LEVELING)) {
+        rotation.debug(PSTR("rotation matrix:\n"));
+        DEBUG_ECHOPAIR_F("LSF Results A=", lsf_results.A, 7);
+        DEBUG_ECHOPAIR_F("  B=", lsf_results.B, 7);
+        DEBUG_ECHOLNPAIR_F("  D=", lsf_results.D, 7);
+        DEBUG_DELAY(55);
+
+        DEBUG_ECHOPAIR_F("bed plane normal = [", normal.x, 7);
+        DEBUG_CHAR(',');
+        DEBUG_ECHO_F(normal.y, 7);
+        DEBUG_CHAR(',');
+        DEBUG_ECHO_F(normal.z, 7);
+        DEBUG_ECHOLNPGM("]");
+        DEBUG_EOL();
+
+        /**
+         * Use the code below to check the validity of the mesh tilting algorithm.
+         * 3-Point Mesh Tilt uses the same algorithm as grid-based tilting, but only
+         * three points are used in the calculation. This guarantees that each probed point
+         * has an exact match when get_z_correction() for that location is calculated.
+         * The Z error between the probed point locations and the get_z_correction()
+         * numbers for those locations should be 0.
+         */
+        #ifdef VALIDATE_MESH_TILT
+          auto d_from = []{ DEBUG_ECHOPGM("D from "); };
+          auto normed = [&](const xy_pos_t &pos, const float &zadd) {
+            return normal.x * pos.x + normal.y * pos.y + zadd;
+          };
+          auto debug_pt = [](PGM_P const pre, const xy_pos_t &pos, const float &zadd) {
+            d_from(); serialprintPGM(pre);
+            DEBUG_ECHO_F(normed(pos, zadd), 6);
+            DEBUG_ECHOLNPAIR_F("   Z error = ", zadd - get_z_correction(pos), 6);
+          };
+          debug_pt(PSTR("1st point: "), probe_pt[0], normal.z * z1);
+          debug_pt(PSTR("2nd point: "), probe_pt[1], normal.z * z2);
+          debug_pt(PSTR("3rd point: "), probe_pt[2], normal.z * z3);
+          d_from(); DEBUG_ECHOPGM("safe home with Z=");
+          DEBUG_ECHOLNPAIR_F("0 : ", normed(safe_homing_xy, 0), 6);
+          d_from(); DEBUG_ECHOPGM("safe home with Z=");
+          DEBUG_ECHOLNPAIR_F("mesh value ", normed(safe_homing_xy, get_z_correction(safe_homing_xy)), 6);
+          DEBUG_ECHOPAIR("   Z error = (", Z_SAFE_HOMING_X_POINT, ",", Z_SAFE_HOMING_Y_POINT);
+          DEBUG_ECHOLNPAIR_F(") = ", get_z_correction(safe_homing_xy), 6);
+        #endif
+      } // DEBUGGING(LEVELING)
+
+    }
+
+  #endif // HAS_BED_PROBE
+
+  #if ENABLED(UBL_G29_P31)
+    void unified_bed_leveling::smart_fill_wlsf(const float &weight_factor) {
+
+      // For each undefined mesh point, compute a distance-weighted least squares fit
+      // from all the originally populated mesh points, weighted toward the point
+      // being extrapolated so that nearby points will have greater influence on
+      // the point being extrapolated.  Then extrapolate the mesh point from WLSF.
+
+      static_assert((GRID_MAX_POINTS_Y) <= 16, "GRID_MAX_POINTS_Y too big");
+      uint16_t bitmap[GRID_MAX_POINTS_X] = { 0 };
+      struct linear_fit_data lsf_results;
+
+      SERIAL_ECHOPGM("Extrapolating mesh...");
+
+      const float weight_scaled = weight_factor * _MAX(MESH_X_DIST, MESH_Y_DIST);
+
+      GRID_LOOP(jx, jy) if (!isnan(z_values[jx][jy])) SBI(bitmap[jx], jy);
+
+      xy_pos_t ppos;
+      LOOP_L_N(ix, GRID_MAX_POINTS_X) {
+        ppos.x = mesh_index_to_xpos(ix);
+        LOOP_L_N(iy, GRID_MAX_POINTS_Y) {
+          ppos.y = mesh_index_to_ypos(iy);
+          if (isnan(z_values[ix][iy])) {
+            // undefined mesh point at (ppos.x,ppos.y), compute weighted LSF from original valid mesh points.
+            incremental_LSF_reset(&lsf_results);
+            xy_pos_t rpos;
+            LOOP_L_N(jx, GRID_MAX_POINTS_X) {
+              rpos.x = mesh_index_to_xpos(jx);
+              LOOP_L_N(jy, GRID_MAX_POINTS_Y) {
+                if (TEST(bitmap[jx], jy)) {
+                  rpos.y = mesh_index_to_ypos(jy);
+                  const float rz = z_values[jx][jy],
+                               w = 1.0f + weight_scaled / (rpos - ppos).magnitude();
+                  incremental_WLSF(&lsf_results, rpos, rz, w);
+                }
+              }
+            }
+            if (finish_incremental_LSF(&lsf_results)) {
+              SERIAL_ECHOLNPGM("Insufficient data");
+              return;
+            }
+            const float ez = -lsf_results.D - lsf_results.A * ppos.x - lsf_results.B * ppos.y;
+            z_values[ix][iy] = ez;
+            TERN_(EXTENSIBLE_UI, ExtUI::onMeshUpdate(ix, iy, z_values[ix][iy]));
+            idle(); // housekeeping
+          }
+        }
+      }
+
+      SERIAL_ECHOLNPGM("done");
+    }
+  #endif // UBL_G29_P31
+
+  #if ENABLED(UBL_DEVEL_DEBUGGING)
+    /**
+     * Much of the 'What?' command can be eliminated. But until we are fully debugged, it is
+     * good to have the extra information. Soon... we prune this to just a few items
+     */
+    void unified_bed_leveling::g29_what_command() {
+      report_state();
+
+      if (storage_slot == -1)
+        SERIAL_ECHOPGM("No Mesh Loaded.");
+      else
+        SERIAL_ECHOPAIR("Mesh ", storage_slot, " Loaded.");
+      SERIAL_EOL();
+      serial_delay(50);
+
+      #if ENABLED(ENABLE_LEVELING_FADE_HEIGHT)
+        SERIAL_ECHOLNPAIR_F("Fade Height M420 Z", planner.z_fade_height, 4);
+      #endif
+
+      adjust_mesh_to_mean(g29_c_flag, g29_constant);
+
+      #if HAS_BED_PROBE
+        SERIAL_ECHOLNPAIR_F("Probe Offset M851 Z", probe.offset.z, 7);
+      #endif
+
+      SERIAL_ECHOLNPAIR("MESH_MIN_X  " STRINGIFY(MESH_MIN_X) "=", MESH_MIN_X); serial_delay(50);
+      SERIAL_ECHOLNPAIR("MESH_MIN_Y  " STRINGIFY(MESH_MIN_Y) "=", MESH_MIN_Y); serial_delay(50);
+      SERIAL_ECHOLNPAIR("MESH_MAX_X  " STRINGIFY(MESH_MAX_X) "=", MESH_MAX_X); serial_delay(50);
+      SERIAL_ECHOLNPAIR("MESH_MAX_Y  " STRINGIFY(MESH_MAX_Y) "=", MESH_MAX_Y); serial_delay(50);
+      SERIAL_ECHOLNPAIR("GRID_MAX_POINTS_X  ", GRID_MAX_POINTS_X);             serial_delay(50);
+      SERIAL_ECHOLNPAIR("GRID_MAX_POINTS_Y  ", GRID_MAX_POINTS_Y);             serial_delay(50);
+      SERIAL_ECHOLNPAIR("MESH_X_DIST  ", MESH_X_DIST);
+      SERIAL_ECHOLNPAIR("MESH_Y_DIST  ", MESH_Y_DIST);                         serial_delay(50);
+
+      SERIAL_ECHOPGM("X-Axis Mesh Points at: ");
+      LOOP_L_N(i, GRID_MAX_POINTS_X) {
+        SERIAL_ECHO_F(LOGICAL_X_POSITION(mesh_index_to_xpos(i)), 3);
+        SERIAL_ECHOPGM("  ");
+        serial_delay(25);
+      }
+      SERIAL_EOL();
+
+      SERIAL_ECHOPGM("Y-Axis Mesh Points at: ");
+      LOOP_L_N(i, GRID_MAX_POINTS_Y) {
+        SERIAL_ECHO_F(LOGICAL_Y_POSITION(mesh_index_to_ypos(i)), 3);
+        SERIAL_ECHOPGM("  ");
+        serial_delay(25);
+      }
+      SERIAL_EOL();
+
+      #if HAS_KILL
+        SERIAL_ECHOLNPAIR("Kill pin on :", int(KILL_PIN), "  state:", int(kill_state()));
+      #endif
+
+      SERIAL_EOL();
+      serial_delay(50);
+
+      #if ENABLED(UBL_DEVEL_DEBUGGING)
+        SERIAL_ECHOLNPAIR("ubl_state_at_invocation :", ubl_state_at_invocation, "\nubl_state_recursion_chk :", ubl_state_recursion_chk);
+        serial_delay(50);
+
+        SERIAL_ECHOLNPAIR("Meshes go from ", hex_address((void*)settings.meshes_start_index()), " to ", hex_address((void*)settings.meshes_end_index()));
+        serial_delay(50);
+
+        SERIAL_ECHOLNPAIR("sizeof(ubl) :  ", (int)sizeof(ubl));         SERIAL_EOL();
+        SERIAL_ECHOLNPAIR("z_value[][] size: ", (int)sizeof(z_values)); SERIAL_EOL();
+        serial_delay(25);
+
+        SERIAL_ECHOLNPAIR("EEPROM free for UBL: ", hex_address((void*)(settings.meshes_end_index() - settings.meshes_start_index())));
+        serial_delay(50);
+
+        SERIAL_ECHOLNPAIR("EEPROM can hold ", settings.calc_num_meshes(), " meshes.\n");
+        serial_delay(25);
+      #endif // UBL_DEVEL_DEBUGGING
+
+      if (!sanity_check()) {
+        echo_name();
+        SERIAL_ECHOLNPGM(" sanity checks passed.");
+      }
+    }
+
+    /**
+     * When we are fully debugged, the EEPROM dump command will get deleted also. But
+     * right now, it is good to have the extra information. Soon... we prune this.
+     */
+    void unified_bed_leveling::g29_eeprom_dump() {
+      uint8_t cccc;
+
+      SERIAL_ECHO_MSG("EEPROM Dump:");
+      persistentStore.access_start();
+      for (uint16_t i = 0; i < persistentStore.capacity(); i += 16) {
+        if (!(i & 0x3)) idle();
+        print_hex_word(i);
+        SERIAL_ECHOPGM(": ");
+        for (uint16_t j = 0; j < 16; j++) {
+          persistentStore.read_data(i + j, &cccc, sizeof(uint8_t));
+          print_hex_byte(cccc);
+          SERIAL_CHAR(' ');
+        }
+        SERIAL_EOL();
+      }
+      SERIAL_EOL();
+      persistentStore.access_finish();
+    }
+
+    /**
+     * When we are fully debugged, this may go away. But there are some valid
+     * use cases for the users. So we can wait and see what to do with it.
+     */
+    void unified_bed_leveling::g29_compare_current_mesh_to_stored_mesh() {
+      const int16_t a = settings.calc_num_meshes();
+
+      if (!a) {
+        SERIAL_ECHOLNPGM("?EEPROM storage not available.");
+        return;
+      }
+
+      if (!parser.has_value() || !WITHIN(g29_storage_slot, 0, a - 1)) {
+        SERIAL_ECHOLNPAIR("?Invalid storage slot.\n?Use 0 to ", a - 1);
+        return;
+      }
+
+      g29_storage_slot = parser.value_int();
+
+      float tmp_z_values[GRID_MAX_POINTS_X][GRID_MAX_POINTS_Y];
+      settings.load_mesh(g29_storage_slot, &tmp_z_values);
+
+      SERIAL_ECHOLNPAIR("Subtracting mesh in slot ", g29_storage_slot, " from current mesh.");
+
+      GRID_LOOP(x, y) {
+        z_values[x][y] -= tmp_z_values[x][y];
+        TERN_(EXTENSIBLE_UI, ExtUI::onMeshUpdate(x, y, z_values[x][y]));
+      }
+    }
+
+  #endif // UBL_DEVEL_DEBUGGING
+
+#endif // AUTO_BED_LEVELING_UBL
diff --git a/Marlin/src/feature/bedlevel/ubl/ubl_motion.cpp b/Marlin/src/feature/bedlevel/ubl/ubl_motion.cpp
new file mode 100644
index 0000000..8b7cd15
--- /dev/null
+++ b/Marlin/src/feature/bedlevel/ubl/ubl_motion.cpp
@@ -0,0 +1,474 @@
+/**
+ * Marlin 3D Printer Firmware
+ * Copyright (c) 2020 MarlinFirmware [https://github.com/MarlinFirmware/Marlin]
+ *
+ * Based on Sprinter and grbl.
+ * Copyright (c) 2011 Camiel Gubbels / Erik van der Zalm
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <https://www.gnu.org/licenses/>.
+ *
+ */
+#include "../../../inc/MarlinConfig.h"
+
+#if ENABLED(AUTO_BED_LEVELING_UBL)
+
+#include "../bedlevel.h"
+#include "../../../module/planner.h"
+#include "../../../module/stepper.h"
+#include "../../../module/motion.h"
+
+#if ENABLED(DELTA)
+  #include "../../../module/delta.h"
+#endif
+
+#include "../../../MarlinCore.h"
+#include <math.h>
+
+#if !UBL_SEGMENTED
+
+  void unified_bed_leveling::line_to_destination_cartesian(const feedRate_t &scaled_fr_mm_s, const uint8_t extruder) {
+    /**
+     * Much of the nozzle movement will be within the same cell. So we will do as little computation
+     * as possible to determine if this is the case. If this move is within the same cell, we will
+     * just do the required Z-Height correction, call the Planner's buffer_line() routine, and leave
+     */
+    #if HAS_POSITION_MODIFIERS
+      xyze_pos_t start = current_position, end = destination;
+      planner.apply_modifiers(start);
+      planner.apply_modifiers(end);
+    #else
+      const xyze_pos_t &start = current_position, &end = destination;
+    #endif
+
+    const xy_int8_t istart = cell_indexes(start), iend = cell_indexes(end);
+
+    // A move within the same cell needs no splitting
+    if (istart == iend) {
+
+      FINAL_MOVE:
+
+      // When UBL_Z_RAISE_WHEN_OFF_MESH is disabled Z correction is extrapolated from the edge of the mesh
+      #ifdef UBL_Z_RAISE_WHEN_OFF_MESH
+        // For a move off the UBL mesh, use a constant Z raise
+        if (!cell_index_x_valid(end.x) || !cell_index_y_valid(end.y)) {
+
+          // Note: There is no Z Correction in this case. We are off the mesh and don't know what
+          // a reasonable correction would be, UBL_Z_RAISE_WHEN_OFF_MESH will be used instead of
+          // a calculated (Bi-Linear interpolation) correction.
+
+          end.z += UBL_Z_RAISE_WHEN_OFF_MESH;
+          planner.buffer_segment(end, scaled_fr_mm_s, extruder);
+          current_position = destination;
+          return;
+        }
+      #endif
+
+      // The distance is always MESH_X_DIST so multiply by the constant reciprocal.
+      const float xratio = (end.x - mesh_index_to_xpos(iend.x)) * RECIPROCAL(MESH_X_DIST),
+                  yratio = (end.y - mesh_index_to_ypos(iend.y)) * RECIPROCAL(MESH_Y_DIST),
+                  z1 = z_values[iend.x][iend.y    ] + xratio * (z_values[iend.x + 1][iend.y    ] - z_values[iend.x][iend.y    ]),
+                  z2 = z_values[iend.x][iend.y + 1] + xratio * (z_values[iend.x + 1][iend.y + 1] - z_values[iend.x][iend.y + 1]);
+
+      // X cell-fraction done. Interpolate the two Z offsets with the Y fraction for the final Z offset.
+      const float z0 = (z1 + (z2 - z1) * yratio) * planner.fade_scaling_factor_for_z(end.z);
+
+      // Undefined parts of the Mesh in z_values[][] are NAN.
+      // Replace NAN corrections with 0.0 to prevent NAN propagation.
+      if (!isnan(z0)) end.z += z0;
+      planner.buffer_segment(end, scaled_fr_mm_s, extruder);
+      current_position = destination;
+      return;
+    }
+
+    /**
+     * Past this point the move is known to cross one or more mesh lines. Check for the most common
+     * case - crossing only one X or Y line - after details are worked out to reduce computation.
+     */
+
+    const xy_float_t dist = end - start;
+    const xy_bool_t neg { dist.x < 0, dist.y < 0 };
+    const xy_int8_t ineg { int8_t(neg.x), int8_t(neg.y) };
+    const xy_float_t sign { neg.x ? -1.0f : 1.0f, neg.y ? -1.0f : 1.0f };
+    const xy_int8_t iadd { int8_t(iend.x == istart.x ? 0 : sign.x), int8_t(iend.y == istart.y ? 0 : sign.y) };
+
+    /**
+     * Compute the extruder scaling factor for each partial move, checking for
+     * zero-length moves that would result in an infinite scaling factor.
+     * A float divide is required for this, but then it just multiplies.
+     * Also select a scaling factor based on the larger of the X and Y
+     * components. The larger of the two is used to preserve precision.
+     */
+
+    const xy_float_t ad = sign * dist;
+    const bool use_x_dist = ad.x > ad.y;
+
+    float on_axis_distance = use_x_dist ? dist.x : dist.y,
+          e_position = end.e - start.e,
+          z_position = end.z - start.z;
+
+    const float e_normalized_dist = e_position / on_axis_distance, // Allow divide by zero
+                z_normalized_dist = z_position / on_axis_distance;
+
+    xy_int8_t icell = istart;
+
+    const float ratio = dist.y / dist.x,        // Allow divide by zero
+                c = start.y - ratio * start.x;
+
+    const bool inf_normalized_flag = isinf(e_normalized_dist),
+               inf_ratio_flag = isinf(ratio);
+
+    /**
+     * Handle vertical lines that stay within one column.
+     * These need not be perfectly vertical.
+     */
+    if (iadd.x == 0) {        // Vertical line?
+      icell.y += ineg.y;      // Line going down? Just go to the bottom.
+      while (icell.y != iend.y + ineg.y) {
+        icell.y += iadd.y;
+        const float next_mesh_line_y = mesh_index_to_ypos(icell.y);
+
+        /**
+         * Skip the calculations for an infinite slope.
+         * For others the next X is the same so this can continue.
+         * Calculate X at the next Y mesh line.
+         */
+        const float rx = inf_ratio_flag ? start.x : (next_mesh_line_y - c) / ratio;
+
+        float z0 = z_correction_for_x_on_horizontal_mesh_line(rx, icell.x, icell.y)
+                   * planner.fade_scaling_factor_for_z(end.z);
+
+        // Undefined parts of the Mesh in z_values[][] are NAN.
+        // Replace NAN corrections with 0.0 to prevent NAN propagation.
+        if (isnan(z0)) z0 = 0.0;
+
+        const float ry = mesh_index_to_ypos(icell.y);
+
+        /**
+         * Without this check, it's possible to generate a zero length move, as in the case where
+         * the line is heading down, starting exactly on a mesh line boundary. Since this is rare
+         * it might be fine to remove this check and let planner.buffer_segment() filter it out.
+         */
+        if (ry != start.y) {
+          if (!inf_normalized_flag) { // fall-through faster than branch
+            on_axis_distance = use_x_dist ? rx - start.x : ry - start.y;
+            e_position = start.e + on_axis_distance * e_normalized_dist;
+            z_position = start.z + on_axis_distance * z_normalized_dist;
+          }
+          else {
+            e_position = end.e;
+            z_position = end.z;
+          }
+
+          planner.buffer_segment(rx, ry, z_position + z0, e_position, scaled_fr_mm_s, extruder);
+        } //else printf("FIRST MOVE PRUNED  ");
+      }
+
+      // At the final destination? Usually not, but when on a Y Mesh Line it's completed.
+      if (xy_pos_t(current_position) != xy_pos_t(end))
+        goto FINAL_MOVE;
+
+      current_position = destination;
+      return;
+    }
+
+    /**
+     * Handle horizontal lines that stay within one row.
+     * These need not be perfectly horizontal.
+     */
+    if (iadd.y == 0) {      // Horizontal line?
+      icell.x += ineg.x;     // Heading left? Just go to the left edge of the cell for the first move.
+      while (icell.x != iend.x + ineg.x) {
+        icell.x += iadd.x;
+        const float rx = mesh_index_to_xpos(icell.x);
+        const float ry = ratio * rx + c;    // Calculate Y at the next X mesh line
+
+        float z0 = z_correction_for_y_on_vertical_mesh_line(ry, icell.x, icell.y)
+                     * planner.fade_scaling_factor_for_z(end.z);
+
+        // Undefined parts of the Mesh in z_values[][] are NAN.
+        // Replace NAN corrections with 0.0 to prevent NAN propagation.
+        if (isnan(z0)) z0 = 0.0;
+
+        /**
+         * Without this check, it's possible to generate a zero length move, as in the case where
+         * the line is heading left, starting exactly on a mesh line boundary. Since this is rare
+         * it might be fine to remove this check and let planner.buffer_segment() filter it out.
+         */
+        if (rx != start.x) {
+          if (!inf_normalized_flag) {
+            on_axis_distance = use_x_dist ? rx - start.x : ry - start.y;
+            e_position = start.e + on_axis_distance * e_normalized_dist;  // is based on X or Y because this is a horizontal move
+            z_position = start.z + on_axis_distance * z_normalized_dist;
+          }
+          else {
+            e_position = end.e;
+            z_position = end.z;
+          }
+
+          if (!planner.buffer_segment(rx, ry, z_position + z0, e_position, scaled_fr_mm_s, extruder))
+            break;
+        } //else printf("FIRST MOVE PRUNED  ");
+      }
+
+      if (xy_pos_t(current_position) != xy_pos_t(end))
+        goto FINAL_MOVE;
+
+      current_position = destination;
+      return;
+    }
+
+    /**
+     * Generic case of a line crossing both X and Y Mesh lines.
+     */
+
+    xy_int8_t cnt = (istart - iend).ABS();
+
+    icell += ineg;
+
+    while (cnt) {
+
+      const float next_mesh_line_x = mesh_index_to_xpos(icell.x + iadd.x),
+                  next_mesh_line_y = mesh_index_to_ypos(icell.y + iadd.y),
+                  ry = ratio * next_mesh_line_x + c,    // Calculate Y at the next X mesh line
+                  rx = (next_mesh_line_y - c) / ratio;  // Calculate X at the next Y mesh line
+                                                        // (No need to worry about ratio == 0.
+                                                        //  In that case, it was already detected
+                                                        //  as a vertical line move above.)
+
+      if (neg.x == (rx > next_mesh_line_x)) { // Check if we hit the Y line first
+        // Yes!  Crossing a Y Mesh Line next
+        float z0 = z_correction_for_x_on_horizontal_mesh_line(rx, icell.x - ineg.x, icell.y + iadd.y)
+                   * planner.fade_scaling_factor_for_z(end.z);
+
+        // Undefined parts of the Mesh in z_values[][] are NAN.
+        // Replace NAN corrections with 0.0 to prevent NAN propagation.
+        if (isnan(z0)) z0 = 0.0;
+
+        if (!inf_normalized_flag) {
+          on_axis_distance = use_x_dist ? rx - start.x : next_mesh_line_y - start.y;
+          e_position = start.e + on_axis_distance * e_normalized_dist;
+          z_position = start.z + on_axis_distance * z_normalized_dist;
+        }
+        else {
+          e_position = end.e;
+          z_position = end.z;
+        }
+        if (!planner.buffer_segment(rx, next_mesh_line_y, z_position + z0, e_position, scaled_fr_mm_s, extruder))
+          break;
+        icell.y += iadd.y;
+        cnt.y--;
+      }
+      else {
+        // Yes!  Crossing a X Mesh Line next
+        float z0 = z_correction_for_y_on_vertical_mesh_line(ry, icell.x + iadd.x, icell.y - ineg.y)
+                   * planner.fade_scaling_factor_for_z(end.z);
+
+        // Undefined parts of the Mesh in z_values[][] are NAN.
+        // Replace NAN corrections with 0.0 to prevent NAN propagation.
+        if (isnan(z0)) z0 = 0.0;
+
+        if (!inf_normalized_flag) {
+          on_axis_distance = use_x_dist ? next_mesh_line_x - start.x : ry - start.y;
+          e_position = start.e + on_axis_distance * e_normalized_dist;
+          z_position = start.z + on_axis_distance * z_normalized_dist;
+        }
+        else {
+          e_position = end.e;
+          z_position = end.z;
+        }
+
+        if (!planner.buffer_segment(next_mesh_line_x, ry, z_position + z0, e_position, scaled_fr_mm_s, extruder))
+          break;
+        icell.x += iadd.x;
+        cnt.x--;
+      }
+
+      if (cnt.x < 0 || cnt.y < 0) break; // Too far! Exit the loop and go to FINAL_MOVE
+    }
+
+    if (xy_pos_t(current_position) != xy_pos_t(end))
+      goto FINAL_MOVE;
+
+    current_position = destination;
+  }
+
+#else // UBL_SEGMENTED
+
+  #if IS_SCARA
+    #define DELTA_SEGMENT_MIN_LENGTH 0.25 // SCARA minimum segment size is 0.25mm
+  #elif ENABLED(DELTA)
+    #define DELTA_SEGMENT_MIN_LENGTH 0.10 // mm (still subject to DELTA_SEGMENTS_PER_SECOND)
+  #else // CARTESIAN
+    #ifdef LEVELED_SEGMENT_LENGTH
+      #define DELTA_SEGMENT_MIN_LENGTH LEVELED_SEGMENT_LENGTH
+    #else
+      #define DELTA_SEGMENT_MIN_LENGTH 1.00 // mm (similar to G2/G3 arc segmentation)
+    #endif
+  #endif
+
+  /**
+   * Prepare a segmented linear move for DELTA/SCARA/CARTESIAN with UBL and FADE semantics.
+   * This calls planner.buffer_segment multiple times for small incremental moves.
+   * Returns true if did NOT move, false if moved (requires current_position update).
+   */
+
+  bool _O2 unified_bed_leveling::line_to_destination_segmented(const feedRate_t &scaled_fr_mm_s) {
+
+    if (!position_is_reachable(destination))  // fail if moving outside reachable boundary
+      return true;                            // did not move, so current_position still accurate
+
+    const xyze_pos_t total = destination - current_position;
+
+    const float cart_xy_mm_2 = HYPOT2(total.x, total.y),
+                cart_xy_mm = SQRT(cart_xy_mm_2);                                     // Total XY distance
+
+    #if IS_KINEMATIC
+      const float seconds = cart_xy_mm / scaled_fr_mm_s;                             // Duration of XY move at requested rate
+      uint16_t segments = LROUND(delta_segments_per_second * seconds),               // Preferred number of segments for distance @ feedrate
+               seglimit = LROUND(cart_xy_mm * RECIPROCAL(DELTA_SEGMENT_MIN_LENGTH)); // Number of segments at minimum segment length
+      NOMORE(segments, seglimit);                                                    // Limit to minimum segment length (fewer segments)
+    #else
+      uint16_t segments = LROUND(cart_xy_mm * RECIPROCAL(DELTA_SEGMENT_MIN_LENGTH)); // Cartesian fixed segment length
+    #endif
+
+    NOLESS(segments, 1U);                                                            // Must have at least one segment
+    const float inv_segments = 1.0f / segments,                                      // Reciprocal to save calculation
+                segment_xyz_mm = SQRT(cart_xy_mm_2 + sq(total.z)) * inv_segments;    // Length of each segment
+
+    #if ENABLED(SCARA_FEEDRATE_SCALING)
+      const float inv_duration = scaled_fr_mm_s / segment_xyz_mm;
+    #endif
+
+    xyze_float_t diff = total * inv_segments;
+
+    // Note that E segment distance could vary slightly as z mesh height
+    // changes for each segment, but small enough to ignore.
+
+    xyze_pos_t raw = current_position;
+
+    // Just do plain segmentation if UBL is inactive or the target is above the fade height
+    if (!planner.leveling_active || !planner.leveling_active_at_z(destination.z)) {
+      while (--segments) {
+        raw += diff;
+        planner.buffer_line(raw, scaled_fr_mm_s, active_extruder, segment_xyz_mm
+          #if ENABLED(SCARA_FEEDRATE_SCALING)
+            , inv_duration
+          #endif
+        );
+      }
+      planner.buffer_line(destination, scaled_fr_mm_s, active_extruder, segment_xyz_mm
+        #if ENABLED(SCARA_FEEDRATE_SCALING)
+          , inv_duration
+        #endif
+      );
+      return false; // Did not set current from destination
+    }
+
+    // Otherwise perform per-segment leveling
+
+    #if ENABLED(ENABLE_LEVELING_FADE_HEIGHT)
+      const float fade_scaling_factor = planner.fade_scaling_factor_for_z(destination.z);
+    #endif
+
+    // Move to first segment destination
+    raw += diff;
+
+    for (;;) {  // for each mesh cell encountered during the move
+
+      // Compute mesh cell invariants that remain constant for all segments within cell.
+      // Note for cell index, if point is outside the mesh grid (in MESH_INSET perimeter)
+      // the bilinear interpolation from the adjacent cell within the mesh will still work.
+      // Inner loop will exit each time (because out of cell bounds) but will come back
+      // in top of loop and again re-find same adjacent cell and use it, just less efficient
+      // for mesh inset area.
+
+      xy_int8_t icell = {
+        int8_t((raw.x - (MESH_MIN_X)) * RECIPROCAL(MESH_X_DIST)),
+        int8_t((raw.y - (MESH_MIN_Y)) * RECIPROCAL(MESH_Y_DIST))
+      };
+      LIMIT(icell.x, 0, (GRID_MAX_POINTS_X) - 1);
+      LIMIT(icell.y, 0, (GRID_MAX_POINTS_Y) - 1);
+
+      float z_x0y0 = z_values[icell.x  ][icell.y  ],  // z at lower left corner
+            z_x1y0 = z_values[icell.x+1][icell.y  ],  // z at upper left corner
+            z_x0y1 = z_values[icell.x  ][icell.y+1],  // z at lower right corner
+            z_x1y1 = z_values[icell.x+1][icell.y+1];  // z at upper right corner
+
+      if (isnan(z_x0y0)) z_x0y0 = 0;              // ideally activating planner.leveling_active (G29 A)
+      if (isnan(z_x1y0)) z_x1y0 = 0;              //   should refuse if any invalid mesh points
+      if (isnan(z_x0y1)) z_x0y1 = 0;              //   in order to avoid isnan tests per cell,
+      if (isnan(z_x1y1)) z_x1y1 = 0;              //   thus guessing zero for undefined points
+
+      const xy_pos_t pos = { mesh_index_to_xpos(icell.x), mesh_index_to_ypos(icell.y) };
+      xy_pos_t cell = raw - pos;
+
+      const float z_xmy0 = (z_x1y0 - z_x0y0) * RECIPROCAL(MESH_X_DIST),   // z slope per x along y0 (lower left to lower right)
+                  z_xmy1 = (z_x1y1 - z_x0y1) * RECIPROCAL(MESH_X_DIST);   // z slope per x along y1 (upper left to upper right)
+
+            float z_cxy0 = z_x0y0 + z_xmy0 * cell.x;        // z height along y0 at cell.x (changes for each cell.x in cell)
+
+      const float z_cxy1 = z_x0y1 + z_xmy1 * cell.x,        // z height along y1 at cell.x
+                  z_cxyd = z_cxy1 - z_cxy0;                 // z height difference along cell.x from y0 to y1
+
+            float z_cxym = z_cxyd * RECIPROCAL(MESH_Y_DIST); // z slope per y along cell.x from pos.y to y1 (changes for each cell.x in cell)
+
+      //    float z_cxcy = z_cxy0 + z_cxym * cell.y;        // interpolated mesh z height along cell.x at cell.y (do inside the segment loop)
+
+      // As subsequent segments step through this cell, the z_cxy0 intercept will change
+      // and the z_cxym slope will change, both as a function of cell.x within the cell, and
+      // each change by a constant for fixed segment lengths.
+
+      const float z_sxy0 = z_xmy0 * diff.x,                                       // per-segment adjustment to z_cxy0
+                  z_sxym = (z_xmy1 - z_xmy0) * RECIPROCAL(MESH_Y_DIST) * diff.x;  // per-segment adjustment to z_cxym
+
+      for (;;) {  // for all segments within this mesh cell
+
+        if (--segments == 0) raw = destination;     // if this is last segment, use destination for exact
+
+        const float z_cxcy = (z_cxy0 + z_cxym * cell.y) // interpolated mesh z height along cell.x at cell.y
+          #if ENABLED(ENABLE_LEVELING_FADE_HEIGHT)
+            * fade_scaling_factor                   // apply fade factor to interpolated mesh height
+          #endif
+        ;
+
+        planner.buffer_line(raw.x, raw.y, raw.z + z_cxcy, raw.e, scaled_fr_mm_s, active_extruder, segment_xyz_mm
+          #if ENABLED(SCARA_FEEDRATE_SCALING)
+            , inv_duration
+          #endif
+        );
+
+        if (segments == 0)                        // done with last segment
+          return false;                           // didn't set current from destination
+
+        raw += diff;
+        cell += diff;
+
+        if (!WITHIN(cell.x, 0, MESH_X_DIST) || !WITHIN(cell.y, 0, MESH_Y_DIST))    // done within this cell, break to next
+          break;
+
+        // Next segment still within same mesh cell, adjust the per-segment
+        // slope and intercept to compute next z height.
+
+        z_cxy0 += z_sxy0;   // adjust z_cxy0 by per-segment z_sxy0
+        z_cxym += z_sxym;   // adjust z_cxym by per-segment z_sxym
+
+      } // segment loop
+    } // cell loop
+
+    return false; // caller will update current_position
+  }
+
+#endif // UBL_SEGMENTED
+
+#endif // AUTO_BED_LEVELING_UBL