[SOLVED] Get 3D affine transformation matrix from mri DICOM files

Get 3D affine transformation matrix from mri DICOM files

I have both RTDose and MRI information in DICOM format. Unfortunately raw MRI and RTDose files share not only different dimensions but incompatible actual spacing as well. Looking Here I've been trying to understand how to build an affine transformation matrix in order to resample RTDose information to the same spacing of the MRI information.

Here you can find my final attempt:

    #read mri
    mr_file_paths = [os.path.join(path_MR, f) for f in os.listdir(path_MR) if f.endswith('.dcm')]        
    mr_data = [pydicom.dcmread(p) for p in mr_file_paths] 
    mr_dtype = mr_data[0].pixel_array.dtype
    mr = np.array([p.pixel_array for p in mr_data], dtype=mr_dtype)
    
    #read rtdose
    rtd_file_path  = [os.path.join(path_RTD, f) for f in os.listdir(path_RTD) if f.endswith('.dcm')][0]
    rtd_data = pydicom.dcmread(rtd_file_path)
    rtd = np.float64(rtd_data.pixel_array)
    
    #get patient's position, orientation and mri pixel spacing
    mr_ippf, mr_ippl = np.array(mr_data[0].ImagePositionPatient), np.array(mr_data[-1].ImagePositionPatient)
    mr_iop = np.array(mr_data[0].ImageOrientationPatient)
    mr_spacing = np.array(mr_data[0].PixelSpacing)
    
    Δr, Δc = mr_spacing[0], mr_spacing[1]
    F11, F21, F31 = mr_iop[3:6]
    F12, F22, F32 = mr_iop[0:3]
    
    k1, k2, k3 = (mr_ippl - mr_ippf) / (len(mr)-1)
    Dx, Dy, Dz = mr_ippf
    #build affine transformation matrix
    affine = np.array([
        [F11*Δr,    F12*Δc, k1, Dx  ],
        [F21*Δr,    F22*Δc, k2, Dy  ],
        [F31*Δr,    F32*Δc, k3, Dz  ],
        [0,         0,      0,  1   ],
    ], dtype=np.float64)
    #apply affine transformation
    rtd = ndimage.affine_transform(rtd, affine, output_shape=mr.shape)

Unfortunately when i plot some slices of both RTDose and MRI I can see that this transforamtion is not working. Does someone have any suggestion? Thanks in advance

Solution

You can use SimpleITK to resample one DICOM series onto another. The code would look something like this.

import SimpleITK as sick

reader1 = sitk.ImageSeriesReader()
reader1.SetFileNames(mr_file_paths)
mr_image = reader1.Execute()

reader2 = sitk.ImageSeriesReader()
reader2.SetFileNames(rtd_file_path)
rtd_image = reader2.Execute()

rtd_resampled_image = sitk.Resample(rtd_image, mr_image)

SimpleITK's Resample function will resample your RTD image into the space of the MR image.

You can read the docs of the function here: https://simpleitk.org/doxygen/latest/html/namespaceitk_1_1simple.html#aadbb243c10d1aedea8e955e8beda4df0