diff --git a/sdcflows/transform.py b/sdcflows/transform.py
index ebf23c7f6f..7247b97c2f 100644
--- a/sdcflows/transform.py
+++ b/sdcflows/transform.py
@@ -296,11 +296,11 @@ def fit(
         projected_reference, _ = ensure_positive_cosines(projected_reference)
 
         # Approximate only if the coordinate systems are not aligned
-        finest_coeffs = listify(self.coeffs)[-1]
+        coeffs = listify(self.coeffs)
         approx &= not np.allclose(
             np.linalg.norm(
                 np.cross(
-                    finest_coeffs.affine[:-1, :-1].T,
+                    coeffs[-1].affine[:-1, :-1].T,
                     target_reference.affine[:-1, :-1].T,
                 ),
                 axis=1,
@@ -310,26 +310,26 @@ def fit(
         )
 
         weights = []
-        coeffs = []
         if approx:
             from sdcflows.utils.tools import deoblique_and_zooms
 
             # Generate a sampling reference on the fieldmap's space that fully covers
             # the target_reference's grid.
             projected_reference = deoblique_and_zooms(
-                listify(self.coeffs)[-1],
+                coeffs[-1],
                 target_reference,
             )
 
         # Generate tensor-product B-Spline weights
-        for level in listify(self.coeffs):
+        coeffs_data = []
+        for level in coeffs:
             wmat = grid_bspline_weights(target_reference, level)
             weights.append(wmat)
-            coeffs.append(level.get_fdata(dtype="float32").reshape(-1))
+            coeffs_data.append(level.get_fdata(dtype="float32").reshape(-1))
 
         # Reconstruct the fieldmap (in Hz) from coefficients
         fmap = np.zeros(projected_reference.shape[:3], dtype="float32")
-        fmap = (np.squeeze(np.hstack(coeffs).T) @ sparse_vstack(weights)).reshape(
+        fmap = (np.squeeze(np.hstack(coeffs_data).T) @ sparse_vstack(weights)).reshape(
             fmap.shape
         )
 
diff --git a/sdcflows/workflows/apply/correction.py b/sdcflows/workflows/apply/correction.py
index 532f123c32..1c5b53db1f 100644
--- a/sdcflows/workflows/apply/correction.py
+++ b/sdcflows/workflows/apply/correction.py
@@ -72,16 +72,6 @@ def init_unwarp_wf(*, free_mem=None, omp_nthreads=1, debug=False, name="unwarp_w
         the target EPI reference image, after applying unwarping.
     corrected_mask
         a fast mask calculated from the corrected EPI reference.
-    fieldmap_ref
-        the actual B\ :sub:`0` inhomogeneity map (also called *fieldmap*)
-        interpolated from the B-Spline coefficients into the reference EPI's
-        grid, given in Hz units.
-        No motion is taken into account.
-    fieldwarp_ref
-        the displacements field interpolated from the B-Spline coefficients
-        and scaled by the appropriate parameters (readout time of the EPI
-        target and voxel size along PE).
-        No motion is taken into account.
 
     """
     from niworkflows.interfaces.images import RobustAverage
@@ -113,8 +103,6 @@ def init_unwarp_wf(*, free_mem=None, omp_nthreads=1, debug=False, name="unwarp_w
                 "corrected",
                 "corrected_ref",
                 "corrected_mask",
-                "fieldwarp_ref",
-                "fieldmap_ref",
             ]
         ),
         name="outputnode",
@@ -161,8 +149,6 @@ def init_unwarp_wf(*, free_mem=None, omp_nthreads=1, debug=False, name="unwarp_w
         (average, brainextraction_wf, [("out_file", "inputnode.in_file")]),
         (merge, outputnode, [("out_file", "corrected")]),
         (resample, outputnode, [("out_field", "fieldmap")]),
-        # (resample_ref, outputnode, [("out_field", "fieldmap_ref")]),
-        # TODO - take reference from brainextraction workflow
         (brainextraction_wf, outputnode, [
             ("outputnode.out_file", "corrected_ref"),
             ("outputnode.out_mask", "corrected_mask"),