16.1.2. Datasets and Volumes¶

Everything volumetric in AFNI is a THD_3dim_dataset. This chapter covers its structure, how to read and load one, how to reach into its voxel data (sub-bricks, brick factors, datum types), how to move between voxel indices and millimeter coordinates, how to restrict work to a mask, and how to build and write a new dataset. The running example is again src/3dToyProg.c.

The two core structures ¶

`MRI_IMAGE` – one image / brick / vector ¶

MRI_IMAGE (declared in mrilib.h) is AFNI’s fundamental n-dimensional array. A 1D column, a 2D slice, a 3D volume, and a single dataset sub-brick are all MRI_IMAGEs. The fields you touch most:

typedef struct {
   int nx, ny, nz, nt, ... ;   /* dimensions (unused ones are 1) */
   int nxy ;                   /* nx*ny */
   int64_t nvox ;              /* total voxels */
   MRI_TYPE kind ;             /* datum: MRI_byte, MRI_short, MRI_float, ... */
   void *im ;                  /* the actual pixel data (cast by kind) */
   float dx, dy, dz, dt ;      /* physical voxel sizes */
   char *name ;                /* optional label / filename */
} MRI_IMAGE ;

The datum kind is one of the MRI_TYPE enum values (mrilib.h):

typedef enum MRI_TYPE {
   MRI_byte, MRI_short, MRI_int,
   MRI_float, MRI_double, MRI_complex, MRI_rgb, MRI_rgba, ...
} MRI_TYPE ;

There are helper accessors – MRI_BYTE_PTR(im), MRI_SHORT_PTR(im), MRI_FLOAT_PTR(im), MRI_COMPLEX_PTR(im) – that cast im->im to the right pointer type. MRI_TYPE_maxval[kind] gives the largest representable value for a type, which is handy when choosing a scale factor for short storage (3dToyProg uses 10.0/MRI_TYPE_maxval[MRI_short]).

A group of images is an MRI_IMARR (image array); the IMARR_SUBIM(imar,i) and IMARR_COUNT(imar) macros walk it.

`THD_3dim_dataset` – a volume plus its header ¶

A THD_3dim_dataset (declared in 3ddata.h) bundles the voxel data (its dblk datablock, a stack of sub-brick MRI_IMAGEs) with all of the header information: the spatial axes (daxes), the time axis (taxis if any), attributes, labels, and statistics. You almost never poke its fields directly – there is a DSET_* macro for nearly everything. The essentials:

Macro	Returns
`DSET_NX(ds)`, `DSET_NY(ds)`, `DSET_NZ(ds)`	number of voxels along i / j / k
`DSET_NVOX(ds)`	voxels per sub-brick (`nxnynz`)
`DSET_NVALS(ds)`	number of sub-bricks (values per voxel)
`DSET_NUM_TIMES(ds)`	number of time points (1 if not a time series)
`DSET_TR(ds)`	TR / time step
`DSET_PREFIX(ds)`	the output prefix string
`DSET_HEADNAME(ds)`	full path of the `.HEAD` file
`DSET_BRICK(ds,iv)`	the `iv`-th sub-brick as an `MRI_IMAGE *`
`DSET_ARRAY(ds,iv)`	raw `void *` data pointer of sub-brick `iv`
`DSET_BRICK_TYPE(ds,iv)`	`MRI_TYPE` of sub-brick `iv`
`DSET_BRICK_FACTOR(ds,iv)`	scale factor of sub-brick `iv` (0 means 1.0)

Note the axis ordering convention: datasets are stored space-first (x fastest, then y, then z, then sub-brick), whereas the MRI_vectim structure discussed in The Function Cookbook is time-first. Keep that straight when you convert between them.

Reading a dataset from disk ¶

Two steps: open the header, then load the voxel data.

THD_3dim_dataset *dset = THD_open_dataset( fname ) ;   /* header only */
if( dset == NULL ) ERROR_exit("Cannot open %s", fname) ;
DSET_mallocize(dset) ;    /* force malloc (vs mmap) so you can modify */
DSET_load(dset) ;         /* actually read the bricks into memory */

THD_open_dataset() accepts every AFNI-recognized name form: dset+orig, dset.nii.gz, dset.HEAD, a NIFTI file, and even sub-brick / index selectors like dset+orig'[0..3]' or dset+orig'[2,5]'. It reads only the header.
DSET_load(dset) is a thin macro over THD_load_datablock that reads (or memory-maps) the actual voxel values. Until you call it, DSET_ARRAY returns NULL.
DSET_mallocize(dset) requests that the data be malloc-ed rather than mmap-ed; do this if you intend to modify the voxels in place.
When you are done, DSET_unload(dset) frees the voxel data but keeps the header; DSET_delete(dset) frees everything.

This exact open/mallocize/load sequence is what src/3dToyProg.c uses for its -input dataset, and what essentially every 3d* program uses.

Reaching into the voxel data ¶

Sub-bricks, datum types, and brick factors ¶

A dataset holds DSET_NVALS(ds) sub-bricks, each a separately typed and separately scaled volume. A sub-brick can be MRI_byte, MRI_short, MRI_float (etc.), and it may carry a brick factor: the true value of a voxel is stored_value * DSET_BRICK_FACTOR(ds,iv) (a factor of 0 is treated as 1.0). This is how AFNI stores, say, correlation coefficients compactly as scaled shorts.

Because of the per-brick type and factor, the safe way to read voxel values is almost never to walk DSET_ARRAY yourself. Instead, ask for the data as scaled floats:

/* one whole sub-brick as float[nvox], factor already applied: */
float *sb = THD_extract_to_float( iv, dset ) ;
...
free(sb) ;

THD_extract_to_float() allocates DSET_NVOX(dset) floats, applies the brick factor, and hands you a clean float volume regardless of on-disk type. This is the pattern in Volumewise_Operations() of src/3dToyProg.c, which accumulates a sum and a transform across all sub-bricks.

If you really do want the raw pointer (for speed, when you know the type), use DSET_ARRAY(ds,iv) and remember to honor DSET_BRICK_FACTOR yourself.

Extracting a voxel time series ¶

The complement of “one brick across all voxels” is “one voxel across all bricks” – i.e. a time series. Given a 1D voxel index nijk:

float *far = (float *)calloc( DSET_NVALS(dset), sizeof(float) ) ;
if( THD_extract_float_array( nijk, dset, far ) == -1 )
   ERROR_message("bad voxel %d", nijk) ;

or, allocating an MRI_IMAGE for you:

MRI_IMAGE *tsim = THD_extract_series( nijk, dset, 0 ) ;   /* 0 = scaled */
float     *ts   = MRI_FLOAT_PTR(tsim) ;

THD_extract_series() (the third argument is a “raw/unscaled” flag) is one of the most widely used routines in AFNI – 3dTcorrelate, 3dDetrend, 3dDeconvolve, 3dTshift, 3dToutcount, 3ddelay and many others pull their per-voxel time courses this way. For pushing a time series into a dataset, use THD_insert_series() (see below).

To read the same set of voxels from many indices efficiently, there is THD_extract_many_arrays; to convert an entire masked dataset to a time-first MRI_vectim in one call, use THD_dset_to_vectim (see The Function Cookbook).

Navigating: voxel indices and mm coordinates ¶

AFNI datasets know their real-world geometry. The Dataset_Navigation() function in src/3dToyProg.c is a compact tour; the key ideas:

1D index <-> 3D index. A whole volume is one flat array of length DSET_NVOX. Convert with:

int ni  = DSET_NX(dset) ;
int nij = DSET_NX(dset) * DSET_NY(dset) ;
int nijk = AFNI_3D_to_1D_index(i, j, k, ni, nij) ;   /* (i,j,k) -> flat */
AFNI_1D_to_3D_index(nijk, i, j, k, ni, nij) ;         /* flat -> (i,j,k) */

Voxel (i,j,k) <-> millimeter (x,y,z). The dataset’s daxes->ijk_to_dicom_real is the 3x4 affine that maps voxel indices to RAI-DICOM millimeters. Load it into an augmented 4x4 and multiply:

float A[4][4], Ai[4][4], I[3], X[3] ;
MAT44_TO_AFF44(A, dset->daxes->ijk_to_dicom_real) ;
I[0]=DSET_NX(dset)/2; I[1]=DSET_NY(dset)/2; I[2]=DSET_NZ(dset)/2 ;
AFF44_MULT_I(X, A, I) ;    /* X = A I  : voxel -> mm RAI */

AFF44_INV(Ai, A) ;          /* invert */
AFF44_MULT_I(I, Ai, X) ;    /* I = Ai X : mm -> voxel */

RAI-DICOM orientation means: x increases Right->Left, y Anterior->Posterior, z Inferior->Superior (so Right, Anterior, and Inferior are the negative directions). Datasets may also be oblique; the affine handles that transparently, which is why you should prefer the matrix over hand-rolled origin + index*delta arithmetic.

Masks ¶

Most analyses restrict computation to brain voxels. A mask is just a byte * of length DSET_NVOX with 1 inside and 0 outside.

From a user-supplied mask dataset (the -mask idiom in src/3dToyProg.c):

THD_3dim_dataset *mset = THD_open_dataset( argv[++iarg] ) ;
if( THD_dataset_mismatch(mset, iset) )     /* same grid as data? */
   ERROR_exit("grid mismatch between input and mask") ;
byte *mask  = THD_makemask( mset, 0, 1.0, -1.0 ) ;  /* brick 0, nonzero */
int   ncount = THD_countmask( DSET_NVOX(mset), mask ) ;
if( ncount <= 0 ) ERROR_exit("empty mask") ;
DSET_delete(mset) ;                                  /* done with it */

THD_makemask(dset, ivol, bot, top) marks voxels whose value is in the [bot,top] range (the bot > top form used above means “any nonzero value”). THD_countmask reports how many survived.

Auto-generated from the data with THD_automask(dset) – the same brain-extraction logic as the 3dAutomask program, reused by dozens of tools (3dAllineate, 3dBandpass, 3dBlurToFWHM, 3dTcorrelate …).

Creating a dataset from scratch ¶

New_Dataset_From_Scratch() in src/3dToyProg.c shows the full pattern. Build the grid, declare the axes, allocate the bricks, fill them.

1. The grid. EDIT_geometry_constructor() builds an empty dataset from a compact geometry string:

/* RAI:D: nx,xorg,dx, ny,yorg,dy, nz,zorg,dz  */
THD_3dim_dataset *oset =
   EDIT_geometry_constructor("RAI:D:32,-12.7,6.0, 32,2.3,6.0, 10,-4.5,5.2",
                             prefix) ;

For an oblique grid you can instead pass a full 3x4 index-to-mm matrix: "MATRIX(a11,...,a34):nx,ny,nz".

2. Declare what the dataset holds with the variadic EDIT_dset_items(). It takes ADN_* key / value pairs and must end with ADN_none:

EDIT_dset_items( oset ,
   ADN_datum_all , MRI_short ,      /* storage type of all bricks */
   ADN_nvals     , 100 ,            /* number of sub-bricks       */
   ADN_ntt       , 100 ,            /* number of time points      */
   ADN_ttdel     , 2.0 ,            /* TR                         */
   ADN_tunits    , UNITS_SEC_TYPE ,
   ADN_prefix    , "my_output" ,
   ADN_none ) ;

Common keys: ADN_prefix, ADN_datum_all, ADN_nvals, ADN_ntt (time points), ADN_ttdel (TR), ADN_tunits, ADN_type / ADN_func_type. If you set ADN_ntt to 0 the dataset is a bucket (no time axis) rather than a time series.

3. Allocate and fill the bricks.

for( tt=0 ; tt < nvals ; tt++ ){
   EDIT_substitute_brick(oset, tt, MRI_short, NULL) ;  /* NULL => allocate */
   EDIT_BRICK_FACTOR(oset, tt, 10.0/MRI_TYPE_maxval[MRI_short]) ;
}
/* push a float time series into flat voxel nijk, converting to storage type */
THD_insert_series(nijk, oset, nvals, MRI_float, ts, 0) ;

EDIT_substitute_brick(ds, iv, type, ptr) installs a data pointer as sub-brick iv (pass NULL to have it allocated for you). The final argument of THD_insert_series() is a “already-scaled” flag: 0 means “scale my floats to the brick’s storage type”, which is what you want when the target is a scaled short.

Deriving a new dataset from an existing one ¶

More often you already have an input dataset and want an output on the same grid. EDIT_empty_copy clones the header (geometry, orientation) but no data, then you attach your results. From Volumewise_Operations() in src/3dToyProg.c:

THD_3dim_dataset *oset = EDIT_empty_copy( dset ) ;
EDIT_dset_items( oset ,
   ADN_prefix , outname ,
   ADN_datum_all , MRI_float ,
   ADN_nvals , 2 , ADN_ntt , 0 ,      /* 2-brick bucket */
   ADN_none ) ;

/* attach already-computed float volumes directly (no copy, no scaling): */
EDIT_substitute_brick( oset, 0, MRI_float, vout ) ; vout = NULL ;
EDIT_substitute_brick( oset, 1, MRI_float, yout ) ; yout = NULL ;

/* or store as scaled shorts, letting AFNI pick the factor (-1.0): */
EDIT_substscale_brick( oset, 0, MRI_float, vout, MRI_short, -1.0 ) ;
free(vout) ;

Two important subtleties, both called out in the toy program’s comments:

EDIT_substitute_brick takes ownership of your pointer when the types match – do not free it afterward (set your copy to NULL).
EDIT_substscale_brick copies the data (converting/scaling into the storage type), so you must free your original buffer. Passing a scale of -1.0 asks AFNI to compute a sensible auto-scale factor.

Label your outputs and, if useful, attach custom header attributes:

EDIT_BRICK_LABEL(oset, 0, "sum") ;
EDIT_BRICK_LABEL(oset, 1, "sisqrt") ;

THD_set_string_atr( oset->dblk, "Pinot_Noir", "sideways" ) ;
THD_set_int_atr   ( oset->dblk, "Lottery_Numbers", 6, iv6 ) ;
THD_set_float_atr ( oset->dblk, "PI", 1, &fv ) ;

(Custom “private” attributes are not preserved when the dataset is later copied by other programs – use them for your own round-trips, not as a long-term storage contract.)

The voxel-wise convenience engine ¶

For the very common “run a function on every voxel’s time series and write the result” pattern, you do not need to write the triple loop yourself. MAKER_4D_to_typed_fbuc (and its siblings) walk the dataset, hand your callback each voxel’s (optionally detrended) time series, and assemble the output bucket. From Voxelwise_Operations() in src/3dToyProg.c:

oset = MAKER_4D_to_typed_fbuc(
         dset,               /* input dataset                       */
         outprefix,          /* output prefix                       */
         MRI_short,          /* output datum                        */
         0,                  /* time points to skip at start        */
         1,                  /* 1 => pass a linearly detrended series*/
         2,                  /* number of output values per voxel   */
         (generic_func *)toy_tsfunc,  /* your per-voxel callback     */
         (void *)&ud,        /* user data handed to the callback    */
         voxmask,            /* byte mask (NULL = all voxels)       */
         0 ) ;               /* allow auto-scaling of output        */

Your callback has a fixed prototype and is also called once at the start and once at the end (with val == NULL) so you can allocate and free scratch space:

static void toy_tsfunc( double tzero, double tdelta,
                        int npts, float ts[],
                        double ts_mean, double ts_slope,
                        void *ud, int nbriks, float *val )
{
   if( val == NULL ){                 /* start (npts>0) / end (npts<=0) */
      ...allocate or free scratch...
      return ;
   }
   /* ts[0..npts-1] is this voxel's series; write nbriks results to val[] */
}

This is the same engine behind 3dTfilter, 3dTnorm, 3dPeriodogram, 3dDWItoDT, 3dRetinoPhase and others – a good thing to reach for before hand-writing voxel loops.

Writing the result ¶

Finally, stamp history and write (with the overwrite guard from Anatomy of an AFNI Program):

tross_Copy_History( iset, oset ) ;
tross_Make_History( "3dProgName", argc, argv, oset ) ;
if( !THD_ok_overwrite() && THD_is_ondisk(DSET_HEADNAME(oset)) )
   ERROR_message("Output %s already exists", DSET_HEADNAME(oset)) ;
else
   DSET_write(oset) ;
DSET_delete(oset) ;

The output format follows the prefix: foo+orig writes BRIK/HEAD, foo.nii.gz writes gzipped NIFTI, and so on – you do not special-case it.