Skip to content

Code/pipeline review #9

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.

Sign up for GitHub

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

Open
wants to merge 1 commit into
base: main
Choose a base branch
from
Open

Code/pipeline review #9

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
12 changes: 12 additions & 0 deletions processing/process_ihmt.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -124,6 +124,7 @@ def process_run(name_source, layout, run_data, out_dir, temp_dir):
name_base = os.path.basename(name_source)

# ihMTRAGE parameters
# These will eventually be read from the BIDS metadata?
interpulse_delay = 100 # ms
n_pulses_per_burst = 1
n_bursts = 10
Expand All @@ -146,6 +147,7 @@ def process_run(name_source, layout, run_data, out_dir, temp_dir):
concat_ihmt_file = os.path.join(temp_dir, f'concat_{name_base}')
concat_ihmt_img.to_filename(concat_ihmt_file)

# How many volumes are in ihmt?
Copy link
Member

Choose a reason for hiding this comment

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

There are five volumes with different contrasts. I took the dwidenoise step from ihmt_proc.

dwidenoise_extent = '3,3,3'
denoised_ihmt_file = os.path.join(temp_dir, f'dwidenoise_{name_base}')
cmd = f'dwidenoise {concat_ihmt_file} {denoised_ihmt_file} -force --extent {dwidenoise_extent}'
Expand Down Expand Up @@ -219,7 +221,11 @@ def process_run(name_source, layout, run_data, out_dir, temp_dir):
ihmtw_img.to_filename(ihmtw_file)

# Calculate ihMTR
# Is there an advantage to keeping ihmt_files_t1space as a list?
# It's less intuitive than using variables
Comment on lines +224 to +225
Copy link
Member

Choose a reason for hiding this comment

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

I wanted to loop over the images, but there are more intuitive ways like you said. Maybe a dictionary so it's still easier to organize the image types in a single variable I can index?

ihmtr_img = image.math_img('ihmt / m0', ihmt=ihmtw_img, m0=ihmt_files_t1space[0])
# What do you think about making a project-wide setting for what to do with calculations that
# involve division? It's an unstable operation and we should handle it consistently.
Comment on lines +227 to +228
Copy link
Member

Choose a reason for hiding this comment

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

Could you expand on this idea?

Copy link
Author

Choose a reason for hiding this comment

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

One of the problems with T1/T2 "myelin maps" is that if the denominator gets close to zero you get extreme numbers. I was imagining something like a config setting for what to do when there is a division and you see extreme values. For example, if some voxels are inf or extreme outliers, should they be capped? Or set to nan? Or should we leave them as is? It makes sense that if we offer this option for one of the calculations that we should offer it for all of them.

ihmtr_file = get_filename(
name_source=name_source,
layout=layout,
Expand Down Expand Up @@ -332,6 +338,7 @@ def process_run(name_source, layout, run_data, out_dir, temp_dir):
fixed=ants.image_read(run_data['t1w_mni']),
moving=ants.image_read(file_),
transformlist=[run_data['t1w2mni_xfm']],
# Add interpolation setting here? Default is linear if you don't.
)
ants.image_write(reg_img, out_file)

Expand Down Expand Up @@ -371,6 +378,7 @@ def iterative_motion_correction(name_sources, layout, in_files, filetypes, out_d
n4_img = ants.n4_bias_field_correction(in_img)

# Skull-stripping
# How well has this been working? It could be worth using synthstrip if it's flaky.
dseg_img = antspynet.utilities.brain_extraction(n4_img, modality='t1threetissue')
dseg_img = dseg_img['segmentation_image']
mask_img = ants.threshold_image(
Expand Down Expand Up @@ -452,6 +460,10 @@ def iterative_motion_correction(name_sources, layout, in_files, filetypes, out_d
)
ants.image_write(out_img, out_file)

# There is a final step here that has to happen, which is you average together
# all the transforms and apply the inverse to the template.
# This prevents the template from drifting.
# See https://github.com/PennLINC/qsiprep/blob/c74ff2f041a9ff6121042d30b5467af6082be85f/qsiprep/workflows/dwi/hmc.py#L273
Comment on lines +463 to +466
Copy link
Member

Choose a reason for hiding this comment

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

Could you expand on this? I think this might be one spot where QSIPrep's approach and ihmt_proc's differ. The template is only recalculated with a voxel-wise average, so I don't think there's any opportunity for it to drift.

Copy link
Author

Choose a reason for hiding this comment

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

It currently won't be an issue here because you're only doing a single iteration. If there were multiple iterations, it's possible that the template brain will move around in the FOV. The workflow in qsiprep is based on the antsMultivariateTemplateConstruction.sh script, where a couple iterations are always used. If there is a lot of motion the first average won't look very good, so the registrations to it will still have some error. Over the multiple iterations the template starts to converge to something nice.

I'm not sure if multiple iterations is overkill here. It would depend on how much motion you expect over these volumes.

Copy link
Member

Choose a reason for hiding this comment

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

The method in ihmt_proc only runs it once. I think it should be fine, but we can revisit and apply a multi-iteration method if it doesn't work well.

return template_file, transforms


Expand Down