Skip to content

Commit

Permalink
RF: Always simply rescale to [0,2pi]
Browse files Browse the repository at this point in the history
  • Loading branch information
effigies committed Feb 13, 2020
1 parent 420a410 commit 543786e
Showing 1 changed file with 6 additions and 19 deletions.
25 changes: 6 additions & 19 deletions sdcflows/interfaces/fmap.py
Original file line number Diff line number Diff line change
Expand Up @@ -663,32 +663,19 @@ def au2rads(in_file, newpath=None):
"""Convert the input phase difference map in arbitrary units (a.u.) to rads."""
from scipy.stats import mode
im = nb.load(in_file)
data = im.get_fdata(dtype='float32')
data = im.get_fdata(caching='unchanged') # Read as float64 for safety
hdr = im.header.copy()

if np.allclose((data.min(), data.max()), (-np.pi, np.pi), atol=0.01):
# Already in rads, but wrap to 0 - 2pi
data[data < 0] += 2 * np.pi
else:
# First center data around 0.0.
data -= mode(data, axis=None)[0][0]
# Rescale to [0, 2*pi]
data = (data - data.min()) * (2 * np.pi / (data.max() - data.min()))

# Scale lower tail
data[data < 0] = - np.pi * data[data < 0] / data[data < 0].min()

# Scale upper tail
data[data > 0] = np.pi * data[data > 0] / data[data > 0].max()

# Offset to 0 - 2pi
data += np.pi

# Clip
data = np.clip(data, 0.0, 2 * np.pi)
# Round to float32 and clip
data = 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(in_file, suffix='_rads', newpath=newpath)
nb.Nifti1Image(data, im.affine, hdr).to_filename(out_file)
nb.Nifti1Image(data, None, hdr).to_filename(out_file)
return out_file


Expand Down

0 comments on commit 543786e

Please sign in to comment.