Project-MONAI / MONAI

AI Toolkit for Healthcare Imaging
https://monai.io/
Apache License 2.0
5.85k stars 1.08k forks source link

LoadImage default reader improper affine handling on nrrd files with non-standard ijk_to_ras #7912

Closed AFanthomme closed 4 months ago

AFanthomme commented 4 months ago

Describe the bug When loading multiple nrrd files with different ijk_to_ras matrices and resampling / concatenating them, results have improper relative positions / orientations

To Reproduce We start with 3 nrrd files, which when opened in 3dSlicer are perfectly aligned but have different ijk_to_ras matrices:

They also have (a priori) different spatial sizes/extents, so I want to resample everyone into a single (3, H, W, D) nrrd file to have everything normalized later on.

To do so, I use the following:

case_dict = {'t2_ax': path_to_t2_sag.nrrd, 't2_sag': path_to_t2_sag.nrrd,  't1_3d': path_to_t1_3d.nrrd} 
all_keys = list(case_dict.keys())

loading_pipe = Compose([
        LoadImaged(keys=all_keys),
        EnsureChannelFirstd(keys=all_keys),
        ToTensord(keys=all_keys), 

        # This block is my (unsuccessful) attempt at solving the issue with transposes
        # # ResizeWithPadOrCropd(keys=['t2_sag'], spatial_size=(224,224,224), mode="constant", value=0),

        # Neither of those yields expectex results
        # # Orientationd(keys=all_keys, as_closest_canonical=True),
        # # Transposed(keys=['t2_sag'], indices=(0, 2, 3, 1)),

        # Pad with zeros since that's background value  
        ResizeWithPadOrCropd(keys=['t2_ax'], spatial_size=(224,224,224), mode="constant", value=0),
        ResampleToMatchd(keys=all_keys, key_dst="t2_ax", padding_mode="zeros", mode='bilinear'),
        ConcatItemsd(keys=all_keys, name='image'),
    ])

images_dict = loading_pipe(case_dict)

Expected behavior At this stage, my understanding is that images_dict["image"] should be a nrrd object with data a (3,224,224,224) tensor, and affine a (common between modalities) 4x4 tensor describing the relationship between position in the data tensor and position in real space.

Therefore, I expect the following:

from monai.data import ITKWriter
writer = ITKWriter()
image = images_dict["image"]
for channel_idx in range(3):
                            writer.set_data_array(image .data[channel_idx], channel_dim=None)
                            writer.set_metadata({"affine": image.affine})
                            writer.write(f"{k}_{channel_idx}.nrrd")

to yield 3 nrrd files that are identical to the ones used as inputs, except with common shapes.

However, loading back into 3dSlicer, I find that to be true for the 2 inputs that were in "standard" orientation", but the "t2_sag" one is completely messed up (aligned neither with the original, nor with the newly saved versions) (looking at the ijk_to_ras of these images, they all have the "standard" diagonal one)

Am I missing a step in the loading pipeline to force every data into a common affine representation? I expected this to be taken care of by ResampleToMatchd or ConcatItemsd, since all necessary information is available in the loaded nrrd files, but apparently this is not the case.

I tried using the transposition commented in the pipeline, but even if I managed to get it working this does not seem like a good solution since it hardcodes a behavior which should be adapted from information available in the nrrds themselves.

I have also tried to replace the last block in the loading pipeline with the same thing and "sag" replacing "ax":

# Pad with zeros since that's background value  
ResizeWithPadOrCropd(keys=['t2_ax'], spatial_size=(224,224,224), mode="constant", value=0),
ResampleToMatchd(keys=all_keys, key_dst="t2_ax", padding_mode="zeros", mode='bilinear'),
ConcatItemsd(keys=all_keys, name='image'),

In that case, all three resulting images were messed up, including the "t2_sag" one which now has to be viewed in coronal view rather than sagittal in 3DSlicer...

Let me know if you need any additional information

Environment

================================ Printing MONAI config...

MONAI version: 1.3.0 Numpy version: 1.26.4 Pytorch version: 2.2.0+cu121 MONAI flags: HAS_EXT = False, USE_COMPILED = False, USE_META_DICT = False MONAI rev id: 865972f7a791bf7b42efbcd87c8402bd865b329e MONAI file: /opt/conda/lib/python3.10/site-packages/monai/init.py

Optional dependencies: Pytorch Ignite version: 0.5.0.post2 ITK version: 5.4.0 Nibabel version: 5.2.1 scikit-image version: NOT INSTALLED or UNKNOWN VERSION. scipy version: 1.14.0 Pillow version: 10.2.0 Tensorboard version: 2.17.0 gdown version: NOT INSTALLED or UNKNOWN VERSION. TorchVision version: NOT INSTALLED or UNKNOWN VERSION. tqdm version: 4.65.0 lmdb version: NOT INSTALLED or UNKNOWN VERSION. psutil version: 5.9.0 pandas version: 2.2.2 einops version: 0.8.0 transformers version: NOT INSTALLED or UNKNOWN VERSION. mlflow version: NOT INSTALLED or UNKNOWN VERSION. pynrrd version: 1.0.0 clearml version: NOT INSTALLED or UNKNOWN VERSION.

For details about installing the optional dependencies, please visit: https://docs.monai.io/en/latest/installation.html#installing-the-recommended-dependencies

================================ Printing system config...

System: Linux Linux version: Ubuntu 22.04.4 LTS Platform: Linux-5.15.146.1-microsoft-standard-WSL2-x86_64-with-glibc2.35 Processor: x86_64 Machine: x86_64 Python version: 3.10.14 Process name: python Command: ['python', 'src/datasets/v2_monai.py'] Open files: [] Num physical CPUs: 16 Num logical CPUs: 32 Num usable CPUs: 32 CPU usage (%): [5.5, 4.9, 4.4, 4.9, 4.9, 4.4, 4.4, 4.9, 4.4, 4.9, 4.4, 4.9, 4.4, 4.9, 4.4, 4.9, 4.9, 4.9, 4.4, 4.9, 4.4, 4.9, 4.9, 4.4, 4.9, 4.9, 4.9, 4.4, 4.4, 4.9, 4.9, 91.3] CPU freq. (MHz): 3187 Load avg. in last 1, 5, 15 mins (%): [0.3, 0.2, 0.6] Disk usage (%): 19.3 Avg. sensor temp. (Celsius): UNKNOWN for given OS Total physical memory (GB): 31.2 Available memory (GB): 28.0 Used memory (GB): 2.8

================================ Printing GPU config...

Num GPUs: 1 Has CUDA: True CUDA version: 12.1 cuDNN enabled: True NVIDIA_TF32_OVERRIDE: None TORCH_ALLOW_TF32_CUBLAS_OVERRIDE: None cuDNN version: 8902 Current device: 0 Library compiled for CUDA architectures: ['sm_50', 'sm_60', 'sm_70', 'sm_75', 'sm_80', 'sm_86', 'sm_90'] GPU 0 Name: NVIDIA GeForce RTX 4090 GPU 0 Is integrated: False GPU 0 Is multi GPU board: False GPU 0 Multi processor count: 128 GPU 0 Total memory (GB): 24.0 GPU 0 CUDA capability (maj.min): 8.9

AFanthomme commented 4 months ago

After a bunch more digging, it turns out the problem was not from resample or any of the following, but rather from LoadImaged() since it already gave errrors even without any resampling or reorientation.

Switching to LoadImaged(reader='ITKReader') rather than default solves this issue, so I assume that when no reader is specified it must default to NrrdReader in that case, which for some reason handles affine incorrectly.

I updated the issue title and closed it, but will keep it here it helps someone else.