# 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 <nipreps@gmail.com>## 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/#"""Utilities to manipulate phase and phase difference maps."""
[docs]defau2rads(in_file,newpath=None):"""Convert the input phase map in arbitrary units (a.u.) to rads."""importnumpyasnpimportnibabelasnbfromnipype.utils.filemanipimportfname_presuffixim=nb.load(in_file)data=im.get_fdata(caching="unchanged")# Read as float64 for safetyhdr=im.header.copy()# Rescale to [0, 2*pi]data=(data-data.min())*(2*np.pi/(data.max()-data.min()))# Round to float32 and clipdata=np.clip(np.float32(data),0.0,2*np.pi)hdr.set_data_dtype(np.float32)hdr.set_xyzt_units("mm")out_file=fname_presuffix(str(in_file),suffix="_rads",newpath=newpath)nb.Nifti1Image(data,None,hdr).to_filename(out_file)returnout_file
[docs]defsubtract_phases(in_phases,in_meta,newpath=None):"""Calculate the phase-difference map, given two input phase maps."""importnumpyasnpimportnibabelasnbfromnipype.utils.filemanipimportfname_presuffixecho_times=tuple([m.pop("EchoTime",None)forminin_meta])ifecho_times[0]>echo_times[1]:in_phases=(in_phases[1],in_phases[0])in_meta=(in_meta[1],in_meta[0])echo_times=(echo_times[1],echo_times[0])in_phases_nii=[nb.load(ph)forphinin_phases]sub_data=in_phases_nii[1].get_fdata(dtype="float32")-in_phases_nii[0].get_fdata(dtype="float32")# wrap negative radians back to [0, 2pi]sub_data[sub_data<0]+=2*np.pisub_data=np.clip(sub_data,0.0,2*np.pi)new_meta=in_meta[1].copy()new_meta.update(in_meta[0])new_meta["EchoTime1"]=echo_times[0]new_meta["EchoTime2"]=echo_times[1]hdr=in_phases_nii[0].header.copy()hdr.set_data_dtype(np.float32)hdr.set_xyzt_units("mm")nii=nb.Nifti1Image(sub_data,in_phases_nii[0].affine,hdr)out_phdiff=fname_presuffix(in_phases[0],suffix="_phdiff",newpath=newpath)nii.to_filename(out_phdiff)returnout_phdiff,new_meta
[docs]defphdiff2fmap(in_file,delta_te,newpath=None):"""Convert the input phase-difference map into a *fieldmap* in Hz."""importnumpyasnpimportnibabelasnbfromnipype.utils.filemanipimportfname_presuffixout_file=fname_presuffix(str(in_file),suffix="_fmap",newpath=newpath)image=nb.load(in_file)data=image.get_fdata(dtype="float32")/(2.0*np.pi*delta_te)nii=nb.Nifti1Image(data,image.affine,image.header)nii.set_data_dtype(np.float32)nii.to_filename(out_file)returnout_file