|
| 1 | +# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- |
| 2 | +# vi: set ft=python sts=4 ts=4 sw=4 et: |
| 3 | +# |
| 4 | +# Copyright 2021 The NiPreps Developers <nipreps@gmail.com> |
| 5 | +# |
| 6 | +# Licensed under the Apache License, Version 2.0 (the "License"); |
| 7 | +# you may not use this file except in compliance with the License. |
| 8 | +# You may obtain a copy of the License at |
| 9 | +# |
| 10 | +# http://www.apache.org/licenses/LICENSE-2.0 |
| 11 | +# |
| 12 | +# Unless required by applicable law or agreed to in writing, software |
| 13 | +# distributed under the License is distributed on an "AS IS" BASIS, |
| 14 | +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
| 15 | +# See the License for the specific language governing permissions and |
| 16 | +# limitations under the License. |
| 17 | +# |
| 18 | +# We support and encourage derived works from this project, please read |
| 19 | +# about our expectations at |
| 20 | +# |
| 21 | +# https://www.nipreps.org/community/licensing/ |
| 22 | +# |
| 23 | +"""Unit tests of the transform object.""" |
| 24 | +from subprocess import check_call |
| 25 | +from itertools import product |
| 26 | +import pytest |
| 27 | +import nibabel as nb |
| 28 | +import numpy as np |
| 29 | +from skimage.morphology import ball |
| 30 | +import scipy.ndimage as nd |
| 31 | + |
| 32 | +from sdcflows import transform as tf |
| 33 | + |
| 34 | + |
| 35 | +def generate_oracle( |
| 36 | + coeff_file, |
| 37 | + rotation=(None, None, None), |
| 38 | + zooms=(2.0, 2.2, 1.5), |
| 39 | + flip=(False, False, False), |
| 40 | +): |
| 41 | + """Generate an in-silico phantom, and a corresponding (aligned) B-Spline field.""" |
| 42 | + data = ball(20) |
| 43 | + data[19:22, ...] = 0 |
| 44 | + data = np.pad(data + nd.binary_erosion(data, ball(3)), 8) |
| 45 | + |
| 46 | + zooms = [z if not f else -z for z, f in zip(zooms, flip)] |
| 47 | + affine = np.diag(zooms + [1]) |
| 48 | + affine[:3, 3] = -affine[:3, :3] @ ((np.array(data.shape) - 1) * 0.5) |
| 49 | + |
| 50 | + coeff_nii = nb.load(coeff_file) |
| 51 | + coeff_aff = np.diag([5.0 if not f else -5.0 for f in flip] + [1]) |
| 52 | + |
| 53 | + if any(rotation): |
| 54 | + R = nb.affines.from_matvec( |
| 55 | + nb.eulerangles.euler2mat( |
| 56 | + x=rotation[0], |
| 57 | + y=rotation[1], |
| 58 | + z=rotation[2], |
| 59 | + ) |
| 60 | + ) |
| 61 | + affine = R @ affine |
| 62 | + coeff_aff = R @ coeff_aff |
| 63 | + |
| 64 | + phantom_nii = nb.Nifti1Image( |
| 65 | + data.astype(np.uint8), |
| 66 | + affine, |
| 67 | + None, |
| 68 | + ) |
| 69 | + coeff_nii = nb.Nifti1Image( |
| 70 | + coeff_nii.dataobj, |
| 71 | + coeff_aff, |
| 72 | + coeff_nii.header, |
| 73 | + ) |
| 74 | + return phantom_nii, coeff_nii |
| 75 | + |
| 76 | + |
| 77 | +@pytest.mark.parametrize("pe_dir", ["j", "j-", "i", "i-", "k", "k-"]) |
| 78 | +@pytest.mark.parametrize("rotation", [(None, None, None), (0.2, None, None)]) |
| 79 | +@pytest.mark.parametrize("flip", list(product(*[(False, True)] * 3))) |
| 80 | +def test_displacements_field(tmpdir, testdata_dir, outdir, pe_dir, rotation, flip): |
| 81 | + """Check the generated displacements fields.""" |
| 82 | + tmpdir.chdir() |
| 83 | + |
| 84 | + # Generate test oracle |
| 85 | + phantom_nii, coeff_nii = generate_oracle( |
| 86 | + testdata_dir / "topup-coeff-fixed.nii.gz", |
| 87 | + rotation=rotation, |
| 88 | + ) |
| 89 | + |
| 90 | + b0 = tf.B0FieldTransform(coeffs=coeff_nii) |
| 91 | + assert b0.fit(phantom_nii) is True |
| 92 | + assert b0.fit(phantom_nii) is False |
| 93 | + |
| 94 | + b0.apply( |
| 95 | + phantom_nii, |
| 96 | + pe_dir=pe_dir, |
| 97 | + ro_time=0.2, |
| 98 | + output_dtype="float32", |
| 99 | + ).to_filename("warped-sdcflows.nii.gz") |
| 100 | + b0.to_displacements( |
| 101 | + ro_time=0.2, |
| 102 | + pe_dir=pe_dir, |
| 103 | + ).to_filename("itk-displacements.nii.gz") |
| 104 | + |
| 105 | + phantom_nii.to_filename("phantom.nii.gz") |
| 106 | + # Run antsApplyTransform |
| 107 | + exit_code = check_call( |
| 108 | + [ |
| 109 | + "antsApplyTransforms -d 3 -r phantom.nii.gz -i phantom.nii.gz " |
| 110 | + "-o warped-ants.nii.gz -n BSpline -t itk-displacements.nii.gz" |
| 111 | + ], |
| 112 | + shell=True, |
| 113 | + ) |
| 114 | + assert exit_code == 0 |
| 115 | + |
| 116 | + ours = np.asanyarray(nb.load("warped-sdcflows.nii.gz").dataobj) |
| 117 | + theirs = np.asanyarray(nb.load("warped-ants.nii.gz").dataobj) |
| 118 | + assert np.all((np.sqrt(((ours - theirs) ** 2).sum()) / ours.size) < 1e-1) |
| 119 | + |
| 120 | + if outdir: |
| 121 | + from niworkflows.interfaces.reportlets.registration import ( |
| 122 | + SimpleBeforeAfterRPT as SimpleBeforeAfter, |
| 123 | + ) |
| 124 | + |
| 125 | + orientation = "".join([ax[bool(f)] for ax, f in zip(("RL", "AP", "SI"), flip)]) |
| 126 | + |
| 127 | + SimpleBeforeAfter( |
| 128 | + after_label="Theirs (ANTs)", |
| 129 | + before_label="Ours (SDCFlows)", |
| 130 | + after="warped-ants.nii.gz", |
| 131 | + before="warped-sdcflows.nii.gz", |
| 132 | + out_report=str( |
| 133 | + outdir / f"xfm_pe-{pe_dir}_flip-{orientation}_x-{rotation[0] or 0}" |
| 134 | + f"_y-{rotation[1] or 0}_z-{rotation[2] or 0}.svg" |
| 135 | + ), |
| 136 | + ).run() |
0 commit comments