How to resample anisotropic voxel size (CT) into isotropic voxel size

Question

I am working on a code to standardize voxel sizes of 3D Dicom series from multiple machines (as part of a pipeline, but for now we'll focus on voxel sizes). It needs to be able to handle input of both isotropic and anisotropic voxels, and also retaining original size of objects because I will be performing measurements.

The code is based on Resample Scalar Volume from Slicer using mostly simple ITK, which works well for all of my dicom inputs. The original code is in C++ but I'm modifying it into Python. However, my code resulted in some bugs. I am testing all of my outputs in Slicer, as I will be using that software for analysis. Please be gentle as I am new to coding.

I tried two different approaches that both resulted in dead end, I am happy to receive advice on any one of them. In both approaches, actual number of voxel size and slice thickness of the output are already correct, but visual appearance is false.

Approach 1: Output as DICOM series. The output only produces .dcm files, but for some reason Slicer detected a transform even though there is no transform file. The actual CTvolume is correct, but the transform resulted in objects measured smaller than it should in xy plane. I want to get rid of the transform effect as I worry it will carry on to the next part in the pipeline. This happens to both isotropic and anisotropic inputs. I already tried to get rid of SetTransform entirely from the code, but the transform effect still persisted. Tried adjusting image direction and position but not sure if it works.

With transform activated (smaller objects in xy plane): With transform activated (smaller objects in xy plane)

Transform deactivated (correct appearance): Transform deactivated (correct appearance)

# Function to resize voxel spacing using SimpleITK
def resize_voxel(itk_image, original_spacing, original_direction, target_spacing):
    original_size = itk_image.GetSize()
    target_size = [
        int(original_size[i] * (original_spacing[i] / target_spacing[i]) + 0.5)
        for i in range(3)
    ]

    resampler = sitk.ResampleImageFilter()
    resampler.SetOutputSpacing(target_spacing)
    resampler.SetSize(target_size)
    resampler.SetOutputDirection(itk_image.GetDirection())
    resampler.SetOutputOrigin(itk_image.GetOrigin())
    
    # Use an explicit identity transform to prevent unintended scaling or rotation
    transform = sitk.Euler3DTransform()
    transform.SetIdentity()
    resampler.SetTransform(transform)
    resampler.SetInterpolator(sitk.sitkBSpline)

    # Adjust to ensure absolute sizes are preserved with high precision
    resampler.SetDefaultPixelValue(itk_image.GetPixelIDValue())
    resampler.SetOutputPixelType(itk_image.GetPixelID())
    resampler.SetNumberOfThreads(4)

    resized_image = resampler.Execute(itk_image)
    return resized_image

# Function to save the 3D image as a DICOM series
def save_dicom_series(image_3d, dicom_names, output_folder, target_spacing, original_image_position_patient):
    if not os.path.exists(output_folder):
        os.makedirs(output_folder)

    # Convert SimpleITK image back to NumPy array
    image_3d_array = sitk.GetArrayFromImage(image_3d)

    for i in range(len(image_3d_array)):
        slice_array = image_3d_array[i]
        output_dicom_path = os.path.join(output_folder, f"{i:03d}.dcm")

        # Load the original slice's DICOM for metadata and anonymize it
        original_dicom = anonymize_dicom(dicom_names[i])

        # Update the pixel data
        original_dicom.PixelData = slice_array.tobytes()
        original_dicom.Rows, original_dicom.Columns = slice_array.shape

        # Update the spacing information in the DICOM metadata
        original_dicom.PixelSpacing = [target_spacing[0], target_spacing[1]]
        original_dicom.SliceThickness = target_spacing[2]
        original_dicom.SpacingBetweenSlices = target_spacing[2]

        # Update Image Position Patient to reflect new slice location
        if original_image_position_patient is not None:
            new_position = original_image_position_patient.copy()
            new_position[2] = round(original_image_position_patient[2] + float(i) * target_spacing[2], 6)
            original_dicom.ImagePositionPatient = new_position

        # Update Slice Location and Instance Number
        original_dicom.InstanceNumber = i + 1
        original_dicom.SliceLocation = i * float(original_dicom.SliceThickness)

        # Debug output for updated metadata
        print(f"Saving slice {i + 1} with SliceThickness: {original_dicom.SliceThickness}, PixelSpacing: {original_dicom.PixelSpacing}, SpacingBetweenSlices: {original_dicom.SpacingBetweenSlices}, ImagePositionPatient: {original_dicom.ImagePositionPatient}")

        original_dicom.save_as(output_dicom_path)

# Main function to process the entire DICOM series
def process_dicom_series(patient_folder, output_folder):
    # Load the DICOM series
    image_3d, dicom_names, spacing, image_position_patient, image_direction = load_dicom_series(patient_folder)

    # Convert numpy array to SimpleITK image for processing
    itk_image = sitk.GetImageFromArray(image_3d)
    itk_image.SetSpacing(spacing)
    itk_image.SetDirection(image_direction)

    # Resize the voxel image to standardized voxel size of 180x180x180 micrometers
    target_spacing = (0.18, 0.18, 0.18)
    resized_voxel_image = resize_voxel(itk_image, spacing, image_direction, target_spacing)
    
    # Save the processed 3D image as a DICOM series
    save_dicom_series(resized_voxel_image, dicom_names, output_folder, target_spacing, image_position_patient)

Approach 2: I suspect issue 1 caused by saving to dicom process, so in the second approach I tried saving to nifti. It worked well if my input is isotropic, no weird transformation, looks great. However if my input is anisotropic, image in xy plane seems fine but is stretched out in z-axis, upside down, with weird alternating lines. Tried adjusting image direction and position but not sure if it does anything because image origin is 0 (it was not 0 in the original image).

I’m not sure what’s happening with the z-axis, looks almost like it’s spliced together with another image/another part of the image except for the top part of the image (actually inferior part of the patient, because the image is upside down). Absolute size in z-axis and height dimension is larger by almost 3 times than it should be.

Image of z axis: Image of z axis

# Function to resize voxel spacing using SimpleITK
def resize_voxel(itk_image, original_spacing, original_direction, target_spacing):
    original_size = itk_image.GetSize()
    target_size = [
        int(original_size[i] * (original_spacing[i] / target_spacing[i]) + 0.5)
        for i in range(3)
    ]

    resampler = sitk.ResampleImageFilter()
    resampler.SetOutputSpacing(target_spacing)
    resampler.SetSize(target_size)
    resampler.SetOutputDirection(original_direction)
    resampler.SetOutputOrigin(itk_image.GetOrigin())
    
    # Use an explicit identity transform to prevent unintended scaling or rotation
    transform = sitk.Euler3DTransform()
    transform.SetIdentity()
    resampler.SetTransform(transform)
    resampler.SetInterpolator(sitk.sitkBSpline)

    # Adjust to ensure absolute sizes are preserved with high precision
    resampler.SetDefaultPixelValue(itk_image.GetPixelIDValue())
    resampler.SetOutputPixelType(itk_image.GetPixelID())
    resampler.SetNumberOfThreads(4)

    resized_image = resampler.Execute(itk_image)
    return resized_image

# Function to save the 3D image as a NIfTI file using SimpleITK
def save_nifti(itk_image, output_path):
    sitk.WriteImage(itk_image, output_path)
    print(f"NIfTI volume successfully saved to {output_path}")

# Main function to process the entire DICOM series
def process_dicom_series(patient_folder, output_folder):
    # Load the DICOM series
    image_3d, dicom_names, spacing, image_position_patient, image_direction = load_dicom_series(patient_folder)

    # Convert numpy array to SimpleITK image for processing
    itk_image = sitk.GetImageFromArray(image_3d)
    itk_image.SetSpacing(spacing)
    itk_image.SetDirection(image_direction)

    # Resize the voxel image to standardized voxel size of 250x250x250 micrometers
    target_spacing = (0.25, 0.25, 0.25)
    resized_voxel_image = resize_voxel(itk_image, spacing, image_direction, target_spacing)
    
    # Save the processed 3D image as a DICOM series using SimpleITK
    save_nifti(resized_voxel_image, os.path.join(output_folder, 'output_volume.nii'))