RF/FIX: Scale all phase maps to `[0, 2pi]` range #88

effigies · 2020-02-13T14:54:04Z

Final decision: Instead of making [-pi, pi] a special case, all data will be naively scaled to [0, 2pi]. This ensures that the source data intervals are respected, but may result in a central mode slightly off of pi.

We judge that the simplicity and consistency outweighs any small adjustments and potential for bugs resulting from a mode-based transform with separate treatment of the left and right tail.

I have some Philips Achieva 2-echo phase/magnitude fieldmaps that come out of dcm2niix with a range of approximately (-pi, pi).

This adjusts au2rads to detect that case and simply wraps the negative values to (pi, 2pi).

This may resolve #83.

cc @ins0mniac2

oesteban · 2020-02-13T15:25:19Z

sdcflows/interfaces/fmap.py

-    data -= mode(data, axis=None)[0][0]
+    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


Why only negatives?

Suggested change

data[data < 0] += 2 * np.pi

data -= data.min()

Because this function sets the range to (0, 2pi), not (pi, 3pi). And it seems unwise to go precisely counter-phase by adding np.pi to everything.

Sorry, I had the suggestion above sign-flipped. I've fixed it.

This is going to wrap and scale only negatives - you don't want either (and even less only for the lower half of the histogram).

If I understand correctly, the initial range goes from -pi to pi. Adding pi will move the range to 0-2pi

So does fieldmap correction treat (θ±π) the same? Shifting by half phase feels problematic.

To be clear, this is equivalent to:

data = (data + 2 * np.pi) % (2 * np.pi)

But without bothering add, divide and round things in the (0, pi) range. It doesn't feel controversial to me. data += np.pi does.

So does fieldmap correction treat (θ±π) the same? Shifting by half phase feels problematic.

FSL FUGUE theoretically only unwraps phase maps in the (0, 2π) range. It seems to work with (-π, π), though.

I don't know why data -= np.pi is controversial - it just moves all values from (-π, π) to (0, 2π).

The transformation you propose will wrap negative phase values (which you definitely don't want to do).

Normally, a phase maps' histogram looks like a centered pyramid. Doing this, we are selecting the left side of the pyramid, shifting it to overlap the right side of the pyramid and finally checking no values have gone over 2pi. That will introduce another wrapping that FUGUE might or might not pick up.

EDIT: The shifted intensities will not overlap, now that I look into this more carefully, but the final distribution will be like a V, with the notch at pi.

Let's talk in the meeting.

oesteban · 2020-02-13T16:55:47Z

sdcflows/interfaces/fmap.py

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


~~This (the old code) looks problematic - we shift the left side of the pyramid first, hence making it overlap with the right hand side. Operations should be done the other way around.~~

Okay, sorry for the noise. There is no shifting yet.

effigies · 2020-02-13T17:15:41Z

Here's my histogram:

effigies · 2020-02-13T18:32:23Z

Results of live discussion:

Replace mode-centering and separately rescaling positive and negative to simply scaling min/max to 0/2pi. This would affect fieldmaps from all scanners.

@mattcieslak @Aksoo Your input here would be appreciated.

effigies

Just a couple explainers before you get to asking. :-)

effigies · 2020-02-13T22:17:07Z

sdcflows/interfaces/fmap.py

@@ -663,28 +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


I'm extremely comfortable using double memory for one volume to keep things deterministic.

effigies · 2020-02-13T22:18:43Z

sdcflows/interfaces/fmap.py


    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)


If you're keeping the affines in the header, there's no need to pass an affine in to override.

sdcflows/interfaces/fmap.py

kimsin98 · 2020-02-14T02:27:23Z

Replace mode-centering and separately rescaling positive and negative to simply scaling min/max to 0/2pi. This would affect fieldmaps from all scanners.

I also had doubts about mode centering. This looks simpler and cleaner.

effigies · 2020-02-14T21:55:42Z

@mgxd I see the CI has gotten clogged. Naturally.

If it makes life easier, I can cancel this and master and just merge this to only have one job in queue.

mgxd · 2020-02-14T22:00:45Z

Doesn't the current master CI needs to pass to upload the new version to PyPI?

I'm fine with cutting corners, but would like to see a green check here prior to merging just for my peace of mind

effigies · 2020-02-14T22:03:33Z

Okay. I'll leave this one up. Still cancelled master since I rebased this on master.

We'll push to PyPI with the tag: https://circleci.com/workflow-run/52c4f144-141b-41ed-acad-e218e1b65f2a

pull-assistant · 2020-02-15T01:19:17Z

Best reviewed: with all changes

Optimal code review plan (2 commits squashed)

FIX: Adjust au2rads ... ... RF: Always simply re...

Squashed 2 commits:

FIX: Adjust au2rads to handle input data with 2pi range

RF: Always simply rescale to [0,2pi]

Powered by Pull Assistant. Last update ee0fe99 ... 410365a. Read the comment docs.

effigies requested review from mattcieslak and oesteban February 13, 2020 14:54

oesteban reviewed Feb 13, 2020

View reviewed changes

effigies commented Feb 13, 2020

View reviewed changes

oesteban approved these changes Feb 13, 2020

View reviewed changes

sdcflows/interfaces/fmap.py Show resolved Hide resolved

effigies force-pushed the fix/rads2rads branch from 543786e to b16156c Compare February 14, 2020 21:12

effigies added 2 commits February 14, 2020 16:43

FIX: Adjust au2rads to handle input data with 2pi range

ee0fe99

RF: Always simply rescale to [0,2pi]

410365a

effigies force-pushed the fix/rads2rads branch from b16156c to 410365a Compare February 14, 2020 21:45

effigies changed the title ~~FIX: Adjust au2rads to handle input data with 2pi range~~ RF/FIX: Scale all phase maps to [0, 2pi] range Feb 15, 2020

effigies merged commit 28b084f into nipreps:master Feb 15, 2020

effigies deleted the fix/rads2rads branch February 15, 2020 01:22

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

RF/FIX: Scale all phase maps to `[0, 2pi]` range #88

RF/FIX: Scale all phase maps to `[0, 2pi]` range #88

effigies commented Feb 13, 2020 •

edited

Loading

oesteban Feb 13, 2020 •

edited

Loading

effigies Feb 13, 2020

oesteban Feb 13, 2020

effigies Feb 13, 2020

effigies Feb 13, 2020

oesteban Feb 13, 2020 •

edited

Loading

effigies Feb 13, 2020

oesteban Feb 13, 2020 •

edited

Loading

effigies commented Feb 13, 2020

effigies commented Feb 13, 2020

effigies left a comment

effigies Feb 13, 2020

effigies Feb 13, 2020

kimsin98 commented Feb 14, 2020

effigies commented Feb 14, 2020

mgxd commented Feb 14, 2020 •

edited

Loading

effigies commented Feb 14, 2020

pull-assistant bot commented Feb 15, 2020

RF/FIX: Scale all phase maps to [0, 2pi] range #88

RF/FIX: Scale all phase maps to [0, 2pi] range #88

Conversation

effigies commented Feb 13, 2020 • edited Loading

oesteban Feb 13, 2020 • edited Loading

Choose a reason for hiding this comment

effigies Feb 13, 2020

Choose a reason for hiding this comment

oesteban Feb 13, 2020

Choose a reason for hiding this comment

effigies Feb 13, 2020

Choose a reason for hiding this comment

effigies Feb 13, 2020

Choose a reason for hiding this comment

oesteban Feb 13, 2020 • edited Loading

Choose a reason for hiding this comment

effigies Feb 13, 2020

Choose a reason for hiding this comment

oesteban Feb 13, 2020 • edited Loading

Choose a reason for hiding this comment

effigies commented Feb 13, 2020

effigies commented Feb 13, 2020

effigies left a comment

Choose a reason for hiding this comment

effigies Feb 13, 2020

Choose a reason for hiding this comment

effigies Feb 13, 2020

Choose a reason for hiding this comment

kimsin98 commented Feb 14, 2020

effigies commented Feb 14, 2020

mgxd commented Feb 14, 2020 • edited Loading

effigies commented Feb 14, 2020

pull-assistant bot commented Feb 15, 2020

RF/FIX: Scale all phase maps to `[0, 2pi]` range #88

RF/FIX: Scale all phase maps to `[0, 2pi]` range #88

effigies commented Feb 13, 2020 •

edited

Loading

oesteban Feb 13, 2020 •

edited

Loading

oesteban Feb 13, 2020 •

edited

Loading

oesteban Feb 13, 2020 •

edited

Loading

mgxd commented Feb 14, 2020 •

edited

Loading