General questions about memory consumption and ANTs for large dataset

valeryozenne commented 3 years ago

Hi,

I have general questions about memory consumption and ANTs.

Context

I'm setting up a new pipeline for processing uCT cardiac datasets. One step of the pipeline require the downsampling of the tensor images. The simple solution (for me) was to do it with ANTs to use the log space computation. The pipeline is the following:

1) Compute Tensor from anatomical images (using home made code) and save each tensor components in nifti format.
2) Create the Tensor (ANTs compatible): ImageMath 3 $TENSOR ComponentTo3DTensor $COMPONENTS_
3) Resample the tensor : AntsApplyTransforms -e 2 -i $TENSOR -t $IDENTITY_MATRIX.mat -r $REF -o $NEW_TENSOR with REF a reference image at the desired resolution.

ps: I also try with the "RebaseTensorImage" command but it has not effect as a linear transformation is used.

My dataset is quite large "matrice size: 3016 x 3392 x 3993" (18um), so I started with lower matrices. While it is perfectly running using a lower dataset "matrice size: 1005 x 1130 x 1331" (18um*3), step 2 failed when my matrices become too large.

I have tried different workaround.

The first solution was to just crop the data. It allows me to validate few things (see screeshot 1-2) but I now need the full dataset.
I checked that my RAM is sufficient to concatenate my 6 components [xx, xy, xz, yy, yz, zz] (17 Gb each) in a 5D nifti file (98 Gb) using another package but the resulting nifti header is not compatible with the tensor format of ITK (see screenshot 3).
Finally, to my knowledge, ANTs is the only high level toolbox that allow manipulation/registration/averaging of tensor images.

Questions

1) Is there any magic solution to decrease the memory consumption with ANTs ? For instance 'AntsApplyTransforms' have a "--float" option , is there any hidden option for ImageMath ?

2) One solution would be to not use ANTs and use directly the ITK function behind ANTs. I would like to avoid this solution because it could be a solution where I recode badly what you did with no final improvement. But in case, does ANTs (in bash) could create a overcost in memory in comparison with a direct use of ITK in python for instance ?

3) For step 1, we are using a python package named DASK (https://dask.org/) and .zarr format(https://zarr.readthedocs.io/en/stable/), do you know if anyone has already tried to combine DASK and ANTs in Python or know if it would be possible, there are few posts suggesting compatibility with ITK (https://blog.dask.org/2019/08/09/image-itk) ?

Thanks in advance for your help

Valéry

System

I'm using a ubuntu computer with 377Gi

Data

Data is available.

Screenshot

Typical results after cropping the data Example_uCT_1 Example_uCT_2

PrintHeader of a tensorImage versus a 5D nifti Capture

gdevenyi commented 3 years ago

I can't answer all these questions, but I will offer some ideas

The MINC family of tools was designed from the start on low-memory systems in the 90's, and so it has voxelwise/slicewise streaming in all its tools so that memory usage is limited (at the cost of IO and processing time of course) so it's possible you could do some of the operations you want to do with MINC
The ITK folks are investigating the Zarr et al as a general format, but I don't think that work is done, https://discourse.itk.org/t/file-format-for-general-purpose-medical-image-computing/4067/2

cookpa commented 3 years ago

If the other package is MRtrix, you can convert to ANTs format:

https://github.com/ANTsX/ANTs/wiki/Importing-diffusion-tensor-data-from-other-software#mrtrix

cookpa commented 3 years ago

I took a look in the code, and it seems the tensor operations are hard-coded to use 32-bit float. Unfortunately there is no easy way to change this.

valeryozenne commented 3 years ago

Hi

Thank you both for your answers and investigation.

@cookpa, if I may , I have two questions regarding the tensor computation in ANTs,

could you confirm that the ITK framework for Tensor Reorientation behind ANTs is the following:

https://www.insight-journal.org/browse/publication/742

As even the "tensor creation" operation crash with ImageMath [...] ComponentTo3DTensor [...]. I made an attempt to replicate this function using ITK in python and I will investigate the use of zarr. The beginning was relatively easy even if I'm not really familiar with ITK except for the last point: I haven't been able to mimic the fonction WriteTensorImage in Python (simply without computing exp(D)). Could you describe/explain the fonction: NiftiDTICheck() and the following line in https://github.com/ANTsX/ANTs/blob/master/Utilities/ReadWriteData.h line 547 to 560 ?

Thanks a lot in advance

My basic script and corresponding (truncated) data are available at this link in case. https://filesender.renater.fr/?s=download&token=09f9b512-8de0-4371-aa16-d6a920883254

# overview
# 1) load all components
# 2) allocate a 3D tensor
# 3) copy the components in the tensor
# 4) save the tensor

# [...]

# 3) now we need to allocate a 3D tensor, 
# typedef itk::SymmetricSecondRankTensor<float, 3>           TensorType;
# typedef itk::Image<TensorType, ImageDimension>             TensorImageType;

TensorType=[itk.itkSymmetricSecondRankTensorPython.itkSymmetricSecondRankTensorD3, 3]

tensorImage = itk.Image[TensorType]

newTensor = tensorImage.New()

start = itk.Index[Dimension]()
start[0] = 0  # first index on X
start[1] = 0  # first index on Y
start[2] = 0  # first index on Z

new_size = itk.Size[Dimension]()
new_size[0] = size[0]  # size along X
new_size[1] = size[1]  # size along Y
new_size[2] = size[2]  # size along Z

new_tensor_region = itk.ImageRegion[Dimension]()
new_tensor_region.SetSize(new_size)
new_tensor_region.SetIndex(start)

newTensor.SetRegions(new_tensor_region)
newTensor.Allocate()

# finally, copy the data from the all components list_of_itk_image into the new tensor
# ['xx', 'xy', 'xz', 'yy', 'yz', 'zz']

for z in range(0,10): #size[2]):
    for y in range(0,10): #size[1]):
        for x in range(0,10): #size[0]):

            pixelIndex[0] = x
            pixelIndex[1] = y
            pixelIndex[2] = z

            pixel_tensor=newTensor.GetPixel(pixelIndex)

            # pix6[0] = xx->GetPixel( ind );
            # pix6[1] = xy->GetPixel( ind );
            # pix6[2] = xz->GetPixel( ind );
            # pix6[3] = yy->GetPixel( ind );
            # pix6[4] = yz->GetPixel( ind );
            # pix6[5] = zz->GetPixel( ind );
            #new_value=x+y+z

            pixel_tensor[0]=list_of_itk_images[0].GetPixel(pixelIndex)
            pixel_tensor[1]=list_of_itk_images[1].GetPixel(pixelIndex)
            pixel_tensor[2]=list_of_itk_images[2].GetPixel(pixelIndex)
            pixel_tensor[3]=list_of_itk_images[3].GetPixel(pixelIndex)
            pixel_tensor[4]=list_of_itk_images[4].GetPixel(pixelIndex)            
            pixel_tensor[5]=list_of_itk_images[5].GetPixel(pixelIndex) 

            newTensor.SetPixel(pixelIndex, pixel_tensor )

# not yet done
itk.WriteTensorImage(newTensor, output_3DTensor_filename)

cookpa commented 3 years ago

could you confirm that the ITK framework for Tensor Reorientation behind ANTs is the following:

https://www.insight-journal.org/browse/publication/742

This is Slicer's algorithm, ours is another implementation but it's similarly based on PPD by Alexander et al.

As even the "tensor creation" operation crash with ImageMath [...] ComponentTo3DTensor [...]. I made an attempt to replicate this function using ITK in python and I will investigate the use of zarr. The beginning was relatively easy even if I'm not really familiar with ITK except for the last point: I haven't been able to mimic the fonction WriteTensorImage in Python (simply without computing exp(D)). Could you describe/explain the fonction: NiftiDTICheck() and the following line in https://github.com/ANTsX/ANTs/blob/master/Utilities/ReadWriteData.h line 547 to 560 ?

I think that function was to change the ordering of the components. ITK eventually fixed this, so that function currently does nothing. NIFTI specifies symmetric matrix data such as tensors should be in lower-triangular order. VTK, ITK, FSL use upper-triangular internally. ITK now reads and writes lower-triangular tensors.

ANTsX / ANTs