Spacing transform only changes spacing in x- and y-axes, not in z-axis

[JaroKuiken]

Hi all,

I have an MRI dataset where each patient has a scan volume of 720 slices of 100x720 (where each slice is a DICOM image), so a shape of (720, 100, 720). They have an isotropic voxel spacing ranging from around 0.2 to 0.4mm. I want to retain the spacing in the x- and y-axes, but change the spacing in the z-axis to 1.6mm. For this i have been trying to use the Spacing transform. I have been able to change the spacing in the x- and y-axes, which was clearly visible in both the resulting shape change of the images and when plotting them. So i know the transform does work on my images, but the spacing in the z-axis refuses to change somehow.

I think it might have to do with how i load my images. I tried to use the LoadImage by itself, but it kept complaining about the affines of slices within the same volume not being the same. So instead i used the LoadImage but with a custom ImageReader, based on: https://simpleitk.readthedocs.io/en/master/link_DicomSeriesReader_docs.html. This allowed me to still load the images with the LoadImage transform, which (i thought) meant that i could use all the MONAI transforms as well. However it turns out that some transforms, such as the Orientation and ScaleIntensityRangePercentiles do work, while changing the spacing in the z-axis using the Spacing transform does not.

The code im currently running:

import os
import numpy as np
from monai.transforms import LoadImage, Compose, ScaleIntensityRangePercentiles, Spacing
from monai.data import ImageReader
from monai.utils.type_conversion import convert_to_numpy
from tqdm import tqdm
import matplotlib.pyplot as plt
import SimpleITK as sitk

class DICOMSeriesReader(ImageReader):
    def verify_suffix(self, filename: str) -> bool:
        # Check if the file is a DICOM by its extension
        return filename.lower().endswith((".dcm"))

    def read(self, data: str):
        # Create reader object
        reader = sitk.ImageSeriesReader()
        # Get the filenames of the DICOMs in the folder and pass to reader
        dicom_series = reader.GetGDCMSeriesFileNames(data)
        reader.SetFileNames(dicom_series)
        
        # Read the 3D volume
        image = reader.Execute()
        return image

    def get_data(self, img):
        # Convert SimpleITK image to numpy array
        img_data = sitk.GetArrayFromImage(img)
        
        # Get metadata to calculate affine
        origin = np.array(img.GetOrigin())
        spacing = np.array(img.GetSpacing())
        print(f"Spacing: {spacing}")
        direction = np.array(img.GetDirection()).reshape(3, 3)
        
        # Calculate affine
        affine = np.eye(4)
        affine[:3, :3] = direction * spacing
        affine[:3, 3] = origin
        
        # Metadata includes affine matrix and shape information
        metadata = {
            "affine": affine,
            "shape": img_data.shape,
            "spacing": spacing,
            "origin": origin,
            "direction": direction
        }
        
        return img_data, metadata

def loadAndTransformScan(scanPaths):
    transform = Compose([
        LoadImage(image_only=True, reader=DICOMSeriesReader, ensure_channel_first=True),
        ScaleIntensityRangePercentiles(lower=5, upper=95, b_min=0, b_max=1),
        Spacing(pixdim=(-1, -1, 1.6), mode='bilinear')
    ])
    patientScan = transform(scanPaths)
    print(f"\n----- Data, metadata and shapes after loading and transforming -----\n Whole shape: {np.asarray(patientScan).shape}, Type [0]: {type(patientScan[0])}, Type [1]: {type(patientScan[1])}, \n patientScan[0].meta: \n{patientScan[0].meta},\n patientScan[0].affine: \n{patientScan[0].affine},\n patientScan[1].meta: \n{patientScan[1].meta},\n patientScan[1].affine: \n{patientScan[1].affine}")

    patientScan = convert_to_numpy(patientScan, wrap_sequence=True)
    print(f"Shape after numpy conversion: {patientScan.shape}\n")  # Result: (720, 1, 720, 100, 1)

    patientScan = np.transpose(np.squeeze(patientScan), (0, 2, 1))
    print(f"Shape after transposing and squeezing: {patientScan.shape}") # Result: (720, 100, 720)
    return patientScan

def sortableNumberFromFilename(filename):
    res = filename.split('I')[-1].split('.')[0]
    res = str(res).rjust(4, '0')
    return res

def main():
    patients = []
    dirPath = "/mnt/toaster/CarotidMontreal/MICCAI_challenge_dataset/careIIChallenge"
    patientFolders = os.listdir(dirPath)

    # Loop through all patientFolders
    for patient in tqdm(patientFolders):
            # Create path per patient folder
            patientFolder = os.path.join(dirPath, patient)
            # List to add only scan paths
            scanPaths = []
            # List to store all patient slices in
            patientSlices = []

            # Sort all slices in each patient folder based on slice number, then add to a list
            for file in sorted(os.listdir(patientFolder), key=sortableNumberFromFilename):
                if not os.path.isdir(os.path.join(patientFolder, file)):
                    scanPaths.append(os.path.join(patientFolder, file))

            patientScan = loadAndTransformScan(scanPaths)

            # Only do it for one patient to check if Spacing works
            break

if __name__ == '__main__':
    main()

Most of this all seems to go correctly, except for the Spacing transform. I would expect the shape of the 3D volume to become about 4-8x smaller in the z-axis, as im increasing the spacing by that much. However, the shape only changes in the x- and y-axes if i indicate that in the Spacing transform, but it never changes in the z-axis. I have looked at putting the values in the 'pixdim' value in every possible configuration, or the patient specific spacing values instead of -1, but the same issue stays. I am starting to think it might be related to my custom ImageReader, but i am also not sure what the issue there would be. Does anyone know what the issue might be or how to solve it?

[KumoLiu]

Hi @JaroKuiken, Spacing should also work for z-axis. You can try the code below.

from monai.transforms import Spacing
from monai.data import MetaTensor
import torch

spacing = torch.eye(4)
spacing[2, 2] = 2.0
input = MetaTensor(torch.randn(1, 128, 128, 128), affine=spacing)

trans = Spacing([-1, -1, 1])
out = trans(input)
print("input affine:")
print(input.affine)
print("output affine:")
print(out.affine, out.pixdim, input.pixdim)

output:

input affine:
tensor([[1., 0., 0., 0.],
        [0., 1., 0., 0.],
        [0., 0., 2., 0.],
        [0., 0., 0., 1.]], dtype=torch.float64)
output affine:
tensor([[1., 0., 0., 0.],
        [0., 1., 0., 0.],
        [0., 0., 1., 0.],
        [0., 0., 0., 1.]], dtype=torch.float64) tensor([1., 1., 1.], dtype=torch.float64) tensor([1., 1., 2.], dtype=torch.float64)

Hope it helps, thanks.

[JaroKuiken]

Hi @KumoLiu,

Thank you for your suggestion! However my issue is not that the affine does not change. The issue is that the affine does change, in the location of the z-axis (as in your example) the value changes to 1.6 as intended, however nothing else changes. I would expect the shape of the volume to become smaller in the z-axis if i increase the spacing value in the affine, but the shape remains unchanged even if the value in the affine is changed to 1.6. Do you know why this could be?

[KumoLiu]

Hi @JaroKuiken, you can try the example I posted and print the out.shape. The shape also changed as expected.