Hi,
I have segmented the liver from a 3D-CT volume and created a volumetric
mesh. I have successfully generated DRRs by using RTK tool kit. The code is
attached below. Now what I want is that I need to project each and every 3D
coordinates of the volumetric mesh on this DRR image plane. This means that
finding the corresponding 2D projection coordinate (i.e. x, y positions)
of each vertex coordinate on the DRR image projection plane.
How can I obtain the (x, y) coordinate values in the DRR image plane for a
given 3D coordinate position? For example, I have (x,y,z) positions for
each vertex and I want to convert into (x,y) i.e. 2D projection coordinates
on the DRR image plane? Can you please elaborate with an example ?
I have python code to project the 3D-mesh coordinates on DRR images based
on ITK SiddonJacobsRayTracing algorithm. However, I want to implement this
code similar to the RTK based DRR projections. I have attached this ITK
based DRR generation code and the point projection code which is compatible
to this ITK-based DRR code for your further reference. Moreover, I have
attached sample DRR image with projected points overlaid on it that is
based on ITK-based code.
*Additional note:* The reason because I have chosen RTK since it is more
convenient than ITK based ray-tracing algorithm when generating DRRs so
that I’m able to generate DRRs to match the exact field of view of real kV
planar x-ray images.
Thanks,
Surang
import os
import itk
from itk import RTK as rtk
import numpy as np
import pydicom as dicom
def generateDRRsRTK(input_file=None, output_file=None, sid=1000., sdd=1536., gantryAngle=0.,
projOffsetX=0., projOffsetY=0., outOfPlaneAngle=0., inPlaneAngle=0.,
sourceOffsetX=0., sourceOffsetY=0., dx=512, dy=512):
CT = itk.imread(input_file, pixel_type=itk.F)
# The CT volume is not correct orientation compatible to the RTK convention and point (0,0,0) mm is not in the CT image
# and this is the default centre of rotation in RTK. Therefore change the origin and the direction to use
# RTK convention to get the correct DRR as expected.
# the input of the Joseph-Filter needs to be oriented in the y-direction. In RTK, the rotation axis is y.
# changed the direction and the iamge origin of the image volume to have the volume in the xz-layer in the y-direction.
# Three rotation angles are used to define the orientation of the detector.
# The ZXY convention of Euler angles is used for detector orientation where GantryAngle is the rotation around y,
# OutOfPlaneAngle the rotation around x and InPlaneAngle the rotation around z.
# change the direction and origin to align with the RTK convention
CTDirection=np.zeros([3,3])
CTDirection[0,0]=1.
CTDirection[1,2]=1.
CTDirection[2,1]=1.
CT.SetDirection(itk.matrix_from_array(CTDirection))
CT.SetOrigin([CT.GetOrigin()[0],CT.GetOrigin()[2],-CT.GetOrigin()[1]])
CTarray = itk.array_view_from_image(CT)
# add 1000 to CT numbers to put air at 0
CTarray += 1000
# Defines the image type
Dimension_CT = 3
PixelType = itk.F
ImageType = itk.Image[PixelType, Dimension_CT]
# Create a stack of empty projection images
ConstantImageSourceType = rtk.ConstantImageSource[ImageType]
constantImageSource = ConstantImageSourceType.New()
# Set origin and spacing based on the Elekta configuration
constantImageSource.SetOrigin([-204.4,-204.4,0])
constantImageSource.SetSpacing([0.8,0.8,0.8])
constantImageSource.SetSize([dx, dy, 1])
constantImageSource.SetConstant(0.)
# Defines the RTK geometry object
geometry = rtk.ThreeDCircularProjectionGeometry.New()
geometry.AddProjection(sid, sdd, gantryAngle, projOffsetX, projOffsetY, outOfPlaneAngle, inPlaneAngle, sourceOffsetX, sourceOffsetY)
REIType = rtk.JosephForwardProjectionImageFilter[ImageType, ImageType]
rei = REIType.New()
rei.SetGeometry(geometry)
rei.SetInput(0, constantImageSource.GetOutput())
rei.SetInput(1, CT)
rei.Update()
Dimension = 3
OutputPixelType = itk.UC
OutputImageType = itk.Image[OutputPixelType, Dimension]
RescaleType = itk.RescaleIntensityImageFilter[ImageType, OutputImageType]
rescaler = RescaleType.New()
rescaler.SetOutputMinimum(0)
rescaler.SetOutputMaximum(255)
rescaler.SetInput(rei.GetOutput())
rescaler.Update()
# Out of some reason, the computed projection is upsided-down.
# Here we use a FilpImageFilter to flip the images in y direction.
FlipFilterType = itk.FlipImageFilter[OutputImageType]
flipFilter = FlipFilterType.New()
FlipAxesArray = itk.FixedArray[itk.D, 3]()
FlipAxesArray[0] = 0
FlipAxesArray[1] = 1
FlipAxesArray[2] = 0
flipFilter.SetFlipAxes(FlipAxesArray)
flipFilter.SetInput(rescaler.GetOutput())
flipFilter.Update()
WriteType = itk.ImageFileWriter[OutputImageType]
writer = WriteType.New()
writer.SetFileName(output_file)
writer.SetInput(flipFilter.GetOutput())
writer.Update()
def getAndSaveDRR(dataFolder, rawReaKVFolder, outputFolder):
all_img_vol_names = [x for x in os.listdir(dataFolder) if x.endswith(".nii.gz")]
for i, relative_ct_path in enumerate(all_img_vol_names):
ct_vol_path = os.path.join(dataFolder, relative_ct_path)
drr_img_name = relative_ct_path.rsplit('.', 2)[0]
all_real_kv_names = [x for x in os.listdir(rawReaKVFolder) if x.endswith(".dcm")]
for j, relative_kv_path in enumerate(all_real_kv_names):
kVImgPath = os.path.join(rawReaKVFolder, relative_kv_path)
ds = dicom.dcmread(kVImgPath)
angle = float(ds.get("PositionerPrimaryAngle"))
# read FOV origin
fovOriginX = float(ds.get("FieldOfViewOrigin")[0])
fovOriginY = float(ds.get("FieldOfViewOrigin")[1])
drr_img_rel_path = drr_img_name + '_' + str(angle)
output_drr_path = os.path.join(outputFolder, drr_img_rel_path + ".png")
generateDRRsRTK(input_file=ct_vol_path, output_file=output_drr_path, gantryAngle=angle, projOffsetX=fovOriginX, projOffsetY=fovOriginY)
rawReaKVFolder = 'F:/RAW_REAL_KV_DICOM_132215'
ctDataFolder = 'F:/all_deformed_CTs'
outputDRRFolder = 'F:/all_drrs'
getAndSaveDRR(ctDataFolder, rawReaKVFolder, outputDRRFolder)import os
import glob
import numpy as np
import itk
import math
from PIL import Image, ImageDraw
import meshio
import matplotlib.pyplot as plt
# Angle in radians betweeen projection central axis and reference axis. [m_ProjectionAngle]
# Source to isocenter (i.e., 3D image center) distance in mm. [m_FocalPointToIsocenterDistance]
# DRR Pixel spacing in DETECTOR plane in mm [im_sx, im_sy]
# We have calculated this value by assuming that the pixel space in the detector plane is 1.0 mm.
# refer https://github.com/InsightSoftwareConsortium/ITKTwoProjectionRegistration/blob/master/test/GetDRRSiddonJacobsRayTracing.cxx#L112-L113
# Size of the DRR image in number of pixels. [dx, dy]
# This function is only use to calculate DRR projection coordinate for only one 3D vertex cooridnate
def getDRRProjectionCoords(ct_img_path=None, points=None, m_ProjectionAngle=0., m_FocalPointToIsocenterDistance=1000., m_Threshold=0.,
source_to_detector_distance=1536., im_sx=0.65104, im_sy=0.65104, dx=1024, dy=1024):
# Use the image center as the central axis position.
# refer https://github.com/InsightSoftwareConsortium/ITKTwoProjectionRegistration/blob/master/test/GetDRRSiddonJacobsRayTracing.cxx#L506-L510
o2Dx = (dx - 1.) / 2.
o2Dy = (dy - 1.) / 2.
# constant for converting degrees into radians
dtr = (math.atan(1.0) * 4.0) / 180.0
Dimension = 3
PointType = itk.Point[itk.D, Dimension]
VectorType = itk.Vector[itk.D, Dimension]
MatrixType = itk.Matrix[itk.D, Dimension, Dimension]
TransformType = itk.Euler3DTransform[itk.D]
# Compute the origin (in mm) of the 2D Image
m_Origin = PointType()
m_Origin[0] = -im_sx * o2Dx
m_Origin[1] = -im_sy * o2Dy
m_Origin[2] = -m_FocalPointToIsocenterDistance
m_Spacing = VectorType()
m_Spacing[0] = im_sx # pixel spacing along X of the 2D DRR image [mm]
m_Spacing[1] = im_sy # pixel spacing along Y of the 2D DRR image [mm]
m_Spacing[2] = 1.0 # slice thickness of the 2D DRR image [mm]
# inverse trasnform to rotate the vertex into DRR detector coordinate system
# user can set an arbitrary transform for the volume. If not explicitly set, it should be identity.
# Here I'm not setting the isocenter location and assume it is in (0, 0, 0) by default,
# but in the DRR algorithm that we used to generate DRR images, Set the center of the 3D image volume (CT) as the isocenter
# refer https://github.com/InsightSoftwareConsortium/ITKTwoProjectionRegistration/blob/master/test/GetDRRSiddonJacobsRayTracing.cxx#L452-L454
m_Transform = TransformType.New()
m_Transform.SetComputeZYX(True)
m_Transform.SetIdentity()
m_InverseTransform = TransformType.New()
m_InverseTransform.SetComputeZYX(True)
m_ComposedTransform = TransformType.New()
m_ComposedTransform.SetComputeZYX(True)
m_GantryRotTransform = TransformType.New()
m_GantryRotTransform.SetComputeZYX(True)
m_GantryRotTransform.SetIdentity()
m_CamShiftTransform = TransformType.New()
m_CamShiftTransform.SetComputeZYX(True)
m_CamShiftTransform.SetIdentity()
m_CamRotTransform = TransformType.New()
m_CamRotTransform.SetComputeZYX(True)
m_CamRotTransform.SetIdentity()
m_CamRotTransform.SetRotation(dtr * (-90.0), 0.0, 0.0)
# isocenter set to (0, 0, 0)
# but in the DRR algorithm that we used to generate DRR images, Set the center of the 3D image volume (CT) as the isocenter
inputImage = itk.imread(ct_img_path)
imOrigin = inputImage.GetOrigin()
imRes = inputImage.GetSpacing()
imRegion = inputImage.GetBufferedRegion()
imSize = imRegion.GetSize()
isocenter = [(imRes[i]*imSize[i]/2.0) for i in range(3)]
m_ComposedTransform.SetIdentity()
m_Transform.SetCenter(isocenter)
m_Transform.SetTranslation([-1.*imOrigin[i] for i in range(3)])
m_ComposedTransform.Compose(m_Transform, False)
# An Euler 3D transform is used to rotate the volume to simulate the roation of the linac gantry.
# The rotation is about z-axis. After the transform, a AP projection geometry (projecting
# towards positive y direction) is established.
m_GantryRotTransform.SetRotation(0.0, 0.0, -1. * dtr * (m_ProjectionAngle))
m_GantryRotTransform.SetCenter(isocenter)
m_ComposedTransform.Compose(m_GantryRotTransform, False)
# An Euler 3D transfrom is used to shift the camera to the detection plane
camtranslation = VectorType()
camtranslation[0] = -isocenter[0]
camtranslation[1] = m_FocalPointToIsocenterDistance - isocenter[1]
camtranslation[2] = -isocenter[2]
#print(camtranslation)
m_CamShiftTransform.SetTranslation(camtranslation)
m_ComposedTransform.Compose(m_CamShiftTransform, False)
# An Euler 3D transform is used to establish the standard negative z-axis projection geometry. (By
# default, the camera is situated at the origin, points down the negative z-axis, and has an up-
# vector of (0, 1, 0))
m_ComposedTransform.Compose(m_CamRotTransform, False)
# The overall inverse transform is computed. The inverse transform will be used by the interpolation
# procedure.
m_ComposedTransform.GetInverse(m_InverseTransform)
MatrixType = itk.Matrix[itk.D, Dimension, Dimension]
m_Direction = MatrixType()
m_Direction.SetIdentity()
m_IndexToPhysicalPoint = MatrixType()
m_IndexToPhysicalPoint.SetIdentity()
m_PhysicalPointToIndex = MatrixType()
m_PhysicalPointToIndex.SetIdentity()
scale = MatrixType()
# Compute helper matrices used to transform Index coordinates to
# PhysicalPoint coordinates and back. This method is virtual and will be
# overloaded in derived classes in order to provide backward compatibility
# behaviour in classes that did not used to take image orientation into
# account.
scale = itk.GetMatrixFromArray(np.diag(m_Spacing))
m_IndexToPhysicalPoint = m_Direction * scale
m_PhysicalPointToIndex = m_IndexToPhysicalPoint.GetInverse()
# this entire bottom part should be iterated over all the points
point_list = points.tolist()
# hold the pixel cooridnates as separate numpy arrays for x and y coordinates
u = np.array([])
v = np.array([])
for point in point_list:
# coordinate to the coordinate system defined by the detector coordinate system
drrPixelWorld = m_ComposedTransform.TransformPoint(point)
#print(drrPixelWorld)
# 2D coordinates on the detector plane (in physical units, not in pixels)
# calculate the physical cooridante in DRR image plane
# i.e. multiply the x and y coordinates by the ratio source-to-detector distance over z
drrPixelWorld[0] = -drrPixelWorld[0] * m_FocalPointToIsocenterDistance/drrPixelWorld[2]
drrPixelWorld[1] = -drrPixelWorld[1] * m_FocalPointToIsocenterDistance/drrPixelWorld[2]
index = VectorType()
for i in range(Dimension):
sum = itk.NumericTraits[itk.D].ZeroValue()
for j in range(Dimension):
sum += m_PhysicalPointToIndex(i, j) * (drrPixelWorld[j] - m_Origin[j])
index[i] = sum
u = np.append(u, index[0])
v = np.append(v, index[1])
return u, v
def getImageDict(img_paths):
img_dict = {}
for i, path in enumerate(img_paths):
key_with_angle = os.path.basename(path).rsplit('.', 1)[0]
# key hold the volumetric mesh name (i.e. ellipsoid) by removing the gantry angle (i.e. ellipsoid_72)
key = os.path.basename(key_with_angle).rsplit('_', 1)[0]
if key in img_dict:
img_dict[key].append(path)
else:
img_dict[key] = [path]
return img_dict
def alignedWithDRR(ct_img_path, meshPaths, img_dict, outputImgFolder):
for j, path in enumerate(meshPaths):
mesh = meshio.read(path)
img_key = os.path.basename(path).rsplit('.', 1)[0]
drr_img_path_list = img_dict[img_key]
for drr_img_path in drr_img_path_list:
# extract the drr image name with extention from the absolute path (e.g. volume-10_0.png)
img_name_str = os.path.basename(drr_img_path).rsplit('.', 1)[0]
# projection angle in degrees
gantry_angle = float(os.path.basename(img_name_str).rsplit('_', 1)[1])
print('gantry_angle:', str(gantry_angle))
img_name = img_name_str + '_aligned'
print('image_name:', img_name)
# need to pass projection angle to get the 2D coordinates for each case
u, v = getDRRProjectionCoords(ct_img_path=ct_img_path, points=mesh.points, m_ProjectionAngle=gantry_angle)
v = 1023-v # Mimicks flipFilter
u = u.tolist()
v = v.tolist()
output_img_path = os.path.join(outputImgFolder, img_name + ".png")
i = Image.open(drr_img_path).convert("RGBA")
draw = ImageDraw.Draw(i)
w, h = i.size
# plot U, V points
for x, y in zip(u, v):
draw.point((x, y), fill="red")
i.save(output_img_path)
ct_img_path = 'F:/all_CTs_HU_corrected/liver_ct.nii.gz'
mesh_paths = glob.glob('F:/meshes/*.vtk')
drr_img_paths = glob.glob('F:/ITK_DRRs/all_drrs/*.png')
output_img_folder = 'F:/results_proejctions'
img_dict = getImageDict(drr_img_paths)
alignedWithDRR(ct_img_path, mesh_paths, img_dict, output_img_folder)import os
import math
import numpy as np
import itk
# <-res float float> DRR Pixel spacing in isocenter plane in mm {im_sx, im_sy}
# <-size int int> Size of DRR in number of pixels [default: 512x512] {dx, dy}
# <-scd float> Source to isocenter (i.e., 3D image center) distance in mm [default: 1000mm]
# <-t float float float> CT volume translation in x, y, and z direction in mm {tx, ty, tz}
# <-rx float> CT volume rotation about x axis in degrees
# <-ry float> CT volume rotation about y axis in degrees
# <-rz float> CT volume rotation about z axis in degrees
# <-2dcx float float> Central axis position of DRR in continuous indices {o2Dx, o2Dy}
# <-iso float float float> Continous voxel indices of CT isocenter (center of rotation and projection center) {cx, cy, cz}
# <-rp float> Projection angle in degrees {rprojection}
# <-threshold float> CT intensity threshold, below which are ignored [default: 0]
def siddonJacobRayTracingDRR(input_name=None, output_name=None, rprojection=0., rx=0., ry=0., rz=0., tx=0., ty=0., tz=0.,
cx=0., cy=0., cz=0., threshold=0., scd=1000, im_sx=0.65104, im_sy=0.65104,
dx=512, dy=512, o2Dx=0.0, o2Dy=0.0):
Dimension = 3
InputPixelType = itk.F
OutputPixelType = itk.UC
InputImageType = itk.Image[InputPixelType, Dimension]
OutputImageType = itk.Image[OutputPixelType, Dimension]
ReaderType = itk.ImageFileReader[InputImageType]
reader = ReaderType.New()
reader.SetFileName(input_name)
reader.Update()
inputImage = reader.GetOutput()
# To simply Siddon-Jacob's fast ray-tracing algorithm, we force the origin of the CT image
# to be (0,0,0). Because we align the CT isocenter with the central axis, the projection
# geometry is fully defined. The origin of the CT image becomes irrelavent.
ctOrigin = [0.0, 0.0, 0.0]
inputImage.SetOrigin(ctOrigin)
FilterType = itk.ResampleImageFilter[InputImageType, InputImageType]
filter = FilterType.New()
filter.SetInput(inputImage)
filter.SetDefaultPixelValue(0)
# An Euler transformation is defined to position the input volume.
TransformType = itk.Euler3DTransform[itk.D]
transform = TransformType.New()
transform.SetComputeZYX(True)
translation = [tx, ty, tz]
# converting degrees into radians old val 1./180.*math.pi
dtr = (math.atan(1.0) * 4.0) / 180.0
transform.SetTranslation(translation)
transform.SetRotation(dtr * rx, dtr * ry, dtr * rz)
imOrigin = inputImage.GetOrigin()
imRes = inputImage.GetSpacing()
imRegion = inputImage.GetBufferedRegion()
imSize = imRegion.GetSize()
isocenter = [imOrigin[i] + (imRes[i]*imSize[i]/2.0) for i in range(3)]
transform.SetCenter(isocenter)
InterpolatorType = itk.SiddonJacobsRayCastInterpolateImageFunction[InputImageType, itk.D]
interpolator = InterpolatorType.New()
interpolator.SetProjectionAngle(dtr * rprojection) # Set angle between projection central axis and -z axis
interpolator.SetFocalPointToIsocenterDistance(scd) # Set source to isocenter distance
interpolator.SetThreshold(threshold) # Set intensity threshold, below which are ignored.
interpolator.SetTransform(transform)
interpolator.Initialize()
filter.SetInterpolator(interpolator)
# number of pixels of the 2D DRR image
size = [dx, dy, 1]
# pixel spacing of the 2D DRR image [mm]
spacing = [im_sx, im_sy, 1]
filter.SetSize(size)
filter.SetOutputSpacing(spacing)
# Use the image centers as the central axis position.
o2Dx = (dx - 1.) / 2.
o2Dy = (dy - 1.) / 2.
# Compute the origin (in mm) of the 2D Image
origin = [-im_sx * o2Dx, -im_sy * o2Dy, -scd]
filter.SetOutputOrigin(origin)
filter.Update()
RescaleType = itk.RescaleIntensityImageFilter[InputImageType, OutputImageType]
rescaler = RescaleType.New()
rescaler.SetOutputMinimum(0)
rescaler.SetOutputMaximum(255)
rescaler.SetInput(filter.GetOutput())
rescaler.Update()
# Out of some reason, the computed projection is upsided-down.
# Here we use a FilpImageFilter to flip the images in y direction.
FlipFilterType = itk.FlipImageFilter[OutputImageType]
flipFilter = FlipFilterType.New()
FlipAxesArray = itk.FixedArray[itk.D, 3]()
FlipAxesArray[0] = 0
FlipAxesArray[1] = 1
FlipAxesArray[2] = 0
flipFilter.SetFlipAxes(FlipAxesArray)
flipFilter.SetInput(rescaler.GetOutput())
flipFilter.Update()
WriteType = itk.ImageFileWriter[OutputImageType]
writer = WriteType.New()
writer.SetFileName(output_name)
writer.SetInput(flipFilter.GetOutput())
writer.Update()
def generateAndSaveDRR(dataFolder, outputFolder, numberOfProjections):
firstAngle = 0
angularArc = 360
all_img_vol_names = [x for x in os.listdir(dataFolder) if x.endswith(".nii.gz")]
for i, relative_path in enumerate(all_img_vol_names):
img_vol_path = os.path.join(dataFolder, relative_path)
drr_img_name = relative_path.rsplit('.', 2)[0]
for noProj in range(numberOfProjections):
angle = firstAngle + noProj * angularArc / numberOfProjections
drr_img_rel_path = drr_img_name + '_' + str(int(angle))
output_drr_Path = os.path.join(outputFolder, drr_img_rel_path + ".png")
siddonJacobRayTracingDRR(input_name=img_vol_path, output_name=output_drr_Path, rprojection=angle)
dataFolder = 'F:/all_CTs_HU_corrected'
outputFolder = 'F:/ITK_DRRs/all_drrs'
generateAndSaveDRR(dataFolder, outputFolder, 10)_______________________________________________
Rtk-users mailing list
rtk-us...@openrtk.org
https://www.creatis.insa-lyon.fr/mailman/listinfo/rtk-users
