I have a patient’s cardiac MRI scan (Dicom file) and its corresponding LA label in a VTK format. I wrote a code to convert the VTK file to binary mask images (like nrrd, nii.gz,etc.) I want the generated binary mask map MRIs. How to set the origin, spacing, and orientation in the code?
My code is as follows:
import numpy as np
import vtk
import cv2
import pydicom
import os
vtk_path =XX
dcms_path=XX
output_path = XXXX + ‘niigz/’ + patient + ‘_LA.nii.gz’
def get_matrix_patienttodicom(dicom_file):
# Load DICOM file
ds = pydicom.dcmread(dicom_file)
## Extract image orientation and position
iop = ds.ImageOrientationPatient
ipp = ds.ImagePositionPatient
print(iop)
image_spacing = ds.PixelSpacing
slice_thickness = ds.SliceThickness
origin_image_spacing = [image_spacing[0], image_spacing[1], slice_thickness]
# Calculate the normal vector of the image plane
normal = np.cross(iop[:3], iop[3:])
w=normal / np.linalg.norm(normal)
# Create a 4x4 transformation matrix
matrix = np.zeros((4, 4))
a=np.reshape(iop, (2, 3))
matrix[:3, :2] = a.transpose()
matrix[:3, 2] = np.array(w)
matrix[:3, 3] = np.array([0,0,0])
matrix[3, 3] = 1.0
#Return the matrix
dicom_to_patient = matrix
patient_to_dicom = np.linalg.inv(dicom_to_patient)
patient_to_dicom[0, :]/=origin_image_spacing[0]
patient_to_dicom[1, :]/=origin_image_spacing[1]
patient_to_dicom[2, :]/=origin_image_spacing[2]
return patient_to_dicom
patient_to_dicom = get_matrix_patienttodicom(dcms_path + “/LGE3D_Slice1.dcm”)
print(patient_to_dicom)
‘’’
[ 0.46748185 -0.46422961 0. 0. ]
[-0. -0. -0.65882351 -0. ]
[ 0.32028829 0.32253213 0. 0. ]
[ 0. 0. 0. 1. ]
‘’’
#vtk_to_patient=??
vtk_dicom = np.dot(patient_to_dicom,vtk_to_patient)
#transform vtk to dicom
reader = vtk.vtkPolyDataReader()
reader.SetFileName(vtk_path)
reader.Update()
#create a vtkMatrix4x4 object
vtk_matrix = vtk.vtkMatrix4x4()
for i in range(4):
for j in range(4):
vtk_matrix.SetElement(i, j, vtk_dicom[i, j])
#set the matrix to the transform
transform.SetMatrix(vtk_matrix)
transform_filter = vtk.vtkTransformPolyDataFilter()
transform_filter.SetInputData(reader.GetOutput())
transform_filter.SetTransform(transform)
transform_filter.Update()
transformed_poly = transform_filter.GetOutput()#
#obtain MRIs size
dcm_path = dcms_path + “/LGE3D_Slice1.dcm”
dcm_file = pydicom.dcmread(dcm_path)
image_position = dcm_file.ImagePositionPatient
h = dcm_file.Rows
w = dcm_file.Columns
z = len(os.listdir(dcms_path))
print("size: %d*%d*%d " %(h,w,z))
#set origin
original_origin_array = [image_position[0], image_position[1], image_position[2],1.0]
transformed_origin_array = np.dot(patient_to_dicom,original_origin_array)
print(transformed_origin_array)
polydata = transformed_poly
origin = np.zeros(3)
spacing = np.zeros(3)
spacing = [1,1,1]
dx=h
dy=z
dz=w
origin[0] = transformed_origin_array[0]
origin[1] = transformed_origin_array[1]
origin[2] = transformed_origin_array[2]#
sx, sy, sz = spacing
ox, oy, oz = origin # Origin
image = vtk.vtkImageData()
image.SetSpacing((sx, sy, sz))
image.SetDimensions((dx, dy, dz))
image.SetExtent(0, dx - 1, 0, dy - 1, 0, dz - 1)
image.SetOrigin((ox, oy, oz))
image.AllocateScalars(vtk.VTK_UNSIGNED_CHAR, 1)
inval = 1
outval = 0
for i in range(image.GetNumberOfPoints()):
image.GetPointData().GetScalars().SetTuple1(i, inval)
pol2stenc = vtk.vtkPolyDataToImageStencil()
pol2stenc.SetInputData(polydata)
pol2stenc.SetOutputOrigin((ox, oy, oz))
pol2stenc.SetOutputSpacing((sx, sy, sz))
pol2stenc.SetOutputWholeExtent(image.GetExtent())
pol2stenc.Update()
imgstenc = vtk.vtkImageStencil()
imgstenc.SetInputData(image)
imgstenc.SetStencilConnection(pol2stenc.GetOutputPort())
imgstenc.ReverseStencilOff()
imgstenc.SetBackgroundValue(outval)
imgstenc.Update()
writer = vtk.vtkNIFTIImageWriter()
writer.SetFileName(output_path)
writer.SetInputConnection(imgstenc.GetOutputPort())
writer.Write()