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: Adds diffusion-related IQMs. #1131

Merged
merged 56 commits into from
Apr 1, 2024
Merged

Conversation

arokem
Copy link
Collaborator

@arokem arokem commented Jul 20, 2023

I am putting up this branch so that folks at OHBM hackathon have an upstream branch to make PRs against, to add the different IQMs that we are all working on.

@arokem
Copy link
Collaborator Author

arokem commented Feb 5, 2024

We would like to be able to group b-values by orders of magnitude. Is there already a function in the code-base that does something like that? Alternatively, we could copy over some code from DIPY:

https://github.com/dipy/dipy/blob/master/dipy/core/gradients.py#L816C1-L842C1
https://github.com/dipy/dipy/blob/master/dipy/core/gradients.py#L923-L950

The idea would be to calculate IQMs by groups, based on the b-value "tiers".

@arokem
Copy link
Collaborator Author

arokem commented Feb 5, 2024

NVM - the previous comment is silly, because obviously we're going to need DIPY as a dependency for so much of this that it doesn't make sense to start copying over code.

@arokem arokem force-pushed the start_on_diffusion branch 2 times, most recently from 68cb9b0 to caf966c Compare February 7, 2024 23:33
@oesteban
Copy link
Member

We would like to be able to group b-values by orders of magnitude. Is there already a function in the code-base that does something like that? Alternatively, we could copy over some code from DIPY:

https://github.com/dipy/dipy/blob/master/dipy/core/gradients.py#L816C1-L842C1 https://github.com/dipy/dipy/blob/master/dipy/core/gradients.py#L923-L950

The idea would be to calculate IQMs by groups, based on the b-value "tiers".

Sorry, I missed your comment. We actually do have this already in the workflow. It uses DIPY underneath but also is able to create b-value "clusters" for DSI, thanks to an algorithm @edickie first proposed in dMRIPrep.

@arokem
Copy link
Collaborator Author

arokem commented Mar 18, 2024

Thanks @oesteban : can you point to this code? Might be useful here as well for DSI-like acquisitions?

@oesteban
Copy link
Member

Thanks @oesteban : can you point to this code? Might be useful here as well for DSI-like acquisitions?

The interface definition is here:

class NumberOfShells(SimpleInterface):

It provides several outputs:

class _NumberOfShellsOutputSpec(_TraitedSpec):
models = traits.List(traits.Int, minlen=1, desc='number of shells ordered by model fit')
n_shells = traits.Int(desc='number of shels')
out_data = traits.List(
traits.Float,
minlen=1,
desc="new b-values table (after 'shell-fying' DSI)",
)
b_values = traits.List(traits.Float, minlen=1, desc='estimated values of b')
b_masks = traits.List(
traits.List(traits.Bool, minlen=1),
minlen=1,
desc='b-value-wise masks')
b_indices = traits.List(
traits.List(traits.Int, minlen=1),
minlen=1,
desc='b-value-wise masks')
b_dict = traits.Dict(
traits.Int, traits.List(traits.Int), desc='b-values dictionary'
)

I think you want a combination of the b_masks and b_values outputs (the former gives you binary masks over the b-vecs/b-vals table corresponding to each entry in the latter, which are ordered by growing magnitude).

This is in the workflow used here:

shells = pe.Node(NumberOfShells(), name='shells')

@oesteban
Copy link
Member

@arokem I will be pushing this over the finish line during the next two days - apologies in advance for a potential deluge in commits :)

The plan is to connect your QC functions into the pipeline, so after merging, DWI will be generating and storing at the output - at least - the metrics coded in this PR.

@oesteban
Copy link
Member

Thanks @oesteban : can you point to this code? Might be useful here as well for DSI-like acquisitions?

After a first dive into the code I can totally see where you will need this. I'll continue working on this for a little longer and let you know when it is ready for your review (see my last commit is marked WIP, I'd suggest waiting for me to finish this tomorrow before adding more commits from your side).

@@ -545,3 +654,108 @@ def _extract_b0(in_file, b0_ixs, out_path=None):

def _exp_func(t, A, K, C):
return A * np.exp(K * t) + C


def get_spike_mask(
Copy link
Member

Choose a reason for hiding this comment

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

@arokem In the process of bringing this from the qc module.

I see the generation of this mask as a "preprocessing" step. The idea is to integrate it at an earlier step of the pipeline and feed it into the node that calculates IQMs. That way, we can define future IQMs using these masks without moving the calculation or re-doing it.

b0 = data[..., ~gtab.b0s_mask]
return np.percentile(np.var(b0[mask], -1), (25, 50, 75))

# OE: Create interface to mask the CC
Copy link
Member

Choose a reason for hiding this comment

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

@arokem Same as above. We already fit a DTI model, so we should use the FA from that step (unless this one has something that makes it special). I think, again, that the CC mask can be really helpful for future IQM definitions beyond those we are adding at this time. Also worth plotting, I think.

return spike_mask


def noise_b0(data, gtab, mask=None):
Copy link
Member

Choose a reason for hiding this comment

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

I think we can integrate the percentile calculation in the ExtractB0 interface (or similar place).

Copy link
Member

Choose a reason for hiding this comment

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

I haven't been able to look throgh this -- will do asap

Copy link
Member

Choose a reason for hiding this comment

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

I refrain from having "preprocessing" here because there is this lingering idea of writing code to make IQMs plug and play. Basically, if we can standardize the arguments passed into the functions, we can iterate over and execute all functions, setting the right inputs into them.

All these functions could be defined with **kwargs to allow passing all the available objects within the IQM executor.

@oesteban oesteban linked an issue Mar 22, 2024 that may be closed by this pull request
@oesteban oesteban force-pushed the start_on_diffusion branch 3 times, most recently from b23e9bd to 9f127d5 Compare March 22, 2024 06:12
@oesteban oesteban force-pushed the start_on_diffusion branch 3 times, most recently from bc5510d to c33b99f Compare March 22, 2024 14:35
arokem and others added 6 commits March 25, 2024 15:40
* Some "preprocessing" functions will be moved to previous steps in the
  pipeline (as compute nodes in the workflow). Therefore, they are moved
  into the interfaces module.
* Clarified and re-worked docstrings
* Added type annotations.
* Added to documentation rendering
@oesteban oesteban force-pushed the start_on_diffusion branch 2 times, most recently from af779bd to 00c4faa Compare March 25, 2024 14:56
@mattcieslak
Copy link

Is there any interest in adding neighboring dwi correlation or fiber coherence from dsi studio?

@arokem
Copy link
Collaborator Author

arokem commented Mar 26, 2024

Yes - I think so. These measures would be quite meaningful even without extensive processing, right?

@mattcieslak
Copy link

I think the raw_neighbor_corr was one of the best scoring features in hbn-pod2. It's calculated on the unprocessed data.

Also, I think the computation is simple enough that we could implement it in python without needing to add DSI Studio to the mriqc image. It would be a great way to confirm the actual formula of NDC too.

@oesteban
Copy link
Member

@mattcieslak @arokem would reimplementing those metrics be a hefty lift? I say this because DSI Studio has commercial use limitations and we got rid of FSL just for that reason.

@mattcieslak
Copy link

@oesteban great point. The neighbor correlation method doesn't have a patent but might still have a non-commercial restriction on the use of the dsi studio implementation. I could implement the method from scratch based on the description in the paper, would that be ok?

@oesteban
Copy link
Member

would that be ok?

That would be an outstanding contribution beyond MRIQC, useful for many. Perhaps DIPY would be a better home for the implementation, but I'd be happy to test it out here and port whenever it is finished -- that as @arokem feels it's best.

@arokem
Copy link
Collaborator Author

arokem commented Mar 27, 2024

Yeah - agreed - I think that it would be really good to have here and also in general (and didn't realize it was calculated on unprocessed data! but even more useful here as such). And DIPY would be a good place for it. Maybe in a new sub-module in the denoise module? https://github.com/dipy/dipy/tree/master/dipy/denoise

- [x] Start with interface for IQMs calculation.
- [x] Accommodate calculation of spikes masks, WM mask, etc.
- [x] Start with workflow connections

Resolves: nipreps#1215.
@oesteban oesteban changed the title WIP: Adds diffusion-related IQMs. ENH: Adds diffusion-related IQMs. Mar 31, 2024
@oesteban
Copy link
Member

@arokem @teresamg this is at a point where it works locally on one DWI example. Considering the limited impact of the changes, I think this should work generally. Please review if you have the time and let's merge it so that we can start thinking of more metrics.

@arokem
Copy link
Collaborator Author

arokem commented Mar 31, 2024

+1 to merge (I can't hit the "approve" button, because GitHub considers me the author of the PR).

@oesteban oesteban merged commit d62756f into nipreps:master Apr 1, 2024
15 checks passed
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.

Do not clean up FA for degenerate tensors within DTI node
4 participants