diff --git a/sdcflows/interfaces/bspline.py b/sdcflows/interfaces/bspline.py
index 2b4baf72d6..e8e1cb7598 100644
--- a/sdcflows/interfaces/bspline.py
+++ b/sdcflows/interfaces/bspline.py
@@ -374,16 +374,11 @@ def _run_interface(self, runtime):
         unwarp = None
 
         # Pre-cached interpolator object
-        unwarp = B0FieldTransform(
-            coeffs=[nb.load(cname) for cname in self.inputs.in_coeff],
-            num_threads=(
-                None if not isdefined(self.inputs.num_threads) else self.inputs.num_threads
-            ),
-        )
+        unwarp = B0FieldTransform(coeffs=[nb.load(cname) for cname in self.inputs.in_coeff])
 
         # Reconstruct the field from the coefficients, on the target dataset's grid.
         unwarp.fit(
-            self.inputs.data,
+            self.inputs.in_data,
             affine=(
                 None if not isdefined(self.inputs.fmap2data_xfm) else self.inputs.fmap2data_xfm
             ),
@@ -397,10 +392,13 @@ def _run_interface(self, runtime):
         # HMC matrices are only necessary when reslicing the data (i.e., apply())
         hmc_mats = None
         self._results["out_corrected"] = unwarp.apply(
-            self.inputs.data,
+            self.inputs.in_data,
             pe_dir,
             ro_time,
             xfms=hmc_mats,
+            # num_threads=(
+            #     None if not isdefined(self.inputs.num_threads) else self.inputs.num_threads
+            # ),
         )
         return runtime
 
diff --git a/sdcflows/transform.py b/sdcflows/transform.py
index 5a25614d21..efbee9412d 100644
--- a/sdcflows/transform.py
+++ b/sdcflows/transform.py
@@ -164,7 +164,7 @@ def apply(
         cval: float = 0.0,
         prefilter: bool = True,
         output_dtype: Union[str, np.dtype] = None,
-        num_threads: int = None,
+        # num_threads: int = None,
     ):
         """
         Apply a transformation to an image, resampling on the reference spatial object.
@@ -202,8 +202,6 @@ def apply(
         output_dtype : :obj:`str` or :obj:`~numpy.dtype`
             Override the output data type, instead of propagating it from the
             moving image.
-        num_threads : :obj:`int`
-            Number of CPUs resampling can be parallelized on.
 
         Returns
         -------
@@ -236,13 +234,14 @@ def apply(
             np.meshgrid(*ijk_axis, indexing="ij"), dtype="float32"
         ).reshape(3, -1)
 
-        if xfms is None:
+        if xfms is not None:
+            mov_ras2vox = np.linalg.inv(moving.affine)
             # Map coordinates from reference to time-step
             xfms.reference = moving
             hmc_xyz = xfms.map(xfms.reference.ndcoords.T)
             # Convert from RAS to voxel coordinates
             voxcoords = (
-                np.linalg.inv(xfms.reference.affine)
+                mov_ras2vox
                 @ _as_homogeneous(np.vstack(hmc_xyz), dim=xfms.reference.ndim).T
             )[:3, ...]
 
diff --git a/sdcflows/workflows/apply/correction.py b/sdcflows/workflows/apply/correction.py
index 7891cac5fa..537ac6dd09 100644
--- a/sdcflows/workflows/apply/correction.py
+++ b/sdcflows/workflows/apply/correction.py
@@ -159,10 +159,8 @@ def init_unwarp_wf(*, free_mem=None, omp_nthreads=1, debug=False, name="unwarp_w
         (merge, average, [("out_file", "in_file")]),
         (average, brainextraction_wf, [("out_file", "inputnode.in_file")]),
         (merge, outputnode, [("out_file", "corrected")]),
-        (resample, outputnode, [("out_field", "fieldmap"),
-                                ("out_warp", "fieldwarp")]),
-        (resample_ref, outputnode, [("out_field", "fieldmap_ref"),
-                                    ("out_warp", "fieldwarp_ref")]),
+        (resample, outputnode, [("out_field", "fieldmap")]),
+        (resample_ref, outputnode, [("out_field", "fieldmap_ref")]),
         (brainextraction_wf, outputnode, [
             ("outputnode.out_file", "corrected_ref"),
             ("outputnode.out_mask", "corrected_mask"),
diff --git a/sdcflows/workflows/fit/pepolar.py b/sdcflows/workflows/fit/pepolar.py
index a50b9d1f51..8d79bd7fef 100644
--- a/sdcflows/workflows/fit/pepolar.py
+++ b/sdcflows/workflows/fit/pepolar.py
@@ -229,8 +229,7 @@ def init_topup_wf(
         (split_blips, unwarp, [("out_files", "in_data")]),
         (sort_pe_blips, unwarp, [("readout_times", "ro_time"),
                                  ("pe_dirs", "pe_dir")]),
-        (unwarp, outputnode, [("out_warp", "out_warps"),
-                              ("out_field", "fmap")]),
+        (unwarp, outputnode, [("out_field", "fmap")]),
         (unwarp, concat_corrected, [("out_corrected", "in_files")]),
         (concat_corrected, ref_average, [("out_file", "in_file")]),
     ])
diff --git a/sdcflows/workflows/fit/syn.py b/sdcflows/workflows/fit/syn.py
index b463bcc229..f595a8ca8a 100644
--- a/sdcflows/workflows/fit/syn.py
+++ b/sdcflows/workflows/fit/syn.py
@@ -328,8 +328,7 @@ def init_syn_sdc_wf(
         (zooms_bmask, outputnode, [("output_image", "fmap_mask")]),
         (bs_filter, outputnode, [("out_coeff", "fmap_coeff")]),
         (unwarp, outputnode, [("out_corrected", "fmap_ref"),
-                              ("out_field", "fmap"),
-                              ("out_warp", "out_warp")]),
+                              ("out_field", "fmap")]),
     ])
     # fmt: on