Nibabel: Cifti XML error when creating CIFTI with multiple subcortical structures

Software Support
,
1

Summary of what happened:

I’m trying to create a CIFTI (dscalar.nii) with left+right vertex cortex and multiple subcortical voxel structures to save processed output from FSL’s probtrackx as CIFTI.

I can successfully create CIFTIs with left+right surfaces and a single voxel structure (e.g. thalamus_left) but when I try to create with many subcortical structures (e.g. thalamus_left, amygdala_left, etc…), it saves but workbench gives errors (copied below) when trying to read.

I’m trying to generate the CIFTI brain model using the .from_mask function, providing the subcortical mask each time.

The error suggests that there is overlap in voxels across the subcortical structures. The subcortical masks (i.e. entries of 1) do not overlap - I’ve checked this by comparing both the coordinates listed by probtrackx and visually.

I’m guessing there is a neater (and correct) way of creating the brain model. Any help would be appreciated.

Command used (and if a helper script was used, a link to the helper script or the command generated):

import nibabel as nib
from nibabel import cifti2
from fsl.data.image import Image
from fsl.data.cifti import cifti2_axes
from fsl.data.cifti import Cifti


# set up cifti brain model axes - seed contains boolean arrays for each hemisphere
bm_l      = cifti2_axes.BrainModelAxis.from_mask(seed[0], name=f'CortexLeft')
bm_r      = cifti2_axes.BrainModelAxis.from_mask(seed[1], name=f'CortexRight')
bm        = bm_l + bm_r
# I know these structures above are not the problem

# loop through subcortical
coords = np.loadtxt(os.path.join(ptx_folder[2], 'coords_for_fdt_matrix2'), dtype='int') # coordinates from probtrackx containing coords for all subcortical voxels
for i,struct in enumerate(np.loadtxt(seed_path[2], dtype="str")):
    # load the seed
    subimg = nib.load(struct) # this is the subcortical seed mask used in probtrackx
    p = subimg.get_fdata().astype(bool)
    struct_name = re.search(r'CIFTI_STRUCTURE_(.*?).nii.gz', struct).group(1)

    # get the coords for the current structure and append to the brain model
    subcort = cifti2_axes.BrainModelAxis.from_mask(p, name=struct_name, affine=subimg.affine)
    if i == 0:
        ida = 0
    else:
        ida = subseed_voxnum_list[i-1] # this is the number of voxels in the mask

    subcort.voxel = coords[ida:ida+subseed_voxnum_list[i],:3]
    bm += subcort # append the current subcortical structure to the brain model

# save cifti
new_fname = os.path.join(out_folder, f'{prefix}BP.LR.dscalar.nii')
sc        = cifti2_axes.ScalarAxis(tracts)
hdr       = cifti2.Cifti2Header.from_axes((sc, bm))
img       = cifti2.Cifti2Image(BP.T, hdr)
nib.save(img, new_fname)

Version:

nibabel v4.0.2
Python v3.10.8

Relevant log outputs (up to 20 lines):

wb_command -file-information ~/Downloads/bp_distthr/BP.LR.allsubcort.dscalar.nii

While running:
/Applications/workbench/bin_macosx64/../macosx64_apps/wb_command.app/Contents/MacOS/wb_command -file-information BP.LR.allsubcort.dscalar.nii

ERROR: Cifti XML error: volume models may not reuse voxels, either internally or from other structures
2

I don’t do a ton with CIFTI, but I’m a bit surprised by

subcort.voxel = coords[ida:ida+subseed_voxnum_list[i],:3]

after constructing subcort = BrainModelAxis.from_mask(...). Shouldn’t that be enough?

3

Unfortunately not. Although skipping that line creates a wb_view/wb_command readable CIFTI file, when viewed the voxels are shuffled. I’m not entirely sure but I think this is because of voxel indexing differences between CIFTI and probtrackx.

image
image790×444 34.2 KB

4

I see. So you shouldn’t be constructing the models from masks, then, as that calculation of indices from nonzero values is pure waste. I think the issue you’re getting with overlapping masks is going to come from a bad calculation of indices, and you could use:

ida += subseed_voxnum_list[i-1]

Assuming that’s right, I think you should be able to do the following:

ref_img = nib.load('any_volume.nii.gz')
start_idx = 0
for struct, nvox in zip(np.loadtxt(seed_path[2], dtype="str"), subseed_voxnum_list):
    struct_name = re.search(r'CIFTI_STRUCTURE_(.*?).nii.gz', struct).group(1)
    bm += nib.cifti2.BrainModelAxis(
        name=struct_name,
        voxel=coords[start_idx : start_idx + nvox, :3],
        affine=ref_img.affine,
        volume_shape=ref_img.shape,
    )
    start_idx += nvox
1 Like
5

Brilliant, this works! Thank you very much.