3dresample on fmriprep output moves dset around

mrinmayi · September 17, 2024, 12:31am

AFNI version info (afni -ver): Precompiled binary linux_centos_7_64: Jun 21 2023 (Version AFNI_23.1.08 'Publius Helvius Pertinax')

Hello AFNI gurus!

I'm trying to resample a manually-segmented hippocampal mask to functional space to extract beta values from it. But for some reason, 3dresample seems to shift the dataset a lot. In the attached image, the reddish dataset is the original ROI, and the green is the ROI after resampling.

This is resample command (nothing fancy!):

3dresample -input sub-1004_space-T1w_desc-HandSeg_mask.nii.gz \
    -master sub-1004_task-aim_run-01_space-T1w_desc-preproc_bold.nii.gz  \
    -prefix sub-1004_space-T1w_desc-HandSeg_mask_resampled.nii.gz

The functional and T1w data was preprocessed using fmriprep (v23.0.2). The hippocampal ROI was created by manually segmenting the hippocampus on the individual's T2w scan. The T2 was then aligned to the preprocessed T1 using ANTs. And this transformation matrix was applied to the ROI to bring it to T1 space.

Could this be happening because, somehow, fmriprep puts the T1 and functional data in different spaces, views, etc?

Here are the outputs from 3dinfo on these datasets:

Preprocessed T1:

Identifier Code: NII_lq875bLvWFd_7NgQ7XaCiA  Creation Date: Mon Sep 16 20:20:51 2024
Template Space:  ORIG
Dataset Type:    Anat Bucket (-abuc)
Byte Order:      LSB_FIRST {assumed} [this CPU native = LSB_FIRST]
Storage Mode:    NIFTI
Storage Space:   46,137,344 (46 million) bytes
Geometry String: "MATRIX(-0.999367,-0.018268,-0.030602,94.54953,0.018787,-0.999683,-0.016754,105.0525,-0.030287,-0.017319,0.999391,-122.564):176,256,256"
Data Axes Tilt:  Oblique (2.042 deg. from plumb)
Data Axes Approximate Orientation:
  first  (x) = Left-to-Right
  second (y) = Posterior-to-Anterior
  third  (z) = Inferior-to-Superior   [-orient LPI]
R-to-L extent:   -80.451 [R] -to-    94.550 [L] -step-     1.000 mm [176 voxels]
A-to-P extent:  -149.947 [A] -to-   105.052 [P] -step-     1.000 mm [256 voxels]
I-to-S extent:  -122.564 [I] -to-   132.436 [S] -step-     1.000 mm [256 voxels]
Number of values stored at each pixel = 1
  -- At sub-brick #0 '?' datum type is float

Preprocessed functional:

Identifier Code: NII_FfAmt0ywfFlY1QiUGQBs7Q  Creation Date: Mon Sep 16 20:21:42 2024
Template Space:  TLRC
Dataset Type:    Echo Planar (-epan)
Byte Order:      LSB_FIRST {assumed} [this CPU native = LSB_FIRST]
Storage Mode:    NIFTI
Storage Space:   557,062,800 (557 million) bytes
Geometry String: "MATRIX(-1.705,0,0,72.38453,0,-1.705,0,75.94521,0,0,1.7,-43.88039):85,105,83"
Data Axes Tilt:  Plumb
Data Axes Orientation:
  first  (x) = Left-to-Right
  second (y) = Posterior-to-Anterior
  third  (z) = Inferior-to-Superior   [-orient LPI]
R-to-L extent:   -70.835 [R] -to-    72.385 [L] -step-     1.705 mm [ 85 voxels]
A-to-P extent:  -101.375 [A] -to-    75.945 [P] -step-     1.705 mm [105 voxels]
I-to-S extent:   -43.880 [I] -to-    95.520 [S] -step-     1.700 mm [ 83 voxels]
Number of time steps = 188  Time step = 1.87000s  Origin = 0.00000s
  -- At sub-brick #0 '?' datum type is float
  -- At sub-brick #1 '?' datum type is float
  -- At sub-brick #2 '?' datum type is float
** For info on all 188 sub-bricks, use '3dinfo -verb' **

Manually segmented hippocampal ROI, aligned to T1 using ANTs:

Dataset File:    sub-1005_space-T1w_desc-HandSeg_mask.nii.gz
Identifier Code: NII_AaC_TgbYphvAalPxsUF1Kw  Creation Date: Mon Sep 16 20:22:11 2024
Template Space:  ORIG
Dataset Type:    Anat Bucket (-abuc)
Byte Order:      LSB_FIRST {assumed} [this CPU native = LSB_FIRST]
Storage Mode:    NIFTI
Storage Space:   46,137,344 (46 million) bytes
Geometry String: "MATRIX(-0.999367,-0.018268,-0.030602,94.54953,0.018787,-0.999683,-0.016754,105.0525,-0.030287,-0.017319,0.999391,-122.564):176,256,256"
Data Axes Tilt:  Oblique (2.042 deg. from plumb)
Data Axes Approximate Orientation:
  first  (x) = Left-to-Right
  second (y) = Posterior-to-Anterior
  third  (z) = Inferior-to-Superior   [-orient LPI]
R-to-L extent:   -80.451 [R] -to-    94.550 [L] -step-     1.000 mm [176 voxels]
A-to-P extent:  -149.948 [A] -to-   105.052 [P] -step-     1.000 mm [256 voxels]
I-to-S extent:  -122.564 [I] -to-   132.436 [S] -step-     1.000 mm [256 voxels]
Number of values stored at each pixel = 1
  -- At sub-brick #0 '?' datum type is float

Finally, here is the comparison of the 3 datasets:

3dinfo -same_all_grid sub-1005_desc-preproc_T1w.nii.gz sub-1005_task-aim_run-01_space-T1w_desc-preproc_bold.nii.gz sub-1005_space-T1w_desc-HandSeg_mask.nii.gz 
0	0	1	0	0
0	0	1	0	0
1	1	1	1	1

Strangely, the same thing happens even if I resample the ROI by explicitly setting -dxyz instead of providing a master dataset. So it may not even be the space/view issue?

Thank you, and sorry for the lengthy message!!
Mrinmayi

ptaylor · September 17, 2024, 11:08am

Hi-

I think the issue is the obliquity that the processed T1 and manually segmented ROI have (like in line 8 of the header output: Data Axes Tilt: Oblique (2.042 deg. from plumb)), but that the functional EPI does not have (line 8 of header: Data Axes Tilt: Plumb). It is odd to see obliquity still present in processed data, but the obliquity means that simple resampling isn't enough here.

I believe fMRIprep outputs data in both "subject" and a "template" space. Is this processing mixing data in one (the T1 and ROI have "ORIG" in their header) but the functional appears to be in the latter (has "TLRC" in the header). The question is if the latter is TLRC, even though the FMRI data is in subject space, because fMRIprep has put a NIFTI header information of qform_code =2 or sform_code=2. That is an ambiguous header code that can be interpreted as either being in "native" or "template" space, confusingly.

In either case, I don't think simple resampling is appropriate to move the ROI between the spaces. But let's see if they are really both in "subject" view (just with the FMRI having confusing NIFTI header info) or whether the FMRI is in "template" space and the T1+ROI are in subject space.

What is the output of:

nifti_tool -disp_hdr -field qform_code -field sform_code -infiles DSET_FMRI DSET_T1 DSET_ROI

?

--pt

mrinmayi · September 17, 2024, 12:00pm

Hi Paul,

Thanks a lot for your quick response!

Sorry, I should've mentioned that - I'm having fmriprep run the analysis in native space. So the T1s, ROIs and functionals are meant to be in ORIG space. I'm not sure why fmriprep changes that (I suppose it is the qform_code and sform_code in the NIFTI header). I had tried 3drefit to put the functional dataset in orig space and ORIG view before using it as a master dataset. But that didn't solve the problem. :(

Here's the output from the nifti_tool:

N-1 header file 'sub-1005_task-aim_run-01_space-T1w_desc-preproc_bold.nii.gz', num_fields = 2
  name                offset  nvals  values
  ------------------- ------  -----  ------
  qform_code           252      1    2
  sform_code           254      1    2

N-1 header file 'sub-1005_desc-preproc_T1w.nii.gz', num_fields = 2
  name                offset  nvals  values
  ------------------- ------  -----  ------
  qform_code           252      1    1
  sform_code           254      1    1

N-1 header file 'sub-1005_space-T1w_desc-HandSeg_mask.nii.gz', num_fields = 2
  name                offset  nvals  values
  ------------------- ------  -----  ------
  qform_code           252      1    1
  sform_code           254      1    1

By the way, the T2s had coronal acquisition with slices perpendicular to the long axis of the hippocampus to enable manual segmentation of the hippocampus - that's why it is oblique? I'm not sure why the T1 is oblique. But the obliquity is also present in the raw, un-preprocessed T1.

Thanks!
Mrinmayi

ptaylor · September 17, 2024, 3:20pm

Hi, Mrinmayi-

OK, that makes sense. The issue with [qs]form_code values---I don't understand why the output, native space EPI is given a value of 2. It would be better to not have an ambiguous code value---here the anatomicals still have values of 1. Sigh.

In terms of how that is translated to either being interpreted as ORIG or TLRC, users can actually control this with an AFNI environment variable,

AFNI_NIFTI_VIEW

    The default view extension used for output when creating AFNI format
    datasets from NIFTI datasets. This variable is only applicable for
    sform and qform codes that do not have clearly defined views
    (sform/qform code = 2). Set to "tlrc" or "orig". See also
    AFNI_DEFAULT_STD_SPACE and AFNI_NIFTI_PRIORITY. Note sform/qform code=5
    can be used for spaces other than MNI or TLRC including MNI_ANAT or D99
    spaces.
    If unset, the default is 'orig'.

We have actually changed our default to be effectively AFNI_NIFTI_VIEW = ORIG because of this happening a lot from a couple different software tools. I think your version of AFNI is old enough that it defaults to -> TLRC.

However, while all of this is kind of inconvenient, it is a bit tangential then to the problem at hand. The main thing is that your FMRI and T1 datasets should be in the same space, however not only are they not on the same grid (which could be overcome or mapped between by simply resampling, as you started with) but there is an obliquity difference between them. It is a bit confusing to have a "processed", final subject-space dataset still having obliquity, but I guess that is what it is. (In the future, it might make sense to deoblique the anatomical before processing; leaving the EPI with obliquity is typically fine and recommended, but having obliquity in the anatomical can lead to minor annoyances.)

So, I think the following is the way to proceed. If the obliquity information were applied to the anatomical, then it should overlap well with the EPI dataset, and the only difference between the two would be the grid sampling. So, to get the ROI from the oblique-anatomical grid to the non-oblique EPI one, you can first estimate what that obliquity transform is (with 3dWarp), and then apply it to the ROI (with 3dAllineate):

# get obliquity difference matrix
3dWarp \
    -oblique_parent DSET_T1W \
    -disp_obl_xform_only \
    DSET_EPI \
    > mat_3dwarp_t1w_obl.1D

# apply the obliquity difference to the T1w-ROI, 
# and use EPI dset as master grid
3dAllineate \
    -1Dmatrix_apply  mat_3dwarp_t1w_obl.1D \
    -source          DSET_T1W_ROI \
    -master          DSET_EPI \
    -prefix          ROI_in_EPI.nii.gz

After than, the ROI should be on the EPI grid, which can be verified with:

3dinfo -same_all_grid DSET_EPI ROI_in_EPI.nii.gz

... and you can visualize the ROI_in_EPI on the DSET_EPI.

--pt

mrinmayi · September 18, 2024, 7:00pm

Hi Paul,

Thanks again for working through this! This works perfectly, and the mask lands up in the right place.

Just one note for anyone else who might have to do this - in my case, the values within the mask dataset are meaningful (e.g., voxel value 1 = left hippocampal head, 2 = left CA1, and so on). Running 3dAllineate interpolates the data so that the voxel values change. Adding -interp NN to the 3dAllineate command fixes this issue, such that the voxel values in the Allineate-d mask still have whole numbers!

Thanks,
Mrinmayi

ptaylor · September 18, 2024, 7:30pm

Hi, Mrinmayi-

Thanks for mentioning that, and glad it worked.

thanks,
pt

mrinmayi · September 30, 2024, 5:34pm

Hi Paul,

Just another question regarding using the mask. I have to extract data from "sub-ROIs" from the entire mask, based on voxel values in the mask. I've tried doing this by using the <> voxel selector in 3dmaskdump. I've also tried using 3dcalc. But in each case, some voxels that fit the criteria seem to be left out.

In the attached image, the values at the selected voxel can be seen in the purple box in left panel. The 2nd dataset here (sub-1004_space...) is the resampled dataset shown in red.
The 3rd dataset is created using this command:

# To "simulate" what is happening the in the 3dmaskdump command
3dcalc  \
      -a "sub-1004_space-T1w_desc-HandSeg_mask_Resampled.nii.gz<1>"  \
      -expr "a" \
      -prefix Area1Mask.nii.gz

The 4th dataset is created using:

3dcalc  \
      -a "sub-1004_space-T1w_desc-HandSeg_mask_Resampled.nii.gz"  \
      -expr "equals(a, 1)" \
      -prefix Area1Mask_Equals.nii.gz

But in both cases, voxels that have a value of 1 in the resampled dataset are being left out. Do you have ideas for why that might be happening? Could it have anything to do with the 3dAllineate command I'm using to de-oblique the mask to functional space?

Thanks!
Mrinmayi

ptaylor · September 30, 2024, 9:28pm

Hi, Mrinmayi-

I'm not sure how to interpret what that GUI interface is doing.

I would expect the two 3dcalc commands to give identical results for integer-valued datasets, which you do indeed have here. When I ran those on your data, and then checked the difference between the results with 3dDiff, I got that there is no difference:

$ 3dDiff -a Area1Mask_Equals.nii.gz -b Area1Mask.nii.gz 
++ Images do NOT differ

Similarly, if I subtract the two outputs, I get a dataset of all zeros:

# subtract the two dsets
$ 3dcalc -a Area1Mask.nii.gz -b Area1Mask_Equals.nii.gz -expr 'a-b' -prefix diff_areas.nii.gz
++ 3dcalc: AFNI version=AFNI_24.2.02 (Aug 16 2024) [64-bit]
++ Authored by: A cast of thousands
++ Output dataset ./diff_areas.nii.gz

# check max value of output: both min and max are 0, hence dset is all zeros
$ 3dinfo -min -max diff_areas.nii.gz 
0	0

When I visualize the results in the AFNI GUI, everything is consistent and matching, as well.

So, given that there is relative rotation between your overlay and underlays, I'm going to guess that the obliquity difference between those over/under datasets is responsible for some interpolation/smoothing happening. I'm not sure how to interpret all the smooth values in red, which lines and boundaries that don't appear to be locked to any voxel.

--pt

mrinmayi · October 1, 2024, 2:16pm

Hi Paul,

Sorry, my explanation wasn't clear at all!

The dataset I'm using in the 3dcalc commands (sub-1004_space-T1w_desc-HandSeg_mask_Resampled.nii.gz), is the one with the obliquity fixed, i.e., the output of the 3dAllineate command you had given me.

Now, I'm trying to extract data from the sub-ROIs within this mask. You're right that both 3dcalc commands yield the same result. But in both cases it's missing some voxels that fit the criteria (i.e., have a voxel value of 1).

In the ITK-SNAP screenshot, in the purple box in the left panel, the values indicate the value at the voxel under the crosshairs. It shows that the second dataset sub-1004_space-T1w_desc-HandSeg_mask_Resampled.nii.gz has a value of 1 in the selected voxel. So this voxel should've been included in the two 3dcalc outputs. But it is being left out, as suggested by the 0 voxel value for the third and fourth dataset. The smoothed red colour scheme is not meaningful - ITK-SNAP just doesn't have a discrete colour scheme suitable for an integer-valued, ROI type datasets!

So my question is, why do you think 3dcalc is dropping certain voxels even though they fit the criteria?

Thanks!
Mrinmayi

mrinmayi · October 1, 2024, 2:29pm

OH MY GOD! It seems the issue might be with ITK-SNAP visualisation, and that 3dcalc was working all along!!

I verified this by 3dmaskdump-ing all the values in the sub-1004_space-T1w_desc-HandSeg_mask_Resampled.nii.gz and counted how many voxels have a value of 1. And then I did the same thing with the 3dcalc output, and both datasets have 83 voxels.

This is what I get for straying from AFNI tools! So sorry to waste your time!
Mrinmayi

rickr · October 1, 2024, 2:30pm

Perhaps the ROI resampled via 3dAllineate above needs something like '-final NN', to preserve the ROI values. Without that, you might not really have values of 1 anymore.

-rick

ptaylor · October 1, 2024, 3:59pm

Howdy-

It's certainly not a waste of time, and no worries. But indeed, dealing with obliquity and visualization is tricky. Different software deal with it differently, and one has to try to make sure that things are being dealt with consistently. Visualizing data with differing obliquity is going to be particularly tricky: software have to choose whether (and how) to interpolate data, or to ignore obliquity.

Visualizing data is help for checking, and so when the visualization is not intuitive, that can make things tricky.

I'm glad things are making more sense now, or at least that there is confidence that the dataset calculations are consistent and accurate. When something unexpected comes up, one has to go back and double (or triple) check each step along the way, verifying it independently. One thing that can help, in addition to things above is this, to verify that grids really are the same if they are expected to be (or to highlight differences):

3dinfo -same_all_grid -prefix DSET1 DSET2 ...

As Rick noted, verifying that integer-valued datasets are still integer-valued is also helpful, particularly when transforming and regridding data. That can be done in the GUI perhaps most easily.

--pt

dglen · October 1, 2024, 5:55pm

I am going to add one or two more things here. From my interaction with ITK-snap, it assumes the same grid for the underlay/background as the overlay/segmentation. That can lead to obvious problems or more subtle problems like left-right flipping problems if you save segmentation from the GUI. In AFNI, the overlay is resampled to the underlay grid, and orientation and grid resolutions can be very different. I'm not sure how obliquity is handled in ITK-snap. It may require the same obliquity for all the involved datasets, but I can't be sure of that.

mrinmayi · October 2, 2024, 6:30pm

Yes, I think that's definitely what is happening - when I load the T1w image as the underlay, the overlay is resampled to match the grid of the T1, which leads to some voxels being interpolated in the mask, but are not captured in the 3dcalc output.

ITK-SNAP allows you to load dsets with different obliquity, but it always resamples the overlay to match the underlay. For instance, if an oblique image is the underlay, a non-oblique image added as an overlay will also look oblique, because it is resampled.

I suppose it's a good reminder to also use non-visualisation tools, where possible, to verify that steps are doing what you think they are!

Thanks so much to you all for your help
Mrinmayi