[Python] NIfTI and DICOM dimensions mismatch despite matching orientations in python using...

Description:

I am working on aligning a NIfTI image to a DICOM coordinate system pixel by pixel. I followed the blog Aligning a DICOM to a NIFTI, pixel-by-pixel, to implement the transformation. However, after reorienting the NIfTI data, the shape of the NIfTI image remains the same as the original, which is unexpected. Follow this blog I write this code for nifti -> Dicom aligning:


import nibabel as nib
import numpy as np
import pydicom
import os

def generate_coordinate_matrix(coord_system):
"""
Generate a transformation matrix from the given coordinate system
to the world coordinate system.
aram coord_system: String representing the coordinate system (e.g., 'LSP', 'RAS', etc.)
:return: 4x4 transformation matrix
"""
if len(coord_system) != 3:
raise ValueError("Coordinate system must have exactly 3 characters (e.g., 'RAS', 'LSP').")

# Mapping from coordinate labels to axis directions
axis_mapping = {
'R': [1, 0, 0], # Right (positive X)
'L': [-1, 0, 0], # Left (negative X)
'A': [0, 1, 0], # Anterior (positive Y)
'P': [0, -1, 0], # Posterior (negative Y)
'S': [0, 0, 1], # Superior (positive Z)
'I': [0, 0, -1] # Inferior (negative Z)
}

# Extract the directions for X, Y, Z
x_axis = axis_mapping[coord_system[0]]
y_axis = axis_mapping[coord_system[1]]
z_axis = axis_mapping[coord_system[2]]

# Ensure the coordinate system is valid
if not np.isclose(np.dot(x_axis, y_axis), 0) or not np.isclose(np.dot(x_axis, z_axis), 0) or not np.isclose(np.dot(y_axis, z_axis), 0):
raise ValueError(f"Invalid coordinate system: {coord_system}. Axes must be orthogonal.")

# Create the transformation matrix
transform_matrix = np.eye(4)
transform_matrix[:3, 0] = x_axis # X axis
transform_matrix[:3, 1] = y_axis # Y axis
transform_matrix[:3, 2] = z_axis # Z axis

return transform_matrix

def calculate_transform_matrix(source_coord, target_coord):
"""
Calculate the transformation matrix to convert coordinates
from source_coord to target_coord.
aram source_coord: Source coordinate system (e.g., 'LSP')
aram target_coord: Target coordinate system (e.g., 'RAS')
:return: 4x4 transformation matrix
"""
# Generate matrices for both systems
source_to_world = generate_coordinate_matrix(source_coord)
target_to_world = generate_coordinate_matrix(target_coord)

# Calculate the transform from source to target
transform_matrix = np.linalg.inv(target_to_world) @ source_to_world
return transform_matrix



Follow the transformation，

M_{DICOM <- NIfTI} = [M_{nifit axes <- Dicom axes} * M_{Dicom axes <- DICOM}]^-1 * M_{nifit axes <- NIfTI}


I computed the transformation matrix from nifti to dicom:

dcm_folder = ''
nii_file = ''

nii_img = nib.load(nii_file)
nii_data = nii_img.get_fdata()

#load dicom files
dicom_files = sorted(
[os.path.join(dcm_folder, f) for f in os.listdir(dcm_folder) if f.endswith('.dcm')],
key=lambda x: pydicom.dcmread(x).InstanceNumber
)

# Compute the affine matrix for the DICOM files
ds1 = pydicom.dcmread(dicom_files[0])
ds2 = pydicom.dcmread(dicom_files[1])
slice_diff = (np.array(ds2.ImagePositionPatient) - np.array(ds1.ImagePositionPatient)) / (
ds2.InstanceNumber - ds1.InstanceNumber)
c = slice_diff / np.sum(slice_diff**2) ** 0.5
matrix = np.zeros((4, 4))
matrix[3, 3] = 1
matrix[:3, 0] = b
matrix[:3, 1] = a
matrix[:3, 2] = c

dicom_axes = nib.aff2axcodes(matrix)
affine_matrix = matrix
print(f'dicom affine_matrix: \n {affine_matrix}')
print(f'dicom affine_matrix_axes: \n {dicom_axes}')

nii_img = nib.load(nii_file)
nii_affine = nii_img.affine
nii_axes = nib.aff2axcodes(nii_affine)
print(f'nii affine_matrix: \n {nii_affine}')
print(f'nii affine_matrix_axes: \n {nii_axes}')

# M_{Dicom -< nifti} = [M_{nifti_space <- Dicom_space} @ M_{Dicom_space<-Dicom} ]^-1 @ M_{nifti_space<-nifti}
M_axes = calculate_transform_matrix(nii_axes ,dicom_axes)
M_nifit2dicom = np.linalg.inv(M_axes @ affine_matrix) @ nii_affine

#transform nifti to Dicom
ornt = nib.orientations.io_orientation(M_nifit2dicom)
reoriented_nifti = nib.orientations.apply_orientation(nii_data, ornt)
print(f'reoriented_nifti_shape: {reoriented_nifti.shape}')


and the output shows:

The affine matrix of dicom

dicom affine_matrix:
[[ 0.02 1. -0.09 0. ]
[ 0. 0.09 1. 0. ]
[-1. 0.02 0. 0. ]
[ 0. 0. 0. 1. ]]


Axes and shape:

dicom affine_matrix_axes:
('I', 'R', 'A')
dicom_shape: (512, 496)


The affine matrix of nifti

nii affine_matrix:
[[ -0.58 -0.19 0.01 149.29]
[ -0.06 1.99 0. 15.83]
[ 0.01 -0. 0.58 180.59]
[ 0. 0. 0. 1. ]]


Axes and shape:

nii affine_matrix_axes:
('L', 'A', 'S')
nii_shape: (496, 114, 512)



The transformation matrix from Nifti to Dicom

array([[ -0.58, -0.22, 0.01, 149.38],
[ 0.06, -1.98, 0.05, -0.95],
[ 0. , 0.19, 0.58, 181.21],
[ 0. , 0. , 0. , 1. ]])



The transformation matrix from IRA to LAS

array([[ 0., 0., -1., 0.],
[-1., 0., 0., 0.],
[ 0., 1., 0., 0.],
[ 0., 0., 0., 1.]])



The reoriented nifti array shape:

reoriented_nifti_shape: (496, 114, 512)



After computing the transformation matrix M_nifit2dicom and applying the orientation using nibabel.orientations.apply_orientation, the shape of the NIfTI data remains unchanged: (496, 114, 512). I know the araay of Dicom files are flipped follow this document, but the axis of 114,512 haven't changed. This is unexpected as I thought the reorientation would adjust the shape to match the DICOM axes.

Questions:

• Why does the shape of the NIfTI data remain unchanged after applying the orientation?
• Did I make any mistakes in the computation of the transformation matrix or the application of apply_orientation?

Any help or guidance would be greatly appreciated!

Continue reading...