ENH: Stop using siemens2rads from old nipype workflows #50
Conversation
A new interface is proposed, to be used in phase-difference based workflows. @Aksoo identified a potential issue with siemens2rads (nipreps#30 (comment)), although the re-centering did not seem to make a big difference. This implementation of the rescaling uses percentiles to robustly calculate the range of the input image (as recommended in the original paper about FSL PRELUDE), and makes sure the output data type is adequate (float32).
Co-Authored-By: Chris Markiewicz <[email protected]>
|
Okay, now I remember why we were not using the 0-6.28 range: it needs recentering (which is not the case for the 2 phases case). I'll add a recentering node to this PR. |
FSL PRELUDE will return some values unwrapped below zero, making it hard to differenciate these pixels from the mask. This commit addresses that particular issue.
|
Trimming wrapped phase images at 2nd and 98th percentiles seems odd, since the values outside that range are not really outliers in a wrapped image. They are the connecting values necessary for smooth unwrapping. In the original PRELUDE paper, the percentile trimming is done on the magnitude image to derive a mask if no mask is specified by the user. |
The intent of this was to deal correctly with phase images cast to float datatypes. Since it doesn't seem to solve any problem, and the particular change goes beyond the scope of this PR, I think you are right I should remove it 👍(done in the latest commit). |
| """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') |
There was a problem hiding this comment.
Now I'm thinking about this again, I wonder if it makes more sense to make it float64 for the calculation, given that it's only one volume. This would ensure that the rounding error of casting to float32 on write (the set_data_dtype below will still take care of that) dominates the arithmetic rounding errors.
There was a problem hiding this comment.
Feels like splitting hairs... especially for images that originally were int16, DYT?
There was a problem hiding this comment.
Yeah, probably. But it's cheap and reduces accumulated error. That accumulated error is almost certainly in the noise, but the extra guard bits also give you lots of room to have suboptimal order of operations without changing the output relative to an optimal (smallest induced error) algorithm, which feels useful for reproducibility/comparability with other tools.
| This interface is equivalent to running the following steps: | ||
| #. Convert from rad to rad/s | ||
| (``niflow.nipype1.workflows.dmri.fsl.utils.rads2radsec``) | ||
| #. FUGUE execution: fsl.FUGUE(save_fmap=True) |
There was a problem hiding this comment.
| #. FUGUE execution: fsl.FUGUE(save_fmap=True) | |
| #. FUGUE execution: fsl.FUGUE(save_fmap=True). |
Putting together the lessons learned in nipreps#30, leveraging nipreps#52 and nipreps#53 (unfolded from nipreps#30 too), and utilizing nipreps#50 and nipreps#51, this workflow adds the phase difference map calculation, considering it one use-case of the general phase-difference fieldmap workflow. On top of this PR, we can continue the discussions held in nipreps#30. Probably, we will want to address nipreps#23 the first - the magnitude segmentation is sometimes really bad (e.g. see the phase1/2 unit test). Another discussion arisen in nipreps#30 is the spatial smoothing of the fieldmap (nipreps#22). Finally, the plan is to revise this implementation and determine whether the subtraction should happen before or after PRELUDE, and whether the arctan2 route is more interesting.
Putting together the lessons learned in nipreps#30, leveraging nipreps#52 and nipreps#53 (unfolded from nipreps#30 too), and utilizing nipreps#50 and nipreps#51, this workflow adds the phase difference map calculation, considering it one use-case of the general phase-difference fieldmap workflow. On top of this PR, we can continue the discussions held in nipreps#30. Probably, we will want to address nipreps#23 the first - the magnitude segmentation is sometimes really bad (e.g. see the phase1/2 unit test). Another discussion arisen in nipreps#30 is the spatial smoothing of the fieldmap (nipreps#22). Finally, the plan is to revise this implementation and determine whether the subtraction should happen before or after PRELUDE, and whether the arctan2 route is more interesting.
A new interface is proposed, to be used in phase-difference based
workflows.
@Aksoo identified a potential issue with siemens2rads
(#30 (comment)),
although the re-centering did not seem to make a big difference.
This implementation of the rescaling uses percentiles to robustly
calculate the range of the input image (as recommended in the original
paper about FSL PRELUDE), and makes sure the output data type is
adequate (float32).