HEAD/BRIK into Matlab (or how to check ROIs against atlas)

Pawel · April 14, 2022, 12:33am

Hello all,

I have an (NHP) cluster map aligned to NMT, which I would like to cross check against the ARM atlas, to produce a report like:
cluster 1: 90% belongs to atlas area 75, 5% to atlas area 133, 5% to nowhere
cluster 2: 14% belongs to atlas area 22, 86 % to nowhere
and then
atlas area 55: 45% in cluster 1, 55% in cluster 3
atlas area 72: 12% in cluster 1, 88% in no cluster
not necessarily exactly in this form, of course. Atlas area names might be included as well.

Is there an AFNI program that would do that for me, ideally with a lot of control over output format?

Alternatively (preferably?) I thought I could just code it in Matlab. But is there a better way to load HEAD/BRIK into Matlab than converting to nifti with 3dAFNItoNIFTI, followed by niftiread in Matlab? I googled up this: https://sscc.nimh.nih.gov/pub/dist/doc/OLD/afni_matlab.shtml but it seems to be old (possibly outdated?) and the link does not work anyway.

Or, if you have other suggestions how to achieve my goal, I’d be happy to hear them.

Two caveats:
If own coding is needed, I would rather stay in Matlab: I know very little of Python and none of R, and coding anything slightly complex in tcsh is not my favorite way of spending time
My Matlab and my AFNI live on the same computer but in different OSes, so I would prefer to avoid mixing up these two, like calling AFNI programs from Matlab script.

ptaylor · April 14, 2022, 1:08pm

Pawel-

I think ‘whereami’ is what you want:
https://afni.nimh.nih.gov/pub/dist/doc/htmldoc/programs/whereami_sphx.html#ahelp-whereami
Don’t think there will be a lot of control over output format.

I don’t know much about Matlab, but there should be a function called “BrikLoad()” for reading in briks. This path seems a bit less ideal, to me.

–pt

Pawel · April 14, 2022, 3:11pm

Thank you!

I was able to get the first answer (what part of the ROI belongs to this and that atlas area) bu using whoami with -bmask or -omask, even if the output is not very machine readable.


Intersection of ROI (all non-zero values) with atlas ARM (sb3):
   31.0 % overlap with CR_area_24, code 506
   29.7 % overlap with CR_rostral_dlPFC
[CR_rostral_dorsolateral_prefrontal_cortex], code 559
   12.1 % overlap with CR_area_8B, code 555
   9.9  % overlap with CR_SMA/preSMA
[CR_medial_supplementary_motor_areas], code 587
   4.2  % overlap with CR_MCC
[CR_midcingulate_cortex], code 511
   3.2  % overlap with CR_rostral_vlPFC
[CR_rostral_ventrolateral_prefrontal_cortex], code 565
   1.3  % overlap with CR_area_32, code 504
   0.1  % overlap with CR_area_10
[CR_frontal_pole], code 518

I can’t find a way though to do the opposite, i.e. for each atlas area that intersects with the ROI, what percentage of that atlas area is covered by the ROI. Can this be done with whereami?

Also am I right that there is no way to limit this to specific subbricks of a multi-brick atlas like ARM?

Thank you for telling me a BrikLoad(), I found an old copy online and I’ll experiment with it in case I decide to go the Matlab way after all.

ptaylor · April 14, 2022, 3:17pm

Hi, Pawel-

Re. subbricks: You should be able to use AFNI subbrick selectors:


INPUT DATASET NAMES
-------------------
 An input dataset is specified using one of these forms:
    'prefix+view', 'prefix+view.HEAD', or 'prefix+view.BRIK'.
 You can also add a sub-brick selection list after the end of the
 dataset name.  This allows only a subset of the sub-bricks to be
 read in (by default, all of a dataset's sub-bricks are input).
 A sub-brick selection list looks like one of the following forms:
   fred+orig[5]                     ==> use only sub-brick #5
   fred+orig[5,9,17]                ==> use #5, #9, and #17
   fred+orig[5..8]     or [5-8]     ==> use #5, #6, #7, and #8
   fred+orig[5..13(2)] or [5-13(2)] ==> use #5, #7, #9, #11, and #13
 Sub-brick indexes start at 0.  You can use the character '$'
 to indicate the last sub-brick in a dataset; for example, you
 can select every third sub-brick by using the selection list
   fred+orig[0..$(3)]

 N.B.: The sub-bricks are read in the order specified, which may
 not be the order in the original dataset.  For example, using
   fred+orig[0..$(2),1..$(2)]
 will cause the sub-bricks in fred+orig to be input into memory
 in an interleaved fashion.  Using
   fred+orig[$..0]
 will reverse the order of the sub-bricks.

 N.B.: You may also use the syntax <a..b> after the name of an input 
 dataset to restrict the range of values read in to the numerical
 values in a..b, inclusive.  For example,
    fred+orig[5..7]<100..200>
 creates a 3 sub-brick dataset with values less than 100 or
 greater than 200 from the original set to zero.
 If you use the <> sub-range selection without the [] sub-brick
 selection, it is the same as if you had put [0..$] in front of
 the sub-range selection.

 N.B.: Datasets using sub-brick/sub-range selectors are treated as:
  - 3D+time if the dataset is 3D+time and more than 1 brick is chosen
  - otherwise, as bucket datasets (-abuc or -fbuc)
    (in particular, fico, fitt, etc datasets are converted to fbuc!)

 N.B.: The characters '$ ( ) [ ] < >'  are special to the shell,
 so you will have to escape them.  This is most easily done by
 putting the entire dataset plus selection list inside forward
 single quotes, as in 'fred+orig[5..7,9]', or double quotes "x".

Pawel · April 14, 2022, 3:23pm

I tried standard subbrick selector and it did not work

whereami -atlas ARM[3] -bmask Clust_mask_0001+tlrc -tab
++ Input coordinates orientation set by default rules to RAI
** ERROR: Atlas ARM[3] not found in list.

But this did!

whereami -atlas ‘ARM Level 3’ -bmask Clust_mask_0001+tlrc -tab

Sorry I did not figure it out myself

ptaylor · April 14, 2022, 5:50pm

Did you put quotes around the selector? In many shells, you have to put quotes around the brackets so the shell doesn’t try to interpret them in the way it normally would:


whereami -atlas ARM"[3]" -bmask Clust_mask_0001+tlrc -tab

–pt

Pawel · April 14, 2022, 5:56pm

That does not work either for me:


whereami -atlas ARM"[3]" -bmask Clust_mask_0001+tlrc
++ Input coordinates orientation set by default rules to RAI
** ERROR: Atlas ARM[3] not found in list.

but again ARM seems to register each level separately in afni so


whereami -atlas 'ARM Level 3' -bmask Clust_mask_0001+tlrc

works fine.

I am now seeing some discrepancies between what whereami shows vs. my Matlab script - which I wrote anyway, to get the percentage of atlas regions covered by each ROI. Most likely I made a mistake in the script, but if I cannot reconcile the results, I may come back for help…

Thanks again,

Pawel

ptaylor · April 14, 2022, 6:07pm

Fascinating. I wonder if the atlas points info is just for the whole dset, and not subbrick by subbrick? That is my suspicion…

I would really be careful with matlab scripting—definitely check orientation/grid matching between dsets.

–pt

Pawel · April 14, 2022, 6:32pm

Per 3dinfo, both the atlas file and the ROI file are of the same geometry and orientation. Which is expected, as alignment happened on the way here, but I checked.

Both report:

Geometry String: "MATRIX(-0.5,0,0,31.75,0,-0.5,0,27.625,0,0,0.5,-7.875):128,156,100"

and

[-orient LPI]

But it really looks to me as if whereami was slightly off… I tried to replicate what (I think) whereami is doing using individual commands, and it produces the same result as my Matlab script, and different from whereami. Which might actually mean that I am missing something about how whereami works. Could you please have a look?

I am focusing on one ROI and two areas of ARM level 3 (=subbrick 2)

Here is what my Matlab code says:


LesionAtlas(cd(), 3, 1)
ROI 1, region 503, overlap 1923, per ROI (size: 5519) = 0.348, per region (size: 3319) = 0.579
ROI 1, region 554, overlap 2496, per ROI (size: 5519) = 0.452, per region (size: 7020) = 0.356
[...]

So 34.8% of ROI 1 is in atlas area 503, and 45.2% of ROI 1 is in area 554. cd() is where to look for files (=pwd), 3 is the atlas level, 1 means use individual ROI files (basically = -bmask)

whereami says:


whereami -atlas 'ARM Level 3' -bmask Clust_mask_0001+tlrc
++ Input coordinates orientation set by default rules to RAI
++ In binary mode ...
++    5519 voxels in ROI
++    43959 voxels in atlas-resampled mask
Intersection of ROI (all non-zero values) with atlas ARM Level 3 (sb0):
   41.8 % overlap with CR_dlPFC
[CR_dorsolateral_prefrontal_cortex], code 554
   32.3 % overlap with CR_ACC
[CR_anterior_cingulate_cortex], code 503
[...]
   -----
   91.5 % of cluster accounted for.

I tried to replicate by:

using 3dcalc to extract areas 503 and 554 from the atlas’ subbrick 2


3dcalc -a ARM_in_NMT_v2.1_sym_05mm.nii.gz[2] -exp 'equals(a, 554)' -prefix area554 -overwrite
3dcalc -a ARM_in_NMT_v2.1_sym_05mm.nii.gz[2] -exp 'equals(a, 503)' -prefix area503 -overwrite

calculating intersection of ROI 1 with each of the extracted areas


3dcalc -a Clust_mask_0001+tlrc -b area554+tlrc -exp 'and(a,b)' -overwrite -prefix intersect554
3dcalc -a Clust_mask_0001+tlrc -b area503+tlrc -exp 'and(a,b)' -overwrite -prefix intersect503

counting the number of voxels in each intersection and in the ROI


3dmaskdump -mask intersect554+tlrc -o /dev/null intersect554+tlrc
++ 3dmaskdump: AFNI version=AFNI_21.3.11 (Dec  9 2021) [64-bit]
++ 2496 voxels in the mask
3dmaskdump -mask intersect503+tlrc -o /dev/null intersect503+tlrc
++ 3dmaskdump: AFNI version=AFNI_21.3.11 (Dec  9 2021) [64-bit]
++ 1923 voxels in the mask
3dmaskdump  -mask Clust_mask_0001+tlrc  -o /dev/null Clust_mask_0001+tlrc
++ 3dmaskdump: AFNI version=AFNI_21.3.11 (Dec  9 2021) [64-bit]
++ 5519 voxels in the mask

dividing them accordingly


awk "BEGIN {print 2496/5519}"
0.452256
awk "BEGIN {print 1923/5519}"
0.348433

So that’s what Matlab says not what whereami says. Am I missing a step, both in Matlab and in my replication, that whereami does?

One hint is that whereami reports “43959 voxels in atlas-resampled mask”. My ROIs come from thresholding of data that had been already masked using GM masks (subcortical and cortical) that come with the ARM and NMT2.1 (NMT_v2.1_sym_05mm_segmentation.nii.gz) but maybe there is some discrepancy there? Like the extent of the segmentation masks differs from the extent of the atlas?

ptaylor · April 14, 2022, 7:51pm

++ In binary mode … You are in “binary” mode—all non-zero voxels are treated like one giant ROI.


-bmask MASK: Report on the overlap of all non-zero voxels in MASK dataset
              with various atlas regions. NOTE: The mask itself is not binary,
              the masking operation results in a binary mask.
 -omask ORDERED_MASK:Report on the overlap of each ROI formed by an integral 
                     value in ORDERED_MASK. For example, if ORDERED_MASK has 
                     ROIs with values 1, 2, and 3, then you'll get three 
                     reports, one for each ROI value. Note that -omask and
                     -bmask are mutually exclusive.

Pawel · April 14, 2022, 8:09pm

I know, but the ROI file I am using (Clust_mask_0001+tlrc) contains only one ROI, it was made by [Write]'ing one specific cluster from GUI’s clustering report.
Either way, if I use the ordered mask (-omask) mode on this single-ROI file, or ordered mask on a multi-ROI file, I get the same numbers on this ROI and regions: 41.8% and 32.3%, vs. 45.2% and 34.8% from doing it step by step or in Matlab


whereami -atlas 'ARM Level 3' -omask Clust_mask_0001+tlrc
++ Input coordinates orientation set by default rules to RAI
++ In ordered mode ...
++ Have 2 unique values of:
   0   1
++ Skipping unique value of 0
++ ========================================================================
++ Processing unique value of 1
++    5519 voxels in ROI
++    43959 voxels in atlas-resampled mask
Intersection of ROI (valued 1) with atlas ARM Level 3 (sb0):
   41.8 % overlap with CR_dlPFC
[CR_dorsolateral_prefrontal_cortex], code 554
   32.3 % overlap with CR_ACC
[CR_anterior_cingulate_cortex], code 503
[...]
   91.5 % of cluster accounted for.


whereami -atlas 'ARM Level 3' -omask Clust_mask+tlrc
++ Input coordinates orientation set by default rules to RAI
++ In ordered mode ...
++ Have 4 unique values of:
   0   1   2   3
++ Skipping unique value of 0
++ ========================================================================
++ Processing unique value of 1
++    5519 voxels in ROI
++    43959 voxels in atlas-resampled mask
Intersection of ROI (valued 1) with atlas ARM Level 3 (sb0):
   41.8 % overlap with CR_dlPFC
[CR_dorsolateral_prefrontal_cortex], code 554
   32.3 % overlap with CR_ACC
[CR_anterior_cingulate_cortex], code 503
 [...]
   -----
   91.5 % of cluster accounted for.

dglen · April 16, 2022, 2:19am

This may have to do with differences in the resampling/interpolation to the template. Try resampling first with 3dresample to the grid of the template and nearest neighbor interpolation.

Pawel · April 16, 2022, 3:01am

Daniel Glen Wrote:

This may have to do with differences in the
resampling/interpolation to the template. Try
resampling first with 3dresample to the grid of
the template and nearest neighbor interpolation.

Hi Daniel,

resampling what to the grid of the template? I think that all my files at this point are on the same grid.(or is there anything more to the grid than the Geometry String, Data Axes Orientation, and Data Axes Tilt?)

I start with a structural T2, which I align to the template. After this point, as far as I can tell, the aligned T2 is on the same grid as the template, which is the same grid as the atlas. (by the way, there is no intermediate T1 anymore, but that’s a different story)

Next I mask the T2 (same grid), and then I threshold and clusterize it in GUI, and save cluster maps (same grid).

Then I calculate the overlaps between cluster maps and atlas ROIs, and I tried to do it in three ways: via whereami, using afni programs to isolate “my” and atlas ROIs and find intersections, and with a Matlab script that does the same (after cluster maps are converted to NIfTI via AFNI program, as Matlab easily loads NIfTI).

And whereami gives slightly different results than the other two. I agree that the magnitude of the difference suggests something like a side effect of resampling, but I do not understand what and why would be resampled, if the intersected files are on the same grid. Or - again - I think they are.

I should also say that yesterday Paul asked me offline to provide the files, which I did. So maybe he will be able to figure what’s going on.

ptaylor · April 16, 2022, 3:50am

Hi, Pawel-

I see you are selecting “ARM Level 3” in your selection, but the parenthetical is “(sb0)”, which means “sub-brick 0”. Did you splice the dataset?

When I run:


whereami -atlas 'ARM' -omask Clust_mask_0001+tlrc > asdf

… with the data you provided, this is what I get for sub-brick 2 (sb2) information, which corresponds to Level 3 (because subbricks are zero-based counting):


Intersection of ROI (valued 1) with atlas ARM (sb2):
   45.2 % overlap with CR_dlPFC
[CR_dorsolateral_prefrontal_cortex], code 554
   34.8 % overlap with CR_ACC
[CR_anterior_cingulate_cortex], code 503
   11.1 % overlap with CR_SMA/preSMA
[CR_medial_supplementary_motor_areas], code 587
   4.6  % overlap with CR_MCC
[CR_midcingulate_cortex], code 511
   3.1  % overlap with CR_vlPFC
[CR_ventrolateral_prefrontal_cortex], code 564
   0.1  % overlap with CR_med_OFC
[CR_medial_orbital_frontal_cortex], code 517
   -----
   98.9 % of cluster accounted for.

This whereami output seems consistent with what you did in Matlab.

So, perhaps the question is what happened to the ARM dataset before you ran whereami?

–pt

Pawel · April 16, 2022, 5:14am

Hi Paul,

I think I know what happened.

I did not split the file, it had been split by ARM/NMT creators. The split versions are registered with AFNI as single-subbrick “ARM Level x” atlases, along with the full 6-subbrick “ARM”. Probably to facilitate subbrick selection which - as we have found out - does not seem to work in standard way with whereami.

But here is the actual reason which, as Daniel suggested, was caused by resampling.

NMT2.1 and ARM come in three versions which differ in resolution and/or FOV. The version that registers with afni by running a script provided with NMT/ARM and is consequently used by whereami, is different from the one I was using as the basis of the alignment and then as the atlas. The script registers the large-FOV/0.25 mm version while I am using the tight-FOV/0.5 mm version.

So “my” whereami was resampling my cluster to match the high-res/large FOV version, while the Matlab script directly used the low-res/tight FOV version, also the version whose grid had been imposed on the data. My step-by-step process also explicitly used the low-res version. And you probably registered the 0.5 mm version that I sent you with your AFNI, and your whereami was using it.

Thank you for your help and patience! The good news is that my Matlab script is accurate - ensuring that was the only reason why I compared it to whereami. I will have to use it instead of whereami because I need the percentages of atlas ROIs intersecting with my ROIs, not the other way around.

ptaylor · April 16, 2022, 11:38am

Great, that makes sense. Thanks for checking to verify, and glad we got to the bottom of the mystery (esp. without any bugs lurking there…).

–pt

Pawel · April 16, 2022, 3:22pm

In hindsight, I had a very good hint here:


++    5519 voxels in ROI
++    43959 voxels in atlas-resampled mask

I think I did not realize that ROI and mask refer to basically the same thing, and that the ratio of the two numbers is almost exactly 8, which strongly suggests resampling to twice the resolution.

I think I am going to register all three ARM atlases to my AFNI and, to avoid confusion, give them meaningful names like ARMfh025, ARM05, and ARM025. Am I correct that this (including unregistering the ARM name for the large-FOV/0.25 mm atlas) can be done by editing the CustomAtlases.niml file?

dglen · April 18, 2022, 2:43am

I’m glad you got this all figured out. The names for the supplementary atlases come exclusively from that CustomAtlases.niml file, so editing there should work. SessionAtlases.niml does something similar for programs to read the current directory. The atlases that are distributed are defined in the similar AFNI_atlas_spaces.niml.