diff --git a/packages/scratch-core/src/conversion/leveling/core.py b/packages/scratch-core/src/conversion/leveling/core.py
index 4f6b5ec3..5cbc9c17 100644
--- a/packages/scratch-core/src/conversion/leveling/core.py
+++ b/packages/scratch-core/src/conversion/leveling/core.py
@@ -1,19 +1,12 @@
 import numpy as np
 from conversion.leveling import SurfaceTerms, LevelingResult
-from conversion.leveling.solver import (
-    fit_surface,
-    get_2d_grid,
-    compute_root_mean_square,
-)
-from conversion.leveling.solver.utils import compute_image_center
+
 from container_models.scan_image import ScanImage
 
+from surfalize import Surface
+
 
-def level_map(
-    scan_image: ScanImage,
-    terms: SurfaceTerms,
-    reference_point: tuple[float, float] | None = None,
-) -> LevelingResult:
+def level_map(scan_image: ScanImage, terms: SurfaceTerms) -> LevelingResult:
     """
     Compute the leveled map by fitting polynomial terms and subtracting them from the image data.
 
@@ -21,41 +14,33 @@ def level_map(
 
     :param scan_image: The scan image containing the image data to level.
     :param terms: The surface terms to use in the fitting. Note: terms can be combined using bit-operators.
-    :param reference_point: A tuple representing a reference point (X, Y) in physical coordinate space.
-        If provided, then the coordinates will be translated such that (X, Y) lies in the origin after translation.
-        If `None`, then the coordinates will be translated such that the center of the image lies in the origin.
     :returns: An instance of `LevelingResult` containing the leveled scan data and estimated physical parameters.
     """
-    if not reference_point:
-        reference_point = compute_image_center(scan_image)
-
-    # Build the 2D grids and translate in the opposite direction of `reference_point`
-    x_grid, y_grid = get_2d_grid(
-        scan_image, offset=(-reference_point[0], -reference_point[1])
-    )
-
-    # Get the point cloud (xs, ys, zs) for the numerical data
-    xs, ys, zs = (
-        x_grid[scan_image.valid_mask],
-        y_grid[scan_image.valid_mask],
-        scan_image.valid_data,
+    if terms == SurfaceTerms.NONE:
+        return LevelingResult(
+            leveled_map=scan_image.data,
+            fitted_surface=np.full_like(scan_image.data, 0.0),
+        )
+
+    surface = Surface(
+        height_data=scan_image.data,
+        step_x=scan_image.scale_x,
+        step_y=scan_image.scale_y,
     )
-
-    # Fit surface by solving the least-squares solution to a linear matrix equation
-    fitted_surface, physical_params = fit_surface(xs, ys, zs, terms)
-    fitted_surface_2d = np.full_like(scan_image.data, np.nan)
-    fitted_surface_2d[scan_image.valid_mask] = fitted_surface
-
-    # Compute the leveled map
-    leveled_map_2d = np.full_like(scan_image.data, np.nan)
-    leveled_map_2d[scan_image.valid_mask] = zs - fitted_surface
-
-    # Calculate RMS of residuals
-    residual_rms = compute_root_mean_square(leveled_map_2d)
-
-    return LevelingResult(
-        leveled_map=leveled_map_2d,
-        parameters=physical_params,
-        residual_rms=residual_rms,
-        fitted_surface=fitted_surface_2d,
+    match terms:
+        case (
+            SurfaceTerms.TILT_X
+            | SurfaceTerms.TILT_Y
+            | SurfaceTerms.ASTIG_45
+            | SurfaceTerms.PLANE
+        ):
+            degree = 1
+        case SurfaceTerms.DEFOCUS | SurfaceTerms.ASTIG_0 | SurfaceTerms.SPHERE:
+            degree = 2
+        case _:
+            degree = 0
+
+    leveled, trend = surface.detrend_polynomial(
+        degree=degree, inplace=False, return_trend=True
     )
+    return LevelingResult(leveled_map=leveled.data, fitted_surface=trend.data)
diff --git a/packages/scratch-core/src/conversion/leveling/data_types.py b/packages/scratch-core/src/conversion/leveling/data_types.py
index e78244c7..9ba82b0d 100644
--- a/packages/scratch-core/src/conversion/leveling/data_types.py
+++ b/packages/scratch-core/src/conversion/leveling/data_types.py
@@ -1,9 +1,7 @@
 from enum import Flag, auto
-import numpy as np
-from typing import Callable
 from pydantic import BaseModel
 
-from container_models.base import FloatArray1D, FloatArray2D
+from container_models.base import FloatArray2D
 
 
 class SurfaceTerms(Flag):
@@ -26,48 +24,13 @@ class SurfaceTerms(Flag):
     SPHERE = OFFSET | TILT_X | TILT_Y | ASTIG_45 | DEFOCUS | ASTIG_0
 
 
-# Mapping mathematical term indices to lambda functions for design matrix generation
-TERM_FUNCTIONS: dict[
-    SurfaceTerms, Callable[[FloatArray1D, FloatArray1D], FloatArray1D]
-] = {
-    SurfaceTerms.OFFSET: lambda xs, ys: np.ones_like(xs),
-    SurfaceTerms.TILT_X: lambda xs, ys: xs,
-    SurfaceTerms.TILT_Y: lambda xs, ys: ys,
-    SurfaceTerms.ASTIG_45: lambda xs, ys: xs * ys,
-    SurfaceTerms.DEFOCUS: lambda xs, ys: xs**2 + ys**2,
-    SurfaceTerms.ASTIG_0: lambda xs, ys: xs**2 - ys**2,
-}
-
-
 class LevelingResult(BaseModel, arbitrary_types_allowed=True):
     """
     Result of a leveling operation.
 
     :param leveled_map: 2D array with the leveled height data
-    :param parameters: Dictionary mapping SurfaceTerms to fitted coefficient values
-    :param residual_rms: Root mean square of residuals after leveling
     :param fitted_surface: 2D array of the fitted surface (same shape as `leveled_map`)
     """
 
     leveled_map: FloatArray2D
-    parameters: dict[SurfaceTerms, float]
-    residual_rms: float
     fitted_surface: FloatArray2D
-
-
-class NormalizedCoordinates(BaseModel, arbitrary_types_allowed=True):
-    """
-    Container class for storing centered and rescaled coordinates.
-
-    :param xs: The rescaled X-coordinates.
-    :param ys: The rescaled Y-coordinates.
-    :param x_mean: The mean of the X-coordinates before rescaling.
-    :param y_mean: The mean of the Y-coordinates before rescaling.
-    :param scale: The multiplier used for rescaling.
-    """
-
-    xs: FloatArray1D
-    ys: FloatArray1D
-    x_mean: float
-    y_mean: float
-    scale: float
diff --git a/packages/scratch-core/src/conversion/leveling/solver/__init__.py b/packages/scratch-core/src/conversion/leveling/solver/__init__.py
deleted file mode 100644
index f6cf3536..00000000
--- a/packages/scratch-core/src/conversion/leveling/solver/__init__.py
+++ /dev/null
@@ -1,17 +0,0 @@
-from conversion.leveling.solver.grid import get_2d_grid
-from conversion.leveling.solver.design import build_design_matrix
-from conversion.leveling.solver.transforms import (
-    denormalize_parameters,
-    normalize_coordinates,
-)
-from conversion.leveling.solver.utils import compute_root_mean_square
-from conversion.leveling.solver.core import fit_surface
-
-__all__ = (
-    "get_2d_grid",
-    "build_design_matrix",
-    "denormalize_parameters",
-    "fit_surface",
-    "normalize_coordinates",
-    "compute_root_mean_square",
-)
diff --git a/packages/scratch-core/src/conversion/leveling/solver/core.py b/packages/scratch-core/src/conversion/leveling/solver/core.py
deleted file mode 100644
index aea0d663..00000000
--- a/packages/scratch-core/src/conversion/leveling/solver/core.py
+++ /dev/null
@@ -1,46 +0,0 @@
-import numpy as np
-
-from container_models.base import FloatArray1D
-from conversion.leveling import SurfaceTerms
-from conversion.leveling.solver import (
-    normalize_coordinates,
-    build_design_matrix,
-    denormalize_parameters,
-)
-
-
-def fit_surface(
-    xs: FloatArray1D, ys: FloatArray1D, zs: FloatArray1D, terms: SurfaceTerms
-) -> tuple[FloatArray1D, dict[SurfaceTerms, float]]:
-    """
-    Core solver: fits a surface to the point cloud (xs, ys, zs).
-
-    :param xs: The X-coordinates.
-    :param ys: The Y-coordinates.
-    :param zs: The Z-values.
-    :param terms: The terms to use in the fitting
-    :return: A tuple containing the fitted surface (z̃s) and the estimated physical parameters.
-    """
-    # 1. Normalize the grid coordinates by centering and rescaling them
-    normalized = normalize_coordinates(xs, ys)
-
-    # 2. Build the design matrix for the least-squares solver
-    design_matrix = build_design_matrix(normalized.xs, normalized.ys, terms)
-
-    # 3. Solve (Least Squares)
-    (
-        coefficients,
-        *_,
-    ) = np.linalg.lstsq(design_matrix, zs, rcond=None)
-
-    # 4. Compute the surface (z̃s-values) from the fitted coefficients
-    fitted_surface = design_matrix @ coefficients
-
-    # 5. Recover physical parameters (optional usage, but part of original spec)
-    physical_params = denormalize_parameters(
-        dict(zip(terms, map(float, coefficients))),
-        normalized.x_mean,
-        normalized.y_mean,
-        normalized.scale,
-    )
-    return fitted_surface, physical_params
diff --git a/packages/scratch-core/src/conversion/leveling/solver/design.py b/packages/scratch-core/src/conversion/leveling/solver/design.py
deleted file mode 100644
index 76c353de..00000000
--- a/packages/scratch-core/src/conversion/leveling/solver/design.py
+++ /dev/null
@@ -1,26 +0,0 @@
-import numpy as np
-
-from container_models.base import FloatArray1D, FloatArray2D
-from conversion.leveling import SurfaceTerms
-from conversion.leveling.data_types import TERM_FUNCTIONS
-
-
-def build_design_matrix(
-    xs: FloatArray1D, ys: FloatArray1D, terms: SurfaceTerms
-) -> FloatArray2D:
-    """
-    Constructs the Least Squares design matrix based on grid coordinates (xs, ys) and requested terms.
-
-    :param xs: The X-coordinates.
-    :param ys: The Y-coordinates.
-    :param terms: The surface terms to use in the design matrix.
-    :returns: The design matrix as a numpy array with shape [n_points, n_terms].
-    """
-    num_points = xs.size
-    matrix = np.zeros((num_points, len(terms)), dtype=np.float64)
-
-    for column_index, term in enumerate(terms):
-        if func := TERM_FUNCTIONS.get(term):
-            matrix[:, column_index] = func(xs, ys)
-
-    return matrix
diff --git a/packages/scratch-core/src/conversion/leveling/solver/grid.py b/packages/scratch-core/src/conversion/leveling/solver/grid.py
deleted file mode 100644
index 8328cbd5..00000000
--- a/packages/scratch-core/src/conversion/leveling/solver/grid.py
+++ /dev/null
@@ -1,27 +0,0 @@
-import numpy as np
-
-from container_models.base import FloatArray2D
-from container_models.scan_image import ScanImage
-
-
-def get_2d_grid(
-    scan_image: ScanImage, offset: tuple[float, float] = (0, 0)
-) -> tuple[FloatArray2D, FloatArray2D]:
-    """
-    Return a 2D grid containing the physical coordinates of the scan data.
-
-    :param scan_image: An instance of `ScanImage` containing the recorded depth data.
-    :param offset: A tuple containing the physical coordinates of the image offset (in meters)
-        relative to the origin by which the grid coordinates need to be translated. The first element
-        corresponds to offset in the X-dimension, and the second element to the offset in the Y-dimension.
-    :returns: A tuple containing the grid coordinates for the X-direction and Y-direction.
-    """
-    # Generate Grid (ij indexing to match matrix coordinates)
-    x_indices, y_indices = np.meshgrid(
-        np.arange(scan_image.width), np.arange(scan_image.height), indexing="xy"
-    )
-    # Translate the grid by `offset`
-    x_grid = (x_indices * scan_image.scale_x) + offset[0]
-    y_grid = (y_indices * scan_image.scale_y) + offset[1]
-
-    return x_grid, y_grid
diff --git a/packages/scratch-core/src/conversion/leveling/solver/transforms.py b/packages/scratch-core/src/conversion/leveling/solver/transforms.py
deleted file mode 100644
index 9fb7f5bc..00000000
--- a/packages/scratch-core/src/conversion/leveling/solver/transforms.py
+++ /dev/null
@@ -1,85 +0,0 @@
-import numpy as np
-from collections.abc import Mapping
-
-from container_models.base import FloatArray1D
-from conversion.leveling import SurfaceTerms
-from conversion.leveling.data_types import NormalizedCoordinates
-
-
-def normalize_coordinates(xs: FloatArray1D, ys: FloatArray1D) -> NormalizedCoordinates:
-    """
-    Normalize grid coordinates by centering and rescaling.
-
-    This method is used to improve numerical stability during fitting.
-
-    :param xs: The X-coordinates to normalize.
-    :param ys: The Y-coordinates to normalize.
-    :returns: An instance of `NormalizedCoordinatesResult` containing the rescaled grid coordinates.
-    """
-    x_mean, y_mean = np.mean(xs), np.mean(ys)
-    vx_norm = xs - x_mean
-    vy_norm = ys - y_mean
-
-    span_x = np.max(vx_norm) - np.min(vx_norm)
-    span_y = np.max(vy_norm) - np.min(vy_norm)
-    # Avoid division by zero
-    max_span = max(span_x, span_y)
-    scale = 1.0 if np.isclose(max_span, 0.0) else 1 / max_span
-
-    return NormalizedCoordinates(
-        xs=vx_norm * scale,  # rescale X-coordinates
-        ys=vy_norm * scale,  # rescale Y-coordinates
-        x_mean=float(x_mean),
-        y_mean=float(y_mean),
-        scale=float(scale),
-    )
-
-
-def denormalize_parameters(
-    coefficients: Mapping[SurfaceTerms, float],
-    x_mean: float,
-    y_mean: float,
-    scale: float,
-) -> dict[SurfaceTerms, float]:
-    """
-    Converts normalized fit parameters back to real-world physical units.
-
-    The computation matches the specific numerical corrections from the original MATLAB script.
-
-    :param coefficients: A dictionary containing the normalized fit parameters.
-    :param x_mean: The mean of the x-coordinates.
-    :param y_mean: The mean of the y-coordinates.
-    :param scale: The scale factor.
-    :returns: A dictionary containing the denormalized fit parameters for all surface terms.
-    """
-    params = np.array(
-        [coefficients.get(term, 0.0) for term in SurfaceTerms], dtype=np.float64
-    )
-
-    # Un-normalize scaling
-    params[1:3] *= scale  # Tilts
-    params[3:] *= scale**2  # Quadratic terms
-
-    # Algebraic corrections for centering (x_mean, y_mean)
-    # Note: These formulas correspond exactly to the MATLAB implementation
-    # P[0] = Offset, P[1] = TiltX, P[2] = TiltY, etc.
-
-    # Adjust Offset (p0)
-    params[0] = (
-        params[0]
-        - params[1] * x_mean
-        - params[2] * y_mean
-        + params[3] * x_mean * y_mean
-        + params[4] * (x_mean**2 + y_mean**2)
-        + params[5] * (x_mean**2 - y_mean**2)
-    )
-    # Adjust Tilt X (p1)
-    params[1] = (
-        params[1] - params[3] * y_mean - 2 * params[4] * x_mean - 2 * params[5] * x_mean
-    )
-    # Adjust Tilt Y (p2)
-    params[2] = (
-        params[2] - params[3] * x_mean - 2 * params[4] * y_mean + 2 * params[5] * y_mean
-    )
-
-    return dict(zip(SurfaceTerms, map(float, params)))
diff --git a/packages/scratch-core/src/conversion/leveling/solver/utils.py b/packages/scratch-core/src/conversion/leveling/solver/utils.py
deleted file mode 100644
index 1cb7c52e..00000000
--- a/packages/scratch-core/src/conversion/leveling/solver/utils.py
+++ /dev/null
@@ -1,16 +0,0 @@
-import numpy as np
-
-from container_models.base import FloatArray
-from container_models.scan_image import ScanImage
-
-
-def compute_root_mean_square(data: FloatArray) -> float:
-    """Compute the root-mean-square from a data array and return as Python float."""
-    return float(np.sqrt(np.nanmean(data**2)))
-
-
-def compute_image_center(scan_image: ScanImage) -> tuple[float, float]:
-    """Compute the centerpoint (X, Y) of a scan image in physical coordinate space."""
-    center_x = (scan_image.width - 1) * scan_image.scale_x * 0.5
-    center_y = (scan_image.height - 1) * scan_image.scale_y * 0.5
-    return center_x, center_y
diff --git a/packages/scratch-core/src/conversion/preprocess_impression/preprocess_impression.py b/packages/scratch-core/src/conversion/preprocess_impression/preprocess_impression.py
index 7272661f..ecce2dd2 100644
--- a/packages/scratch-core/src/conversion/preprocess_impression/preprocess_impression.py
+++ b/packages/scratch-core/src/conversion/preprocess_impression/preprocess_impression.py
@@ -89,7 +89,7 @@ def preprocess_impression_mark(
         )
 
     # Stage 8: Final leveling
-    mark_filtered, _ = _level_mark(mark_filtered, params.surface_terms, mark.center)
+    mark_filtered, _ = _level_mark(mark_filtered, params.surface_terms)
 
     # Prepare leveled-only output
     mark_leveled_final = _finalize_leveled_output(
@@ -109,11 +109,8 @@ def preprocess_impression_mark(
 def _level_mark(
     mark: Mark,
     terms: SurfaceTerms,
-    reference_point: Point2D | None = None,
 ) -> tuple[Mark, DepthData]:
-    result = level_map(
-        mark.scan_image, terms=terms, reference_point=reference_point or mark.center
-    )
+    result = level_map(mark.scan_image, terms=terms)
     leveled_mark = update_mark_data(mark, result.leveled_map)
     return leveled_mark, result.fitted_surface
 
@@ -156,6 +153,5 @@ def _finalize_leveled_output(
 
     # Apply PLANE-only leveling (after resampling, like MATLAB)
     rigid_terms = surface_terms & SurfaceTerms.PLANE
-    leveled_mark, _ = _level_mark(mark_restored, rigid_terms, reference_point)
-
+    leveled_mark, _ = _level_mark(mark_restored, rigid_terms)
     return leveled_mark
diff --git a/packages/scratch-core/tests/conversion/leveling/resources/baseline_design_matrix.npy b/packages/scratch-core/tests/conversion/leveling/resources/baseline_design_matrix.npy
deleted file mode 100644
index 76fbf21b..00000000
Binary files a/packages/scratch-core/tests/conversion/leveling/resources/baseline_design_matrix.npy and /dev/null differ
diff --git a/packages/scratch-core/tests/conversion/leveling/solver/test_core.py b/packages/scratch-core/tests/conversion/leveling/solver/test_core.py
deleted file mode 100644
index 7d3e31fe..00000000
--- a/packages/scratch-core/tests/conversion/leveling/solver/test_core.py
+++ /dev/null
@@ -1,110 +0,0 @@
-from container_models.base import FloatArray1D
-from conversion.leveling import SurfaceTerms
-from conversion.leveling.solver import fit_surface
-import numpy as np
-import pytest
-
-from ..constants import SINGLE_AND_COMBINED_TERMS, SINGLE_TERMS
-
-
-@pytest.mark.parametrize("terms", SINGLE_AND_COMBINED_TERMS)
-def test_fit_surface_reduces_variance(
-    xs: FloatArray1D,
-    ys: FloatArray1D,
-    zs: FloatArray1D,
-    terms: SurfaceTerms,
-):
-    fitted_surface, _ = fit_surface(xs=xs, ys=ys, zs=zs, terms=terms)
-    leveled_map = zs - fitted_surface
-    assert np.var(fitted_surface) < np.var(zs)
-    assert np.var(fitted_surface) < np.var(leveled_map)
-
-
-@pytest.mark.parametrize("terms", list(SurfaceTerms.PLANE))
-def test_fit_surface_plane_reduces_variance(
-    xs: FloatArray1D,
-    ys: FloatArray1D,
-    zs: FloatArray1D,
-    terms: SurfaceTerms,
-):
-    leveled_map_terms = zs - fit_surface(xs=xs, ys=ys, zs=zs, terms=terms)[0]
-    leveled_map_plane = (
-        zs - fit_surface(xs=xs, ys=ys, zs=zs, terms=SurfaceTerms.PLANE)[0]
-    )
-
-    single_term_var = np.var(leveled_map_terms)
-    plane_var = np.var(leveled_map_plane)
-
-    assert plane_var < single_term_var or np.isclose(plane_var, single_term_var)
-
-
-@pytest.mark.parametrize("terms", SINGLE_TERMS + [SurfaceTerms.PLANE])
-def test_fit_surface_sphere_reduces_variance(
-    xs: FloatArray1D,
-    ys: FloatArray1D,
-    zs: FloatArray1D,
-    terms: SurfaceTerms,
-):
-    leveled_map_terms = zs - fit_surface(xs=xs, ys=ys, zs=zs, terms=terms)[0]
-    leveled_map_sphere = (
-        zs - fit_surface(xs=xs, ys=ys, zs=zs, terms=SurfaceTerms.SPHERE)[0]
-    )
-
-    assert np.var(leveled_map_sphere) < np.var(leveled_map_terms)
-
-
-@pytest.mark.parametrize("terms", SINGLE_AND_COMBINED_TERMS)
-def test_fit_surface_fits_terms(
-    xs: FloatArray1D,
-    ys: FloatArray1D,
-    zs: FloatArray1D,
-    terms: SurfaceTerms,
-):
-    _, physical_params = fit_surface(xs=xs, ys=ys, zs=zs, terms=terms)
-    assert all(not np.isclose(v, 0.0) for p, v in physical_params.items() if p in terms)
-    assert all(np.isclose(v, 0.0) for p, v in physical_params.items() if p not in terms)
-
-
-def test_fit_surface_offset_is_invariant_for_tilted_surfaces(
-    xs: FloatArray1D, ys: FloatArray1D, zs: FloatArray1D
-):
-    _, physical_params_tilt_x = fit_surface(
-        xs=xs, ys=ys, zs=zs, terms=SurfaceTerms.TILT_X
-    )
-    _, physical_params_tilt_y = fit_surface(
-        xs=xs, ys=ys, zs=zs, terms=SurfaceTerms.TILT_Y
-    )
-    assert np.isclose(
-        physical_params_tilt_x[SurfaceTerms.OFFSET],
-        physical_params_tilt_y[SurfaceTerms.OFFSET],
-    )
-
-
-def test_fit_surface_offset_equals_mean(
-    xs: FloatArray1D, ys: FloatArray1D, zs: FloatArray1D
-):
-    fitted_surface_offset, physical_params_offset = fit_surface(
-        xs=xs, ys=ys, zs=zs, terms=SurfaceTerms.OFFSET
-    )
-    param = physical_params_offset[SurfaceTerms.OFFSET]
-    assert np.isclose(np.mean(zs), param)
-    assert np.allclose(fitted_surface_offset, param)
-
-
-def test_fit_surface_defocus_is_positive(
-    xs: FloatArray1D, ys: FloatArray1D, zs: FloatArray1D
-):
-    fitted_surface_defocus, _ = fit_surface(
-        xs=xs, ys=ys, zs=zs, terms=SurfaceTerms.DEFOCUS
-    )
-    assert np.all(fitted_surface_defocus > 0.0)
-
-
-def test_fit_surface_none_has_no_effect(
-    xs: FloatArray1D, ys: FloatArray1D, zs: FloatArray1D
-):
-    fitted_surface_none, physical_params_none = fit_surface(
-        xs=xs, ys=ys, zs=zs, terms=SurfaceTerms.NONE
-    )
-    assert np.allclose(fitted_surface_none, 0.0)
-    assert all(np.isclose(v, 0.0) for v in physical_params_none.values())
diff --git a/packages/scratch-core/tests/conversion/leveling/solver/test_design.py b/packages/scratch-core/tests/conversion/leveling/solver/test_design.py
deleted file mode 100644
index b532d45b..00000000
--- a/packages/scratch-core/tests/conversion/leveling/solver/test_design.py
+++ /dev/null
@@ -1,21 +0,0 @@
-from container_models.base import FloatArray1D
-from conversion.leveling import SurfaceTerms
-from conversion.leveling.solver import build_design_matrix
-import numpy as np
-import pytest
-
-from ..constants import SINGLE_AND_COMBINED_TERMS, RESOURCES_DIR
-
-
-@pytest.mark.parametrize("terms", SINGLE_AND_COMBINED_TERMS)
-def test_design_matrix_shape_matches_number_of_terms(
-    xs: FloatArray1D, ys: FloatArray1D, terms: SurfaceTerms
-):
-    design_matrix = build_design_matrix(xs=xs, ys=ys, terms=terms)
-    assert design_matrix.shape == (len(xs), len(terms))
-
-
-def test_design_matrix_matches_baseline_output(xs: FloatArray1D, ys: FloatArray1D):
-    design_matrix = build_design_matrix(xs=xs, ys=ys, terms=SurfaceTerms.SPHERE)
-    verified = np.load(RESOURCES_DIR / "baseline_design_matrix.npy")
-    assert np.allclose(design_matrix, verified)
diff --git a/packages/scratch-core/tests/conversion/leveling/solver/test_grid.py b/packages/scratch-core/tests/conversion/leveling/solver/test_grid.py
deleted file mode 100644
index fbf63087..00000000
--- a/packages/scratch-core/tests/conversion/leveling/solver/test_grid.py
+++ /dev/null
@@ -1,62 +0,0 @@
-from conversion.leveling.solver import get_2d_grid
-from conversion.leveling.solver.utils import compute_image_center
-from container_models.scan_image import ScanImage
-import numpy as np
-import pytest
-
-
-def test_grid_matches_scan_image_shape(scan_image_rectangular_with_nans: ScanImage):
-    x_grid, y_grid = get_2d_grid(scan_image_rectangular_with_nans)
-    assert x_grid.shape == (
-        scan_image_rectangular_with_nans.height,
-        scan_image_rectangular_with_nans.width,
-    )
-    assert y_grid.shape == (
-        scan_image_rectangular_with_nans.height,
-        scan_image_rectangular_with_nans.width,
-    )
-
-
-def test_grid_is_a_meshgrid(scan_image_rectangular_with_nans: ScanImage):
-    x_grid, y_grid = get_2d_grid(scan_image_rectangular_with_nans)
-
-    assert all(
-        np.array_equal(x_grid[i, :], x_grid[i + 1, :])
-        for i in range(x_grid.shape[0] - 1)
-    )
-    assert all(
-        np.array_equal(y_grid[:, i], y_grid[:, i + 1])
-        for i in range(y_grid.shape[1] - 1)
-    )
-    assert np.all(x_grid[0, :-1] < x_grid[0, 1:])
-    assert np.all(y_grid[:-1, 0] < y_grid[1:, 0])
-
-
-def test_grid_is_translated_by_offset(scan_image_rectangular_with_nans: ScanImage):
-    si = scan_image_rectangular_with_nans
-    offset_x, offset_y = 0.1, 0.5
-    x_grid, y_grid = get_2d_grid(si, offset=(offset_x, offset_y))
-
-    xs, ys = x_grid[0, :], y_grid[:, 0]
-
-    assert np.isclose(xs[0], offset_x)
-    assert np.isclose(xs[-1], offset_x + (si.width - 1) * si.scale_x)
-    assert np.isclose(ys[0], offset_y)
-    assert np.isclose(ys[-1], offset_y + (si.height - 1) * si.scale_y)
-
-
-@pytest.mark.integration
-def test_grid_is_centered_around_origin(
-    scan_image_rectangular_with_nans: ScanImage,
-):
-    center_x, center_y = compute_image_center(scan_image_rectangular_with_nans)
-    x_grid, y_grid = get_2d_grid(
-        scan_image_rectangular_with_nans, offset=(-center_x, -center_y)
-    )
-
-    # With xy indexing: x from row 0, y from column 0
-    mid_x, mid_y = x_grid.shape[1] // 2, y_grid.shape[0] // 2
-    xs, ys = x_grid[0, :], y_grid[:, 0]
-
-    assert xs[mid_x - 1] < 0.0 < xs[mid_x]
-    assert ys[mid_y - 1] < 0.0 < ys[mid_y]
diff --git a/packages/scratch-core/tests/conversion/leveling/solver/test_transforms.py b/packages/scratch-core/tests/conversion/leveling/solver/test_transforms.py
deleted file mode 100644
index 3ce7e3e5..00000000
--- a/packages/scratch-core/tests/conversion/leveling/solver/test_transforms.py
+++ /dev/null
@@ -1,74 +0,0 @@
-import numpy as np
-import pytest
-
-from container_models.base import FloatArray1D
-from conversion.leveling import SurfaceTerms
-from conversion.leveling.solver import normalize_coordinates, denormalize_parameters
-from ..constants import SINGLE_TERMS
-
-
-class TestNormalizeCoordinates:
-    def test_normalized_coordinates_have_zero_mean(
-        self, xs: FloatArray1D, ys: FloatArray1D
-    ):
-        normalized = normalize_coordinates(xs=xs, ys=ys)
-        assert np.isclose(np.mean(normalized.xs), 0.0)
-        assert np.isclose(np.mean(normalized.ys), 0.0)
-        assert normalized.xs[0] < 0 < normalized.xs[-1]
-        assert normalized.ys[0] < 0 < normalized.ys[-1]
-
-    def test_normalized_coordinates_are_strictly_increasing(
-        self, xs: FloatArray1D, ys: FloatArray1D
-    ):
-        normalized = normalize_coordinates(xs=xs, ys=ys)
-        assert np.all(normalized.xs[:-1] < normalized.xs[1:])
-        assert np.all(normalized.ys[:-1] < normalized.ys[1:])
-
-    def test_normalize_coordinates_returns_means(
-        self, xs: FloatArray1D, ys: FloatArray1D
-    ):
-        normalized = normalize_coordinates(xs=xs, ys=ys)
-        assert np.isclose(normalized.x_mean, np.mean(xs))
-        assert np.isclose(normalized.y_mean, np.mean(ys))
-
-    def test_normalized_coordinates_are_bounded_by_unit_disk(
-        self, xs: FloatArray1D, ys: FloatArray1D
-    ):
-        normalized = normalize_coordinates(xs=xs, ys=ys)
-        assert np.isclose(normalized.xs[-1] - normalized.xs[0], 1.0)
-        assert np.isclose(normalized.ys[-1] - normalized.ys[0], 0.5)
-
-
-class TestDenormalizeParameters:
-    @pytest.mark.parametrize("terms", SINGLE_TERMS)
-    def test_denormalize_parameters_returns_all_terms(self, terms: SurfaceTerms):
-        params = denormalize_parameters(
-            coefficients={terms: 1.5}, x_mean=0.1, y_mean=0.5, scale=0.21
-        )
-        assert all(term in params for term in SurfaceTerms)
-
-    @pytest.mark.parametrize("terms", SINGLE_TERMS)
-    def test_denormalize_parameters_adjusts_terms(self, terms: SurfaceTerms):
-        initial_value = 1.5
-        params = denormalize_parameters(
-            coefficients={terms: initial_value}, x_mean=0.1, y_mean=0.5, scale=0.21
-        )
-        assert all(
-            not np.isclose(params[t], initial_value)
-            for t in SurfaceTerms
-            if t in terms and t != SurfaceTerms.OFFSET
-        )
-
-    def test_denormalize_parameters_matches_baseline_output(self):
-        params = denormalize_parameters(
-            coefficients={term: 1.0 for term in SurfaceTerms},
-            x_mean=0.1,
-            y_mean=0.5,
-            scale=0.21,
-        )
-        assert np.isclose(params[SurfaceTerms.OFFSET], 0.877087)
-        assert np.isclose(params[SurfaceTerms.TILT_X], 0.17031)
-        assert np.isclose(params[SurfaceTerms.TILT_Y], 0.20559)
-        assert np.isclose(params[SurfaceTerms.ASTIG_45], 0.0441)
-        assert np.isclose(params[SurfaceTerms.DEFOCUS], 0.0441)
-        assert np.isclose(params[SurfaceTerms.ASTIG_0], 0.0441)
diff --git a/packages/scratch-core/tests/conversion/leveling/solver/test_utils.py b/packages/scratch-core/tests/conversion/leveling/solver/test_utils.py
deleted file mode 100644
index 8e515259..00000000
--- a/packages/scratch-core/tests/conversion/leveling/solver/test_utils.py
+++ /dev/null
@@ -1,45 +0,0 @@
-import numpy as np
-import pytest
-
-from conversion.leveling.solver import compute_root_mean_square
-from conversion.leveling.solver.utils import compute_image_center
-from container_models.scan_image import ScanImage
-
-
-class TestRootMeanSquare:
-    @pytest.mark.parametrize("value", [0.0, 1.0, 2.0, -3.15, 40.123, -80, 100])
-    def test_rms_is_constant_for_constant_input(self, value: float):
-        result = compute_root_mean_square(np.array([value] * 100))
-        assert np.isclose(result, abs(value))
-
-    @pytest.mark.parametrize("value", list(range(-10, 10)))
-    def test_rms_is_non_negative(self, value: float):
-        result = compute_root_mean_square(np.array([value] * 100))
-        assert result > 0 or np.isclose(result, 0.0)
-
-    def test_rms_can_handle_nans(self):
-        array_with_nans = np.array(
-            [0, 1, 0.15, 2, np.nan, 3, 4, np.nan, -5], dtype=np.float64
-        )
-        array_without_nans = array_with_nans[~np.isnan(array_with_nans)]
-
-        result_with_nans = compute_root_mean_square(array_with_nans)
-        result_without_nans = compute_root_mean_square(array_without_nans)
-        assert np.isclose(result_with_nans, result_without_nans)
-        assert np.isclose(result_with_nans, 2.8036328473709147)
-
-
-class TestComputeImageCenter:
-    def test_compute_image_center_for_square_image(
-        self, scan_image_with_nans: ScanImage
-    ):
-        center_x, center_y = compute_image_center(scan_image_with_nans)
-        assert np.isclose(center_x, 0.0004037151235)
-        assert np.isclose(center_y, center_x)
-
-    def test_compute_image_center_for_rectangular_image(
-        self, scan_image_rectangular_with_nans: ScanImage
-    ):
-        center_x, center_y = compute_image_center(scan_image_rectangular_with_nans)
-        assert np.isclose(center_x, 0.00030245829675)
-        assert np.isclose(center_y, 0.0004037151235)
diff --git a/packages/scratch-core/tests/conversion/leveling/test_level_map.py b/packages/scratch-core/tests/conversion/leveling/test_level_map.py
index d217c55d..f849bccc 100644
--- a/packages/scratch-core/tests/conversion/leveling/test_level_map.py
+++ b/packages/scratch-core/tests/conversion/leveling/test_level_map.py
@@ -11,12 +11,6 @@ def test_map_level_sphere(scan_image_with_nans: ScanImage):
     result = level_map(scan_image_with_nans, SurfaceTerms.SPHERE)
     assert result
     assert np.allclose(result.leveled_map, verified, equal_nan=True)
-    assert all(
-        np.isclose(result.parameters[p], 0.0)
-        for p in SurfaceTerms
-        if p not in SurfaceTerms.SPHERE
-    )
-    assert all(not np.isclose(result.parameters[p], 0.0) for p in SurfaceTerms.SPHERE)
 
 
 @pytest.mark.integration
@@ -25,12 +19,6 @@ def test_map_level_plane(scan_image_with_nans: ScanImage):
     result = level_map(scan_image_with_nans, SurfaceTerms.PLANE)
     assert result
     assert np.allclose(result.leveled_map, verified, equal_nan=True)
-    assert all(
-        np.isclose(result.parameters[p], 0.0)
-        for p in SurfaceTerms
-        if p not in SurfaceTerms.PLANE
-    )
-    assert all(not np.isclose(result.parameters[p], 0.0) for p in SurfaceTerms.PLANE)
 
 
 @pytest.mark.integration
@@ -38,7 +26,6 @@ def test_map_level_none(scan_image_with_nans: ScanImage):
     result = level_map(scan_image_with_nans, SurfaceTerms.NONE)
     assert result
     assert np.allclose(result.leveled_map, scan_image_with_nans.data, equal_nan=True)
-    assert all(np.isclose(result.parameters[p], 0.0) for p in SurfaceTerms)
 
 
 @pytest.mark.integration
@@ -51,35 +38,3 @@ def test_map_level_offset(scan_image_with_nans: ScanImage):
         scan_image_with_nans.data,
         equal_nan=True,
     )
-    assert result.parameters[SurfaceTerms.OFFSET] != 0
-    assert all(
-        np.isclose(result.parameters[p], 0.0)
-        for p in SurfaceTerms
-        if p != SurfaceTerms.OFFSET
-    )
-
-
-@pytest.mark.integration
-def test_map_level_reference_point_has_no_effect_with_none(
-    scan_image_with_nans: ScanImage,
-):
-    result_centered = level_map(scan_image_with_nans, SurfaceTerms.NONE)
-    result_ref = level_map(
-        scan_image_with_nans, SurfaceTerms.NONE, reference_point=(10.5, -5.2)
-    )
-    assert np.allclose(
-        result_centered.leveled_map, result_ref.leveled_map, equal_nan=True
-    )
-
-
-@pytest.mark.integration
-def test_map_level_reference_point_has_effect_with_plane(
-    scan_image_with_nans: ScanImage,
-):
-    result_centered = level_map(scan_image_with_nans, SurfaceTerms.PLANE)
-    result_ref = level_map(
-        scan_image_with_nans, SurfaceTerms.NONE, reference_point=(10.5, -5.2)
-    )
-    assert not np.allclose(
-        result_centered.leveled_map, result_ref.leveled_map, equal_nan=True
-    )