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

Added scripts to filter B1 and B0 maps #2

Open
wants to merge 14 commits into
base: main
Choose a base branch
from
Open

Conversation

jcohenadad
Copy link
Member

@jcohenadad jcohenadad commented Jul 29, 2022

Fixes #1

This filtering depends on fixing astropy/astropy#13511

Testing this PR:

@jcohenadad
Copy link
Member Author

Preliminary results with 3D Gaussian filtering (ignoring masked voxels) with a 9x9x9 kernel:

anim

@@ -1,4 +1,5 @@
# Licensed under a 3-clause BSD style license - see LICENSE.rst
# Author: Mathieu Guay-Paquet

Choose a reason for hiding this comment

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

Most of the code is copy-pasted from astropy, so this should probably be:

# Author: Astropy Developers, Mathieu Guay-Paquet

(See the original LICENSE.rst.)

Copy link
Member Author

Choose a reason for hiding this comment

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

fixed in 0324a34

@mathieuboudreau
Copy link
Member

FYI: we have 4 different types of B1 filtering implemented in qMRLab, it's in MATLAB but here's a blog post about it: https://qmrlab.org/2019/08/08/b1filter.html

@mathieuboudreau
Copy link
Member

In particular, there's an interactive figure at the end that compares the B1 maps for 1 type of artifact. I didn't explore more in the blog post, but I've thought a lot about this over the years so have some more opinions on this.

@jcohenadad
Copy link
Member Author

jcohenadad commented Jul 29, 2022

FYI: we have 4 different types of B1 filtering implemented in qMRLab, it's in MATLAB but here's a blog post about it: https://qmrlab.org/2019/08/08/b1filter.html

Amazing!! Thank you so much for chipping in @mathieuboudreau ! Would you mind giving it a try and see if it improves things compared to our Gaussian filtering? I've made data available in the PR body.

EDIT 2022-07-29 10:51:51: From https://qmrlab.org/2019/08/08/b1filter.html Figure 6, I'd say that median (vs. gaussian) filtering might be more appropriate for our maps (which show some outliers).

@mathieuboudreau
Copy link
Member

mathieuboudreau commented Jul 29, 2022

EDIT 2022-07-29 10:51:51: From https://qmrlab.org/2019/08/08/b1filter.html Figure 6, I'd say that median (vs. gaussian) filtering might be more appropriate for our maps (which show some outliers).

Yeah, something that's always annoyed me is the lack of discussion in the B1 mapping community about when and why they use certain filters**. It's really going to be artefact/image dependent; for example, in AFI a common issue is a "stripping pattern" from Gibbs Ringing, and and median isn't going to work well unless you use a very large kernel. At the end of the day, you know that the true B1 map is smoothly varying and the purpose of the filter should just be to get close to that, but I found certain filters can do more harm than good in certain images. Unless I get a very odd artefact (eg, due to phase poles in a Bloch-Siegert map for example), my first choice would be median.

Another good way to check is by doing an absolute or % difference map between filtered and unfiltered, to see how the filter impacted your final map (eg just denoising, or did it attenuate values in regions, etc).

I'll test out this feature as you implemented this afternoon, and test out a median filtering implementation!

@mathieuboudreau
Copy link
Member

mathieuboudreau commented Jul 30, 2022

Opened a PR into this branch #3

It implements an initial draft of the median filter. It has better edge preservation than the scipy median_filter; it ignores NaNs from the median calculation, so it has better edge preservation. I also made a Jupyter Notebook to help generate some comparisons.

Here are the three orientation views of the original, median (5x5x5 voxel), and gaussian-filtered B1 maps:
sagittal
coronal
axial

To have a better view of what the filters are "removing", difference maps (%) relative to the original map are plotted next for one orientation:

percdiff

Red -> filtered maps have higher value than original map, blue -> filtered maps have lower values

Lastly, a comparison between the two filtered-maps are plotted:

gauss-med

Zoom in to view specific differences (arrow):

Screen Shot 2022-07-30 at 1 30 35 AM

(Note that this would probably smoothed out with a larger Gaussian kernel, or preserved with a 3x3x3 Median one.)

Also note that because the voxels are anisotropic (2x2x4), but the median filter is 5x5x5 voxel, more "filtering" is happening along the slice direction (which might lead to some systemic underestimations in B1 values), but this could be fixed by using the "footprint" option of the generic_filter tool.

Does this improve vs Gaussian filter? Using these parameters, I think so, but I haven't played around with the Gaussian
kernel size any more (I could play around with it a bit more next week). Also, I imagine that only the B1 values exactly in the spinal cord are of interest, so segmenting this ROI and comparing the values would help with the comparison. It'll still likely be subjective, there's no real ground truth here, it's a pull-push balance of 1) knowing the true B1 imap s expected to be a smooth function, but 2) if you try and make it too smooth by increasing the filter/kernel size, it will create systemic bias by typically lowering the B1 values due to including voxels too far away from the calculating pixel during the filtering.

@jcohenadad
Copy link
Member Author

@mathieuboudreau this is amazing! thank you so much for spending time doing this. I like results from the median filter better. I agree the kernel should be adjusted to the effective spatial resolution.

@mathieuboudreau
Copy link
Member

Great! I just pushed the changes to handle anisotropic pixel dimentions (eg slice thickness different from in-plane resolution) in commit f4db728 of my PR to this branch #3

Here's a gif of the difference (5x5x5 voxel kernel size vs 5x5x3 (rounded up the slice size)).
gif

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Filter B1 and B0 maps
3 participants