From 64a814c27143c78d34444a38525144ecda4378fd Mon Sep 17 00:00:00 2001
From: Sebastien Tourbier <sebastien.tourbier1@gmail.com>
Date: Thu, 11 Nov 2021 14:37:40 +0100
Subject: [PATCH 01/11] ENH: Add notebook to generate data to test integration
 of EddyMotionEstimator

---
 ...t integration of EddyMotionEstimator.ipynb | 585 ++++++++++++++++++
 1 file changed, 585 insertions(+)
 create mode 100644 docs/notebooks/Generate data used to test integration of EddyMotionEstimator.ipynb

diff --git a/docs/notebooks/Generate data used to test integration of EddyMotionEstimator.ipynb b/docs/notebooks/Generate data used to test integration of EddyMotionEstimator.ipynb
new file mode 100644
index 00000000..45cb5b81
--- /dev/null
+++ b/docs/notebooks/Generate data used to test integration of EddyMotionEstimator.ipynb	
@@ -0,0 +1,585 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "bacterial-mileage",
+   "metadata": {},
+   "source": [
+    "# Generate data used to test integration of EddyMotionEstimator"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "artificial-albuquerque",
+   "metadata": {},
+   "source": [
+    "This notebook provides all code for the generation of files used in `test_integration.py` to test the `EddyMotionEstimator` using the trivial b0 model. In brief, it consists of:\n",
+    "\n",
+    "1. Estimation of motion parameters from a fMRI recording obtained on OpenNeuro\n",
+    "\n",
+    "2. Extraction of a diffusion-free b0 volume from a diffusion MRI scan (`b0.nii.gz`) and creation of the gradient table in RAS+B format (`gradients.tsv`)\n",
+    "\n",
+    "3. Creation of (1) a 4D Nifti image (`b0.moving.nii.gz`) that contains a serie of 20 diffusion-free b0 volumes transformed using the motion parameters estimated on the last 20 frames of the fMRI recording, and (2) a new gradient table `gradients.moving.tsv`. Each of the 20 transforms are saved in ITK `.tfm` format.\n",
+    "\n",
+    "\n",
+    "## Notes\n",
+    "* The fMRI recording (`sub-01_task-mixedgamblestask_run-01_bold.nii.gz`) is not provided directy in `eddymotion/tests/data` but can be easily get from https://openneuro.org/datasets/ds000005/versions/00001. Please download it before going further in this notebook. but should be s is not included\n",
+    "* The notebook also assumes that a (`sub-01_dwi.nii.gz`, `sub-01_dwi.bval`, `sub-01_dwi.bvec`) file set is included inside `eddymotion/tests/data`."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "small-terrorist",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import os\n",
+    "from glob import glob\n",
+    "from pathlib import Path\n",
+    "import subprocess\n",
+    "import nibabel\n",
+    "import numpy as np\n",
+    "from nibabel.eulerangles import euler2mat\n",
+    "from nibabel.affines import from_matvec\n",
+    "import nitransforms.linear as nitl"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "id": "decreased-information",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "notebookdir = Path(os.path.dirname(os.path.realpath(\"__file__\")))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "determined-replica",
+   "metadata": {},
+   "source": [
+    "## Estimate motion parameters from a BOLD fMRI using AFNI `3dvolreg`"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "id": "specialized-baker",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "/Users/sebastientourbier/Softwares/forks/EddyMotionCorrection/eddymotion/tests/data/sub-01_task-mixedgamblestask_run-01_bold.nii.gz\n"
+     ]
+    }
+   ],
+   "source": [
+    "datadir = Path((notebookdir.parents[1] / 'eddymotion' / 'tests' / 'data'))\n",
+    "bold_file = datadir / 'sub-01_task-mixedgamblestask_run-01_bold.nii.gz'\n",
+    "print(bold_file)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "id": "rural-tension",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def getNumberOfTimeFrames(bold_file):\n",
+    "    \"\"\"Get the number of time frames in the fMRI recording.\"\"\"\n",
+    "    _img = nibabel.load(str(bold_file))\n",
+    "    return _img.get_fdata().shape[3]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "id": "controlled-booth",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def afniEstimateMotionParameters(bold_file, out_file, n):\n",
+    "    \"\"\"Call AFNI 3dvolreg to register each 3D sub-brick from the input fMRI to the n-th base brick.\"\"\"\n",
+    "    \n",
+    "    # Get output filename without .nii.gz extension and make path for \n",
+    "    # output motion parameters file.\n",
+    "    # (Path .stem was only able to remove the first extension (.gz))\n",
+    "    out_filename = os.path.basename(str(out_file))\n",
+    "    index_of_dot = out_filename.index('.')\n",
+    "    out_params_filename = f'{out_filename[:index_of_dot]}.motionparams.1d'\n",
+    "    out_params_file = os.path.join(out_file.parent, out_params_filename)\n",
+    "    \n",
+    "    cmd = ['3dvolreg']\n",
+    "    cmd.append('-prefix')\n",
+    "    cmd.append(f'{str(out_file)}')\n",
+    "    cmd.append('-base')\n",
+    "    cmd.append(str(n))\n",
+    "    cmd.append('-twopass')  # Make the alignment more precise\n",
+    "    cmd.append('-dfile')\n",
+    "    cmd.append(f'{out_params_file}')\n",
+    "    cmd.append(str(bold_file))\n",
+    "\n",
+    "    print(f'Command to register each 3D sub-brick from the input fMRI to the n-th base brick.:\\n {\" \".join(cmd)}')\n",
+    "    try:\n",
+    "        p = subprocess.Popen(cmd, stdout=subprocess.PIPE, stderr=subprocess.PIPE)\n",
+    "        output, err = p.communicate()\n",
+    "        rc = p.returncode\n",
+    "    except Exception as e:\n",
+    "        print(e)\n",
+    "        return 1\n",
+    "    return rc, output, err, out_params_file"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "id": "young-margin",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Command to register each 3D sub-brick from the input fMRI to the n-th base brick.:\n",
+      " 3dvolreg -prefix /Users/sebastientourbier/Softwares/forks/EddyMotionCorrection/eddymotion/tests/data/sub-01_task-mixedgamblestask_run-01_desc-motioncorr_bold.nii.gz -base 120 -twopass -dfile /Users/sebastientourbier/Softwares/forks/EddyMotionCorrection/eddymotion/tests/data/sub-01_task-mixedgamblestask_run-01_desc-motioncorr_bold.motionparams.1d /Users/sebastientourbier/Softwares/forks/EddyMotionCorrection/eddymotion/tests/data/sub-01_task-mixedgamblestask_run-01_bold.nii.gz\n"
+     ]
+    }
+   ],
+   "source": [
+    "number_timepoints = getNumberOfTimeFrames(bold_file)\n",
+    "n = int(0.5 * number_timepoints) # Take the middle frame as reference\n",
+    "out_prefix = datadir / 'sub-01_task-mixedgamblestask_run-01_desc-motioncorr_bold.nii.gz'\n",
+    "rc, output, err, afni_params_file = afniEstimateMotionParameters(bold_file, out_prefix, n)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "id": "burning-netherlands",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Output:\n",
+      "++ 3dvolreg: AFNI version=AFNI_21.3.05 (Nov  7 2021) [64-bit]\n",
+      "++ Authored by: RW Cox\n",
+      "++ Coarse del was 10, replaced with 3\n",
+      "++ Max displacement in automask = 0.63 (mm) at sub-brick 0\n",
+      "++ Max delta displ  in automask = 0.15 (mm) at sub-brick 215\n",
+      "\n"
+     ]
+    }
+   ],
+   "source": [
+    "print(f'Output:\\n{err.decode(\"utf-8\")}')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "perceived-battle",
+   "metadata": {},
+   "source": [
+    "## Extract B0 from diffusion and create RAS+B gradient table"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "id": "established-realtor",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "dwi_file = datadir / 'sub-01_dwi.nii.gz'\n",
+    "fbval = datadir / 'sub-01_dwi.bval'\n",
+    "fbvec = datadir / 'sub-01_dwi.bvec'\n",
+    "\n",
+    "out_rasbn_file = datadir / 'gradients.tsv'\n",
+    "out_b0_file = datadir / 'b0.nii.gz'"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "id": "provincial-eating",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def extract_diffusion_b0(dwi_file, out_b0_file):\n",
+    "    _img = nibabel.load(str(dwi_file))\n",
+    "    b0_data = _img.get_fdata()[..., 0]\n",
+    "    _b0_img = nibabel.Nifti1Image(b0_data, affine=_img.affine)\n",
+    "    print(_b0_img)\n",
+    "    nibabel.save(_b0_img, out_b0_file)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 10,
+   "id": "young-allowance",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "<class 'nibabel.nifti1.Nifti1Image'>\n",
+      "data shape (128, 128, 66)\n",
+      "affine: \n",
+      "[[ -2.           0.          -0.         126.84400177]\n",
+      " [ -0.           2.          -0.         -94.82400513]\n",
+      " [  0.           0.           2.         -26.26519966]\n",
+      " [  0.           0.           0.           1.        ]]\n",
+      "metadata:\n",
+      "<class 'nibabel.nifti1.Nifti1Header'> object, endian='<'\n",
+      "sizeof_hdr      : 348\n",
+      "data_type       : b''\n",
+      "db_name         : b''\n",
+      "extents         : 0\n",
+      "session_error   : 0\n",
+      "regular         : b''\n",
+      "dim_info        : 0\n",
+      "dim             : [  3 128 128  66   1   1   1   1]\n",
+      "intent_p1       : 0.0\n",
+      "intent_p2       : 0.0\n",
+      "intent_p3       : 0.0\n",
+      "intent_code     : none\n",
+      "datatype        : float64\n",
+      "bitpix          : 64\n",
+      "slice_start     : 0\n",
+      "pixdim          : [-1.  2.  2.  2.  1.  1.  1.  1.]\n",
+      "vox_offset      : 0.0\n",
+      "scl_slope       : nan\n",
+      "scl_inter       : nan\n",
+      "slice_end       : 0\n",
+      "slice_code      : unknown\n",
+      "xyzt_units      : 0\n",
+      "cal_max         : 0.0\n",
+      "cal_min         : 0.0\n",
+      "slice_duration  : 0.0\n",
+      "toffset         : 0.0\n",
+      "glmax           : 0\n",
+      "glmin           : 0\n",
+      "descrip         : b''\n",
+      "aux_file        : b''\n",
+      "qform_code      : unknown\n",
+      "sform_code      : aligned\n",
+      "quatern_b       : 0.0\n",
+      "quatern_c       : 1.0\n",
+      "quatern_d       : 0.0\n",
+      "qoffset_x       : 126.844\n",
+      "qoffset_y       : -94.824005\n",
+      "qoffset_z       : -26.2652\n",
+      "srow_x          : [ -2.      0.     -0.    126.844]\n",
+      "srow_y          : [ -0.         2.        -0.       -94.824005]\n",
+      "srow_z          : [  0.       0.       2.     -26.2652]\n",
+      "intent_name     : b''\n",
+      "magic           : b'n+1'\n"
+     ]
+    }
+   ],
+   "source": [
+    "extract_diffusion_b0(dwi_file, out_b0_file)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 11,
+   "id": "neutral-melissa",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def fslgrad2rasb(dwi_file, fbval, fbvec, out_rasbn_file):\n",
+    "    \"\"\"Save gradient table in RAS+B format taking as input the DWI with FSL `.bval` and `.bvec`.\"\"\"\n",
+    "    import numpy as np\n",
+    "    from nibabel import load\n",
+    "    from dipy.io import read_bvals_bvecs\n",
+    "    \n",
+    "    # Read / Load\n",
+    "    img = load(str(dwi_file))\n",
+    "    bvals, bvecs = read_bvals_bvecs(str(fbval), str(fbvec))\n",
+    "    \n",
+    "    # Apply the affine transform to bvecs\n",
+    "    bvecs_tr = np.matmul(img.affine[:3,:3], bvecs.T).T\n",
+    "    \n",
+    "    # Normalize the bvecs\n",
+    "    norm = np.sum(bvecs_tr**2, axis=1)\n",
+    "    bvecs_tr_norm = np.zeros_like(bvecs_tr)\n",
+    "    for i in range(bvecs_tr.shape[0]):\n",
+    "        bvecs_tr_norm[i, :] = bvecs_tr[i, :] / norm[i] \n",
+    "    # Handles bzeros\n",
+    "    bvecs_tr_norm = np.nan_to_num(bvecs_tr_norm)\n",
+    "    \n",
+    "    rasbn = np.c_[bvecs_tr_norm, bvals]\n",
+    "    \n",
+    "    # Save Nx4 numpy matrix in TSV text file\n",
+    "    np.savetxt(fname=str(out_rasbn_file),\n",
+    "               delimiter=\"\\t\",\n",
+    "               X=rasbn)\n",
+    "    \n",
+    "    return rasbn\n",
+    "    "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 12,
+   "id": "daily-liverpool",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "<ipython-input-11-cb78b93f219c>:18: RuntimeWarning: invalid value encountered in true_divide\n",
+      "  bvecs_tr_norm[i, :] = bvecs_tr[i, :] / norm[i]\n"
+     ]
+    }
+   ],
+   "source": [
+    "rasbn = fslgrad2rasb(dwi_file, fbval, fbvec, out_rasbn_file)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 13,
+   "id": "secure-showcase",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def create_moving_b0s_rasbn_gradient(rasbn_file, start_index, nb_of_gradient):\n",
+    "    rasbn = np.genfromtxt(fname=out_rasbn_file, delimiter=\"\\t\", skip_header=0)\n",
+    "    moving_b0s_rasbn = np.zeros((nb_of_gradient,4))\n",
+    "    moving_b0s_rasbn[0:nb_of_gradient, ...] = rasbn[start_index:start_index+nb_of_gradient, ...]\n",
+    "    # Save Nx4 numpy matrix in TSV text file\n",
+    "    np.savetxt(fname=str(rasbn_file.parent / 'gradients.moving.tsv'),\n",
+    "               delimiter=\"\\t\",\n",
+    "               X=moving_b0s_rasbn)\n",
+    "    return moving_b0s_rasbn"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 14,
+   "id": "north-stone",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "[[ 1.80353913e-01  4.66122441e-01  0.00000000e+00  1.00000000e+03]\n",
+      " [-1.27439466e-01  2.82365876e-01  3.92313651e-01  1.00000000e+03]\n",
+      " [ 4.30191553e-01 -2.31833775e-01  1.04924769e-01  1.00000000e+03]\n",
+      " [-1.53638121e-01 -3.83344627e-01  2.82253746e-01  1.00000000e+03]\n",
+      " [-3.67928622e-01  6.49873924e-03  3.38434344e-01  1.00000000e+03]\n",
+      " [ 2.66066783e-01  1.71543057e-01  3.87097161e-01  1.00000000e+03]\n",
+      " [ 8.84487882e-02  4.82220794e-01  9.74435802e-02  1.00000000e+03]\n",
+      " [ 3.85795133e-01  8.15623952e-02  3.07735418e-01  1.00000000e+03]\n",
+      " [ 3.95176644e-02 -4.98222706e-01 -1.80080496e-02  1.00000000e+03]\n",
+      " [ 5.45186999e-02 -4.60157834e-01 -1.88064506e-01  1.00000000e+03]\n",
+      " [ 1.51083700e-01 -3.89715903e-01 -2.74652157e-01  1.00000000e+03]\n",
+      " [-2.32012761e-01 -2.30012651e-01  3.78520819e-01  1.00000000e+03]\n",
+      " [-2.32029468e-01  4.19553283e-01  1.42018036e-01  1.00000000e+03]\n",
+      " [-2.64278535e-01  6.99413891e-03 -4.24144567e-01  1.00000000e+03]\n",
+      " [-4.12452980e-01 -2.69969224e-01 -8.34904821e-02  1.00000000e+03]\n",
+      " [ 3.48560301e-01 -1.43524830e-01 -3.28556840e-01  1.00000000e+03]\n",
+      " [ 2.39551504e-01 -1.69036343e-01 -4.05087094e-01  1.00000000e+03]\n",
+      " [-1.06503195e-01 -4.25512765e-01 -2.40007200e-01  1.00000000e+03]\n",
+      " [-3.24121546e-01  3.68138052e-01 -9.75365762e-02  1.00000000e+03]\n",
+      " [ 2.89905781e-01  3.72378977e-01  1.64946392e-01  1.00000000e+03]]\n",
+      "(20, 4)\n"
+     ]
+    }
+   ],
+   "source": [
+    "moving_b0s_rasbn = create_moving_b0s_rasbn_gradient(out_rasbn_file, 6, 20)\n",
+    "print(moving_b0s_rasbn)\n",
+    "print(moving_b0s_rasbn.shape)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "suffering-equation",
+   "metadata": {},
+   "source": [
+    "## Create moving b0 volumes using the motion parameters estimated from fMRI"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 15,
+   "id": "marine-custom",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def create_moving_b0s_and_transforms(b0_file, params_file, indices):\n",
+    "    \"\"\"Move the b0 volume according to motion parameters estimated from a list of time frame indices.\"\"\"\n",
+    "    \n",
+    "    # Get b0 filename without .nii.gz extension and make path from it\n",
+    "    b0_filename = os.path.basename(str(b0_file))\n",
+    "    index_of_dot = b0_filename.index('.')\n",
+    "    \n",
+    "    # Get b0 affine transform\n",
+    "    b0_affine = nibabel.load(str(b0_file)).affine\n",
+    "    \n",
+    "    # Initialize the matrix containing all b0 volumes\n",
+    "    nb_of_moving_b0s = len(indices) + 1\n",
+    "    x_shape, y_shape, z_shape = nibabel.load(str(b0_file)).get_fdata().shape\n",
+    "    moving_b0s_data = np.zeros((x_shape, y_shape, z_shape, nb_of_moving_b0s))\n",
+    "    \n",
+    "    # Add first b0\n",
+    "    moving_b0s_data[..., 0] = nibabel.load(str(b0_file)).get_fdata()\n",
+    "    \n",
+    "    # Read the motion parameter file generated by 3dvolreg\n",
+    "    with open(params_file) as file:\n",
+    "        lines = file.readlines()\n",
+    "        xfm_files = []\n",
+    "        xforms = []\n",
+    "        cnt=1\n",
+    "        for i in indices:\n",
+    "            # Extract the parameters for time frame i\n",
+    "            formatted_line = lines[i].split()\n",
+    "            transform_params = np.array(formatted_line).astype(float)\n",
+    "            print(transform_params)\n",
+    "            \n",
+    "            # Create a nitransforms object\n",
+    "            T = from_matvec(\n",
+    "                euler2mat(\n",
+    "                    x=transform_params[1],\n",
+    "                    y=transform_params[2],\n",
+    "                    z=transform_params[3],\n",
+    "                ),\n",
+    "                [\n",
+    "                    transform_params[4],\n",
+    "                    transform_params[5],\n",
+    "                    transform_params[6]\n",
+    "                ]\n",
+    "            )\n",
+    "            xfm = nitl.Affine(T, reference=b0_file)\n",
+    "            \n",
+    "            # Apply the transform and save into a nifti file\n",
+    "            moving_b0_filename = f'{b0_filename[:index_of_dot]}.motion-{i}.nii.gz'\n",
+    "            moving_b0_file = b0_file.parent / moving_b0_filename\n",
+    "            (~xfm).apply(b0_file, reference=b0_file).to_filename(moving_b0_file)\n",
+    "            \n",
+    "            # Save the transform\n",
+    "            xfm_filename = f'{b0_filename[:index_of_dot]}.motion-{i}.tfm'\n",
+    "            xfm_file = b0_file.parent / xfm_filename\n",
+    "            xfm.to_filename(xfm_file, fmt=\"itk\")\n",
+    "            \n",
+    "            moving_b0s_data[..., cnt] = nibabel.load(str(moving_b0_file)).get_fdata().copy()\n",
+    "            os.remove(moving_b0_file)\n",
+    "\n",
+    "            xfm_files.append(xfm_file)\n",
+    "            xforms.append(xfm)\n",
+    "            cnt += 1\n",
+    "            \n",
+    "        moving_b0s_img = nibabel.Nifti1Image(moving_b0s_data, affine=b0_affine)\n",
+    "        moving_b0s_filename = f'{b0_filename[:index_of_dot]}.moving.nii.gz'\n",
+    "        moving_b0s_file = b0_file.parent / moving_b0s_filename\n",
+    "        moving_b0s_img.to_filename(str(moving_b0s_file))\n",
+    "\n",
+    "        return xforms, xfm_files, moving_b0s_file"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 16,
+   "id": "imported-thong",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "[220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237\n",
+      " 238 239]\n",
+      "[ 2.200e+02 -4.410e-02  1.570e-02  4.780e-02  2.840e-02 -6.820e-02\n",
+      "  1.140e-02  8.900e+00  8.682e+00]\n",
+      "[ 2.210e+02 -3.690e-02  7.600e-03  4.750e-02  2.610e-02 -5.710e-02\n",
+      " -2.720e-02  8.782e+00  8.595e+00]\n",
+      "[ 2.220e+02 -4.050e-02  1.660e-02  5.230e-02  1.960e-02 -6.230e-02\n",
+      "  1.140e-02  8.784e+00  8.599e+00]\n",
+      "[ 2.230e+02 -4.140e-02  4.300e-03  3.870e-02  2.070e-02 -5.670e-02\n",
+      " -3.700e-02  8.663e+00  8.481e+00]\n",
+      "[ 2.240e+02 -3.870e-02  2.570e-02  4.570e-02  2.050e-02 -5.860e-02\n",
+      "  2.450e-02  8.830e+00  8.636e+00]\n",
+      "[ 2.250e+02 -4.170e-02  5.700e-03  3.330e-02  3.490e-02 -5.530e-02\n",
+      " -4.410e-02  8.415e+00  8.215e+00]\n",
+      "[ 2.260e+02 -4.650e-02  3.980e-02  3.780e-02  1.450e-02 -5.600e-02\n",
+      "  2.370e-02  8.328e+00  8.131e+00]\n",
+      "[ 2.270e+02 -5.420e-02  5.400e-03  3.400e-02  3.390e-02 -6.010e-02\n",
+      " -5.500e-02  7.585e+00  7.300e+00]\n",
+      "[ 2.280e+02 -5.270e-02  2.790e-02  3.140e-02  2.180e-02 -6.360e-02\n",
+      " -9.700e-03  8.115e+00  7.929e+00]\n",
+      "[ 2.290e+02 -5.300e-02 -2.200e-03  3.760e-02  3.380e-02 -6.450e-02\n",
+      " -6.480e-02  8.393e+00  8.065e+00]\n",
+      "[ 2.300e+02 -4.770e-02  6.700e-03  3.980e-02  1.290e-02 -5.750e-02\n",
+      " -9.100e-03  8.625e+00  8.480e+00]\n",
+      "[ 2.310e+02 -5.050e-02 -9.200e-03  3.880e-02  4.350e-02 -5.690e-02\n",
+      " -5.860e-02  8.495e+00  8.221e+00]\n",
+      "[ 2.320e+02 -4.190e-02  1.740e-02  3.940e-02 -1.340e-02 -5.250e-02\n",
+      " -1.200e-03  7.845e+00  7.712e+00]\n",
+      "[ 2.330e+02 -4.080e-02 -4.500e-03  4.110e-02  3.220e-02 -5.220e-02\n",
+      " -5.640e-02  8.123e+00  7.873e+00]\n",
+      "[ 2.340e+02 -4.400e-02  2.470e-02  3.420e-02 -1.190e-02 -5.110e-02\n",
+      "  1.580e-02  8.452e+00  8.311e+00]\n",
+      "[ 2.350e+02 -4.720e-02  7.000e-03  3.410e-02  2.910e-02 -5.230e-02\n",
+      " -5.420e-02  8.309e+00  8.080e+00]\n",
+      "[ 2.360e+02 -4.680e-02  1.790e-02  3.190e-02  5.700e-03 -5.650e-02\n",
+      "  1.430e-02  8.100e+00  7.946e+00]\n",
+      "[ 2.370e+02 -5.500e-02  7.400e-03  3.690e-02  3.770e-02 -5.840e-02\n",
+      " -4.430e-02  8.286e+00  8.050e+00]\n",
+      "[ 2.380e+02 -4.840e-02  1.000e-02  3.620e-02  1.900e-02 -5.940e-02\n",
+      " -3.780e-02  8.568e+00  8.395e+00]\n",
+      "[ 2.390e+02 -5.350e-02  9.300e-03  4.950e-02  4.350e-02 -6.640e-02\n",
+      " -5.300e-02  8.919e+00  8.632e+00]\n"
+     ]
+    }
+   ],
+   "source": [
+    "indices=np.arange(220,240)\n",
+    "print(indices)\n",
+    "xforms, xfm_files, moving_b0s_file = create_moving_b0s_and_transforms(\n",
+    "    out_b0_file, afni_params_file, indices\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "advanced-belgium",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.8.3"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

From ad02acad1bc76b39192bba6c0a133962e9809faf Mon Sep 17 00:00:00 2001
From: Sebastien Tourbier <sebastien.tourbier1@gmail.com>
Date: Thu, 11 Nov 2021 14:44:33 +0100
Subject: [PATCH 02/11] ENH: Add `test_proximity_estimator_trivial_model` to
 `test_integration.py`

---
 eddymotion/tests/test_integration.py | 48 ++++++++++++++++++++++++++++
 1 file changed, 48 insertions(+)

diff --git a/eddymotion/tests/test_integration.py b/eddymotion/tests/test_integration.py
index 60918591..d6a097be 100644
--- a/eddymotion/tests/test_integration.py
+++ b/eddymotion/tests/test_integration.py
@@ -21,3 +21,51 @@
 #     https://www.nipreps.org/community/licensing/
 #
 """Integration tests."""
+
+import pytest  # noqa
+import numpy as np
+import nibabel as nb
+from eddymotion.dmri import DWI
+from eddymotion.estimator import EddyMotionEstimator
+import nitransforms as nit
+
+
+def test_proximity_estimator_trivial_model(pkg_datadir):
+    """Check the proximity of the transforms estimated by :obj:`~eddymotion.estimator.EddyMotionEstimator` with a trivial B0 model."""
+    _img = nb.load((pkg_datadir / 'b0.moving.nii.gz'))
+    _moving_b0s_data = _img.get_fdata()[..., 1:]
+    _b0 = _img.get_fdata()[..., 0]
+    _affine = _img.affine.copy()
+    _gradients = np.genfromtxt(
+        fname=str(pkg_datadir / 'gradients.moving.tsv'),
+        delimiter="\t",
+        skip_header=0
+    ).T
+    _DWI = DWI(
+            dataobj=_moving_b0s_data,
+            affine=_affine,
+            bzero=_b0,
+            gradients=_gradients,
+    )
+
+    estimator = EddyMotionEstimator()
+    em_affines = estimator.fit(
+        dwdata=_DWI,
+        n_iter=1,
+        model="b0",
+        align_kwargs=None,
+        seed=None
+    )
+
+    for i, xfm in enumerate(em_affines):
+        fixed_b0_img = nb.Nifti1Image(_b0, affine=_affine)
+        moving_b0_img = nb.Nifti1Image(_moving_b0s_data[..., i], affine=_affine)
+        xfm2 = nit.linear.Affine(
+            nit.io.itk.ITKLinearTransform.from_filename(
+                str(pkg_datadir / f'b0.motion-{i + 220}.tfm')
+            ).to_ras(reference=fixed_b0_img, moving=moving_b0_img),
+            reference=(pkg_datadir / 'b0.nii.gz')
+        )
+        assert np.all(
+            abs(xfm.map(xfm.reference.ndcoords.T) - xfm2.map(xfm.reference.ndcoords.T)) < 0.3
+        )

From 910f2f0a173961b9160873c82d17dd7653515b22 Mon Sep 17 00:00:00 2001
From: Sebastien Tourbier <sebastien.tourbier1@gmail.com>
Date: Thu, 11 Nov 2021 14:51:33 +0100
Subject: [PATCH 03/11] ENH: Remove use of 'b0.nii.gz' in
 `test_proximity_estimator_trivial_model`

---
 eddymotion/tests/test_integration.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/eddymotion/tests/test_integration.py b/eddymotion/tests/test_integration.py
index d6a097be..fdf6cba9 100644
--- a/eddymotion/tests/test_integration.py
+++ b/eddymotion/tests/test_integration.py
@@ -64,7 +64,7 @@ def test_proximity_estimator_trivial_model(pkg_datadir):
             nit.io.itk.ITKLinearTransform.from_filename(
                 str(pkg_datadir / f'b0.motion-{i + 220}.tfm')
             ).to_ras(reference=fixed_b0_img, moving=moving_b0_img),
-            reference=(pkg_datadir / 'b0.nii.gz')
+            reference=fixed_b0_img
         )
         assert np.all(
             abs(xfm.map(xfm.reference.ndcoords.T) - xfm2.map(xfm.reference.ndcoords.T)) < 0.3

From 19febdd0b9513f629c932be3c98ee7a6eab6ce9a Mon Sep 17 00:00:00 2001
From: Sebastien Tourbier <sebastien.tourbier1@gmail.com>
Date: Thu, 11 Nov 2021 14:59:27 +0100
Subject: [PATCH 04/11] MISC: Increase upper bound of
 `test_proximity_estimator_trivial_model`

---
 eddymotion/tests/test_integration.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/eddymotion/tests/test_integration.py b/eddymotion/tests/test_integration.py
index fdf6cba9..5ea64394 100644
--- a/eddymotion/tests/test_integration.py
+++ b/eddymotion/tests/test_integration.py
@@ -67,5 +67,5 @@ def test_proximity_estimator_trivial_model(pkg_datadir):
             reference=fixed_b0_img
         )
         assert np.all(
-            abs(xfm.map(xfm.reference.ndcoords.T) - xfm2.map(xfm.reference.ndcoords.T)) < 0.3
+            abs(xfm.map(xfm.reference.ndcoords.T) - xfm2.map(xfm.reference.ndcoords.T)) < 0.4
         )

From 7e2f3834ade9fc937e9bc41955ecff3988d33663 Mon Sep 17 00:00:00 2001
From: Sebastien Tourbier <sebastien.tourbier1@gmail.com>
Date: Thu, 11 Nov 2021 15:36:46 +0100
Subject: [PATCH 05/11] ENH: Save all transform matrices in one `.npy` file
 that are loaded by `test_proximity_estimator_trivial_model`

---
 ...t integration of EddyMotionEstimator.ipynb | 131 ++++++++++++------
 eddymotion/tests/test_integration.py          |   7 +-
 2 files changed, 91 insertions(+), 47 deletions(-)

diff --git a/docs/notebooks/Generate data used to test integration of EddyMotionEstimator.ipynb b/docs/notebooks/Generate data used to test integration of EddyMotionEstimator.ipynb
index 45cb5b81..71baea84 100644
--- a/docs/notebooks/Generate data used to test integration of EddyMotionEstimator.ipynb	
+++ b/docs/notebooks/Generate data used to test integration of EddyMotionEstimator.ipynb	
@@ -10,16 +10,16 @@
   },
   {
    "cell_type": "markdown",
-   "id": "artificial-albuquerque",
+   "id": "exterior-magazine",
    "metadata": {},
    "source": [
-    "This notebook provides all code for the generation of files used in `test_integration.py` to test the `EddyMotionEstimator` using the trivial b0 model. In brief, it consists of:\n",
+    "This notebook provides all code for the generation of `b0.moving.nii.gz`, `gradients.moving.tsv`, and `b0.moving.transforms.npy`, used in `test_integration.py` to test the `EddyMotionEstimator` using the trivial b0 model. In brief, it consists of:\n",
     "\n",
     "1. Estimation of motion parameters from a fMRI recording obtained on OpenNeuro\n",
     "\n",
     "2. Extraction of a diffusion-free b0 volume from a diffusion MRI scan (`b0.nii.gz`) and creation of the gradient table in RAS+B format (`gradients.tsv`)\n",
     "\n",
-    "3. Creation of (1) a 4D Nifti image (`b0.moving.nii.gz`) that contains a serie of 20 diffusion-free b0 volumes transformed using the motion parameters estimated on the last 20 frames of the fMRI recording, and (2) a new gradient table `gradients.moving.tsv`. Each of the 20 transforms are saved in ITK `.tfm` format.\n",
+    "3. Creation of (1) a 4D Nifti image (`b0.moving.nii.gz`) that contains a serie of 10 diffusion-free b0 volumes transformed using the motion parameters estimated on the last 10 frames of the fMRI recording, and (2) a new gradient table `gradients.moving.tsv`. Each of the 10 transforms are temporarily saved in ITK `.tfm` format before having their matrix extracted, concatenated in a numpy array, saved as `b0.moving.transforms.npy`.\n",
     "\n",
     "\n",
     "## Notes\n",
@@ -48,7 +48,7 @@
   {
    "cell_type": "code",
    "execution_count": 2,
-   "id": "decreased-information",
+   "id": "pointed-apache",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -366,7 +366,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 14,
+   "execution_count": 27,
    "id": "north-stone",
    "metadata": {},
    "outputs": [
@@ -383,23 +383,13 @@
       " [ 8.84487882e-02  4.82220794e-01  9.74435802e-02  1.00000000e+03]\n",
       " [ 3.85795133e-01  8.15623952e-02  3.07735418e-01  1.00000000e+03]\n",
       " [ 3.95176644e-02 -4.98222706e-01 -1.80080496e-02  1.00000000e+03]\n",
-      " [ 5.45186999e-02 -4.60157834e-01 -1.88064506e-01  1.00000000e+03]\n",
-      " [ 1.51083700e-01 -3.89715903e-01 -2.74652157e-01  1.00000000e+03]\n",
-      " [-2.32012761e-01 -2.30012651e-01  3.78520819e-01  1.00000000e+03]\n",
-      " [-2.32029468e-01  4.19553283e-01  1.42018036e-01  1.00000000e+03]\n",
-      " [-2.64278535e-01  6.99413891e-03 -4.24144567e-01  1.00000000e+03]\n",
-      " [-4.12452980e-01 -2.69969224e-01 -8.34904821e-02  1.00000000e+03]\n",
-      " [ 3.48560301e-01 -1.43524830e-01 -3.28556840e-01  1.00000000e+03]\n",
-      " [ 2.39551504e-01 -1.69036343e-01 -4.05087094e-01  1.00000000e+03]\n",
-      " [-1.06503195e-01 -4.25512765e-01 -2.40007200e-01  1.00000000e+03]\n",
-      " [-3.24121546e-01  3.68138052e-01 -9.75365762e-02  1.00000000e+03]\n",
-      " [ 2.89905781e-01  3.72378977e-01  1.64946392e-01  1.00000000e+03]]\n",
-      "(20, 4)\n"
+      " [ 5.45186999e-02 -4.60157834e-01 -1.88064506e-01  1.00000000e+03]]\n",
+      "(10, 4)\n"
      ]
     }
    ],
    "source": [
-    "moving_b0s_rasbn = create_moving_b0s_rasbn_gradient(out_rasbn_file, 6, 20)\n",
+    "moving_b0s_rasbn = create_moving_b0s_rasbn_gradient(out_rasbn_file, 6, 10)\n",
     "print(moving_b0s_rasbn)\n",
     "print(moving_b0s_rasbn.shape)"
    ]
@@ -491,7 +481,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 16,
+   "execution_count": 23,
    "id": "imported-thong",
    "metadata": {},
    "outputs": [
@@ -499,28 +489,7 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "[220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237\n",
-      " 238 239]\n",
-      "[ 2.200e+02 -4.410e-02  1.570e-02  4.780e-02  2.840e-02 -6.820e-02\n",
-      "  1.140e-02  8.900e+00  8.682e+00]\n",
-      "[ 2.210e+02 -3.690e-02  7.600e-03  4.750e-02  2.610e-02 -5.710e-02\n",
-      " -2.720e-02  8.782e+00  8.595e+00]\n",
-      "[ 2.220e+02 -4.050e-02  1.660e-02  5.230e-02  1.960e-02 -6.230e-02\n",
-      "  1.140e-02  8.784e+00  8.599e+00]\n",
-      "[ 2.230e+02 -4.140e-02  4.300e-03  3.870e-02  2.070e-02 -5.670e-02\n",
-      " -3.700e-02  8.663e+00  8.481e+00]\n",
-      "[ 2.240e+02 -3.870e-02  2.570e-02  4.570e-02  2.050e-02 -5.860e-02\n",
-      "  2.450e-02  8.830e+00  8.636e+00]\n",
-      "[ 2.250e+02 -4.170e-02  5.700e-03  3.330e-02  3.490e-02 -5.530e-02\n",
-      " -4.410e-02  8.415e+00  8.215e+00]\n",
-      "[ 2.260e+02 -4.650e-02  3.980e-02  3.780e-02  1.450e-02 -5.600e-02\n",
-      "  2.370e-02  8.328e+00  8.131e+00]\n",
-      "[ 2.270e+02 -5.420e-02  5.400e-03  3.400e-02  3.390e-02 -6.010e-02\n",
-      " -5.500e-02  7.585e+00  7.300e+00]\n",
-      "[ 2.280e+02 -5.270e-02  2.790e-02  3.140e-02  2.180e-02 -6.360e-02\n",
-      " -9.700e-03  8.115e+00  7.929e+00]\n",
-      "[ 2.290e+02 -5.300e-02 -2.200e-03  3.760e-02  3.380e-02 -6.450e-02\n",
-      " -6.480e-02  8.393e+00  8.065e+00]\n",
+      "[230 231 232 233 234 235 236 237 238 239]\n",
       "[ 2.300e+02 -4.770e-02  6.700e-03  3.980e-02  1.290e-02 -5.750e-02\n",
       " -9.100e-03  8.625e+00  8.480e+00]\n",
       "[ 2.310e+02 -5.050e-02 -9.200e-03  3.880e-02  4.350e-02 -5.690e-02\n",
@@ -545,7 +514,7 @@
     }
    ],
    "source": [
-    "indices=np.arange(220,240)\n",
+    "indices=np.arange(230,240)\n",
     "print(indices)\n",
     "xforms, xfm_files, moving_b0s_file = create_moving_b0s_and_transforms(\n",
     "    out_b0_file, afni_params_file, indices\n",
@@ -554,9 +523,85 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 28,
    "id": "advanced-belgium",
    "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "[[ 0.99918598 -0.0397886   0.00669995  0.0129    ]\n",
+      " [ 0.039425    0.99808401  0.0476808  -0.0575    ]\n",
+      " [-0.00858427 -0.0473779   0.99883997 -0.0091    ]\n",
+      " [ 0.          0.          0.          1.        ]]\n",
+      "[[ 0.99920499 -0.0387886  -0.00919987  0.0435    ]\n",
+      " [ 0.0392049   0.99795502  0.0504764  -0.0569    ]\n",
+      " [ 0.00722315 -0.050797    0.99868298 -0.0586    ]\n",
+      " [ 0.          0.          0.          1.        ]]\n",
+      "[[ 0.99907303 -0.0393838   0.0173991  -0.0134    ]\n",
+      " [ 0.038627    0.99837601  0.0418814  -0.0525    ]\n",
+      " [-0.0190203  -0.0411705   0.99897099 -0.0012    ]\n",
+      " [ 0.          0.          0.          1.        ]]\n",
+      "[[ 0.99914497 -0.041088   -0.00449998  0.0322    ]\n",
+      " [ 0.0412376   0.99831599  0.0407883  -0.0522    ]\n",
+      " [ 0.0028165  -0.040939    0.99915802 -0.0564    ]\n",
+      " [ 0.          0.          0.          1.        ]]\n",
+      "[[ 0.99910998 -0.0341829   0.0246975  -0.0119    ]\n",
+      " [ 0.0330745   0.99848503  0.0439724  -0.0511    ]\n",
+      " [-0.0261632  -0.0431164   0.99872702  0.0158    ]\n",
+      " [ 0.          0.          0.          1.        ]]\n",
+      "[[ 0.999394   -0.0340926   0.00699994  0.0291    ]\n",
+      " [ 0.0337253   0.998317    0.0471813  -0.0523    ]\n",
+      " [-0.00859669 -0.0469167   0.99886203 -0.0542    ]\n",
+      " [ 0.          0.          0.          1.        ]]\n",
+      "[[ 0.999331   -0.0318895   0.017899    0.0057    ]\n",
+      " [ 0.0310227   0.99842399  0.0467754  -0.0565    ]\n",
+      " [-0.0193625  -0.0461889   0.99874502  0.0143    ]\n",
+      " [ 0.          0.          0.          1.        ]]\n",
+      "[[ 0.99929202 -0.0368906   0.00739993  0.0377    ]\n",
+      " [ 0.0364293   0.997823    0.0549708  -0.0584    ]\n",
+      " [-0.00941173 -0.0546623   0.99846101 -0.0443    ]\n",
+      " [ 0.          0.          0.          1.        ]]\n",
+      "[[ 0.999295   -0.0361903   0.00999983  0.019     ]\n",
+      " [ 0.0356662   0.99819201  0.0483787  -0.0594    ]\n",
+      " [-0.0117326  -0.0479879   0.998779   -0.0378    ]\n",
+      " [ 0.          0.          0.          1.        ]]\n",
+      "[[ 0.99873197 -0.0494776   0.00929987  0.0435    ]\n",
+      " [ 0.0489123   0.99737102  0.0534722  -0.0664    ]\n",
+      " [-0.0119211  -0.0529495   0.99852598 -0.053     ]\n",
+      " [ 0.          0.          0.          1.        ]]\n"
+     ]
+    }
+   ],
+   "source": [
+    "import nitransforms as nit\n",
+    "\n",
+    "_img = nibabel.load((datadir / 'b0.moving.nii.gz'))\n",
+    "_moving_b0s_data = _img.get_fdata()[..., 1:]\n",
+    "_b0 = _img.get_fdata()[..., 0]\n",
+    "_affine = _img.affine.copy()\n",
+    "\n",
+    "xfms_array = np.zeros((4,4,len(xfm_files)))\n",
+    "\n",
+    "for i, xfm in enumerate(xfm_files):\n",
+    "    fixed_b0_img = nibabel.Nifti1Image(_b0, affine=_affine)\n",
+    "    moving_b0_img = nibabel.Nifti1Image(_moving_b0s_data[..., i], affine=_affine)\n",
+    "    T = nit.io.itk.ITKLinearTransform.from_filename(\n",
+    "        str(datadir / f'b0.motion-{i + 230}.tfm')\n",
+    "    ).to_ras(reference=fixed_b0_img, moving=moving_b0_img)\n",
+    "    xfms_array[..., i] = T\n",
+    "    print(xfms_array[..., i])\n",
+    "\n",
+    "with open((datadir / 'b0.moving.transforms.npy'), 'wb') as f:\n",
+    "    np.save(f, xfms_array)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "stunning-ministry",
+   "metadata": {},
    "outputs": [],
    "source": []
   }
diff --git a/eddymotion/tests/test_integration.py b/eddymotion/tests/test_integration.py
index 5ea64394..60be1760 100644
--- a/eddymotion/tests/test_integration.py
+++ b/eddymotion/tests/test_integration.py
@@ -57,13 +57,12 @@ def test_proximity_estimator_trivial_model(pkg_datadir):
         seed=None
     )
 
+    ref_xfms = np.load((pkg_datadir / "b0.moving.transforms.npy"))
+
     for i, xfm in enumerate(em_affines):
         fixed_b0_img = nb.Nifti1Image(_b0, affine=_affine)
-        moving_b0_img = nb.Nifti1Image(_moving_b0s_data[..., i], affine=_affine)
         xfm2 = nit.linear.Affine(
-            nit.io.itk.ITKLinearTransform.from_filename(
-                str(pkg_datadir / f'b0.motion-{i + 220}.tfm')
-            ).to_ras(reference=fixed_b0_img, moving=moving_b0_img),
+            ref_xfms[..., i],
             reference=fixed_b0_img
         )
         assert np.all(

From f5124fb216bcade7bfe541d01de5f497403b51b6 Mon Sep 17 00:00:00 2001
From: Sebastien Tourbier <sebastien.tourbier1@gmail.com>
Date: Thu, 11 Nov 2021 15:45:04 +0100
Subject: [PATCH 06/11] DOC: Add comments in `test_integration.py`

---
 eddymotion/tests/test_integration.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/eddymotion/tests/test_integration.py b/eddymotion/tests/test_integration.py
index 60be1760..4df97f5a 100644
--- a/eddymotion/tests/test_integration.py
+++ b/eddymotion/tests/test_integration.py
@@ -57,8 +57,10 @@ def test_proximity_estimator_trivial_model(pkg_datadir):
         seed=None
     )
 
+    # Load reference transforms originally applied
     ref_xfms = np.load((pkg_datadir / "b0.moving.transforms.npy"))
 
+    # For each moved b0 volume
     for i, xfm in enumerate(em_affines):
         fixed_b0_img = nb.Nifti1Image(_b0, affine=_affine)
         xfm2 = nit.linear.Affine(

From 44dfb5a95dcd4513b51d06ef617dabb8416f1fbe Mon Sep 17 00:00:00 2001
From: Sebastien Tourbier <sebastien.tourbier1@gmail.com>
Date: Thu, 11 Nov 2021 16:03:55 +0100
Subject: [PATCH 07/11] STY: Remove extra indent in `test_integration.py`

---
 eddymotion/tests/test_integration.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/eddymotion/tests/test_integration.py b/eddymotion/tests/test_integration.py
index 4df97f5a..0dd7d562 100644
--- a/eddymotion/tests/test_integration.py
+++ b/eddymotion/tests/test_integration.py
@@ -42,10 +42,10 @@ def test_proximity_estimator_trivial_model(pkg_datadir):
         skip_header=0
     ).T
     _DWI = DWI(
-            dataobj=_moving_b0s_data,
-            affine=_affine,
-            bzero=_b0,
-            gradients=_gradients,
+        dataobj=_moving_b0s_data,
+        affine=_affine,
+        bzero=_b0,
+        gradients=_gradients,
     )
 
     estimator = EddyMotionEstimator()

From 6c9240d4b3b1eb8f10620a72816151dc807d814a Mon Sep 17 00:00:00 2001
From: Sebastien Tourbier <sebastien.tourbier1@gmail.com>
Date: Thu, 11 Nov 2021 16:28:10 +0100
Subject: [PATCH 08/11] MISC: Remove unused `import pytest`

---
 eddymotion/tests/test_integration.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/eddymotion/tests/test_integration.py b/eddymotion/tests/test_integration.py
index 0dd7d562..6a8efacf 100644
--- a/eddymotion/tests/test_integration.py
+++ b/eddymotion/tests/test_integration.py
@@ -22,7 +22,6 @@
 #
 """Integration tests."""
 
-import pytest  # noqa
 import numpy as np
 import nibabel as nb
 from eddymotion.dmri import DWI

From 2789f0cb87a0e3f22607bfed19f81c7756631aa7 Mon Sep 17 00:00:00 2001
From: Sebastien Tourbier <sebastien.tourbier1@gmail.com>
Date: Thu, 11 Nov 2021 16:29:26 +0100
Subject: [PATCH 09/11] STY: Reduce docstring line length for
 `test_proximity_estimator_trivial_model`

---
 eddymotion/tests/test_integration.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/eddymotion/tests/test_integration.py b/eddymotion/tests/test_integration.py
index 6a8efacf..c168768f 100644
--- a/eddymotion/tests/test_integration.py
+++ b/eddymotion/tests/test_integration.py
@@ -30,7 +30,7 @@
 
 
 def test_proximity_estimator_trivial_model(pkg_datadir):
-    """Check the proximity of the transforms estimated by :obj:`~eddymotion.estimator.EddyMotionEstimator` with a trivial B0 model."""
+    """Check the proximity of transforms estimated by `EddyMotionEstimator` with a trivial B0 model."""
     _img = nb.load((pkg_datadir / 'b0.moving.nii.gz'))
     _moving_b0s_data = _img.get_fdata()[..., 1:]
     _b0 = _img.get_fdata()[..., 0]

From 94120464fc96de42f67b6b24230eeaf8e0f83a4a Mon Sep 17 00:00:00 2001
From: Sebastien Tourbier <sebastien.tourbier1@gmail.com>
Date: Thu, 11 Nov 2021 16:34:16 +0100
Subject: [PATCH 10/11] STY: Reduce docstring line length even more for
 `test_proximity_estimator_trivial_model`

---
 eddymotion/tests/test_integration.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/eddymotion/tests/test_integration.py b/eddymotion/tests/test_integration.py
index c168768f..33243ad9 100644
--- a/eddymotion/tests/test_integration.py
+++ b/eddymotion/tests/test_integration.py
@@ -30,7 +30,7 @@
 
 
 def test_proximity_estimator_trivial_model(pkg_datadir):
-    """Check the proximity of transforms estimated by `EddyMotionEstimator` with a trivial B0 model."""
+    """Check the proximity of transforms estimated by the estimator with a trivial B0 model."""
     _img = nb.load((pkg_datadir / 'b0.moving.nii.gz'))
     _moving_b0s_data = _img.get_fdata()[..., 1:]
     _b0 = _img.get_fdata()[..., 0]

From a9e67cc9503c5bc77da6af75791e39dce6df2645 Mon Sep 17 00:00:00 2001
From: Oscar Esteban <code@oscaresteban.es>
Date: Thu, 7 Apr 2022 20:21:30 +0200
Subject: [PATCH 11/11] enh: minimize test towards the testing goal

---
 ...t integration of EddyMotionEstimator.ipynb | 630 ------------------
 eddymotion/tests/data/README.md               |  13 +
 .../data/head-motion-parameters.aff12.1D      |  11 +
 eddymotion/tests/test_integration.py          |  55 +-
 4 files changed, 50 insertions(+), 659 deletions(-)
 delete mode 100644 docs/notebooks/Generate data used to test integration of EddyMotionEstimator.ipynb
 create mode 100644 eddymotion/tests/data/README.md
 create mode 100644 eddymotion/tests/data/head-motion-parameters.aff12.1D

diff --git a/docs/notebooks/Generate data used to test integration of EddyMotionEstimator.ipynb b/docs/notebooks/Generate data used to test integration of EddyMotionEstimator.ipynb
deleted file mode 100644
index 71baea84..00000000
--- a/docs/notebooks/Generate data used to test integration of EddyMotionEstimator.ipynb	
+++ /dev/null
@@ -1,630 +0,0 @@
-{
- "cells": [
-  {
-   "cell_type": "markdown",
-   "id": "bacterial-mileage",
-   "metadata": {},
-   "source": [
-    "# Generate data used to test integration of EddyMotionEstimator"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "id": "exterior-magazine",
-   "metadata": {},
-   "source": [
-    "This notebook provides all code for the generation of `b0.moving.nii.gz`, `gradients.moving.tsv`, and `b0.moving.transforms.npy`, used in `test_integration.py` to test the `EddyMotionEstimator` using the trivial b0 model. In brief, it consists of:\n",
-    "\n",
-    "1. Estimation of motion parameters from a fMRI recording obtained on OpenNeuro\n",
-    "\n",
-    "2. Extraction of a diffusion-free b0 volume from a diffusion MRI scan (`b0.nii.gz`) and creation of the gradient table in RAS+B format (`gradients.tsv`)\n",
-    "\n",
-    "3. Creation of (1) a 4D Nifti image (`b0.moving.nii.gz`) that contains a serie of 10 diffusion-free b0 volumes transformed using the motion parameters estimated on the last 10 frames of the fMRI recording, and (2) a new gradient table `gradients.moving.tsv`. Each of the 10 transforms are temporarily saved in ITK `.tfm` format before having their matrix extracted, concatenated in a numpy array, saved as `b0.moving.transforms.npy`.\n",
-    "\n",
-    "\n",
-    "## Notes\n",
-    "* The fMRI recording (`sub-01_task-mixedgamblestask_run-01_bold.nii.gz`) is not provided directy in `eddymotion/tests/data` but can be easily get from https://openneuro.org/datasets/ds000005/versions/00001. Please download it before going further in this notebook. but should be s is not included\n",
-    "* The notebook also assumes that a (`sub-01_dwi.nii.gz`, `sub-01_dwi.bval`, `sub-01_dwi.bvec`) file set is included inside `eddymotion/tests/data`."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 1,
-   "id": "small-terrorist",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "import os\n",
-    "from glob import glob\n",
-    "from pathlib import Path\n",
-    "import subprocess\n",
-    "import nibabel\n",
-    "import numpy as np\n",
-    "from nibabel.eulerangles import euler2mat\n",
-    "from nibabel.affines import from_matvec\n",
-    "import nitransforms.linear as nitl"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 2,
-   "id": "pointed-apache",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "notebookdir = Path(os.path.dirname(os.path.realpath(\"__file__\")))"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "id": "determined-replica",
-   "metadata": {},
-   "source": [
-    "## Estimate motion parameters from a BOLD fMRI using AFNI `3dvolreg`"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 3,
-   "id": "specialized-baker",
-   "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "/Users/sebastientourbier/Softwares/forks/EddyMotionCorrection/eddymotion/tests/data/sub-01_task-mixedgamblestask_run-01_bold.nii.gz\n"
-     ]
-    }
-   ],
-   "source": [
-    "datadir = Path((notebookdir.parents[1] / 'eddymotion' / 'tests' / 'data'))\n",
-    "bold_file = datadir / 'sub-01_task-mixedgamblestask_run-01_bold.nii.gz'\n",
-    "print(bold_file)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 4,
-   "id": "rural-tension",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "def getNumberOfTimeFrames(bold_file):\n",
-    "    \"\"\"Get the number of time frames in the fMRI recording.\"\"\"\n",
-    "    _img = nibabel.load(str(bold_file))\n",
-    "    return _img.get_fdata().shape[3]"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 5,
-   "id": "controlled-booth",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "def afniEstimateMotionParameters(bold_file, out_file, n):\n",
-    "    \"\"\"Call AFNI 3dvolreg to register each 3D sub-brick from the input fMRI to the n-th base brick.\"\"\"\n",
-    "    \n",
-    "    # Get output filename without .nii.gz extension and make path for \n",
-    "    # output motion parameters file.\n",
-    "    # (Path .stem was only able to remove the first extension (.gz))\n",
-    "    out_filename = os.path.basename(str(out_file))\n",
-    "    index_of_dot = out_filename.index('.')\n",
-    "    out_params_filename = f'{out_filename[:index_of_dot]}.motionparams.1d'\n",
-    "    out_params_file = os.path.join(out_file.parent, out_params_filename)\n",
-    "    \n",
-    "    cmd = ['3dvolreg']\n",
-    "    cmd.append('-prefix')\n",
-    "    cmd.append(f'{str(out_file)}')\n",
-    "    cmd.append('-base')\n",
-    "    cmd.append(str(n))\n",
-    "    cmd.append('-twopass')  # Make the alignment more precise\n",
-    "    cmd.append('-dfile')\n",
-    "    cmd.append(f'{out_params_file}')\n",
-    "    cmd.append(str(bold_file))\n",
-    "\n",
-    "    print(f'Command to register each 3D sub-brick from the input fMRI to the n-th base brick.:\\n {\" \".join(cmd)}')\n",
-    "    try:\n",
-    "        p = subprocess.Popen(cmd, stdout=subprocess.PIPE, stderr=subprocess.PIPE)\n",
-    "        output, err = p.communicate()\n",
-    "        rc = p.returncode\n",
-    "    except Exception as e:\n",
-    "        print(e)\n",
-    "        return 1\n",
-    "    return rc, output, err, out_params_file"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 6,
-   "id": "young-margin",
-   "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "Command to register each 3D sub-brick from the input fMRI to the n-th base brick.:\n",
-      " 3dvolreg -prefix /Users/sebastientourbier/Softwares/forks/EddyMotionCorrection/eddymotion/tests/data/sub-01_task-mixedgamblestask_run-01_desc-motioncorr_bold.nii.gz -base 120 -twopass -dfile /Users/sebastientourbier/Softwares/forks/EddyMotionCorrection/eddymotion/tests/data/sub-01_task-mixedgamblestask_run-01_desc-motioncorr_bold.motionparams.1d /Users/sebastientourbier/Softwares/forks/EddyMotionCorrection/eddymotion/tests/data/sub-01_task-mixedgamblestask_run-01_bold.nii.gz\n"
-     ]
-    }
-   ],
-   "source": [
-    "number_timepoints = getNumberOfTimeFrames(bold_file)\n",
-    "n = int(0.5 * number_timepoints) # Take the middle frame as reference\n",
-    "out_prefix = datadir / 'sub-01_task-mixedgamblestask_run-01_desc-motioncorr_bold.nii.gz'\n",
-    "rc, output, err, afni_params_file = afniEstimateMotionParameters(bold_file, out_prefix, n)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 7,
-   "id": "burning-netherlands",
-   "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "Output:\n",
-      "++ 3dvolreg: AFNI version=AFNI_21.3.05 (Nov  7 2021) [64-bit]\n",
-      "++ Authored by: RW Cox\n",
-      "++ Coarse del was 10, replaced with 3\n",
-      "++ Max displacement in automask = 0.63 (mm) at sub-brick 0\n",
-      "++ Max delta displ  in automask = 0.15 (mm) at sub-brick 215\n",
-      "\n"
-     ]
-    }
-   ],
-   "source": [
-    "print(f'Output:\\n{err.decode(\"utf-8\")}')"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "id": "perceived-battle",
-   "metadata": {},
-   "source": [
-    "## Extract B0 from diffusion and create RAS+B gradient table"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 8,
-   "id": "established-realtor",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "dwi_file = datadir / 'sub-01_dwi.nii.gz'\n",
-    "fbval = datadir / 'sub-01_dwi.bval'\n",
-    "fbvec = datadir / 'sub-01_dwi.bvec'\n",
-    "\n",
-    "out_rasbn_file = datadir / 'gradients.tsv'\n",
-    "out_b0_file = datadir / 'b0.nii.gz'"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 9,
-   "id": "provincial-eating",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "def extract_diffusion_b0(dwi_file, out_b0_file):\n",
-    "    _img = nibabel.load(str(dwi_file))\n",
-    "    b0_data = _img.get_fdata()[..., 0]\n",
-    "    _b0_img = nibabel.Nifti1Image(b0_data, affine=_img.affine)\n",
-    "    print(_b0_img)\n",
-    "    nibabel.save(_b0_img, out_b0_file)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 10,
-   "id": "young-allowance",
-   "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "<class 'nibabel.nifti1.Nifti1Image'>\n",
-      "data shape (128, 128, 66)\n",
-      "affine: \n",
-      "[[ -2.           0.          -0.         126.84400177]\n",
-      " [ -0.           2.          -0.         -94.82400513]\n",
-      " [  0.           0.           2.         -26.26519966]\n",
-      " [  0.           0.           0.           1.        ]]\n",
-      "metadata:\n",
-      "<class 'nibabel.nifti1.Nifti1Header'> object, endian='<'\n",
-      "sizeof_hdr      : 348\n",
-      "data_type       : b''\n",
-      "db_name         : b''\n",
-      "extents         : 0\n",
-      "session_error   : 0\n",
-      "regular         : b''\n",
-      "dim_info        : 0\n",
-      "dim             : [  3 128 128  66   1   1   1   1]\n",
-      "intent_p1       : 0.0\n",
-      "intent_p2       : 0.0\n",
-      "intent_p3       : 0.0\n",
-      "intent_code     : none\n",
-      "datatype        : float64\n",
-      "bitpix          : 64\n",
-      "slice_start     : 0\n",
-      "pixdim          : [-1.  2.  2.  2.  1.  1.  1.  1.]\n",
-      "vox_offset      : 0.0\n",
-      "scl_slope       : nan\n",
-      "scl_inter       : nan\n",
-      "slice_end       : 0\n",
-      "slice_code      : unknown\n",
-      "xyzt_units      : 0\n",
-      "cal_max         : 0.0\n",
-      "cal_min         : 0.0\n",
-      "slice_duration  : 0.0\n",
-      "toffset         : 0.0\n",
-      "glmax           : 0\n",
-      "glmin           : 0\n",
-      "descrip         : b''\n",
-      "aux_file        : b''\n",
-      "qform_code      : unknown\n",
-      "sform_code      : aligned\n",
-      "quatern_b       : 0.0\n",
-      "quatern_c       : 1.0\n",
-      "quatern_d       : 0.0\n",
-      "qoffset_x       : 126.844\n",
-      "qoffset_y       : -94.824005\n",
-      "qoffset_z       : -26.2652\n",
-      "srow_x          : [ -2.      0.     -0.    126.844]\n",
-      "srow_y          : [ -0.         2.        -0.       -94.824005]\n",
-      "srow_z          : [  0.       0.       2.     -26.2652]\n",
-      "intent_name     : b''\n",
-      "magic           : b'n+1'\n"
-     ]
-    }
-   ],
-   "source": [
-    "extract_diffusion_b0(dwi_file, out_b0_file)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 11,
-   "id": "neutral-melissa",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "def fslgrad2rasb(dwi_file, fbval, fbvec, out_rasbn_file):\n",
-    "    \"\"\"Save gradient table in RAS+B format taking as input the DWI with FSL `.bval` and `.bvec`.\"\"\"\n",
-    "    import numpy as np\n",
-    "    from nibabel import load\n",
-    "    from dipy.io import read_bvals_bvecs\n",
-    "    \n",
-    "    # Read / Load\n",
-    "    img = load(str(dwi_file))\n",
-    "    bvals, bvecs = read_bvals_bvecs(str(fbval), str(fbvec))\n",
-    "    \n",
-    "    # Apply the affine transform to bvecs\n",
-    "    bvecs_tr = np.matmul(img.affine[:3,:3], bvecs.T).T\n",
-    "    \n",
-    "    # Normalize the bvecs\n",
-    "    norm = np.sum(bvecs_tr**2, axis=1)\n",
-    "    bvecs_tr_norm = np.zeros_like(bvecs_tr)\n",
-    "    for i in range(bvecs_tr.shape[0]):\n",
-    "        bvecs_tr_norm[i, :] = bvecs_tr[i, :] / norm[i] \n",
-    "    # Handles bzeros\n",
-    "    bvecs_tr_norm = np.nan_to_num(bvecs_tr_norm)\n",
-    "    \n",
-    "    rasbn = np.c_[bvecs_tr_norm, bvals]\n",
-    "    \n",
-    "    # Save Nx4 numpy matrix in TSV text file\n",
-    "    np.savetxt(fname=str(out_rasbn_file),\n",
-    "               delimiter=\"\\t\",\n",
-    "               X=rasbn)\n",
-    "    \n",
-    "    return rasbn\n",
-    "    "
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 12,
-   "id": "daily-liverpool",
-   "metadata": {},
-   "outputs": [
-    {
-     "name": "stderr",
-     "output_type": "stream",
-     "text": [
-      "<ipython-input-11-cb78b93f219c>:18: RuntimeWarning: invalid value encountered in true_divide\n",
-      "  bvecs_tr_norm[i, :] = bvecs_tr[i, :] / norm[i]\n"
-     ]
-    }
-   ],
-   "source": [
-    "rasbn = fslgrad2rasb(dwi_file, fbval, fbvec, out_rasbn_file)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 13,
-   "id": "secure-showcase",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "def create_moving_b0s_rasbn_gradient(rasbn_file, start_index, nb_of_gradient):\n",
-    "    rasbn = np.genfromtxt(fname=out_rasbn_file, delimiter=\"\\t\", skip_header=0)\n",
-    "    moving_b0s_rasbn = np.zeros((nb_of_gradient,4))\n",
-    "    moving_b0s_rasbn[0:nb_of_gradient, ...] = rasbn[start_index:start_index+nb_of_gradient, ...]\n",
-    "    # Save Nx4 numpy matrix in TSV text file\n",
-    "    np.savetxt(fname=str(rasbn_file.parent / 'gradients.moving.tsv'),\n",
-    "               delimiter=\"\\t\",\n",
-    "               X=moving_b0s_rasbn)\n",
-    "    return moving_b0s_rasbn"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 27,
-   "id": "north-stone",
-   "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "[[ 1.80353913e-01  4.66122441e-01  0.00000000e+00  1.00000000e+03]\n",
-      " [-1.27439466e-01  2.82365876e-01  3.92313651e-01  1.00000000e+03]\n",
-      " [ 4.30191553e-01 -2.31833775e-01  1.04924769e-01  1.00000000e+03]\n",
-      " [-1.53638121e-01 -3.83344627e-01  2.82253746e-01  1.00000000e+03]\n",
-      " [-3.67928622e-01  6.49873924e-03  3.38434344e-01  1.00000000e+03]\n",
-      " [ 2.66066783e-01  1.71543057e-01  3.87097161e-01  1.00000000e+03]\n",
-      " [ 8.84487882e-02  4.82220794e-01  9.74435802e-02  1.00000000e+03]\n",
-      " [ 3.85795133e-01  8.15623952e-02  3.07735418e-01  1.00000000e+03]\n",
-      " [ 3.95176644e-02 -4.98222706e-01 -1.80080496e-02  1.00000000e+03]\n",
-      " [ 5.45186999e-02 -4.60157834e-01 -1.88064506e-01  1.00000000e+03]]\n",
-      "(10, 4)\n"
-     ]
-    }
-   ],
-   "source": [
-    "moving_b0s_rasbn = create_moving_b0s_rasbn_gradient(out_rasbn_file, 6, 10)\n",
-    "print(moving_b0s_rasbn)\n",
-    "print(moving_b0s_rasbn.shape)"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "id": "suffering-equation",
-   "metadata": {},
-   "source": [
-    "## Create moving b0 volumes using the motion parameters estimated from fMRI"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 15,
-   "id": "marine-custom",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "def create_moving_b0s_and_transforms(b0_file, params_file, indices):\n",
-    "    \"\"\"Move the b0 volume according to motion parameters estimated from a list of time frame indices.\"\"\"\n",
-    "    \n",
-    "    # Get b0 filename without .nii.gz extension and make path from it\n",
-    "    b0_filename = os.path.basename(str(b0_file))\n",
-    "    index_of_dot = b0_filename.index('.')\n",
-    "    \n",
-    "    # Get b0 affine transform\n",
-    "    b0_affine = nibabel.load(str(b0_file)).affine\n",
-    "    \n",
-    "    # Initialize the matrix containing all b0 volumes\n",
-    "    nb_of_moving_b0s = len(indices) + 1\n",
-    "    x_shape, y_shape, z_shape = nibabel.load(str(b0_file)).get_fdata().shape\n",
-    "    moving_b0s_data = np.zeros((x_shape, y_shape, z_shape, nb_of_moving_b0s))\n",
-    "    \n",
-    "    # Add first b0\n",
-    "    moving_b0s_data[..., 0] = nibabel.load(str(b0_file)).get_fdata()\n",
-    "    \n",
-    "    # Read the motion parameter file generated by 3dvolreg\n",
-    "    with open(params_file) as file:\n",
-    "        lines = file.readlines()\n",
-    "        xfm_files = []\n",
-    "        xforms = []\n",
-    "        cnt=1\n",
-    "        for i in indices:\n",
-    "            # Extract the parameters for time frame i\n",
-    "            formatted_line = lines[i].split()\n",
-    "            transform_params = np.array(formatted_line).astype(float)\n",
-    "            print(transform_params)\n",
-    "            \n",
-    "            # Create a nitransforms object\n",
-    "            T = from_matvec(\n",
-    "                euler2mat(\n",
-    "                    x=transform_params[1],\n",
-    "                    y=transform_params[2],\n",
-    "                    z=transform_params[3],\n",
-    "                ),\n",
-    "                [\n",
-    "                    transform_params[4],\n",
-    "                    transform_params[5],\n",
-    "                    transform_params[6]\n",
-    "                ]\n",
-    "            )\n",
-    "            xfm = nitl.Affine(T, reference=b0_file)\n",
-    "            \n",
-    "            # Apply the transform and save into a nifti file\n",
-    "            moving_b0_filename = f'{b0_filename[:index_of_dot]}.motion-{i}.nii.gz'\n",
-    "            moving_b0_file = b0_file.parent / moving_b0_filename\n",
-    "            (~xfm).apply(b0_file, reference=b0_file).to_filename(moving_b0_file)\n",
-    "            \n",
-    "            # Save the transform\n",
-    "            xfm_filename = f'{b0_filename[:index_of_dot]}.motion-{i}.tfm'\n",
-    "            xfm_file = b0_file.parent / xfm_filename\n",
-    "            xfm.to_filename(xfm_file, fmt=\"itk\")\n",
-    "            \n",
-    "            moving_b0s_data[..., cnt] = nibabel.load(str(moving_b0_file)).get_fdata().copy()\n",
-    "            os.remove(moving_b0_file)\n",
-    "\n",
-    "            xfm_files.append(xfm_file)\n",
-    "            xforms.append(xfm)\n",
-    "            cnt += 1\n",
-    "            \n",
-    "        moving_b0s_img = nibabel.Nifti1Image(moving_b0s_data, affine=b0_affine)\n",
-    "        moving_b0s_filename = f'{b0_filename[:index_of_dot]}.moving.nii.gz'\n",
-    "        moving_b0s_file = b0_file.parent / moving_b0s_filename\n",
-    "        moving_b0s_img.to_filename(str(moving_b0s_file))\n",
-    "\n",
-    "        return xforms, xfm_files, moving_b0s_file"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 23,
-   "id": "imported-thong",
-   "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "[230 231 232 233 234 235 236 237 238 239]\n",
-      "[ 2.300e+02 -4.770e-02  6.700e-03  3.980e-02  1.290e-02 -5.750e-02\n",
-      " -9.100e-03  8.625e+00  8.480e+00]\n",
-      "[ 2.310e+02 -5.050e-02 -9.200e-03  3.880e-02  4.350e-02 -5.690e-02\n",
-      " -5.860e-02  8.495e+00  8.221e+00]\n",
-      "[ 2.320e+02 -4.190e-02  1.740e-02  3.940e-02 -1.340e-02 -5.250e-02\n",
-      " -1.200e-03  7.845e+00  7.712e+00]\n",
-      "[ 2.330e+02 -4.080e-02 -4.500e-03  4.110e-02  3.220e-02 -5.220e-02\n",
-      " -5.640e-02  8.123e+00  7.873e+00]\n",
-      "[ 2.340e+02 -4.400e-02  2.470e-02  3.420e-02 -1.190e-02 -5.110e-02\n",
-      "  1.580e-02  8.452e+00  8.311e+00]\n",
-      "[ 2.350e+02 -4.720e-02  7.000e-03  3.410e-02  2.910e-02 -5.230e-02\n",
-      " -5.420e-02  8.309e+00  8.080e+00]\n",
-      "[ 2.360e+02 -4.680e-02  1.790e-02  3.190e-02  5.700e-03 -5.650e-02\n",
-      "  1.430e-02  8.100e+00  7.946e+00]\n",
-      "[ 2.370e+02 -5.500e-02  7.400e-03  3.690e-02  3.770e-02 -5.840e-02\n",
-      " -4.430e-02  8.286e+00  8.050e+00]\n",
-      "[ 2.380e+02 -4.840e-02  1.000e-02  3.620e-02  1.900e-02 -5.940e-02\n",
-      " -3.780e-02  8.568e+00  8.395e+00]\n",
-      "[ 2.390e+02 -5.350e-02  9.300e-03  4.950e-02  4.350e-02 -6.640e-02\n",
-      " -5.300e-02  8.919e+00  8.632e+00]\n"
-     ]
-    }
-   ],
-   "source": [
-    "indices=np.arange(230,240)\n",
-    "print(indices)\n",
-    "xforms, xfm_files, moving_b0s_file = create_moving_b0s_and_transforms(\n",
-    "    out_b0_file, afni_params_file, indices\n",
-    ")"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 28,
-   "id": "advanced-belgium",
-   "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "[[ 0.99918598 -0.0397886   0.00669995  0.0129    ]\n",
-      " [ 0.039425    0.99808401  0.0476808  -0.0575    ]\n",
-      " [-0.00858427 -0.0473779   0.99883997 -0.0091    ]\n",
-      " [ 0.          0.          0.          1.        ]]\n",
-      "[[ 0.99920499 -0.0387886  -0.00919987  0.0435    ]\n",
-      " [ 0.0392049   0.99795502  0.0504764  -0.0569    ]\n",
-      " [ 0.00722315 -0.050797    0.99868298 -0.0586    ]\n",
-      " [ 0.          0.          0.          1.        ]]\n",
-      "[[ 0.99907303 -0.0393838   0.0173991  -0.0134    ]\n",
-      " [ 0.038627    0.99837601  0.0418814  -0.0525    ]\n",
-      " [-0.0190203  -0.0411705   0.99897099 -0.0012    ]\n",
-      " [ 0.          0.          0.          1.        ]]\n",
-      "[[ 0.99914497 -0.041088   -0.00449998  0.0322    ]\n",
-      " [ 0.0412376   0.99831599  0.0407883  -0.0522    ]\n",
-      " [ 0.0028165  -0.040939    0.99915802 -0.0564    ]\n",
-      " [ 0.          0.          0.          1.        ]]\n",
-      "[[ 0.99910998 -0.0341829   0.0246975  -0.0119    ]\n",
-      " [ 0.0330745   0.99848503  0.0439724  -0.0511    ]\n",
-      " [-0.0261632  -0.0431164   0.99872702  0.0158    ]\n",
-      " [ 0.          0.          0.          1.        ]]\n",
-      "[[ 0.999394   -0.0340926   0.00699994  0.0291    ]\n",
-      " [ 0.0337253   0.998317    0.0471813  -0.0523    ]\n",
-      " [-0.00859669 -0.0469167   0.99886203 -0.0542    ]\n",
-      " [ 0.          0.          0.          1.        ]]\n",
-      "[[ 0.999331   -0.0318895   0.017899    0.0057    ]\n",
-      " [ 0.0310227   0.99842399  0.0467754  -0.0565    ]\n",
-      " [-0.0193625  -0.0461889   0.99874502  0.0143    ]\n",
-      " [ 0.          0.          0.          1.        ]]\n",
-      "[[ 0.99929202 -0.0368906   0.00739993  0.0377    ]\n",
-      " [ 0.0364293   0.997823    0.0549708  -0.0584    ]\n",
-      " [-0.00941173 -0.0546623   0.99846101 -0.0443    ]\n",
-      " [ 0.          0.          0.          1.        ]]\n",
-      "[[ 0.999295   -0.0361903   0.00999983  0.019     ]\n",
-      " [ 0.0356662   0.99819201  0.0483787  -0.0594    ]\n",
-      " [-0.0117326  -0.0479879   0.998779   -0.0378    ]\n",
-      " [ 0.          0.          0.          1.        ]]\n",
-      "[[ 0.99873197 -0.0494776   0.00929987  0.0435    ]\n",
-      " [ 0.0489123   0.99737102  0.0534722  -0.0664    ]\n",
-      " [-0.0119211  -0.0529495   0.99852598 -0.053     ]\n",
-      " [ 0.          0.          0.          1.        ]]\n"
-     ]
-    }
-   ],
-   "source": [
-    "import nitransforms as nit\n",
-    "\n",
-    "_img = nibabel.load((datadir / 'b0.moving.nii.gz'))\n",
-    "_moving_b0s_data = _img.get_fdata()[..., 1:]\n",
-    "_b0 = _img.get_fdata()[..., 0]\n",
-    "_affine = _img.affine.copy()\n",
-    "\n",
-    "xfms_array = np.zeros((4,4,len(xfm_files)))\n",
-    "\n",
-    "for i, xfm in enumerate(xfm_files):\n",
-    "    fixed_b0_img = nibabel.Nifti1Image(_b0, affine=_affine)\n",
-    "    moving_b0_img = nibabel.Nifti1Image(_moving_b0s_data[..., i], affine=_affine)\n",
-    "    T = nit.io.itk.ITKLinearTransform.from_filename(\n",
-    "        str(datadir / f'b0.motion-{i + 230}.tfm')\n",
-    "    ).to_ras(reference=fixed_b0_img, moving=moving_b0_img)\n",
-    "    xfms_array[..., i] = T\n",
-    "    print(xfms_array[..., i])\n",
-    "\n",
-    "with open((datadir / 'b0.moving.transforms.npy'), 'wb') as f:\n",
-    "    np.save(f, xfms_array)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "stunning-ministry",
-   "metadata": {},
-   "outputs": [],
-   "source": []
-  }
- ],
- "metadata": {
-  "kernelspec": {
-   "display_name": "Python 3",
-   "language": "python",
-   "name": "python3"
-  },
-  "language_info": {
-   "codemirror_mode": {
-    "name": "ipython",
-    "version": 3
-   },
-   "file_extension": ".py",
-   "mimetype": "text/x-python",
-   "name": "python",
-   "nbconvert_exporter": "python",
-   "pygments_lexer": "ipython3",
-   "version": "3.8.3"
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 5
-}
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 33243ad9..81d134c5 100644
--- a/eddymotion/tests/test_integration.py
+++ b/eddymotion/tests/test_integration.py
@@ -26,46 +26,43 @@
 import nibabel as nb
 from eddymotion.dmri import DWI
 from eddymotion.estimator import EddyMotionEstimator
-import nitransforms as nit
+import nitransforms as nt
 
 
-def test_proximity_estimator_trivial_model(pkg_datadir):
+def test_proximity_estimator_trivial_model(pkg_datadir, tmp_path):
     """Check the proximity of transforms estimated by the estimator with a trivial B0 model."""
-    _img = nb.load((pkg_datadir / 'b0.moving.nii.gz'))
-    _moving_b0s_data = _img.get_fdata()[..., 1:]
-    _b0 = _img.get_fdata()[..., 0]
-    _affine = _img.affine.copy()
-    _gradients = np.genfromtxt(
-        fname=str(pkg_datadir / 'gradients.moving.tsv'),
-        delimiter="\t",
-        skip_header=0
-    ).T
-    _DWI = DWI(
-        dataobj=_moving_b0s_data,
-        affine=_affine,
-        bzero=_b0,
-        gradients=_gradients,
+
+    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,
+        dwdata=dwi_motion,
         n_iter=1,
         model="b0",
         align_kwargs=None,
         seed=None
     )
 
-    # Load reference transforms originally applied
-    ref_xfms = np.load((pkg_datadir / "b0.moving.transforms.npy"))
-
     # For each moved b0 volume
-    for i, xfm in enumerate(em_affines):
-        fixed_b0_img = nb.Nifti1Image(_b0, affine=_affine)
-        xfm2 = nit.linear.Affine(
-            ref_xfms[..., i],
-            reference=fixed_b0_img
-        )
-        assert np.all(
-            abs(xfm.map(xfm.reference.ndcoords.T) - xfm2.map(xfm.reference.ndcoords.T)) < 0.4
-        )
+    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