Misalignment of preprocessed T1w with freesurfer segmentation?

[bwinsto2]

Dear @mattcieslak, we are running an analysis in which we're aligning streamline endpoints to native space surface vertices. There seemed to be some issue of misalignment and I think we've traced it back to the qsiprep preprocessed T1w images.

Below are data from HCP subject 105923. The first screenshot is the native freesurfer pial surface (after hcp preprocessing) overlaid on the T1w_acpc_dc_restore.nii image found in the HCP minimally preprocessed download. Looks good:

The below screenshot is the same surface segmentation overlaid on the qsiprep preprocessed T1w image from HCP unprocessed data. Here was the command I used:

qsiprep-docker \
/data/HCP_Raw/raw_data /data/HCP_Raw/derivatives/ participant \
--output-resolution 1.25 \
--output-space T1w \
--hmc-model eddy \
--fs-license-file /data2/Brian/random/license.txt \
--participant_label 105923 \
--work-dir /data/HCP_Raw/qsiprep/tmp_105923 \
--distortion-group-merge average \
--gpus all \
--eddy-config /data2/Brian/connectome_harmonics/eddy_config.json \
--skip-t1-based-spatial-normalization

I used --skip-t1-based-spatial-normalization in the hopes that this would leave everything in native T1w space. I don't have screenshots of the alignment with the default settings (i.e. using spatial normalization), but we had similar downstream issues, so I assume this misalignment issue remains. If you think this was the issue I can re-run with that setting off.

Do you know why we see this misalignment between the qsiprep preprocessed T1 and the native freesurfer segmentation? Is there a way we can configure qsiprep to output a T1 (and thus a T1 mask and 5tt) that are in the same space as the HCP T1w_acpc_dc_restore.nii, e.g.?

Thank you very much! Brian

[mattcieslak]

I wonder if this is related to gradient distortion correction. QSIPrep doesn't perform GDC at the moment, so it's possible that your HCP preprocessed T1w may have a different shape than the QSIPrep preprocessed T1w.

Another possibility is that there is a slightly different registration to AC-PC alignment in the HCP pipeline than in QSIPrep. Both QSIPrep and HCP perform a rigid transformation of the T1w into AC-PC alignment but they go about it in different ways.

Can you try registering the acpc_dc_restore T1w to the desc-preproc T1w from qsiprep to see if they have the same shape when registered to one another with a rigid transform? An important note is that the data coming out of the HCP pipelines is LAS+ image orientation and QSIPrep is LPS+. This should only cause very obvious problems if not handled properly, but it's worth being aware of.

[bwinsto2]

Thank you, @mattcieslak! After a rigid transform, one of the brains looked a little more squished than the other, so I'm inclined to believe that you're correct that GDC is the difference. Before I test that, two questions, if I may:

Is it possible to configure QSIPrep in such a way to output (and use) a T1w in the same space as fMRIprep preprocessed (or raw T1) space? As I mentioned, we are trying to align streamline endpoints to surface vertices, so we'd need surfaces and tractography in the same space. Native freesurfer surfaces align with fMRIprep/raw T1w, so if QSIPrep could use a T1w in this space, everything downstream would also be in that space:

If this is not possible, I see that QSIPrep also has an option to run recon-all. But it looks like the input to recon-all is the raw T1w image(s):

Theoretically another option for us would be if freesurfer were run on the QSIPrep preprocessed T1w since then the surfaces and QSIPrep tractography would be in the same space.

Do you think either of these options are possible within QSIPrep?

I suppose worst case we could run --dwi-only and then register the dwi to the fMRIprep T1w and run tractograpghy ourselves.

Thank you!

[mattcieslak]

Is it possible to configure QSIPrep in such a way to output (and use) a T1w in the same space as fMRIprep preprocessed (or raw T1) space?

The ACPC-aligned output from qsiprep will always be one rigid transformation away from the raw T1w image. Since fmriprep also doesn't use GDC, qsiprep's ACPC-aligned t1w will also be one rigid transform away from the fmriprep's desc-preproc_space-t1w images. At the moment, the rigid transform from native t1w to acpc isn't saved in the derivatives, but is still available in the working directory (if you keep that around).

I'd suggest performing all your reconstruction in the qsiprep-aligned world so your streamlines would be one tcktransform away from your native t1w surface data - with no image resampling required.

Theoretically another option for us would be if freesurfer were run on the QSIPrep preprocessed T1w since then the surfaces and QSIPrep tractography would be in the same space.

The freesurfer workflow comes from fmriprep 1.26 (very old) but will give you gifti surfaces aligned with the ACPC T1w output. I haven't thoroughly tested this though, and would probably recommend the earlier suggestion. Especially if you've already run freesurfer

One side note - the next phase of qsiprep development is going to focus heavily on working with outputs from fmriprep and the hcp pipelines. You're already a bit ahead of us on this, so this is super helpful

[bwinsto2]

Thanks for the reply! Unfortunately after a lot of trying, we still haven't been able to resolve this. A couple questions:

the rigid transform from native t1w to acpc isn't saved in the derivatives, but is still available in the working directory

I'm assuming this is the transformation inside skullstrip_wf/rigid_acpc_align from highres001_N4Corrected0.nii.gz (t1) -> mni_1mm_lps_brain.nii.gz (acpc). It seems to me that highres001_N4Corrected0.nii.gz is not in the same space as native T1 or fmriprep desc-preproc_space T1, though, since it doesn't align with native freesurfer segmentation:

Here is the skullstripped fmriprep desc-preproc_space T1, for comparison (raw t1 is essentially the same):

I'd suggest performing all your reconstruction in the qsiprep-aligned world so your streamlines would be one tcktransform away from your native t1w surface data

We used the acpc -> native_t1 transform in qsiprep's work dir for a linear tcktransform, but the streamlines and surfaces were not aligned (which makes sense given above). So I replaced highres001_N4Corrected0.nii.gz with the skullstripped fmriprep desc-preproc_space T1 in the antsRegistration call and did the same process. The results were significantly better (the streamlines were more registered to the native surface), but still not very good. I also tried replacing mni_1mm_lps_brain.nii.gz in the antsRegistration call with qsiprep's desc-preproc T1, since I assume the tracks are registered to this space. This yielded another close-ish result, but still not aligned.

To demonstrate, below are downsampled streamline endpoints as vtks overlaid on the surfaces. On the left is our HCP pipeline (using hcp preprocessed native 32k surfaces, with streamlines aligned to T1w_acpc_dc_restore.nii), which fits well. On the right is the native 32k surface with the tracks tcktransform-ed using the fmriprep desc-preproc_space T1 -> qsiprep's desc-preproc T1 warp (tcktransform needs the inverse direction i.e. native -> acpc to get the tracks from acpc -> native). Clear misalignment of endpoints to surface on the right:

Any thoughts why we seemingly can't get these qsiprep tracks to align with native surfaces? If this will be addressed in the next phase of development, that would be wonderful! But if there's some variation of what I'm trying above that can serve as a workaround for now, that would be even more wonderful. We're actually working on a pipeline that we plan to open-source that will hopefully accept the outputs of qsiprep/recon (would be happy to specify more details over email). This is looking like the only major roadblock for now.

For reference, here is the shell script I'm using to do all this:

I've been swapping different things in and out for native_t1 and acpc_template, as described above (e.g. native_t1 = skullstripped fmriprep desc-preproc_space T1; acpc_template = qsiprep's desc-preproc T1; tracks = qsirecon tracks)

native_t1=$1
acpc_template=$2
tracks=$3

#register native t1 to acpc with antsRegistration
antsRegistration --collapse-output-transforms 1 --dimensionality 3 --initial-moving-transform [ ${acpc_template}, ${native_t1}, 1 ] --initialize-transforms-per-stage 0 --interpolation LanczosWindowedSinc --output [ transform, transform_Warped.nii.gz ] --transform Rigid[ 0.2 ] --metric Mattes[ ${acpc_template}, ${native_t1}, 1, 32, Random, 0.25 ] --convergence [ 10000x1000x10000x10000, 1e-06, 10 ] --smoothing-sigmas 7.0x3.0x1.0x0.0vox --shrink-factors 8x4x2x1 --use-histogram-matching 1 --winsorize-image-intensities [ 0.025, 0.975 ]  --write-composite-transform 0

#convert ANTs affine matrix (native t1 to acpc) to ITK txt format, then to mrtrix format
ConvertTransformFile 3 transform0GenericAffine.mat transform.txt
transformconvert transform.txt itk_import transform_mrtrix.txt

#convert native t1 (moving) and acpc template (reference) images to mrtrix format
mrconvert ${native_t1} native_t1.mif
mrconvert ${acpc_template} acpc_template.mif

#initialize identity warp for native t1 (moving)
warpinit native_t1.mif id_warp.mif

#transform identity warp to image warp with t1->acpc matrix
mrtransform id_warp.mif -linear transform_mrtrix.txt native_to_acpc_warp.mif -template acpc_template.mif -interp linear -nan

#transform tracks from acpc space to native space (INVERSE OF ABOVE TRANSFORM)
tcktransform ${tracks} native_to_acpc_warp.mif transformed_tracks.tck

tckresample -downsample 15 transformed_tracks.tck ds_tf_tracks.tck

tckresample -endpoints ds_tf_tracks.tck ds_tf_endp.tck

tckconvert ds_tf_endp.tck ds_tf_endp.vtk

[mattcieslak]

Based on the first image I wonder if some of these problems are coming from different image orientations (there is a big red warning on the bottom right of your first screenshot). Each of the 3+ different software packages involved here take different approaches to how much they care about the image headers.

In my experience, you will be doing yourself a huge favor by converting everything to the same orientation before doing any sort of spatial operations. keep in mind, your HCP results will be in LAS+, qsiprep is in LPS+ and mrtrix is flexible but uses RAS+ coordinates in the tck files.

Also, I'd definitely be interested in hearing more about your pipeline. It sounds like we're going to be working on similar problems in the next few months

[bwinsto2]

Thanks @mattcieslak for your continued support, and please see my email re: details about our pipeline :)

Your suggestion about image orientation is well-taken; however, I’ve confirmed that it wasn’t the issue. It seems the source of this issue is likely the tractography output I am getting from qsiprep/recon (before the tcktransform warp to native space). Apologies if I’m wrong, but I assume the tractography is by default registered to qsiprep’s desc-preproc T1w. I did a clean run of qsiprep/recon with the following options, and I’m getting an unexpected result. These data are HCP unprocessed data BIDS-ified with hcp2bids:

Preprocessing:

qsiprep-docker \
/data/HCP_Raw/raw_data /data/HCP_Raw/derivatives/ participant \
--output-resolution 1.25 \
--output-space T1w \
--hmc-model eddy \
--fs-license-file /data2/Brian/random/license.txt \
--participant_label 105923 \
--work-dir /data/HCP_Raw/qsiprep/tmp_105923 \
--distortion-group-merge average 

Reconstruction:

qsiprep-docker \
/data/HCP_Raw/raw_data /data/HCP_Raw/derivatives/ participant \
--recon_input /data/HCP_Raw/derivatives/qsiprep \
--work-dir /data/HCP_Raw/ignore/qsiprep/tmp_105923_recon \
--fs-license-file /data2/Brian/random/license.txt \
--recon_spec mrtrix_multishell_msmt \
--participant_label 105923 \
--recon-only

Here are downsampled tracks overlaid on the desc-preproc T1w. Upon close inspection, there are quite a few white matter areas with no streamlines, and many places with streamlines that aren't white matter (more screenshots below):

For comparison, here is our HCP pipeline, which has cleaner alignment:

Is there anything you suspect to be the problem with the qsiprep settings I'm using? Am I overlaying the streamlines on the wrong anatomical? Any insight would be very much appreciated!

Here are a couple more screenshots of the qsiprep/recon output for reference:

Thank you so much!

[mattcieslak]

I wonder if something is going on with ACT. Could you try running the mrtrix_multishell_msmt_noACT pipeline? It could also be helpful to do a quick dsi_studio_gqi? The latter is pretty bulletproof and should give us an idea of something really bizarre is happening

[bwinsto2]

Thanks :) after some trial-and-error, it seems as if the streamlines were misaligned with the white matter due to 5tt images with poor boundaries. The default 5ttgen fsl produces:

But this was completely resolved by using the new 5ttgen hsvs option (note these are slightly different slices), and the streamlines fit the freesurfer surfaces beautifully:

I understand that qsiprep is quite customizable using the .json files, and it looks like 5ttgen freesurfer might be supported, which was an improvement for us but not quite as good as hsvs.

[mattcieslak]

oh wow, that is a huge difference in the quality of the segmentation. I haven't heard of the HSV option, but it seems like users may want to avoid ACT in qsiprep until this is further investigated. Thanks!

[ghost]

Dear experts, regarding this topic, I am also not very happy with the default 5ttgen based on fsl segmentation. I ran qsiprep with the --do-reconall option and I would like to use freesurfer segmentation to generate the 5tt file. Is there a way to customize the 5ttgen call using the freesurfer algorithm (or even better the hsvs one) by creating a custom reconstruction .json file?

Thank you, Giuseppe

[mattcieslak]

Currently working on supporting 5ttgen's HSVS (which are a significantly better) #287. A rough sketch of the workflow is this:

# Get 5tt
5ttgen hsvs freesurfer/sub-X 5tthsv.nii -white_stem

# Preprocessed T1w from qsiprep - aligned with DWI
qpt1=/path/to/qsiprep/sub-X_desc-preproc_T1w.nii.gz

# get fs image into LPS+
mrconvert freesurfer/sub-X/mri/orig_nu.mgz fsimg.nii -strides -1,-2,3
fst1=fsimg.nii

# Standard rigid registration with ANTs
antsRegistration --dimensionality 3 --float 0 \
        --output [fs2qp_,fs2qp_Warped.nii.gz] \
        --interpolation LanczosWindowedSinc \
        --initial-moving-transform [$qpt1,$fst1,1] \
        --transform Rigid[0.1] \
        --metric MI[$qpt1,$fst1,1,32,Regular,0.25] \
        --convergence [1000x500x250x100,1e-6,10] \
        --shrink-factors 8x4x2x1 \
        --smoothing-sigmas 3x2x1x0vox

# ANTs in qsiprep makes binary transform matrices - convert to text for mrtrix
ConvertTransformFile 3 fs2qp_0GenericAffine.mat fs2qp_0GenericAffine.txt
# Convert text transform to mrtrix format
transformconvert fs2qp_0GenericAffine.txt itk_import fs2qp_transform.txt
# Apply directly to affine of the 5tt image
mrtransform -linear fs2qp_transform.txt 5tthsv.nii hsd_trf_ants_5tt.nii

Essentially, we're going to calculate the 5tt image in freesurfer space, find the transform from freesurfer space to qsiprep t1w space, then adjust the header of the 5tt image so it's aligned with the qsiprep t1w in world coordinates. I also plan on generating a visual report so you can check how the segmentation looks on your qsiprep t1w.

[ghost]

Dear @mattcieslak, this is great, looking forward to it!

This would be the workflow assuming that freesurfer reconall is already run on the raw T1w (that's why you would have to align the segmentation with the qsiprep-preprocessed T1w). Otherwise, if freesurfer is run from inside qsiprep with the --do-reconall option, the segmentation should be already in the qsiprep-T1w space. Am I correct?

[mattcieslak]

For the big picture, this approach makes as few assumptions about the freesurfer run as possible. I want to only require a freesurfer directory because the next immediate plan is to support hcp pipeline outputs in the reconstruction workflows.

From the detailed perspective, 5ttgen hsvs requires both surfaces and the volume segmentations. the *prep outputs produce gifti surfaces and volume segmentations aligned with the t1ws and the labels re-mapped to bids, which 5ttgen hsvs doesn't know what to do with.

The --do-reconall option isn't particularly useful right now, other than providing a standard freesurfer output.

[ghost]

Ok, seems fair to me.

Thank you for your work!