Slicing of large (compressed) 4D files

[ChristianHinge]

I am working with dynamic 4D PET data (x,y,z,t), and one challenge is that it is not possible to do load the full data onto memory (often +30gb files).
Luckily, most analyses in dynamic PET can be performed independently on subsets of the array by either slicing at a specific time point (t=t1) or spatial location (z=z1). I want to show you some cool results and hear you out on your general approach to this problem - maybe I missed another handy approach. I use .dataobj to access the stored array dynamically.

My example files:
img.nii (30gb) shape=(440,440,645,62)
img.nii.gz (2.5gb) (identical but compressed version)

The compressed version was created by nib.save(nib.load("img.nii"),"img.nii.gz")
Now accessing the data:

img = nib.load("img.nii") 
img.dataobj[..., 3].sum() #(slicing on time) FAST
img.dataobj[...,3,:].sum() #(slicing on axial dim) FAST
img.dataobj[...,3,:,:].sum() #(slicing on cor/sag) SLOW/does not work

I assume this has something to do with how the data is stored on disk (fortran/C-layout), but maybe I am wrong.
For fun, I also tried this with the .gz compressed file:

img = nib.load("img.nii.gz") 
img.dataobj[..., 3].sum() #(slicing on time) FAST
img.dataobj[...,3,:].sum() #(slicing on axial dim) SLOW/does not work
img.dataobj[...,3,:,:].sum() #(slicing on cor/sag) SLOW/does not work

I was very surprised that it is possible to slice on the time dimension. I assumed that the compression would make it impossible to do this kind of dynamic slicing and loading of data subsets from disk. Interestingly, it is not possible to slice in the axial dimension. My guess is that 4D images are compressed in way where the x-y-z 3D volumes are compressed independently, which makes it possible to quickly obtain these arrays for a specified t.

Some algorithms (for instance voxel-wise Patlak modelling) require that I can slice on one of the spatial dimensions to obtain an x-y-t 3D array. This is possible with the .nii file, but I would love to get it to work with the .nii.gz file due to the significantly lower disk footprint. My idea is to somehow transpose either the data array or the save-order on the disk. If I am able to save it as (x,y,t,z) or (t,z,y,x) , then my hope is that the compression algorithm will make it possible to slice on the last spatial dimension!

However, I was unable to find a good way to do this. Any ideas? :)

[effigies]

Do you have indexed-gzip installed? It will dramatically speed up operations on compressed datasets by creating an mapping from gzip block offsets to byte offsets in the underlying data. Nibabel will use it transparently, if found.

NIfTI is Fortran-ordered, so volumes are the slowest-changing index, and hence contiguous. There's nothing special about the compression; it's just gzip, and there's no chunking to improve slicing.

You cannot reorder axes in NIfTI. You can do that in other formats. You could also convert to a generic format like .zarr, which allows chunked compression to improve data access patterns. @balbasty has written a spec called NIfTI-Zarr (neuroscales/nifti-zarr#7) which allows lossless preservation of NIfTI metadata, so you can easily go back and forth.

I am surprised that coronal/sagittal slicing is slow, as I haven't really experienced that. That's worth looking into, but it may just be that loading chunks with such large strides is inefficient over 30GB. If you're having trouble loading 30GB into memory, it could be that swapping and garbage collection are where you're actually spending your time.

[ChristianHinge]

Thank you so much! indexed-gzip works fantastically. Reordering axes felt like a very hacky solution anyways, so I am glad the indexed-gzip fixed it for me.

Before:

img.dataobj[...,i,:].sum() #(slicing on axial dim) 60 seconds per i
img.dataobj[...,i,:].sum() #(slicing on axial dim) 20 seconds first i, <1 seconds for subsequent i's .

My whole department has always been doing dynamic analyses on dicom data. Not only does the niffti data reduce space by a factor 10, it also much cleaner and significantly faster!

The loading of chunks with large stride was also what I expected to be the problem. I am loading the file from a network mounted drive (smb). Maybe this strided access pattern only becomes slow in these cases.

For context, I am able to load the 30GB into memory (I have 64gb RAM), but the dynamic analyses often require me to have running means and residuals (for datafitting) in which case the memory footprint can easily become too large.

Thanks for the help!