Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[ENH] workflow for calculating phasediff from phase images #30

Closed
wants to merge 134 commits into from
Closed

[ENH] workflow for calculating phasediff from phase images #30

wants to merge 134 commits into from

Conversation

mattcieslak
Copy link
Collaborator

@mattcieslak mattcieslak commented Aug 17, 2019

This PR covers

  • Handling phase1/phase2 exactly as described in the FSL documentation
  • Adding subworkflows for processing magnitude images for phdiff, phase1/phase2, fmap
  • Adding subworkflows for post-processing fieldmaps (denoise, demean, bspline, etc)
  • Ensuring that phdiff, phase1/phase2 and fmap uniformly use these subworkflows

A couple questions:

  • What units should the fieldmaps be in when they go through the postprocessing workflow?
  • Should I add a report in the postprocessing workflow that plots the fieldmap in Hz like you did for pepolar?

@mattcieslak
Copy link
Collaborator Author

@oesteban do you have an example of a GE Fieldmap file? What is the typical input for FieldEnhance?

@mattcieslak mattcieslak changed the title [WIP] workflow for calculating phasediff from phase images [ENH] workflow for calculating phasediff from phase images Aug 20, 2019
@mattcieslak mattcieslak requested a review from oesteban October 2, 2019 18:17
Copy link
Member

@oesteban oesteban left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks very much for this! I'm midway through the PR. Left a few comments to get you going, most of them worried with focusing the PR on the target and asking to address other issues in separate PRs.

That said, this is looking good. I'll review the remainder tomorrow.

.gitignore Show resolved Hide resolved
sdcflows/interfaces/fmap.py Outdated Show resolved Hide resolved
sdcflows/interfaces/fmap.py Outdated Show resolved Hide resolved
sdcflows/interfaces/fmap.py Outdated Show resolved Hide resolved
sdcflows/interfaces/fmap.py Outdated Show resolved Hide resolved
sdcflows/workflows/b0.py Outdated Show resolved Hide resolved
sdcflows/workflows/b0.py Outdated Show resolved Hide resolved
sdcflows/workflows/b0.py Outdated Show resolved Hide resolved
sdcflows/workflows/base.py Show resolved Hide resolved
sdcflows/workflows/fmap.py Show resolved Hide resolved
Copy link
Collaborator Author

@mattcieslak mattcieslak left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'll open the bold renaming in a separate PR

@oesteban
Copy link
Member

Getting excited after seeing https://266-189236569-gh.circle-artifacts.com/0/tmp/tests/sdcflows/sub-1/fmap/sub-1_acq-v2_phase1.svg

I guess the outstanding question at this point is #30 (comment)

@mattcieslak
Copy link
Collaborator Author

Something in the changes in the last few months broke the workflow. The phasediff and phase1/phase2 used to produce fieldmaps that were within 0.1Hz of each other in ds001600.

wf.run()


def test_phases_workflow(bids_layouts, tmpdir, output_path):
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

duplicate function

@kimsin98
Copy link
Contributor

kimsin98 commented Nov 15, 2019

@oesteban In my case, there are sudden patches of fieldmap with maximum values (300+) in regions with low inhomogeneity (absolute values < 10). The patches are outside brain regions where susceptibility distortion is usually expected (ventral frontal/temporal).

The bright patches are also in calculated phasediff, but I noticed that the original wrapped phases do not have such large differences. However, the region does have a jump in phase values. I suspected some sort of overaggressive regularization because the location of that jump gets smoothed too much for short TE phase after PRELUDE. This leads to large differences in unwrapped phases. N4 is the only regularizing step between raw phase images and difference in unwrapped phases.

p.s. After writing all that, I realized having such jump in unwrapped phase values is suspicious, since the point of unwrapping is to smooth wrapped phases into a continuous function.

@kimsin98
Copy link
Contributor

Well, this is very unexpected. I thought the blame on N4 was tenuous, but the real culprit was masked PRELUDE. I know this is bad practice, but unmasked PRELUDE with non-extracted magnitude properly preserves the phase value jump location. This results in proper continuous phasediff without sudden patches of high differences.

Of course, no extraction/mask comes with its own slew of problems, so I definitely won't do that. However, this finally does pinpoint the problem as PRELUDE's odd behavior near boundaries. I will try messing with PRELUDE and BET options.

FUGUE user guide mentions brain extraction but not bias correction. I tried following the user guide exactly without N4 and got a more reasonable fieldmap.

Ha, this only happened because the FUGUE guide omits the mask parameter in their PRELUDE example. Sorry, N4. You are innocent.

@kimsin98
Copy link
Contributor

Using a tighter mask (0.7 threshold) fixed the problem. I am guessing some bad outliers near the brain edge threw off PRELUDE.

How was the current threshold (0.6) chosen? If the choice was arbitrary to this workflow, perhaps we could add a parameter for adjusting the magnitude brain extraction threshold.

Also, is there a good method of estimating the ideal threshold based on input? It would suck if we had to determine it case by case.

oesteban added a commit to oesteban/sdcflows that referenced this pull request Nov 18, 2019
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).
@kimsin98
Copy link
Contributor

kimsin98 commented Nov 19, 2019

Perhaps less severe than before, but my fieldmaps still have odd blobs in various places (mostly ventral).

fieldmap_blobs

Most blobs are small enough that I feel like they could be erased using proper regularization, but they are large enough that despiking/median filter is not sufficient.

Are these blobs expected in phase differences? If so, I will just try to look for stronger regularization methods.

@mattcieslak
Copy link
Collaborator Author

@Aksoo The original workflow for processing those PNC fieldmaps didn't use prelude on the phase images, but directly took arctan2 of them to get the phasediff and ran prelude on the result. I'm also finding these suspicious jumps when using the official fsl method of prelude-ing the phase images and then subtracting and running prelude again on that subtraction. Sometimes these will go away if you increase the size of the denoising kernel, but that is not great

@oesteban
Copy link
Member

@mattcieslak sorry I missed the benefits of the arctan2 approach. This PR is just too big for review: I had to juggle between several moving pieces (plus a very unstable period where docs were not generated and circle would fail often). I hope that, tackling #17 separately, and with #50 and #51 merged everything will show more clear.

Please allow me to finish #17 and after we can compare the arctan2 approach.

@mattcieslak
Copy link
Collaborator Author

@oesteban I agree. Should we close this and start fresh with a different PR once the others are ready?

@oesteban
Copy link
Member

No rush with closing this, but that is probably the cleanest route to finish up.

@kimsin98
Copy link
Contributor

kimsin98 commented Nov 20, 2019

@mattcieslak Oh right! I remember seeing arctan2 in the original nipreps/fmriprep#1359 PR before finding sdcflows. I will try getting phasediff first then unwrapping.

Using arctan2 for phase difference seems unnecessarily convoluted though. You are

  1. converting phase angles to a+bi form
  2. calculating phase0 * conj(phase1), which is equivalent to just taking the difference of angles
  3. then applying -arctan2 to get the phase angle difference

The only possible benefit of arctan2 here is ensuring that the output falls in [-pi, pi] range, but you can conditionally add/subtract 2*pi to ensure that easily.

@kimsin98
Copy link
Contributor

kimsin98 commented Nov 20, 2019

After thinking about what PRELUDE does and testing a few fieldmaps, I think I have the best workflow for calculating phase difference from 2 phase images. FUGUE guide's method is not robust.

PRELUDE tries to minimize the differences between regions of similar voxels. The blobs I kept getting in the fieldmap is exactly what PRELUDE tries to eliminate. However, both PRELUDE and real world phase images are not perfect. Given 2 phase images that should be similar, PRELUDE can still produce results with +/- 2n pi differences in corresponding voxels. These unwrapped phase images could make sense individually but produce those sudden blobs together.

Thus, it is simply better to run PRELUDE after subtracting the raw phase images. This way, the final phasediff and fieldmap would be as smooth as possible because we are working directly off of PRELUDE's difference-minimized output. In ideal case, there is no difference between running PRELUDE before or after taking subtracting phases anyway, since all PRELUDE does is adjust voxels by +/- 2n pi.

So the workflow should be

  1. Convert phase images to radians. (Scale to [0, 2 pi])
  2. Subtract short phase from long phase. (Output values in [-2 pi, 2pi])
  3. PRELUDE requires values to be within 2 pi range. Add/subtract 2 pi to values outside [- pi, pi].
  4. PRELUDE
  5. Regularization, edge cleaning

@oesteban
Copy link
Member

@Aksoo please hold off for one or two days - the api is changing quickly, so working off this branch will become very impractical.

@oesteban oesteban force-pushed the master branch 2 times, most recently from b4edfbc to 065a2c4 Compare November 22, 2019 00:14
oesteban added a commit that referenced this pull request Nov 23, 2019
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).
oesteban added a commit to oesteban/sdcflows that referenced this pull request Nov 25, 2019
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.
oesteban added a commit to oesteban/sdcflows that referenced this pull request Nov 25, 2019
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.
@oesteban
Copy link
Member

Closed in favor of #60 for the base implementation, and subthreads #22, #23, and #59 for relevant details to be discussed.

@oesteban oesteban closed this Nov 25, 2019
@nipreps nipreps locked as resolved and limited conversation to collaborators Nov 25, 2019
Sign up for free to subscribe to this conversation on GitHub. Already have an account? Sign in.
Projects
None yet
Development

Successfully merging this pull request may close these issues.

6 participants