diff --git a/Marlin/Marlin_main.cpp b/Marlin/Marlin_main.cpp
index 8593b20581..c02c02c746 100755
--- a/Marlin/Marlin_main.cpp
+++ b/Marlin/Marlin_main.cpp
@@ -2224,11 +2224,15 @@ static void clean_up_after_endstop_or_probe_move() {
   void set_bed_leveling_enabled(bool enable=true) {
     #if ENABLED(MESH_BED_LEVELING)
 
-      if (!enable && mbl.active())
-        current_position[Z_AXIS] +=
-          mbl.get_z(RAW_CURRENT_POSITION(X_AXIS), RAW_CURRENT_POSITION(Y_AXIS)) - (MESH_HOME_SEARCH_Z);
+      if (enable != mbl.active()) {
 
-      mbl.set_active(enable && mbl.has_mesh()); // was set_has_mesh(). Is this not correct?
+        if (!enable)
+          planner.apply_leveling(current_position[X_AXIS], current_position[Y_AXIS], current_position[Z_AXIS]);
+
+        mbl.set_active(enable && mbl.has_mesh());
+
+        if (enable) planner.unapply_leveling(current_position);
+      }
 
     #elif HAS_ABL
 
@@ -3162,8 +3166,10 @@ inline void gcode_G4() {
     #elif ENABLED(MESH_BED_LEVELING)
       SERIAL_ECHOPGM("Mesh Bed Leveling");
       if (mbl.active()) {
+        float lz = current_position[Z_AXIS];
+        planner.apply_leveling(current_position[X_AXIS], current_position[Y_AXIS], lz);
         SERIAL_ECHOLNPGM(" (enabled)");
-        SERIAL_ECHOPAIR("MBL Adjustment Z", mbl.get_z(RAW_CURRENT_POSITION(X_AXIS), RAW_CURRENT_POSITION(Y_AXIS)));
+        SERIAL_ECHOPAIR("MBL Adjustment Z", lz);
       }
       SERIAL_EOL;
     #endif
@@ -3321,13 +3327,15 @@ inline void gcode_G28() {
       #if ENABLED(DEBUG_LEVELING_FEATURE)
         if (DEBUGGING(LEVELING)) SERIAL_ECHOLNPGM("MBL was active");
       #endif
-      // Save known Z position if already homed
+      // Use known Z position if already homed
       if (axis_homed[X_AXIS] && axis_homed[Y_AXIS] && axis_homed[Z_AXIS]) {
+        set_bed_leveling_enabled(false);
         pre_home_z = current_position[Z_AXIS];
-        pre_home_z += mbl.get_z(RAW_CURRENT_POSITION(X_AXIS), RAW_CURRENT_POSITION(Y_AXIS));
       }
-      mbl.set_active(false);
-      current_position[Z_AXIS] = pre_home_z;
+      else {
+        mbl.set_active(false);
+        current_position[Z_AXIS] = pre_home_z;
+      }
       #if ENABLED(DEBUG_LEVELING_FEATURE)
         if (DEBUGGING(LEVELING)) DEBUG_POS("Set Z to pre_home_z", current_position);
       #endif
@@ -3703,8 +3711,8 @@ inline void gcode_G28() {
 
       case MeshReset:
         if (mbl.active()) {
-          current_position[Z_AXIS] +=
-            mbl.get_z(RAW_CURRENT_POSITION(X_AXIS), RAW_CURRENT_POSITION(Y_AXIS)) - MESH_HOME_SEARCH_Z;
+          current_position[Z_AXIS] -= MESH_HOME_SEARCH_Z;
+          planner.apply_leveling(current_position[X_AXIS], current_position[Y_AXIS], current_position[Z_AXIS]);
           mbl.reset();
           SYNC_PLAN_POSITION_KINEMATIC();
         }
@@ -7640,9 +7648,12 @@ void tool_change(const uint8_t tmp_extruder, const float fr_mm_s/*=0.0*/, bool n
                 #if ENABLED(DEBUG_LEVELING_FEATURE)
                   if (DEBUGGING(LEVELING)) SERIAL_ECHOPAIR("Z before MBL: ", current_position[Z_AXIS]);
                 #endif
-                float xpos = RAW_CURRENT_POSITION(X_AXIS),
-                      ypos = RAW_CURRENT_POSITION(Y_AXIS);
-                current_position[Z_AXIS] += mbl.get_z(xpos + xydiff[X_AXIS], ypos + xydiff[Y_AXIS]) - mbl.get_z(xpos, ypos);
+                float x2 = current_position[X_AXIS] + xydiff[X_AXIS],
+                      y2 = current_position[Y_AXIS] + xydiff[Y_AXIS],
+                      z1 = current_position[Z_AXIS], z2 = z1;
+                planner.apply_leveling(current_position[X_AXIS], current_position[Y_AXIS], z1);
+                planner.apply_leveling(x2, y2, z2);
+                current_position[Z_AXIS] += z2 - z1;
                 #if ENABLED(DEBUG_LEVELING_FEATURE)
                   if (DEBUGGING(LEVELING))
                     SERIAL_ECHOLNPAIR(" after: ", current_position[Z_AXIS]);
diff --git a/Marlin/mesh_bed_leveling.h b/Marlin/mesh_bed_leveling.h
index 0e66bb8882..217cd60340 100644
--- a/Marlin/mesh_bed_leveling.h
+++ b/Marlin/mesh_bed_leveling.h
@@ -36,54 +36,58 @@
 
     void reset();
 
-    static FORCE_INLINE float get_probe_x(int8_t i) { return MESH_MIN_X + (MESH_X_DIST) * i; }
-    static FORCE_INLINE float get_probe_y(int8_t i) { return MESH_MIN_Y + (MESH_Y_DIST) * i; }
-    void set_z(const int8_t px, const int8_t py, const float z) { z_values[py][px] = z; }
+    static FORCE_INLINE float get_probe_x(const int8_t i) { return MESH_MIN_X + (MESH_X_DIST) * i; }
+    static FORCE_INLINE float get_probe_y(const int8_t i) { return MESH_MIN_Y + (MESH_Y_DIST) * i; }
+    void set_z(const int8_t px, const int8_t py, const float &z) { z_values[py][px] = z; }
 
-    bool active()                 { return TEST(status, MBL_STATUS_ACTIVE_BIT); }
-    void set_active(bool onOff)   { if (onOff) SBI(status, MBL_STATUS_ACTIVE_BIT); else CBI(status, MBL_STATUS_ACTIVE_BIT); }
-    bool has_mesh()               { return TEST(status, MBL_STATUS_HAS_MESH_BIT); }
-    void set_has_mesh(bool onOff) { if (onOff) SBI(status, MBL_STATUS_HAS_MESH_BIT); else CBI(status, MBL_STATUS_HAS_MESH_BIT); }
+    bool active() const                 { return TEST(status, MBL_STATUS_ACTIVE_BIT); }
+    void set_active(const bool onOff)   { onOff ? SBI(status, MBL_STATUS_ACTIVE_BIT) : CBI(status, MBL_STATUS_ACTIVE_BIT); }
+    bool has_mesh() const               { return TEST(status, MBL_STATUS_HAS_MESH_BIT); }
+    void set_has_mesh(const bool onOff) { onOff ? SBI(status, MBL_STATUS_HAS_MESH_BIT) : CBI(status, MBL_STATUS_HAS_MESH_BIT); }
 
-    inline void zigzag(int8_t index, int8_t &px, int8_t &py) {
+    inline void zigzag(const int8_t index, int8_t &px, int8_t &py) const {
       px = index % (MESH_NUM_X_POINTS);
       py = index / (MESH_NUM_X_POINTS);
       if (py & 1) px = (MESH_NUM_X_POINTS - 1) - px; // Zig zag
     }
 
-    void set_zigzag_z(int8_t index, float z) {
+    void set_zigzag_z(const int8_t index, const float &z) {
       int8_t px, py;
       zigzag(index, px, py);
       set_z(px, py, z);
     }
 
-    int8_t cell_index_x(float x) {
+    int8_t cell_index_x(const float &x) const {
       int8_t cx = (x - (MESH_MIN_X)) * (1.0 / (MESH_X_DIST));
       return constrain(cx, 0, (MESH_NUM_X_POINTS) - 2);
     }
 
-    int8_t cell_index_y(float y) {
+    int8_t cell_index_y(const float &y) const {
       int8_t cy = (y - (MESH_MIN_Y)) * (1.0 / (MESH_Y_DIST));
       return constrain(cy, 0, (MESH_NUM_Y_POINTS) - 2);
     }
 
-    int8_t probe_index_x(float x) {
+    int8_t probe_index_x(const float &x) const {
       int8_t px = (x - (MESH_MIN_X) + (MESH_X_DIST) * 0.5) * (1.0 / (MESH_X_DIST));
       return (px >= 0 && px < (MESH_NUM_X_POINTS)) ? px : -1;
     }
 
-    int8_t probe_index_y(float y) {
+    int8_t probe_index_y(const float &y) const {
       int8_t py = (y - (MESH_MIN_Y) + (MESH_Y_DIST) * 0.5) * (1.0 / (MESH_Y_DIST));
       return (py >= 0 && py < (MESH_NUM_Y_POINTS)) ? py : -1;
     }
 
-    float calc_z0(float a0, float a1, float z1, float a2, float z2) {
-      float delta_z = (z2 - z1) / (a2 - a1);
-      float delta_a = a0 - a1;
+    float calc_z0(const float &a0, const float &a1, const float &z1, const float &a2, const float &z2) const {
+      const float delta_z = (z2 - z1) / (a2 - a1);
+      const float delta_a = a0 - a1;
       return z1 + delta_a * delta_z;
     }
 
-    float get_z(float x0, float y0) {
+    float get_z(const float &x0, const float &y0
+      #if ENABLED(ENABLE_LEVELING_FADE_HEIGHT)
+        , const float &factor
+      #endif
+    ) const {
       int8_t cx = cell_index_x(x0),
              cy = cell_index_y(y0);
       if (cx < 0 || cy < 0) return z_offset;
@@ -96,7 +100,12 @@
       float z0 = calc_z0(y0,
                          get_probe_y(cy), z1,
                          get_probe_y(cy + 1), z2);
-      return z0 + z_offset;
+
+      #if ENABLED(ENABLE_LEVELING_FADE_HEIGHT)
+        return z0 * factor + z_offset;
+      #else
+        return z0 + z_offset;
+      #endif
     }
   };
 
diff --git a/Marlin/planner.cpp b/Marlin/planner.cpp
index ea9dcf284d..80cfb550ff 100644
--- a/Marlin/planner.cpp
+++ b/Marlin/planner.cpp
@@ -559,7 +559,7 @@ void Planner::check_axes_activity() {
     #if ENABLED(MESH_BED_LEVELING)
 
       if (mbl.active())
-        lz += mbl.get_z(RAW_X_POSITION(lx), RAW_Y_POSITION(ly)) * z_fade_factor;
+        lz += mbl.get_z(RAW_X_POSITION(lx), RAW_Y_POSITION(ly), z_fade_factor);
 
     #elif ABL_PLANAR
 
@@ -595,7 +595,7 @@ void Planner::check_axes_activity() {
 
       if (mbl.active()) {
         #if ENABLED(ENABLE_LEVELING_FADE_HEIGHT)
-          const float c = mbl.get_z(RAW_X_POSITION(logical[X_AXIS]), RAW_Y_POSITION(logical[Y_AXIS]));
+          const float c = mbl.get_z(RAW_X_POSITION(logical[X_AXIS]), RAW_Y_POSITION(logical[Y_AXIS]), 1.0);
           logical[Z_AXIS] = (z_fade_height * (RAW_Z_POSITION(logical[Z_AXIS]) - c)) / (z_fade_height - c);
         #else
           logical[Z_AXIS] -= mbl.get_z(RAW_X_POSITION(logical[X_AXIS]), RAW_Y_POSITION(logical[Y_AXIS]));