# -*- coding: utf-8 -*-
# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
# vi: set ft=python sts=4 ts=4 sw=4 et:
"""Miscellaneous algorithms."""
import os
import os.path as op
import nibabel as nb
import numpy as np
from math import floor, ceil
import itertools
import warnings
from .. import logging
from . import metrics as nam
from ..interfaces.base import (
BaseInterface,
traits,
TraitedSpec,
File,
InputMultiPath,
OutputMultiPath,
BaseInterfaceInputSpec,
isdefined,
DynamicTraitedSpec,
Undefined,
)
from ..utils.filemanip import fname_presuffix, split_filename, ensure_list
from . import confounds
iflogger = logging.getLogger("nipype.interface")
class PickAtlasInputSpec(BaseInterfaceInputSpec):
atlas = File(
exists=True, desc="Location of the atlas that will be used.", mandatory=True
)
labels = traits.Either(
traits.Int,
traits.List(traits.Int),
desc=(
"Labels of regions that will be included in the mask. Must be\
compatible with the atlas used."
),
mandatory=True,
)
hemi = traits.Enum(
"both",
"left",
"right",
desc="Restrict the mask to only one hemisphere: left or right",
usedefault=True,
)
dilation_size = traits.Int(
usedefault=True,
desc="Defines how much the mask will be dilated (expanded in 3D).",
)
output_file = File(desc="Where to store the output mask.")
class PickAtlasOutputSpec(TraitedSpec):
mask_file = File(exists=True, desc="output mask file")
class PickAtlas(BaseInterface):
"""Returns ROI masks given an atlas and a list of labels. Supports dilation
and left right masking (assuming the atlas is properly aligned).
"""
input_spec = PickAtlasInputSpec
output_spec = PickAtlasOutputSpec
def _run_interface(self, runtime):
nim = self._get_brodmann_area()
nb.save(nim, self._gen_output_filename())
return runtime
def _gen_output_filename(self):
if not isdefined(self.inputs.output_file):
output = fname_presuffix(
fname=self.inputs.atlas,
suffix="_mask",
newpath=os.getcwd(),
use_ext=True,
)
else:
output = os.path.realpath(self.inputs.output_file)
return output
def _get_brodmann_area(self):
nii = nb.load(self.inputs.atlas)
origdata = np.asanyarray(nii.dataobj)
newdata = np.zeros(origdata.shape)
if not isinstance(self.inputs.labels, list):
labels = [self.inputs.labels]
else:
labels = self.inputs.labels
for lab in labels:
newdata[origdata == lab] = 1
if self.inputs.hemi == "right":
newdata[int(floor(float(origdata.shape[0]) / 2)) :, :, :] = 0
elif self.inputs.hemi == "left":
newdata[: int(ceil(float(origdata.shape[0]) / 2)), :, :] = 0
if self.inputs.dilation_size != 0:
from scipy.ndimage.morphology import grey_dilation
newdata = grey_dilation(
newdata,
(
2 * self.inputs.dilation_size + 1,
2 * self.inputs.dilation_size + 1,
2 * self.inputs.dilation_size + 1,
),
)
return nb.Nifti1Image(newdata, nii.affine, nii.header)
def _list_outputs(self):
outputs = self._outputs().get()
outputs["mask_file"] = self._gen_output_filename()
return outputs
class SimpleThresholdInputSpec(BaseInterfaceInputSpec):
volumes = InputMultiPath(
File(exists=True), desc="volumes to be thresholded", mandatory=True
)
threshold = traits.Float(
desc="volumes to be thresholdedeverything below this value will be set\
to zero",
mandatory=True,
)
class SimpleThresholdOutputSpec(TraitedSpec):
thresholded_volumes = OutputMultiPath(File(exists=True), desc="thresholded volumes")
class SimpleThreshold(BaseInterface):
"""Applies a threshold to input volumes"""
input_spec = SimpleThresholdInputSpec
output_spec = SimpleThresholdOutputSpec
def _run_interface(self, runtime):
for fname in self.inputs.volumes:
img = nb.load(fname)
data = img.get_fdata()
active_map = data > self.inputs.threshold
thresholded_map = np.zeros(data.shape)
thresholded_map[active_map] = data[active_map]
new_img = nb.Nifti1Image(thresholded_map, img.affine, img.header)
_, base, _ = split_filename(fname)
nb.save(new_img, base + "_thresholded.nii")
return runtime
def _list_outputs(self):
outputs = self._outputs().get()
outputs["thresholded_volumes"] = []
for fname in self.inputs.volumes:
_, base, _ = split_filename(fname)
outputs["thresholded_volumes"].append(
os.path.abspath(base + "_thresholded.nii")
)
return outputs
class ModifyAffineInputSpec(BaseInterfaceInputSpec):
volumes = InputMultiPath(
File(exists=True),
desc="volumes which affine matrices will be modified",
mandatory=True,
)
transformation_matrix = traits.Array(
value=np.eye(4),
shape=(4, 4),
desc="transformation matrix that will be left multiplied by the\
affine matrix",
usedefault=True,
)
class ModifyAffineOutputSpec(TraitedSpec):
transformed_volumes = OutputMultiPath(File(exist=True))
class ModifyAffine(BaseInterface):
"""Left multiplies the affine matrix with a specified values. Saves the volume
as a nifti file.
"""
input_spec = ModifyAffineInputSpec
output_spec = ModifyAffineOutputSpec
def _gen_output_filename(self, name):
_, base, _ = split_filename(name)
return os.path.abspath(base + "_transformed.nii")
def _run_interface(self, runtime):
for fname in self.inputs.volumes:
img = nb.load(fname)
affine = img.affine
affine = np.dot(self.inputs.transformation_matrix, affine)
nb.save(
nb.Nifti1Image(img.dataobj, affine, img.header),
self._gen_output_filename(fname),
)
return runtime
def _list_outputs(self):
outputs = self._outputs().get()
outputs["transformed_volumes"] = []
for fname in self.inputs.volumes:
outputs["transformed_volumes"].append(self._gen_output_filename(fname))
return outputs
class CreateNiftiInputSpec(BaseInterfaceInputSpec):
data_file = File(exists=True, mandatory=True, desc="ANALYZE img file")
header_file = File(
exists=True, mandatory=True, desc="corresponding ANALYZE hdr file"
)
affine = traits.Array(desc="affine transformation array")
class CreateNiftiOutputSpec(TraitedSpec):
nifti_file = File(exists=True)
class CreateNifti(BaseInterface):
"""Creates a nifti volume"""
input_spec = CreateNiftiInputSpec
output_spec = CreateNiftiOutputSpec
def _gen_output_file_name(self):
_, base, _ = split_filename(self.inputs.data_file)
return os.path.abspath(base + ".nii")
def _run_interface(self, runtime):
with open(self.inputs.header_file, "rb") as hdr_file:
hdr = nb.AnalyzeHeader.from_fileobj(hdr_file)
if isdefined(self.inputs.affine):
affine = self.inputs.affine
else:
affine = None
with open(self.inputs.data_file, "rb") as data_file:
data = hdr.data_from_fileobj(data_file)
img = nb.Nifti1Image(data, affine, hdr)
nb.save(img, self._gen_output_file_name())
return runtime
def _list_outputs(self):
outputs = self._outputs().get()
outputs["nifti_file"] = self._gen_output_file_name()
return outputs
class GunzipInputSpec(BaseInterfaceInputSpec):
in_file = File(exists=True, mandatory=True)
class GunzipOutputSpec(TraitedSpec):
out_file = File(exists=True)
class Gunzip(BaseInterface):
"""Gunzip wrapper
>>> from nipype.algorithms.misc import Gunzip
>>> gunzip = Gunzip(in_file='tpms_msk.nii.gz')
>>> res = gunzip.run()
>>> res.outputs.out_file # doctest: +ELLIPSIS
'.../tpms_msk.nii'
.. testcleanup::
>>> os.unlink('tpms_msk.nii')
"""
input_spec = GunzipInputSpec
output_spec = GunzipOutputSpec
def _gen_output_file_name(self):
_, base, ext = split_filename(self.inputs.in_file)
if ext[-3:].lower() == ".gz":
ext = ext[:-3]
return os.path.abspath(base + ext)
def _run_interface(self, runtime):
import gzip
import shutil
with gzip.open(self.inputs.in_file, "rb") as in_file:
with open(self._gen_output_file_name(), "wb") as out_file:
shutil.copyfileobj(in_file, out_file)
return runtime
def _list_outputs(self):
outputs = self._outputs().get()
outputs["out_file"] = self._gen_output_file_name()
return outputs
def replaceext(in_list, ext):
out_list = list()
for filename in in_list:
path, name, _ = split_filename(op.abspath(filename))
out_name = op.join(path, name) + ext
out_list.append(out_name)
return out_list
def _matlab2csv(in_array, name, reshape):
output_array = np.asarray(in_array)
if reshape:
if len(np.shape(output_array)) > 1:
output_array = np.reshape(
output_array, (np.shape(output_array)[0] * np.shape(output_array)[1], 1)
)
iflogger.info(np.shape(output_array))
output_name = op.abspath(name + ".csv")
np.savetxt(output_name, output_array, delimiter=",")
return output_name
class Matlab2CSVInputSpec(TraitedSpec):
in_file = File(exists=True, mandatory=True, desc="Input MATLAB .mat file")
reshape_matrix = traits.Bool(
True,
usedefault=True,
desc="The output of this interface is meant for R, so matrices will be\
reshaped to vectors by default.",
)
class Matlab2CSVOutputSpec(TraitedSpec):
csv_files = OutputMultiPath(
File(
desc="Output CSV files for each variable saved in the input .mat\
file"
)
)
class Matlab2CSV(BaseInterface):
"""
Save the components of a MATLAB .mat file as a text file with comma-separated values (CSVs).
CSV files are easily loaded in R, for use in statistical processing.
For further information, see cran.r-project.org/doc/manuals/R-data.pdf
Example
-------
>>> from nipype.algorithms import misc
>>> mat2csv = misc.Matlab2CSV()
>>> mat2csv.inputs.in_file = 'cmatrix.mat'
>>> mat2csv.run() # doctest: +SKIP
"""
input_spec = Matlab2CSVInputSpec
output_spec = Matlab2CSVOutputSpec
def _run_interface(self, runtime):
import scipy.io as sio
in_dict = sio.loadmat(op.abspath(self.inputs.in_file))
# Check if the file has multiple variables in it. If it does, loop
# through them and save them as individual CSV files.
# If not, save the variable as a single CSV file using the input file
# name and a .csv extension.
saved_variables = list()
for key in list(in_dict.keys()):
if not key.startswith("__"):
if isinstance(in_dict[key][0], np.ndarray):
saved_variables.append(key)
else:
iflogger.info(
"One of the keys in the input file, %s, is "
"not a Numpy array",
key,
)
if len(saved_variables) > 1:
iflogger.info("%i variables found:", len(saved_variables))
iflogger.info(saved_variables)
for variable in saved_variables:
iflogger.info(
"...Converting %s - type %s - to CSV",
variable,
type(in_dict[variable]),
)
_matlab2csv(in_dict[variable], variable, self.inputs.reshape_matrix)
elif len(saved_variables) == 1:
_, name, _ = split_filename(self.inputs.in_file)
variable = saved_variables[0]
iflogger.info(
"Single variable found %s, type %s:", variable, type(in_dict[variable])
)
iflogger.info(
"...Converting %s to CSV from %s", variable, self.inputs.in_file
)
_matlab2csv(in_dict[variable], name, self.inputs.reshape_matrix)
else:
iflogger.error("No values in the MATLAB file?!")
return runtime
def _list_outputs(self):
import scipy.io as sio
outputs = self.output_spec().get()
in_dict = sio.loadmat(op.abspath(self.inputs.in_file))
saved_variables = list()
for key in list(in_dict.keys()):
if not key.startswith("__"):
if isinstance(in_dict[key][0], np.ndarray):
saved_variables.append(key)
else:
iflogger.error(
"One of the keys in the input file, %s, is "
"not a Numpy array",
key,
)
if len(saved_variables) > 1:
outputs["csv_files"] = replaceext(saved_variables, ".csv")
elif len(saved_variables) == 1:
_, name, ext = split_filename(self.inputs.in_file)
outputs["csv_files"] = op.abspath(name + ".csv")
else:
iflogger.error("No values in the MATLAB file?!")
return outputs
def merge_csvs(in_list):
for idx, in_file in enumerate(in_list):
try:
in_array = np.loadtxt(in_file, delimiter=",")
except ValueError:
try:
in_array = np.loadtxt(in_file, delimiter=",", skiprows=1)
except ValueError:
with open(in_file, "r") as first:
header_line = first.readline()
header_list = header_line.split(",")
n_cols = len(header_list)
try:
in_array = np.loadtxt(
in_file,
delimiter=",",
skiprows=1,
usecols=list(range(1, n_cols)),
)
except ValueError:
in_array = np.loadtxt(
in_file,
delimiter=",",
skiprows=1,
usecols=list(range(1, n_cols - 1)),
)
if idx == 0:
out_array = in_array
else:
out_array = np.dstack((out_array, in_array))
out_array = np.squeeze(out_array)
iflogger.info("Final output array shape:")
iflogger.info(np.shape(out_array))
return out_array
def remove_identical_paths(in_files):
import os.path as op
from ..utils.filemanip import split_filename
if len(in_files) > 1:
out_names = list()
commonprefix = op.commonprefix(in_files)
lastslash = commonprefix.rfind("/")
commonpath = commonprefix[0 : (lastslash + 1)]
for fileidx, in_file in enumerate(in_files):
path, name, ext = split_filename(in_file)
in_file = op.join(path, name)
name = in_file.replace(commonpath, "")
name = name.replace("_subject_id_", "")
out_names.append(name)
else:
path, name, ext = split_filename(in_files[0])
out_names = [name]
return out_names
def maketypelist(rowheadings, shape, extraheadingBool, extraheading):
typelist = []
if rowheadings:
typelist.append(("heading", "a40"))
if len(shape) > 1:
for idx in range(1, (min(shape) + 1)):
typelist.append((str(idx), float))
else:
for idx in range(1, (shape[0] + 1)):
typelist.append((str(idx), float))
if extraheadingBool:
typelist.append((extraheading, "a40"))
iflogger.info(typelist)
return typelist
def makefmtlist(output_array, typelist, rowheadingsBool, shape, extraheadingBool):
fmtlist = []
if rowheadingsBool:
fmtlist.append("%s")
if len(shape) > 1:
output = np.zeros(max(shape), typelist)
for idx in range(1, min(shape) + 1):
output[str(idx)] = output_array[:, idx - 1]
fmtlist.append("%f")
else:
output = np.zeros(1, typelist)
for idx in range(1, len(output_array) + 1):
output[str(idx)] = output_array[idx - 1]
fmtlist.append("%f")
if extraheadingBool:
fmtlist.append("%s")
fmt = ",".join(fmtlist)
return fmt, output
class MergeCSVFilesInputSpec(TraitedSpec):
in_files = InputMultiPath(
File(exists=True),
mandatory=True,
desc="Input comma-separated value (CSV) files",
)
out_file = File(
"merged.csv", usedefault=True, desc="Output filename for merged CSV file"
)
column_headings = traits.List(
traits.Str,
desc="List of column headings to save in merged CSV file\
(must be equal to number of input files). If left undefined, these\
will be pulled from the input filenames.",
)
row_headings = traits.List(
traits.Str,
desc="List of row headings to save in merged CSV file\
(must be equal to number of rows in the input files).",
)
row_heading_title = traits.Str(
"label",
usedefault=True,
desc="Column heading for the row headings\
added",
)
extra_column_heading = traits.Str(desc="New heading to add for the added field.")
extra_field = traits.Str(
desc="New field to add to each row. This is useful for saving the\
group or subject ID in the file."
)
class MergeCSVFilesOutputSpec(TraitedSpec):
csv_file = File(desc="Output CSV file containing columns ")
class MergeCSVFiles(BaseInterface):
"""
Merge several CSV files into a single CSV file.
This interface is designed to facilitate data loading in the R environment.
If provided, it will also incorporate column heading names into the
resulting CSV file.
CSV files are easily loaded in R, for use in statistical processing.
For further information, see cran.r-project.org/doc/manuals/R-data.pdf
Example
-------
>>> from nipype.algorithms import misc
>>> mat2csv = misc.MergeCSVFiles()
>>> mat2csv.inputs.in_files = ['degree.mat','clustering.mat']
>>> mat2csv.inputs.column_headings = ['degree','clustering']
>>> mat2csv.run() # doctest: +SKIP
"""
input_spec = MergeCSVFilesInputSpec
output_spec = MergeCSVFilesOutputSpec
def _run_interface(self, runtime):
extraheadingBool = False
extraheading = ""
rowheadingsBool = False
"""
This block defines the column headings.
"""
if isdefined(self.inputs.column_headings):
iflogger.info("Column headings have been provided:")
headings = self.inputs.column_headings
else:
iflogger.info("Column headings not provided! Pulled from input filenames:")
headings = remove_identical_paths(self.inputs.in_files)
if isdefined(self.inputs.extra_field):
if isdefined(self.inputs.extra_column_heading):
extraheading = self.inputs.extra_column_heading
iflogger.info("Extra column heading provided: %s", extraheading)
else:
extraheading = "type"
iflogger.info('Extra column heading was not defined. Using "type"')
headings.append(extraheading)
extraheadingBool = True
if len(self.inputs.in_files) == 1:
iflogger.warning("Only one file input!")
if isdefined(self.inputs.row_headings):
iflogger.info(
'Row headings have been provided. Adding "labels"' "column header."
)
prefix = '"{p}","'.format(p=self.inputs.row_heading_title)
csv_headings = prefix + '","'.join(itertools.chain(headings)) + '"\n'
rowheadingsBool = True
else:
iflogger.info("Row headings have not been provided.")
csv_headings = '"' + '","'.join(itertools.chain(headings)) + '"\n'
iflogger.info("Final Headings:")
iflogger.info(csv_headings)
"""
Next we merge the arrays and define the output text file
"""
output_array = merge_csvs(self.inputs.in_files)
_, name, ext = split_filename(self.inputs.out_file)
if not ext == ".csv":
ext = ".csv"
out_file = op.abspath(name + ext)
with open(out_file, "w") as file_handle:
file_handle.write(csv_headings)
shape = np.shape(output_array)
typelist = maketypelist(rowheadingsBool, shape, extraheadingBool, extraheading)
fmt, output = makefmtlist(
output_array, typelist, rowheadingsBool, shape, extraheadingBool
)
if rowheadingsBool:
row_heading_list = self.inputs.row_headings
row_heading_list_with_quotes = []
for row_heading in row_heading_list:
row_heading_with_quotes = '"' + row_heading + '"'
row_heading_list_with_quotes.append(row_heading_with_quotes)
row_headings = np.array(row_heading_list_with_quotes, dtype="|S40")
output["heading"] = row_headings
if isdefined(self.inputs.extra_field):
extrafieldlist = []
if len(shape) > 1:
mx = shape[0]
else:
mx = 1
for idx in range(0, mx):
extrafieldlist.append(self.inputs.extra_field)
iflogger.info(len(extrafieldlist))
output[extraheading] = extrafieldlist
iflogger.info(output)
iflogger.info(fmt)
with open(out_file, "a") as file_handle:
np.savetxt(file_handle, output, fmt, delimiter=",")
return runtime
def _list_outputs(self):
outputs = self.output_spec().get()
_, name, ext = split_filename(self.inputs.out_file)
if not ext == ".csv":
ext = ".csv"
out_file = op.abspath(name + ext)
outputs["csv_file"] = out_file
return outputs
class AddCSVColumnInputSpec(TraitedSpec):
in_file = File(
exists=True, mandatory=True, desc="Input comma-separated value (CSV) files"
)
out_file = File(
"extra_heading.csv", usedefault=True, desc="Output filename for merged CSV file"
)
extra_column_heading = traits.Str(desc="New heading to add for the added field.")
extra_field = traits.Str(
desc="New field to add to each row. This is useful for saving the\
group or subject ID in the file."
)
class AddCSVColumnOutputSpec(TraitedSpec):
csv_file = File(desc="Output CSV file containing columns ")
class AddCSVColumn(BaseInterface):
"""
Short interface to add an extra column and field to a text file.
Example
-------
>>> from nipype.algorithms import misc
>>> addcol = misc.AddCSVColumn()
>>> addcol.inputs.in_file = 'degree.csv'
>>> addcol.inputs.extra_column_heading = 'group'
>>> addcol.inputs.extra_field = 'male'
>>> addcol.run() # doctest: +SKIP
"""
input_spec = AddCSVColumnInputSpec
output_spec = AddCSVColumnOutputSpec
def _run_interface(self, runtime):
in_file = open(self.inputs.in_file, "r")
_, name, ext = split_filename(self.inputs.out_file)
if not ext == ".csv":
ext = ".csv"
out_file = op.abspath(name + ext)
out_file = open(out_file, "w")
firstline = in_file.readline()
firstline = firstline.replace("\n", "")
new_firstline = firstline + ',"' + self.inputs.extra_column_heading + '"\n'
out_file.write(new_firstline)
for line in in_file:
new_line = line.replace("\n", "")
new_line = new_line + "," + self.inputs.extra_field + "\n"
out_file.write(new_line)
in_file.close()
out_file.close()
return runtime
def _list_outputs(self):
outputs = self.output_spec().get()
_, name, ext = split_filename(self.inputs.out_file)
if not ext == ".csv":
ext = ".csv"
out_file = op.abspath(name + ext)
outputs["csv_file"] = out_file
return outputs
class AddCSVRowInputSpec(DynamicTraitedSpec, BaseInterfaceInputSpec):
in_file = File(mandatory=True, desc="Input comma-separated value (CSV) files")
_outputs = traits.Dict(traits.Any, value={}, usedefault=True)
def __setattr__(self, key, value):
if key not in self.copyable_trait_names():
if not isdefined(value):
super(AddCSVRowInputSpec, self).__setattr__(key, value)
self._outputs[key] = value
else:
if key in self._outputs:
self._outputs[key] = value
super(AddCSVRowInputSpec, self).__setattr__(key, value)
class AddCSVRowOutputSpec(TraitedSpec):
csv_file = File(desc="Output CSV file containing rows ")
class AddCSVRow(BaseInterface):
"""
Simple interface to add an extra row to a CSV file.
.. note:: Requires `pandas `_
.. warning:: Multi-platform thread-safe execution is possible with
`lockfile `_. Please
recall that (1) this module is alpha software; and (2) it should be
installed for thread-safe writing.
If lockfile is not installed, then the interface is not thread-safe.
Example
-------
>>> from nipype.algorithms import misc
>>> addrow = misc.AddCSVRow()
>>> addrow.inputs.in_file = 'scores.csv'
>>> addrow.inputs.si = 0.74
>>> addrow.inputs.di = 0.93
>>> addrow.inputs.subject_id = 'S400'
>>> addrow.inputs.list_of_values = [ 0.4, 0.7, 0.3 ]
>>> addrow.run() # doctest: +SKIP
"""
input_spec = AddCSVRowInputSpec
output_spec = AddCSVRowOutputSpec
def __init__(self, infields=None, force_run=True, **kwargs):
super(AddCSVRow, self).__init__(**kwargs)
undefined_traits = {}
self._infields = infields
self._have_lock = False
self._lock = None
if infields:
for key in infields:
self.inputs.add_trait(key, traits.Any)
self.inputs._outputs[key] = Undefined
undefined_traits[key] = Undefined
self.inputs.trait_set(trait_change_notify=False, **undefined_traits)
if force_run:
self._always_run = True
def _run_interface(self, runtime):
try:
import pandas as pd
except ImportError as e:
raise ImportError(
"This interface requires pandas " "(http://pandas.pydata.org/) to run."
) from e
try:
from filelock import SoftFileLock
self._have_lock = True
except ImportError:
from warnings import warn
warn(
(
"Python module filelock was not found: AddCSVRow will not be"
" thread-safe in multi-processor execution"
)
)
input_dict = {}
for key, val in list(self.inputs._outputs.items()):
# expand lists to several columns
if key == "trait_added" and val in self.inputs.copyable_trait_names():
continue
if isinstance(val, list):
for i, v in enumerate(val):
input_dict["%s_%d" % (key, i)] = v
else:
input_dict[key] = val
df = pd.DataFrame([input_dict])
if self._have_lock:
self._lock = SoftFileLock("%s.lock" % self.inputs.in_file)
# Acquire lock
self._lock.acquire()
if op.exists(self.inputs.in_file):
formerdf = pd.read_csv(self.inputs.in_file, index_col=0)
df = pd.concat([formerdf, df], ignore_index=True)
with open(self.inputs.in_file, "w") as f:
df.to_csv(f)
if self._have_lock:
self._lock.release()
return runtime
def _list_outputs(self):
outputs = self.output_spec().get()
outputs["csv_file"] = self.inputs.in_file
return outputs
def _outputs(self):
return self._add_output_traits(super(AddCSVRow, self)._outputs())
def _add_output_traits(self, base):
return base
class CalculateNormalizedMomentsInputSpec(TraitedSpec):
timeseries_file = File(
exists=True,
mandatory=True,
desc="Text file with timeseries in columns and timepoints in rows,\
whitespace separated",
)
moment = traits.Int(
mandatory=True,
desc="Define which moment should be calculated, 3 for skewness, 4 for\
kurtosis.",
)
class CalculateNormalizedMomentsOutputSpec(TraitedSpec):
moments = traits.List(traits.Float(), desc="Moments")
class CalculateNormalizedMoments(BaseInterface):
"""
Calculates moments of timeseries.
Example
-------
>>> from nipype.algorithms import misc
>>> skew = misc.CalculateNormalizedMoments()
>>> skew.inputs.moment = 3
>>> skew.inputs.timeseries_file = 'timeseries.txt'
>>> skew.run() # doctest: +SKIP
"""
input_spec = CalculateNormalizedMomentsInputSpec
output_spec = CalculateNormalizedMomentsOutputSpec
def _run_interface(self, runtime):
self._moments = calc_moments(self.inputs.timeseries_file, self.inputs.moment)
return runtime
def _list_outputs(self):
outputs = self.output_spec().get()
outputs["skewness"] = self._moments
return outputs
def calc_moments(timeseries_file, moment):
"""Returns nth moment (3 for skewness, 4 for kurtosis) of timeseries
(list of values; one per timeseries).
Keyword arguments:
timeseries_file -- text file with white space separated timepoints in rows
"""
import scipy.stats as stats
timeseries = np.genfromtxt(timeseries_file)
m2 = stats.moment(timeseries, 2, axis=0)
m3 = stats.moment(timeseries, moment, axis=0)
zero = m2 == 0
return np.where(zero, 0, m3 / m2 ** (moment / 2.0))
class AddNoiseInputSpec(TraitedSpec):
in_file = File(
exists=True,
mandatory=True,
desc="input image that will be corrupted with noise",
)
in_mask = File(
exists=True,
desc=("input mask, voxels outside this mask " "will be considered background"),
)
snr = traits.Float(10.0, desc="desired output SNR in dB", usedefault=True)
dist = traits.Enum(
"normal",
"rician",
usedefault=True,
mandatory=True,
desc=("desired noise distribution"),
)
bg_dist = traits.Enum(
"normal",
"rayleigh",
usedefault=True,
mandatory=True,
desc=("desired noise distribution, currently " "only normal is implemented"),
)
out_file = File(desc="desired output filename")
class AddNoiseOutputSpec(TraitedSpec):
out_file = File(exists=True, desc="corrupted image")
class AddNoise(BaseInterface):
"""
Corrupts with noise the input image.
Example
-------
>>> from nipype.algorithms.misc import AddNoise
>>> noise = AddNoise()
>>> noise.inputs.in_file = 'T1.nii'
>>> noise.inputs.in_mask = 'mask.nii'
>>> noise.snr = 30.0
>>> noise.run() # doctest: +SKIP
"""
input_spec = AddNoiseInputSpec
output_spec = AddNoiseOutputSpec
def _run_interface(self, runtime):
in_image = nb.load(self.inputs.in_file)
in_data = in_image.get_fdata()
snr = self.inputs.snr
if isdefined(self.inputs.in_mask):
in_mask = np.asanyarray(nb.load(self.inputs.in_mask).dataobj)
else:
in_mask = np.ones_like(in_data)
result = self.gen_noise(
in_data,
mask=in_mask,
snr_db=snr,
dist=self.inputs.dist,
bg_dist=self.inputs.bg_dist,
)
res_im = nb.Nifti1Image(result, in_image.affine, in_image.header)
res_im.to_filename(self._gen_output_filename())
return runtime
def _gen_output_filename(self):
if not isdefined(self.inputs.out_file):
_, base, ext = split_filename(self.inputs.in_file)
out_file = os.path.abspath("%s_SNR%03.2f%s" % (base, self.inputs.snr, ext))
else:
out_file = self.inputs.out_file
return out_file
def _list_outputs(self):
outputs = self.output_spec().get()
outputs["out_file"] = self._gen_output_filename()
return outputs
def gen_noise(self, image, mask=None, snr_db=10.0, dist="normal", bg_dist="normal"):
"""
Generates a copy of an image with a certain amount of
added gaussian noise (rayleigh for background in mask)
"""
from math import sqrt
snr = sqrt(np.power(10.0, snr_db / 10.0))
if mask is None:
mask = np.ones_like(image)
else:
mask[mask > 0] = 1
mask[mask < 1] = 0
if mask.ndim < image.ndim:
mask = np.rollaxis(np.array([mask] * image.shape[3]), 0, 4)
signal = image[mask > 0].reshape(-1)
if dist == "normal":
signal = signal - signal.mean()
sigma_n = sqrt(signal.var() / snr)
noise = np.random.normal(size=image.shape, scale=sigma_n)
if (np.any(mask == 0)) and (bg_dist == "rayleigh"):
bg_noise = np.random.rayleigh(size=image.shape, scale=sigma_n)
noise[mask == 0] = bg_noise[mask == 0]
im_noise = image + noise
elif dist == "rician":
sigma_n = signal.mean() / snr
n_1 = np.random.normal(size=image.shape, scale=sigma_n)
n_2 = np.random.normal(size=image.shape, scale=sigma_n)
stde_1 = n_1 / sqrt(2.0)
stde_2 = n_2 / sqrt(2.0)
im_noise = np.sqrt((image + stde_1) ** 2 + (stde_2) ** 2)
else:
raise NotImplementedError(
("Only normal and rician distributions " "are supported")
)
return im_noise
class NormalizeProbabilityMapSetInputSpec(TraitedSpec):
in_files = InputMultiPath(
File(exists=True, mandatory=True, desc="The tpms to be normalized")
)
in_mask = File(exists=True, desc="Masked voxels must sum up 1.0, 0.0 otherwise.")
class NormalizeProbabilityMapSetOutputSpec(TraitedSpec):
out_files = OutputMultiPath(File(exists=True), desc="normalized maps")
class NormalizeProbabilityMapSet(BaseInterface):
"""
Returns the input tissue probability maps (tpms, aka volume fractions).
The tissue probability maps are normalized to sum up 1.0 at each voxel within the mask.
.. note:: Please recall this is not a spatial normalization algorithm
Example
-------
>>> from nipype.algorithms import misc
>>> normalize = misc.NormalizeProbabilityMapSet()
>>> normalize.inputs.in_files = [ 'tpm_00.nii.gz', 'tpm_01.nii.gz', \
'tpm_02.nii.gz' ]
>>> normalize.inputs.in_mask = 'tpms_msk.nii.gz'
>>> normalize.run() # doctest: +SKIP
"""
input_spec = NormalizeProbabilityMapSetInputSpec
output_spec = NormalizeProbabilityMapSetOutputSpec
def _run_interface(self, runtime):
mask = None
if isdefined(self.inputs.in_mask):
mask = self.inputs.in_mask
self._out_filenames = normalize_tpms(self.inputs.in_files, mask)
return runtime
def _list_outputs(self):
outputs = self.output_spec().get()
outputs["out_files"] = self._out_filenames
return outputs
class SplitROIsInputSpec(TraitedSpec):
in_file = File(exists=True, mandatory=True, desc="file to be splitted")
in_mask = File(exists=True, desc="only process files inside mask")
roi_size = traits.Tuple(traits.Int, traits.Int, traits.Int, desc="desired ROI size")
class SplitROIsOutputSpec(TraitedSpec):
out_files = OutputMultiPath(File(exists=True), desc="the resulting ROIs")
out_masks = OutputMultiPath(
File(exists=True), desc="a mask indicating valid values"
)
out_index = OutputMultiPath(
File(exists=True), desc="arrays keeping original locations"
)
class SplitROIs(BaseInterface):
"""
Splits a 3D image in small chunks to enable parallel processing.
ROIs keep time series structure in 4D images.
Example
-------
>>> from nipype.algorithms import misc
>>> rois = misc.SplitROIs()
>>> rois.inputs.in_file = 'diffusion.nii'
>>> rois.inputs.in_mask = 'mask.nii'
>>> rois.run() # doctest: +SKIP
"""
input_spec = SplitROIsInputSpec
output_spec = SplitROIsOutputSpec
def _run_interface(self, runtime):
mask = None
roisize = None
self._outnames = {}
if isdefined(self.inputs.in_mask):
mask = self.inputs.in_mask
if isdefined(self.inputs.roi_size):
roisize = self.inputs.roi_size
res = split_rois(self.inputs.in_file, mask, roisize)
self._outnames["out_files"] = res[0]
self._outnames["out_masks"] = res[1]
self._outnames["out_index"] = res[2]
return runtime
def _list_outputs(self):
outputs = self.output_spec().get()
for k, v in list(self._outnames.items()):
outputs[k] = v
return outputs
class MergeROIsInputSpec(TraitedSpec):
in_files = InputMultiPath(
File(exists=True, mandatory=True, desc="files to be re-merged")
)
in_index = InputMultiPath(
File(exists=True, mandatory=True), desc="array keeping original locations"
)
in_reference = File(exists=True, desc="reference file")
class MergeROIsOutputSpec(TraitedSpec):
merged_file = File(exists=True, desc="the recomposed file")
class MergeROIs(BaseInterface):
"""
Splits a 3D image in small chunks to enable parallel processing.
ROIs keep time series structure in 4D images.
Example
-------
>>> from nipype.algorithms import misc
>>> rois = misc.MergeROIs()
>>> rois.inputs.in_files = ['roi%02d.nii' % i for i in range(1, 6)]
>>> rois.inputs.in_reference = 'mask.nii'
>>> rois.inputs.in_index = ['roi%02d_idx.npz' % i for i in range(1, 6)]
>>> rois.run() # doctest: +SKIP
"""
input_spec = MergeROIsInputSpec
output_spec = MergeROIsOutputSpec
def _run_interface(self, runtime):
res = merge_rois(
self.inputs.in_files, self.inputs.in_index, self.inputs.in_reference
)
self._merged = res
return runtime
def _list_outputs(self):
outputs = self.output_spec().get()
outputs["merged_file"] = self._merged
return outputs
def normalize_tpms(in_files, in_mask=None, out_files=None):
"""
Returns the input tissue probability maps (tpms, aka volume fractions)
normalized to sum up 1.0 at each voxel within the mask.
"""
import nibabel as nb
import numpy as np
import os.path as op
in_files = np.atleast_1d(in_files).tolist()
if out_files is None:
out_files = []
if len(out_files) != len(in_files):
for i, finname in enumerate(in_files):
fname, fext = op.splitext(op.basename(finname))
if fext == ".gz":
fname, fext2 = op.splitext(fname)
fext = fext2 + fext
out_file = op.abspath("%s_norm_%02d%s" % (fname, i, fext))
out_files += [out_file]
imgs = [nb.load(fim) for fim in in_files]
if len(in_files) == 1:
img_data = imgs[0].get_fdata(dtype=np.float32)
img_data[img_data > 0.0] = 1.0
hdr = imgs[0].header.copy()
hdr.set_data_dtype(np.float32)
nb.save(nb.Nifti1Image(img_data, imgs[0].affine, hdr), out_files[0])
return out_files[0]
img_data = np.stack(
[im.get_fdata(caching="unchanged", dtype=np.float32) for im in imgs]
)
# img_data[img_data>1.0] = 1.0
img_data[img_data < 0.0] = 0.0
weights = np.sum(img_data, axis=0)
msk = np.ones(imgs[0].shape)
msk[weights <= 0] = 0
if in_mask is not None:
msk = np.asanyarray(nb.load(in_mask).dataobj)
msk[msk <= 0] = 0
msk[msk > 0] = 1
msk = np.ma.masked_equal(msk, 0)
for i, out_file in enumerate(out_files):
data = np.ma.masked_equal(img_data[i], 0)
probmap = data / weights
hdr = imgs[i].header.copy()
hdr.set_data_dtype("float32")
nb.save(
nb.Nifti1Image(probmap.astype(np.float32), imgs[i].affine, hdr), out_file
)
return out_files
def split_rois(in_file, mask=None, roishape=None):
"""
Splits an image in ROIs for parallel processing
"""
import nibabel as nb
import numpy as np
from math import sqrt, ceil
import os.path as op
if roishape is None:
roishape = (10, 10, 1)
im = nb.load(in_file)
imshape = im.shape
dshape = imshape[:3]
nvols = imshape[-1]
roisize = roishape[0] * roishape[1] * roishape[2]
droishape = (roishape[0], roishape[1], roishape[2], nvols)
if mask is not None:
mask = np.asanyarray(nb.load(mask).dataobj)
mask[mask > 0] = 1
mask[mask < 1] = 0
else:
mask = np.ones(dshape)
mask = mask.reshape(-1).astype(np.uint8)
nzels = np.nonzero(mask)
els = np.sum(mask)
nrois = int(ceil(els / float(roisize)))
data = np.asanyarray(im.dataobj).reshape((mask.size, -1))
data = np.squeeze(data.take(nzels, axis=0))
nvols = data.shape[-1]
roidefname = op.abspath("onesmask.nii.gz")
nb.Nifti1Image(np.ones(roishape, dtype=np.uint8), None, None).to_filename(
roidefname
)
out_files = []
out_mask = []
out_idxs = []
for i in range(nrois):
first = i * roisize
last = (i + 1) * roisize
fill = 0
if last > els:
fill = last - els
last = els
droi = data[first:last, ...]
iname = op.abspath("roi%010d_idx" % i)
out_idxs.append(iname + ".npz")
np.savez(iname, (nzels[0][first:last],))
if fill > 0:
droi = np.vstack(
(droi, np.zeros((int(fill), int(nvols)), dtype=np.float32))
)
partialmsk = np.ones((roisize,), dtype=np.uint8)
partialmsk[-int(fill) :] = 0
partname = op.abspath("partialmask.nii.gz")
nb.Nifti1Image(partialmsk.reshape(roishape), None, None).to_filename(
partname
)
out_mask.append(partname)
else:
out_mask.append(roidefname)
fname = op.abspath("roi%010d.nii.gz" % i)
nb.Nifti1Image(droi.reshape(droishape), None, None).to_filename(fname)
out_files.append(fname)
return out_files, out_mask, out_idxs
def merge_rois(in_files, in_idxs, in_ref, dtype=None, out_file=None):
"""
Re-builds an image resulting from a parallelized processing
"""
import nibabel as nb
import numpy as np
import os.path as op
import subprocess as sp
if out_file is None:
out_file = op.abspath("merged.nii.gz")
if dtype is None:
dtype = np.float32
# if file is compressed, uncompress using os
# to avoid memory errors
if op.splitext(in_ref)[1] == ".gz":
try:
iflogger.info("uncompress %s", in_ref)
sp.check_call(["gunzip", in_ref], stdout=sp.PIPE, shell=True)
in_ref = op.splitext(in_ref)[0]
except:
pass
ref = nb.load(in_ref)
aff = ref.affine
hdr = ref.header.copy()
rsh = ref.shape
del ref
npix = rsh[0] * rsh[1] * rsh[2]
fcimg = nb.load(in_files[0])
if len(fcimg.shape) == 4:
ndirs = fcimg.shape[-1]
else:
ndirs = 1
newshape = (rsh[0], rsh[1], rsh[2], ndirs)
hdr.set_data_dtype(dtype)
hdr.set_xyzt_units("mm", "sec")
if ndirs < 300:
data = np.zeros((npix, ndirs), dtype=dtype)
for cname, iname in zip(in_files, in_idxs):
f = np.load(iname)
idxs = np.squeeze(f["arr_0"])
cdata = np.asanyarray(nb.load(cname).dataobj).reshape(-1, ndirs)
nels = len(idxs)
idata = (idxs,)
try:
data[idata, ...] = cdata[0:nels, ...]
except:
print(
(
"Consistency between indexes and chunks was "
"lost: data=%s, chunk=%s"
)
% (str(data.shape), str(cdata.shape))
)
raise
nb.Nifti1Image(data.reshape(newshape), aff, hdr).to_filename(out_file)
else:
hdr.set_data_shape(rsh[:3])
nii = []
for d in range(ndirs):
fname = op.abspath("vol%06d.nii" % d)
nb.Nifti1Image(np.zeros(rsh[:3]), aff, hdr).to_filename(fname)
nii.append(fname)
for cname, iname in zip(in_files, in_idxs):
f = np.load(iname)
idxs = np.squeeze(f["arr_0"])
for d, fname in enumerate(nii):
data = np.asanyarray(nb.load(fname).dataobj).reshape(-1)
cdata = nb.load(cname).dataobj[..., d].reshape(-1)
nels = len(idxs)
idata = (idxs,)
data[idata] = cdata[0:nels]
nb.Nifti1Image(data.reshape(rsh[:3]), aff, hdr).to_filename(fname)
imgs = [nb.load(im) for im in nii]
allim = nb.concat_images(imgs)
allim.to_filename(out_file)
return out_file
class CalculateMedianInputSpec(BaseInterfaceInputSpec):
in_files = InputMultiPath(
File(
exists=True,
mandatory=True,
desc="One or more realigned Nifti 4D timeseries",
)
)
median_file = traits.Str(desc="Filename prefix to store median images")
median_per_file = traits.Bool(
False, usedefault=True, desc="Calculate a median file for each Nifti"
)
class CalculateMedianOutputSpec(TraitedSpec):
median_files = OutputMultiPath(File(exists=True), desc="One or more median images")
class CalculateMedian(BaseInterface):
"""
Computes an average of the median across one or more 4D Nifti timeseries
Example
-------
>>> from nipype.algorithms.misc import CalculateMedian
>>> mean = CalculateMedian()
>>> mean.inputs.in_files = 'functional.nii'
>>> mean.run() # doctest: +SKIP
"""
input_spec = CalculateMedianInputSpec
output_spec = CalculateMedianOutputSpec
def __init__(self, *args, **kwargs):
super(CalculateMedian, self).__init__(*args, **kwargs)
self._median_files = []
def _gen_fname(self, suffix, idx=None, ext=None):
if idx:
in_file = self.inputs.in_files[idx]
else:
if isinstance(self.inputs.in_files, list):
in_file = self.inputs.in_files[0]
else:
in_file = self.inputs.in_files
fname, in_ext = op.splitext(op.basename(in_file))
if in_ext == ".gz":
fname, in_ext2 = op.splitext(fname)
in_ext = in_ext2 + in_ext
if ext is None:
ext = in_ext
if ext.startswith("."):
ext = ext[1:]
if self.inputs.median_file:
outname = self.inputs.median_file
else:
outname = "{}_{}".format(fname, suffix)
if idx:
outname += str(idx)
return op.abspath("{}.{}".format(outname, ext))
def _run_interface(self, runtime):
total = None
self._median_files = []
for idx, fname in enumerate(ensure_list(self.inputs.in_files)):
img = nb.load(fname)
data = np.median(img.get_fdata(), axis=3)
if self.inputs.median_per_file:
self._median_files.append(self._write_nifti(img, data, idx))
else:
if total is None:
total = data
else:
total += data
if not self.inputs.median_per_file:
self._median_files.append(self._write_nifti(img, total, idx))
return runtime
def _list_outputs(self):
outputs = self._outputs().get()
outputs["median_files"] = self._median_files
return outputs
def _write_nifti(self, img, data, idx, suffix="median"):
if self.inputs.median_per_file:
median_img = nb.Nifti1Image(data, img.affine, img.header)
filename = self._gen_fname(suffix, idx=idx)
else:
median_img = nb.Nifti1Image(data / (idx + 1), img.affine, img.header)
filename = self._gen_fname(suffix)
median_img.to_filename(filename)
return filename
# Deprecated interfaces ------------------------------------------------------
class Distance(nam.Distance):
"""Calculates distance between two volumes.
.. deprecated:: 0.10.0
Use :py:class:`nipype.algorithms.metrics.Distance` instead.
"""
def __init__(self, **inputs):
super(nam.Distance, self).__init__(**inputs)
warnings.warn(
(
"This interface has been deprecated since 0.10.0,"
" please use nipype.algorithms.metrics.Distance"
),
DeprecationWarning,
)
class Overlap(nam.Overlap):
"""Calculates various overlap measures between two maps.
.. deprecated:: 0.10.0
Use :py:class:`nipype.algorithms.metrics.Overlap` instead.
"""
def __init__(self, **inputs):
super(nam.Overlap, self).__init__(**inputs)
warnings.warn(
(
"This interface has been deprecated since 0.10.0,"
" please use nipype.algorithms.metrics.Overlap"
),
DeprecationWarning,
)
class FuzzyOverlap(nam.FuzzyOverlap):
"""Calculates various overlap measures between two maps, using a fuzzy
definition.
.. deprecated:: 0.10.0
Use :py:class:`nipype.algorithms.metrics.FuzzyOverlap` instead.
"""
def __init__(self, **inputs):
super(nam.FuzzyOverlap, self).__init__(**inputs)
warnings.warn(
(
"This interface has been deprecated since 0.10.0,"
" please use nipype.algorithms.metrics.FuzzyOverlap"
),
DeprecationWarning,
)
class TSNR(confounds.TSNR):
"""
.. deprecated:: 0.12.1
Use :py:class:`nipype.algorithms.confounds.TSNR` instead
"""
def __init__(self, **inputs):
super(confounds.TSNR, self).__init__(**inputs)
warnings.warn(
(
"This interface has been moved since 0.12.0,"
" please use nipype.algorithms.confounds.TSNR"
),
UserWarning,
)