From bffe18e300dd1a8e266cd149b494c7fb4cfcb7ec Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Thu, 30 Sep 2021 10:41:26 +0200 Subject: [PATCH 1/6] FIX: Generation of RAS displacements fields from VSMs After much thought, I have come to the conclusion that the VOX2RAS affine must be applied in all cases, regardless of whether the dataset is oblique or not. After all, the VSM (voxel-shift-map) is no more than a regularly gridded field of voxel coordinates. The transformation between voxels and mm is biyective and univocal, and formalized by the affine matrix of the image. This PR simplifies the calculation, and theoretically should resolve most of the issues we are experiencing when resampling data. Resolves: #218. Resolves: #236. Related: nipreps/fmriprep#2210. --- sdcflows/transform.py | 44 ++++++++++++++++++------------------------- 1 file changed, 18 insertions(+), 26 deletions(-) diff --git a/sdcflows/transform.py b/sdcflows/transform.py index dc5a9f00c4..537781605b 100644 --- a/sdcflows/transform.py +++ b/sdcflows/transform.py @@ -215,38 +215,30 @@ def to_displacements(self, ro_time, pe_dir): A NIfTI 1.0 object containing the distortion. """ - from math import pi - from nibabel.affines import voxel_sizes, obliquity - from nibabel.orientations import io_orientation - - # Generate warp field - data = self.shifts.get_fdata(dtype="float32").copy() + # Set polarity & scale VSM (voxel-shift-map) by readout time + vsm = self.shifts.get_fdata().copy() pe_axis = "ijk".index(pe_dir[0]) - pe_sign = -1.0 if pe_dir.endswith("-") else 1.0 - pe_size = self.shifts.header.get_zooms()[pe_axis] - data *= pe_sign * ro_time * pe_size + vsm *= -1.0 if pe_dir.endswith("-") else 1.0 + vsm *= ro_time - fieldshape = tuple(list(data.shape[:3]) + [3]) + # Shape of displacements field + fieldshape = tuple(list(vsm.shape[:3]) + [3]) - # Compose a vector field - field = np.zeros((data.size, 3), dtype="float32") - field[..., pe_axis] = data.reshape(-1) + # Convert VSM to voxel displacements + ijk_deltas = np.zeros((vsm.size, 3), dtype="float32") + ijk_deltas[..., pe_axis] = vsm.reshape(-1) - # If coordinate system is oblique, project displacements through directions matrix + # To convert from VSM to RAS field we just apply the affine aff = self.shifts.affine - if obliquity(aff).max() * 180 / pi > 0.01: - dirmat = np.eye(4) - dirmat[:3, :3] = aff[:3, :3] / ( - voxel_sizes(aff) * io_orientation(aff)[:, 1] - ) - field = nb.affines.apply_affine(dirmat, field) - - warpnii = nb.Nifti1Image( - field.reshape(fieldshape)[:, :, :, np.newaxis, :], aff, None + xyz_deltas = nb.affines.apply_affine(aff, ijk_deltas) + xyz_nii = nb.Nifti1Image( + # WARNING: ITK NIfTI fields are 5D (have an empty 4th dim). + # Hence, the ``np.newaxis``. + xyz_deltas.reshape(fieldshape)[:, :, :, np.newaxis, :], aff, None ) - warpnii.header.set_intent("vector", (), "") - warpnii.header.set_xyzt_units("mm") - return warpnii + xyz_nii.header.set_intent("vector", (), "") + xyz_nii.header.set_xyzt_units("mm") + return xyz_nii def _cubic_bspline(d): From 5a8a75cede3c40bd66c21739ddc74a12ac2a19a3 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Thu, 30 Sep 2021 11:43:29 +0200 Subject: [PATCH 2/6] fix: do not apply translations of affine --- sdcflows/transform.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/sdcflows/transform.py b/sdcflows/transform.py index 537781605b..9996723335 100644 --- a/sdcflows/transform.py +++ b/sdcflows/transform.py @@ -229,12 +229,15 @@ def to_displacements(self, ro_time, pe_dir): ijk_deltas[..., pe_axis] = vsm.reshape(-1) # To convert from VSM to RAS field we just apply the affine - aff = self.shifts.affine + aff = self.shifts.affine.copy() + aff[:3, 3] = 0 # Translations MUST NOT be applied, though. xyz_deltas = nb.affines.apply_affine(aff, ijk_deltas) xyz_nii = nb.Nifti1Image( # WARNING: ITK NIfTI fields are 5D (have an empty 4th dim). # Hence, the ``np.newaxis``. - xyz_deltas.reshape(fieldshape)[:, :, :, np.newaxis, :], aff, None + xyz_deltas.reshape(fieldshape)[:, :, :, np.newaxis, :], + self.shifts.affine, + None, ) xyz_nii.header.set_intent("vector", (), "") xyz_nii.header.set_xyzt_units("mm") From 5b8a641ff9d6f194147fb5cf63003ff0dac05721 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Thu, 30 Sep 2021 16:27:24 +0200 Subject: [PATCH 3/6] fix: FOUND THE BUG - ITK's coordinate system is LPS! See https://github.com/poldracklab/nitransforms/blob/46004d14b80daf581b055c3f36033b309c043eb3/nitransforms/io/itk.py#L296 --- sdcflows/transform.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/sdcflows/transform.py b/sdcflows/transform.py index 9996723335..d7a10b6970 100644 --- a/sdcflows/transform.py +++ b/sdcflows/transform.py @@ -194,7 +194,7 @@ def apply( moved.header.set_data_dtype(output_dtype) return moved - def to_displacements(self, ro_time, pe_dir): + def to_displacements(self, ro_time, pe_dir, itk_format=True): """ Generate a NIfTI file containing a displacements field transform compatible with ITK/ANTs. @@ -232,6 +232,9 @@ def to_displacements(self, ro_time, pe_dir): aff = self.shifts.affine.copy() aff[:3, 3] = 0 # Translations MUST NOT be applied, though. xyz_deltas = nb.affines.apply_affine(aff, ijk_deltas) + if itk_format: + xyz_deltas[..., (0, 1)] *= -1.0 + xyz_nii = nb.Nifti1Image( # WARNING: ITK NIfTI fields are 5D (have an empty 4th dim). # Hence, the ``np.newaxis``. From 057570d045b6fa0b72eb253c4ea0b8e34327a9ec Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Thu, 30 Sep 2021 16:30:54 +0200 Subject: [PATCH 4/6] Apply suggestions from code review Co-authored-by: Chris Markiewicz --- sdcflows/transform.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/sdcflows/transform.py b/sdcflows/transform.py index d7a10b6970..5afee96641 100644 --- a/sdcflows/transform.py +++ b/sdcflows/transform.py @@ -222,11 +222,12 @@ def to_displacements(self, ro_time, pe_dir, itk_format=True): vsm *= ro_time # Shape of displacements field - fieldshape = tuple(list(vsm.shape[:3]) + [3]) + # Note that ITK NIfTI fields are 5D (have an empty 4th dimension) + fieldshape = tuple(list(vsm.shape[:3]) + [1, 3]) # Convert VSM to voxel displacements ijk_deltas = np.zeros((vsm.size, 3), dtype="float32") - ijk_deltas[..., pe_axis] = vsm.reshape(-1) + ijk_deltas[:, pe_axis] = vsm.reshape(-1) # To convert from VSM to RAS field we just apply the affine aff = self.shifts.affine.copy() @@ -236,9 +237,7 @@ def to_displacements(self, ro_time, pe_dir, itk_format=True): xyz_deltas[..., (0, 1)] *= -1.0 xyz_nii = nb.Nifti1Image( - # WARNING: ITK NIfTI fields are 5D (have an empty 4th dim). - # Hence, the ``np.newaxis``. - xyz_deltas.reshape(fieldshape)[:, :, :, np.newaxis, :], + xyz_deltas.reshape(fieldshape), self.shifts.affine, None, ) From 3642514326d135ea604079d8098b0fd43e78578d Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Thu, 30 Sep 2021 16:36:25 +0200 Subject: [PATCH 5/6] Update sdcflows/transform.py Co-authored-by: Chris Markiewicz --- sdcflows/transform.py | 1 + 1 file changed, 1 insertion(+) diff --git a/sdcflows/transform.py b/sdcflows/transform.py index 5afee96641..93c809650b 100644 --- a/sdcflows/transform.py +++ b/sdcflows/transform.py @@ -234,6 +234,7 @@ def to_displacements(self, ro_time, pe_dir, itk_format=True): aff[:3, 3] = 0 # Translations MUST NOT be applied, though. xyz_deltas = nb.affines.apply_affine(aff, ijk_deltas) if itk_format: + # ITK displacement vectors are in LPS orientation xyz_deltas[..., (0, 1)] *= -1.0 xyz_nii = nb.Nifti1Image( From d9c142b68696bab67de44204eb110e1fc4d9099e Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Thu, 30 Sep 2021 18:48:38 +0200 Subject: [PATCH 6/6] enh: add a test with visual artifacts on circle Resolves: #31. --- sdcflows/tests/test_transform.py | 136 +++++++++++++++++++++++++++++++ 1 file changed, 136 insertions(+) create mode 100644 sdcflows/tests/test_transform.py diff --git a/sdcflows/tests/test_transform.py b/sdcflows/tests/test_transform.py new file mode 100644 index 0000000000..642b991c3e --- /dev/null +++ b/sdcflows/tests/test_transform.py @@ -0,0 +1,136 @@ +# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- +# vi: set ft=python sts=4 ts=4 sw=4 et: +# +# Copyright 2021 The NiPreps Developers +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# +# We support and encourage derived works from this project, please read +# about our expectations at +# +# https://www.nipreps.org/community/licensing/ +# +"""Unit tests of the transform object.""" +from subprocess import check_call +from itertools import product +import pytest +import nibabel as nb +import numpy as np +from skimage.morphology import ball +import scipy.ndimage as nd + +from sdcflows import transform as tf + + +def generate_oracle( + coeff_file, + rotation=(None, None, None), + zooms=(2.0, 2.2, 1.5), + flip=(False, False, False), +): + """Generate an in-silico phantom, and a corresponding (aligned) B-Spline field.""" + data = ball(20) + data[19:22, ...] = 0 + data = np.pad(data + nd.binary_erosion(data, ball(3)), 8) + + zooms = [z if not f else -z for z, f in zip(zooms, flip)] + affine = np.diag(zooms + [1]) + affine[:3, 3] = -affine[:3, :3] @ ((np.array(data.shape) - 1) * 0.5) + + coeff_nii = nb.load(coeff_file) + coeff_aff = np.diag([5.0 if not f else -5.0 for f in flip] + [1]) + + if any(rotation): + R = nb.affines.from_matvec( + nb.eulerangles.euler2mat( + x=rotation[0], + y=rotation[1], + z=rotation[2], + ) + ) + affine = R @ affine + coeff_aff = R @ coeff_aff + + phantom_nii = nb.Nifti1Image( + data.astype(np.uint8), + affine, + None, + ) + coeff_nii = nb.Nifti1Image( + coeff_nii.dataobj, + coeff_aff, + coeff_nii.header, + ) + return phantom_nii, coeff_nii + + +@pytest.mark.parametrize("pe_dir", ["j", "j-", "i", "i-", "k", "k-"]) +@pytest.mark.parametrize("rotation", [(None, None, None), (0.2, None, None)]) +@pytest.mark.parametrize("flip", list(product(*[(False, True)] * 3))) +def test_displacements_field(tmpdir, testdata_dir, outdir, pe_dir, rotation, flip): + """Check the generated displacements fields.""" + tmpdir.chdir() + + # Generate test oracle + phantom_nii, coeff_nii = generate_oracle( + testdata_dir / "topup-coeff-fixed.nii.gz", + rotation=rotation, + ) + + b0 = tf.B0FieldTransform(coeffs=coeff_nii) + assert b0.fit(phantom_nii) is True + assert b0.fit(phantom_nii) is False + + b0.apply( + phantom_nii, + pe_dir=pe_dir, + ro_time=0.2, + output_dtype="float32", + ).to_filename("warped-sdcflows.nii.gz") + b0.to_displacements( + ro_time=0.2, + pe_dir=pe_dir, + ).to_filename("itk-displacements.nii.gz") + + phantom_nii.to_filename("phantom.nii.gz") + # Run antsApplyTransform + exit_code = check_call( + [ + "antsApplyTransforms -d 3 -r phantom.nii.gz -i phantom.nii.gz " + "-o warped-ants.nii.gz -n BSpline -t itk-displacements.nii.gz" + ], + shell=True, + ) + assert exit_code == 0 + + ours = np.asanyarray(nb.load("warped-sdcflows.nii.gz").dataobj) + theirs = np.asanyarray(nb.load("warped-ants.nii.gz").dataobj) + assert np.all((np.sqrt(((ours - theirs) ** 2).sum()) / ours.size) < 1e-1) + + if outdir: + from niworkflows.interfaces.reportlets.registration import ( + SimpleBeforeAfterRPT as SimpleBeforeAfter, + ) + + orientation = "".join([ax[bool(f)] for ax, f in zip(("RL", "AP", "SI"), flip)]) + + SimpleBeforeAfter( + after_label="Theirs (ANTs)", + before_label="Ours (SDCFlows)", + after="warped-ants.nii.gz", + before="warped-sdcflows.nii.gz", + out_report=str( + outdir / f"xfm_pe-{pe_dir}_flip-{orientation}_x-{rotation[0] or 0}" + f"_y-{rotation[1] or 0}_z-{rotation[2] or 0}.svg" + ), + ).run()