Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

ENH: Update image utility output path behaviour #81

Merged
merged 9 commits into from
Mar 20, 2020
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
112 changes: 39 additions & 73 deletions dmriprep/interfaces/images.py
Original file line number Diff line number Diff line change
@@ -1,23 +1,27 @@
"""Image tools interfaces."""
import numpy as np
import nibabel as nb
from nipype.utils.filemanip import fname_presuffix
from pathlib import Path

from nipype import logging
from nipype.interfaces.base import (
traits, TraitedSpec, BaseInterfaceInputSpec, SimpleInterface, File
BaseInterfaceInputSpec,
File,
SimpleInterface,
TraitedSpec,
traits
)

from dmriprep.utils.images import extract_b0, median, rescale_b0

LOGGER = logging.getLogger('nipype.interface')


class _ExtractB0InputSpec(BaseInterfaceInputSpec):
in_file = File(exists=True, mandatory=True, desc='dwi file')
b0_ixs = traits.List(traits.Int, mandatory=True,
desc='Index of b0s')
b0_ixs = traits.List(traits.Int, mandatory=True, desc='Index of b0s')


class _ExtractB0OutputSpec(TraitedSpec):
out_file = File(exists=True, desc='b0 file')
out_file = File(exists=True, desc='output b0 file')


class ExtractB0(SimpleInterface):
Expand All @@ -38,29 +42,19 @@ class ExtractB0(SimpleInterface):
output_spec = _ExtractB0OutputSpec

def _run_interface(self, runtime):
self._results['out_file'] = extract_b0(
self.inputs.in_file,
self.inputs.b0_ixs,
newpath=runtime.cwd)
return runtime


def extract_b0(in_file, b0_ixs, newpath=None):
"""Extract the *b0* volumes from a DWI dataset."""
out_file = fname_presuffix(
in_file, suffix='_b0', newpath=newpath)

img = nb.load(in_file)
data = img.get_fdata(dtype='float32')
from nipype.utils.filemanip import fname_presuffix

b0 = data[..., b0_ixs]
out_file = fname_presuffix(
self.inputs.in_file,
suffix='_b0',
use_ext=True,
newpath=str(Path(runtime.cwd).absolute()),
)

hdr = img.header.copy()
hdr.set_data_shape(b0.shape)
hdr.set_xyzt_units('mm')
hdr.set_data_dtype(np.float32)
nb.Nifti1Image(b0, img.affine, hdr).to_filename(out_file)
return out_file
self._results['out_file'] = extract_b0(
self.inputs.in_file, self.inputs.b0_ixs, out_file
)
return runtime


class _RescaleB0InputSpec(BaseInterfaceInputSpec):
Expand Down Expand Up @@ -91,55 +85,27 @@ class RescaleB0(SimpleInterface):
output_spec = _RescaleB0OutputSpec

def _run_interface(self, runtime):
from nipype.utils.filemanip import fname_presuffix

out_b0s = fname_presuffix(
self.inputs.in_file,
suffix='_rescaled',
use_ext=True,
newpath=str(Path(runtime.cwd).absolute())
)
out_ref = fname_presuffix(
self.inputs.in_file,
suffix='_ref',
use_ext=True,
newpath=str(Path(runtime.cwd).absolute())
)

self._results['out_b0s'] = rescale_b0(
self.inputs.in_file,
self.inputs.mask_file,
newpath=runtime.cwd
self.inputs.mask_file, out_b0s
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this is what introduced #90 - the mask and the b0s are swapped.

Copy link
Collaborator Author

@josephmje josephmje Apr 6, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

it failed at the very first step where the b0s get extracted. the difference between the 2 dwi images for that subject is that one contains multiple b0s and the borked contains a single b0. it looks like this single b0 case is messing things up.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

yup I realized after commenting. I'm still going through this - I might have found a solution.

)
self._results['out_ref'] = median(
self._results['out_b0s'],
newpath=runtime.cwd
out_path=out_ref
)
return runtime


def rescale_b0(in_file, mask_file, newpath=None):
"""Rescale the input volumes using the median signal intensity."""
out_file = fname_presuffix(
in_file, suffix='_rescaled_b0', newpath=newpath)

img = nb.load(in_file)
if img.dataobj.ndim == 3:
return in_file

data = img.get_fdata(dtype='float32')
mask_img = nb.load(mask_file)
mask_data = mask_img.get_fdata(dtype='float32')

median_signal = np.median(data[mask_data > 0, ...], axis=0)
rescaled_data = 1000 * data / median_signal
hdr = img.header.copy()
nb.Nifti1Image(rescaled_data, img.affine, hdr).to_filename(out_file)
return out_file


def median(in_file, newpath=None):
"""Average a 4D dataset across the last dimension using median."""
out_file = fname_presuffix(
in_file, suffix='_b0ref', newpath=newpath)

img = nb.load(in_file)
if img.dataobj.ndim == 3:
return in_file
if img.shape[-1] == 1:
nb.squeeze_image(img).to_filename(out_file)
return out_file

median_data = np.median(img.get_fdata(dtype='float32'),
axis=-1)

hdr = img.header.copy()
hdr.set_xyzt_units('mm')
hdr.set_data_dtype(np.float32)
nb.Nifti1Image(median_data, img.affine, hdr).to_filename(out_file)
return out_file
70 changes: 70 additions & 0 deletions dmriprep/utils/images.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,70 @@
import nibabel as nb
import numpy as np
from nipype.utils.filemanip import fname_presuffix


def extract_b0(in_file, b0_ixs, out_path=None):
"""Extract the *b0* volumes from a DWI dataset."""
if out_path is None:
out_path = fname_presuffix(
in_file, suffix='_b0', use_ext=True)

img = nb.load(in_file)
data = img.get_fdata()

b0 = data[..., b0_ixs]

hdr = img.header.copy()
hdr.set_data_shape(b0.shape)
hdr.set_xyzt_units('mm')
nb.Nifti1Image(b0.astype(hdr.get_data_dtype()), img.affine, hdr).to_filename(out_path)
return out_path


def rescale_b0(in_file, mask_file, out_path=None):
"""Rescale the input volumes using the median signal intensity."""
if out_path is None:
out_path = fname_presuffix(
in_file, suffix='_rescaled', use_ext=True)

img = nb.load(in_file)
if img.dataobj.ndim == 3:
return in_file

data = img.get_fdata()
mask_img = nb.load(mask_file)
mask_data = mask_img.get_fdata()

median_signal = np.median(data[mask_data > 0, ...], axis=0)
rescaled_data = 1000 * data / median_signal
hdr = img.header.copy()
nb.Nifti1Image(
rescaled_data.astype(hdr.get_data_dtype()), img.affine, hdr
).to_filename(out_path)
return out_path


def median(in_file, dtype=None, out_path=None):
"""Average a 4D dataset across the last dimension using median."""
if out_path is None:
out_path = fname_presuffix(
in_file, suffix='_b0ref', use_ext=True)

img = nb.load(in_file)
if img.dataobj.ndim == 3:
return in_file
if img.shape[-1] == 1:
nb.squeeze_image(img).to_filename(out_path)
return out_path

median_data = np.median(img.get_fdata(dtype=dtype),
axis=-1)

hdr = img.header.copy()
hdr.set_xyzt_units('mm')
if dtype is not None:
hdr.set_data_dtype(dtype)
else:
dtype = hdr.get_data_dtype()
nb.Nifti1Image(median_data.astype(dtype), img.affine, hdr).to_filename(out_path)
return out_path