Recommended reader for `.nrrd` files? `NRRDReader` gives wrong output for non-diagonal orientations

[surajpaib]

Describe the bug

The default reader for '.nrrd' file extensions is the NRRDReader as shown in the docs, (https://docs.monai.io/en/stable/transforms.html#loadimage) as well as used by LoadImage when trying to load the file type (all dependencies are installed)

However, NRRDReader does not work accurately for non-diagonal orientation matrices.

ITKReader seems to work reliably for each use case, whereas NRRDReader doesn't.

Is there a reason why NRRDReader is preferred over the ITKReader? Also, since ITKReader reads nifti files, why is it not the default reader selected for all such formats?

To Reproduce

Code using the default NRRDReader:

import monai
import matplotlib.pyplot as plt

nrrd_datalist = ["path_to_nrrd_file"]

transform = monai.transforms.Compose([
                                      monai.transforms.LoadImage(image_only=True, ensure_channel_first=True),
                                      monai.transforms.DataStats(additional_info=lambda x: x.affine, prefix="On load"),
                                      monai.transforms.Orientation(axcodes="LPS"), 
                                      monai.transforms.Transpose(indices=[0, 3, 2, 1]),
                                      monai.transforms.DataStats(additional_info=lambda x: x.affine, prefix="After transform")
                                      ])

out = transform(nrrd_datalist)


# Get the mid indexes for each dimension
img = out[0][0]
mid_dim0 = img.shape[0] // 2
mid_dim1 = img.shape[1] // 2
mid_dim2 = img.shape[2] // 2

# Plot mid-axial, sagittal and coronal views
fig, axs = plt.subplots(1, 3, figsize=(15, 5))
axs[0].imshow(img[mid_dim0, :, :], cmap='gray')
axs[0].set_title('Mid dim0')
axs[1].imshow(img[:, mid_dim1, :], cmap='gray')
axs[1].set_title('Mid dim1')
axs[2].imshow(img[:, :, mid_dim2], cmap='gray')
axs[2].set_title('Mid dim2')
plt.show()

On load statistics:
Type: <class 'monai.data.meta_tensor.MetaTensor'> torch.float32
Shape: torch.Size([1, 512, 512, 78])
Value range: (-1024.0, 3071.0)
Additional info: tensor([[ 0.0000, -0.4941, 0.0000, -262.2529],
[ 0.4941, 0.0000, 0.0000, 56.2529],
[ 0.0000, 0.0000, 1.5000, 59.5000],
[ 0.0000, 0.0000, 0.0000, 1.0000]], dtype=torch.float64)

After transform statistics:
Type: <class 'monai.data.meta_tensor.MetaTensor'> torch.float32
Shape: torch.Size([1, 78, 512, 512])
Value range: (-1024.0, 3071.0)
Additional info: tensor([[ -0.4941, 0.0000, 0.0000, -262.2529],
[ 0.0000, -0.4941, 0.0000, 308.7588],
[ 0.0000, 0.0000, 1.5000, 59.5000],
[ 0.0000, 0.0000, 0.0000, 1.0000]], dtype=torch.float64)

Code using ITKReader:

import monai
import matplotlib.pyplot as plt

nrrd_datalist = ["path_to_nrrd_file"]

transform = monai.transforms.Compose([
                                      monai.transforms.LoadImage(image_only=True, ensure_channel_first=True, reader="ITKReader"),
                                      monai.transforms.DataStats(additional_info=lambda x: x.affine, prefix="On load"),
                                      monai.transforms.Orientation(axcodes="LPS"), 
                                      monai.transforms.Transpose(indices=[0, 3, 2, 1]),
                                      monai.transforms.DataStats(additional_info=lambda x: x.affine, prefix="After transform")
                                      ])

out = transform(nrrd_datalist)


# Get the mid indexes for each dimension
img = out[0][0]
mid_dim0 = img.shape[0] // 2
mid_dim1 = img.shape[1] // 2
mid_dim2 = img.shape[2] // 2

# Plot mid-axial, sagittal and coronal views
fig, axs = plt.subplots(1, 3, figsize=(15, 5))
axs[0].imshow(img[mid_dim0, :, :], cmap='gray')
axs[0].set_title('Mid dim0')
axs[1].imshow(img[:, mid_dim1, :], cmap='gray')
axs[1].set_title('Mid dim1')
axs[2].imshow(img[:, :, mid_dim2], cmap='gray')
axs[2].set_title('Mid dim2')
plt.show()

On load statistics:
Type: <class 'monai.data.meta_tensor.MetaTensor'> torch.float32
Shape: torch.Size([1, 512, 512, 78])
Value range: (-1024.0, 3071.0)
Additional info: tensor([[ 0.0000, 0.4941, 0.0000, -262.2529],
[ -0.4941, 0.0000, 0.0000, 56.2529],
[ 0.0000, 0.0000, 1.5000, 59.5000],
[ 0.0000, 0.0000, 0.0000, 1.0000]], dtype=torch.float64)

After transform statistics:
Type: <class 'monai.data.meta_tensor.MetaTensor'> torch.float32
Shape: torch.Size([1, 78, 512, 512])
Value range: (-1024.0, 3071.0)
Additional info: tensor([[-0.4941, 0.0000, 0.0000, -9.7471],
[ 0.0000, -0.4941, 0.0000, 56.2529],
[ 0.0000, 0.0000, 1.5000, 59.5000],
[ 0.0000, 0.0000, 0.0000, 1.0000]], dtype=torch.float64)

Environment
monai_env_config.txt

[surajpaib]

To add, this is the output from the NRRDReader for a diagonal orientation matrix.

import monai
import matplotlib.pyplot as plt

nrrd_datalist = ["/mnt/data1/RadiomicsFoundationModel/LUNG1/NRRDs/LUNG1-001/CT.nrrd"]


transform = monai.transforms.Compose([
                                      monai.transforms.LoadImage(image_only=True, ensure_channel_first=True),
                                      monai.transforms.DataStats(additional_info=lambda x: x.affine, prefix="On load"),
                                      monai.transforms.Orientation(axcodes="LPS"),
                                      monai.transforms.Transpose(indices=[0, 3, 2, 1]),
                                      monai.transforms.DataStats(additional_info=lambda x: x.affine, prefix="After transform")
                                      ])

out = transform(nrrd_datalist)


# Get the mid indexes for each dimension
img = out[0][0]
mid_dim0 = img.shape[0] // 2
mid_dim1 = img.shape[1] // 2
mid_dim2 = img.shape[2] // 2

# Plot mid-axial, sagittal and coronal views
fig, axs = plt.subplots(1, 3, figsize=(15, 5))
axs[0].imshow(img[mid_dim0, :, :], cmap='gray')
axs[0].set_title('Mid dim0')
axs[1].imshow(img[:, mid_dim1, :], cmap='gray')
axs[1].set_title('Mid dim1')
axs[2].imshow(img[:, :, mid_dim2], cmap='gray')
axs[2].set_title('Mid dim2')
plt.show()

On load statistics:
Type: <class 'monai.data.meta_tensor.MetaTensor'> torch.float32
Shape: torch.Size([1, 512, 512, 134])
Value range: (-1024.0, 3034.0)
Additional info: tensor([[ -0.9766, 0.0000, 0.0000, 249.5117],
[ 0.0000, -0.9766, 0.0000, 460.5117],
[ 0.0000, 0.0000, 3.0000, -681.5000],
[ 0.0000, 0.0000, 0.0000, 1.0000]], dtype=torch.float64)

After transform statistics:
Type: <class 'monai.data.meta_tensor.MetaTensor'> torch.float32
Shape: torch.Size([1, 134, 512, 512])
Value range: (-1024.0, 3034.0)
Additional info: tensor([[ -0.9766, 0.0000, 0.0000, 249.5117],
[ 0.0000, -0.9766, 0.0000, 460.5117],
[ 0.0000, 0.0000, 3.0000, -681.5000],
[ 0.0000, 0.0000, 0.0000, 1.0000]], dtype=torch.float64)

As you can see, the orientations differ from the non-diagonal case while ITKReader is consistent

[KumoLiu]

Hi @surajpaib, would you mind sharing the data then I could reproduce the issue? Thanks.

[maunzzz]

Hi @surajpaib
I've noticed the same problem. The issue is that MONAI interprets the "space direction" variable in .nrrd files incorrectly in NRRDReader. Looking at the definition of the .nrrd format https://teem.sourceforge.net/nrrd/format.html#spacedirections, the vectors in "space direction" should be interpreted as the columns in the upper left 3x3 part of the affine matrix, however MONAI instead inserts them as rows.

The solution is easy, if you want to use NRRDReader you need to transpose the "space direction" matrix read by pynrrd. Ideally, it should of course be fixed in MONAI.

Changes needed:
monai/data/image_reader.py
L1339: direction = img.header["space directions"].T

[KumoLiu]

Free feel to reopen it the PR not solve the issue, thanks.