Can 3dClusterize .HEAD-embedded cluster tables?

Given that afni_proc.py provides a procedure to embed cluster tables directly in images, I was curious if 3dClusterize is able to directly parse this table?

(on a related note, I've also read that the afni gui can read embedded cluster-tables. is there documentation on how to define clusters directly from the embedded table using the afni gui? the clusterize gui still requires one to manually identify cluster sizes)

Hi-

I'm not quite understanding what is being asked, and I might just be unaware of the desired functionality here. But I have a couple questions:

  • What is an example of this: Given that afni_proc.py provides a procedure to embed cluster tables directly in images? Does that refer to something in the APQC HTML file, or to something else?

  • I'm not sure what the goal of "reading" a cluster table into 3dClusterize or into the AFNI GUI would be; each of those things outputs clustermaps, as well as cluster tables (the Clusterize feature in the GUI calls 3dClusterize under the hood). I don't think there's enough information in the table itself to generate a cluster image.

  • Are you picturing the GUI/3dClusterize reading the comments at the top of the table output, to know the parameters to generate a 3dClusterize clustermap+table of its own?

  • I might be unaware of the clusterize functionality you want---could you please point me at the description of the GUI reading an embedded cluster table?

More generally, what is your desired input and desired output for this particular cluster related project? Do you want to regenerate a dataset from existing cluster tables? Or are you trying to generate images?

--pt

Thank you for the response, I will try to do my best to answer.

My apologies, that was admittedly worded badly. This refers to the cluster tables embedded in the stats image header using 3drefit. e.g.:

type = string-attribute
name = AFNI_CLUSTSIM_NN1_bisided
count = 3250
'<3dClustSim_NN1
  ni_type="10*float"
  ni_dimen="29"
  commandline="3dClustSim -both -mask full_mask.ab+tlrc -acf 0.880482 3.2657 15.6258 -cmd 3dClustSim.ACF.cmd -prefix files_ClustSim/ClustSim.ACF"

The idea I was thinking of: if a user provides an 'athr' probability and either provides a cluster-table or an image with a header-embedded table (e.g. as produced by 3drefit / afni_proc.py), then 3dClusterize could internally parse the desired cluster size.

I was thinking more that the clusterize gui could have an option to show the table if it is stored in the stats-image header.

For the other question, in the 3dClustSim documentation, I read:

* To add the cluster simulation C(p,alpha) table to the header of an AFNI
  dataset, something like the following can be done [tcsh syntax]:
     set fx = ( `3dFWHMx -detrend time_series_dataset+orig` )
     3dClustSim -mask mask+orig -acf $fx[5] $fx[6] $fx[7] -niml -prefix CStemp
     3drefit -atrstring AFNI_CLUSTSIM_NN1_1sided file:CStemp.NN1_1sided.niml \
             -atrstring AFNI_CLUSTSIM_MASK file:CStemp.mask    \
             statistics_dataset+orig
     rm -f CStemp.*
  AFNI's Clusterize GUI makes use of these attributes, if stored in a
  statistics dataset (e.g., something from 3dDeconvolve, 3dREMLfit, etc.).

   ** Nota Bene: afni_proc.py will automatically run 3dClustSim,  and **
  *** put the results  into the statistical results  dataset for you. ***
 **** Another reason to use afni_proc.py for single-subject analyses! ****

Looks like I inferred "reads table" from "makes use of". Still curious how I can use the header-embedded 3dClustSim attributes in the gui, though, this is more for my own curiosity!

Oh, now I get it---thanks for those additional details.

I was picturing the cluster table as the thing output by 3dClusterize, which looks like this:

#
#  Cluster report 
#[ Dataset prefix      = stats.FT ]
#[ Threshold vol       = [2] 'vis#0_Tstat' ]
#[ Supplement dat vol  = [1] 'vis#0_Coef' ]
#[ Option summary      = bisided,-2.556,2.556,clust_nvox,150,NN2 ]
#[ Threshold value(s)  = left-tail thr=-2.556000;  right-tail thr=2.556000 ]
#[ Aux. stat. info.    = DEGREES-of-FREEDOM : 412  ]
#[ Nvoxel threshold    = 150;  Volume threshold = 2343.750 ]
#[ Single voxel volume = 15.625 (microliters) ]
#[ Neighbor type, NN   = 2 ]
#[ Voxel datum type    = float ]
#[ Voxel dimensions    = 2.500 mm X 2.500 mm X 2.500 mm ]
#[ Coordinates Order   = RAI ]
#[ Mean and SEM based on signed voxel intensities ]
#
#Volume  CM RL  CM AP  CM IS  minRL  maxRL  minAP  maxAP  minIS  maxIS    Mean     SEM    Max Int  MI RL  MI AP  MI IS
#------  -----  -----  -----  -----  -----  -----  -----  -----  -----  -------  -------  -------  -----  -----  -----
  16378   -4.0   34.6   27.2  -78.8   68.8  -71.2   98.8  -21.2   66.2  -0.4271   0.0026  -3.9898    1.2   68.8   51.2 
  12901    1.1   42.9    2.8  -78.8   76.2  -78.8   98.8  -21.2   66.2   0.9981   0.0093     12.7   43.8  -46.2  -18.8 
    839  -47.2   -2.2   42.5  -66.2  -26.2  -31.2   16.2   16.2   66.2   0.5937   0.0196   5.8923  -36.2    8.8   66.2 
    421   46.1    8.7  -13.2   21.2   66.2  -18.8   48.8  -21.2    3.8  -0.4573   0.0111  -1.3876   53.8   18.8  -11.2 
    221    3.5  -28.5   48.1  -11.2   13.8  -53.8   -1.2   38.8   58.8   0.3041   0.0124   1.1872    8.8  -51.2   48.8 
    185  -73.6  -10.8    7.3  -78.8  -48.8  -18.8   11.2   -3.8   28.8   2.1604    0.122   7.2476  -78.8  -13.8    6.2 
#------  -----  -----  -----  -----  -----  -----  -----  -----  -----  -------  -------  -------  -----  -----  -----
# 30945   -2.7   37.5   11.9                                              0.215   0.0059                             

... not the cluster-simulation results from 3dClustSim, which looks like this:

# 3dClustSim -mask TMP_RES.nii.gz -acf 0.5128 3.8451 11.3839 -nodec -prefix TMP_RES_CLUST
# 2-sided thresholding
# Grid: 91x109x91 2.00x2.00x2.00 mm^3 (178964 voxels in mask)
#
# CLUSTER SIZE THRESHOLD(pthr,alpha) in Voxels
# -NN 3  | alpha = Prob(Cluster >= given size)
#  pthr  | .10000 .05000 .02000 .01000
# ------ | ------ ------ ------ ------
 0.050000     844   1016   1264   1466
 0.020000     317    386    486    574
 0.010000     172    212    269    311
 0.005000     104    126    161    190
 0.002000      57     70     90    107
 0.001000      39     47     61     74
 0.000500      26     33     43     50
 0.000200      16     21     28     33
 0.000100      11     15     20     24

Hence my confusion to the query. I now understand!

However, I think the answer is that there is not a way for 3dClusterize to directly take inputs from the cluster-simulation (clustsim) table. The reason is that the clustsim table contains a lot of outputs, and the user has to choose which pthr+alpha combination they want to use. That is, you must choose whether you want pthr=0.0005 and alpha=0.02, and get the threshold number sitting at that cross-section of the grid, or whether you want pthr=0.002 and alpha=0.10, or ... whatever. Because cluster simulations are so computationally expensive, but multiple pthr+alpha combos can be calculated simultaneously (even for multiple NN and sidedness, which are typically in separate files), it makes sense to output a table to choose from later---but that choice has to be made.

Now, that being said, it is possible to scriptically select your row/column entry from the clustsim table, such as doing this (NB: 1dcat ignores commented lines, and zerobased counting is used) to get the pthr=0.0005 entry from the 6th row and alpha=0.02 from the 3rd column:

1dcat  csim_output.NN3_2sided.1D'{6}[3]'

You could dump that to a variable in a script:

set holy_cluster_threshold = `1dcat csim_output.NN3_2sided.1D'{6}[3]'`

The downside of this kind of selection is that it is not by row/column label, so you have to be sure of that (probably also checking the information, to make sure the person providing the file didn't use different pthr/alpha values).

All of that is for when you have the file, and you asked about getting it from a file header. It's true that that table can be stored in a header, in an attribute typically linked to the particular NN and sidedness used in the 3dClustSim call, as you showed above; it can be dumped out with 3dAttribute:

3dAttribute AFNI_CLUSTSIM_NN1_1sided DSET_NAME

That displays the NIML-formatted version of the file. However, I'm not sure how to easily get at the table that is there. Perhaps someone else will, though.

--pt

Thank you for the information about the 1dcat solution! I figured there may not be a pre-packaged solution.

Guess for now I'll build my own parsers, haha!

(code below for reading an embedded table, in case a passerby is looking for a solution)

import xml.etree.ElementTree as ET

from math import floor

import nibabel as nib
import numpy as np
import pandas as pd

def _read_embedded_table(tbl_str: str) -> pd.DataFrame:
    """
    Read a cluster table embedded in an afni image.
    """
    ct_xml: xml.etree.ElementTree.Element = ET.fromstring(tbl_str.replace('3d', 'threed'))
    vox_probs = [float(n) for n in ct_xml.atrib['pthr'].split(',')]
    clust_probs = [float(n) for n in ct_xml.atrrib['athr'].split(',')]

    return pd.DataFrame(
        np.reshape(
            map(lambda x: floor(float(x)), ct_xml),
            (len(vox_probs), len(clust_probs))
        ),
        index=vox_probs,
        columns=clust_probs
    )
athr = 0.05
pthr = 0.05
NN = 1
clust_tail = 'bisided'
img = nib.load('stats.ab+tlrc.HEAD')
clust_tbl = _read_embedded_table(img.header.info[f'AFNI_CLUSTSIM_NN{NN}_{clust_tail}'])

floor(clust_tbl.loc[pthr, athr])
# e.g. 100

Cool, thanks. It is possible that there is an existing command line tool for extracting NIML file information from that attribute's clustsim table (@rickr ?), but I am not aware of one.

But that code looks like a nice way of doing it, with the benefit of specifying athr and pthr directly, too, thanks.

--pt

We do have 1d_tool.py to parse these tables, and an extracted table could be piped to it. But it might be preferable to take the direct tables as Paul suggests. That would be more clear in a script.

Consider sample 1d_tool.py commands like these. The "-verb 0" consideration allows you to just capture the single value text output, while -csim_pthr and -csim_alpha allow you to control the row and column from the table.

1d_tool.py -infile csim_output.NN3_2sided.1D -csim_show_clustsize
1d_tool.py -infile csim_output.NN3_2sided.1D -csim_show_clustsize -verb 0

1d_tool.py -infile csim_output.NN3_2sided.1D -csim_show_clustsize \
                  -csim_pthr 0.0001 -csim_alpha 0.01
1d_tool.py -infile csim_output.NN3_2sided.1D -csim_show_clustsize \
                  -csim_pthr 0.0001 -csim_alpha 0.01 -verb 0
  • rick