using NIfTI
ni = niread("my.nii")
ni = niread("my.nii.gz") # gzipped NIfTI files are detected automatically

ni = niread("my.nii", mmap=true)

vsize = voxel_size(ni.header)    # In mm
tr = time_step(ni.header)        # In ms

d = vox(ni, x, y, z, t)       # Scaled by slope and intercept given in header, zero-based indexes
d = ni[x, y, z, t]            # Scaled by slope and intercept given in header, one-based indexes
d = ni.raw[x, y, z, t]        # Unscaled, one-based indexes

niwrite("my.nii", ni)