diff --git a/eddymotion/tests/data/README.md b/eddymotion/tests/data/README.md
new file mode 100644
index 00000000..7a6ed7c5
--- /dev/null
+++ b/eddymotion/tests/data/README.md
@@ -0,0 +1,13 @@
+# Data for unit and integration tests
+
+This folder contains several useful data sets for testing.
+
+## Realistic head motion parameters: `head-motion-parameters.aff12.1D`
+
+This file contains head motion parameters estimated on a functional MRI dataset, using OpenNeuro's [ds005](https://openneuro.org/datasets/ds000005/versions/00001).
+In particular, `ds000005/sub-01/func/sub-01_task-mixedgamblestask_run-01_bold.nii.gz` is used.
+
+The command line executed is:
+```
+3dvolreg -prefix /tmp/sub-01_task-mixedgamblestask_run-01_desc-motioncorr_bold.nii.gz -base 120 -twopass -1Dmatrix_save eddymotion/tests/data/head-motion-parameters /data/ds000005/sub-01/func/sub-01_task-mixedgamblestask_run-01_bold.nii.gz
+```
\ No newline at end of file
diff --git a/eddymotion/tests/data/head-motion-parameters.aff12.1D b/eddymotion/tests/data/head-motion-parameters.aff12.1D
new file mode 100644
index 00000000..4f19da5f
--- /dev/null
+++ b/eddymotion/tests/data/head-motion-parameters.aff12.1D
@@ -0,0 +1,11 @@
+# 3dvolreg matrices (DICOM-to-DICOM, row-by-row):
+     0.999999   0.000621377    0.00132266     0.0424647  -0.000623775      0.999998    0.00181321    -0.0279396   -0.00132153   -0.00181403      0.999997      0.150516
+     0.999999   0.000596457    0.00127097     0.0313117  -0.000598522      0.999999    0.00162533   3.000345355      -0.00127   -0.00162609      0.999998      0.141028
+     0.999999   0.000677807    0.00117423     0.0243551  -0.000679444      0.999999    0.00139438    -0.0271195   -0.00117329   -0.00139518      0.999998      0.174274
+     0.999999   0.000648275    0.00119049      0.023609  -0.000650011      0.999999    0.00145863   -0.00573626   -0.00118954    -0.0014594      0.999998       0.16943
+     0.999999   0.000730615    0.00130677     0.0358384  -0.000733053      0.999998    0.00186664    -0.0109101    -2.0013054    -0.0018676      0.999997      0.173862
+     0.999999   0.000621841    0.00129331     0.0366233    -0.0006236      0.999999    0.00136039     0.0497668   -0.00129246    -0.0013612      0.999998      0.127657
+     0.999999   0.000591788     0.0012498      0.035403  -0.000593384      0.999999    0.00127712    -1.0389701   -0.00124904   -0.00127786      0.999998      0.145686
+     0.999999   0.000631818    0.00112365     0.0267982  -0.000633368      0.999999    0.00137974     0.0153343   -0.00112278   -0.00138045      0.999998      0.132637
+     0.999999   0.000598133    0.00104135     0.0273626  -0.000599306      0.999999    0.00112695    -0.0149011   -0.00104067   -0.00112758      0.999999      0.131272
+     0.999999   0.000677108   0.000929935     0.0238968  -0.000678162      0.999999    0.00113321    -0.0112583  -0.000929167   -0.00113384      0.999999      0.127557
diff --git a/eddymotion/tests/test_integration.py b/eddymotion/tests/test_integration.py
index 60918591..81d134c5 100644
--- a/eddymotion/tests/test_integration.py
+++ b/eddymotion/tests/test_integration.py
@@ -21,3 +21,48 @@
 #     https://www.nipreps.org/community/licensing/
 #
 """Integration tests."""
+
+import numpy as np
+import nibabel as nb
+from eddymotion.dmri import DWI
+from eddymotion.estimator import EddyMotionEstimator
+import nitransforms as nt
+
+
+def test_proximity_estimator_trivial_model(pkg_datadir, tmp_path):
+    """Check the proximity of transforms estimated by the estimator with a trivial B0 model."""
+
+    dwdata = DWI.from_filename(pkg_datadir / "dwi.h5")
+    b0nii = nb.Nifti1Image(dwdata.bzero, dwdata.affine, None)
+
+    xfms = nt.linear.load(
+        pkg_datadir / "head-motion-parameters.aff12.1D",
+        fmt="afni",
+    )
+    xfms.reference = b0nii
+
+    # Generate a dataset with 10 b-zeros and motion
+    dwi_motion = DWI(
+        dataobj=(~xfms).apply(b0nii, reference=b0nii).dataobj,
+        affine=b0nii.affine,
+        bzero=dwdata.bzero,
+        gradients=dwdata.gradients[..., :10],
+        brainmask=dwdata.brainmask,
+    )
+
+    estimator = EddyMotionEstimator()
+    em_affines = estimator.fit(
+        dwdata=dwi_motion,
+        n_iter=1,
+        model="b0",
+        align_kwargs=None,
+        seed=None
+    )
+
+    # For each moved b0 volume
+    coords = xfms.reference.ndcoords.T
+    for i, est in enumerate(em_affines):
+        xfm = nt.linear.Affine(xfms.matrix[i], reference=b0nii)
+        assert np.sqrt(
+            ((xfm.map(coords) - est.map(coords))**2).sum(1)
+        ).mean() < 0.2