From 014370ee043464dad7bfa50e69a999019740d48a Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Thu, 7 Apr 2022 20:21:30 +0200 Subject: [PATCH] enh: minimize test towards the testing goal Former-commit-id: a9e67cc9503c5bc77da6af75791e39dce6df2645 --- ...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": [ - "\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", - " 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": [ - ":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