From 2c3e3a3f4998ff5131a0a2a02d047dd78340e41d Mon Sep 17 00:00:00 2001
From: Chris Markiewicz <markiewicz@stanford.edu>
Date: Fri, 14 Feb 2020 12:57:48 -0500
Subject: [PATCH] FIX: Center phase maps around second mode when first is
 minimum

---
 sdcflows/interfaces/fmap.py | 21 ++++++++++++++++-----
 1 file changed, 16 insertions(+), 5 deletions(-)

diff --git a/sdcflows/interfaces/fmap.py b/sdcflows/interfaces/fmap.py
index 3d347366af..f60245ae05 100644
--- a/sdcflows/interfaces/fmap.py
+++ b/sdcflows/interfaces/fmap.py
@@ -661,19 +661,30 @@ def _delta_te(in_values, te1=None, te2=None):
 
 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')
     hdr = im.header.copy()
 
-    # First center data around 0.0.
-    data -= mode(data, axis=None)[0][0]
+    vals, counts = np.unique(data, return_counts=True)
+    modes = vals[np.argsort(counts)[::-1]]
+
+    # Second mode (idx 1) if first node (idx 0) is minimum, else first
+    idx = int(modes[0] == data.min())
+
+    # Center data around 0.0
+    data -= modes[idx]
+
+    # Provide a less opaque error if we still can't figure it out
+    neg = data < 0
+    pos = data > 0
+    if not (neg.any() and pos.any())
+        raise ValueError("Could not find an appropriate mode to center phase values around")
 
     # Scale lower tail
-    data[data < 0] = - np.pi * data[data < 0] / data[data < 0].min()
+    data[neg] = - np.pi * data[neg] / data[neg].min()
 
     # Scale upper tail
-    data[data > 0] = np.pi * data[data > 0] / data[data > 0].max()
+    data[pos] = np.pi * data[pos] / data[pos].max()
 
     # Offset to 0 - 2pi
     data += np.pi