投影到DRR图像平面上的3D点

编程问答 2022-07-30

问题描述

我已经通过使用vtk工具对它进行体素化来创建了体积网格和相应的3D CT图像。然后使用龙门角度为0的SidddonJacobRay跟踪算法获得了3D CT图像的DRR图像。

现在我想获取DRR图像投影平面上3D网格的每个顶点坐标的对应2D投影坐标（即x，y位置）。

我经历过contributed note at the function documentation Algirhtm，但实际上它想输入3D CT图像体积并输出DRR图像。

如何获取给定3D坐标位置在DRR图像平面中的（x，y）坐标值？例如，我每个顶点都有（x，y，z）个位置，我想转换为（x，y），即DRR图像平面上的2D投影坐标吗？您能举例说明一下吗？

我刚刚将C ++代码GetDRRSiddonJacobsRayTracing转换为python版本，以获得给定3D顶点坐标的DRR投影，如下所示。

但是我认为此插值函数用于将当前体素坐标旋转到检测器坐标系定义的坐标系，因为我已经获得了三个坐标，即（x，y，z）作为给定3d的输出顶点坐标。

m_FocalPointToIsocenterdistance = 1000.; # Focal point to isocenter distance in mm.
m_ProjectionAngle = 0.;                  # Angle in radians betweeen projection central axis and reference axis
m_Threshold = 0.;

TransformType = itk.Euler3DTransform[itk.D]

# user can set an arbitrary transform for the volume. If not explicitly set,it should be identity.
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()
# constant for converting degrees into radians
dtr = (math.atan(1.0) * 4.0) / 180.0
m_CamRottransform.SetRotation(dtr * (-90.0),0.0,0.0)

m_ComposedTransform.SetIdentity()
m_ComposedTransform.Compose(m_Transform,False)

isocenter = m_Transform.GetCenter()
# 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,-m_ProjectionAngle)
m_GantryRottransform.SetCenter(isocenter)
m_ComposedTransform.Compose(m_GantryRottransform,False)

# An Euler 3D transfrom is used to shift the source to the origin.
focalpointtranslation = [None] * 3
focalpointtranslation[0] = -isocenter[0]
focalpointtranslation[1] = m_FocalPointToIsocenterdistance - isocenter[1]
focalpointtranslation[2] = -isocenter[2]

m_CamShiftTransform.SetTranslation(focalpointtranslation)
m_ComposedTransform.Compose(m_CamShiftTransform,False)

# A 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)

# You probably need to discard one of these 3 coordinates. 
# it should be easier for to figure out which one is depth,and which two are X and Y.
# Coordinate of a DRR pixel in the world coordinate system
drrPixelWorld = m_InverseTransform.TransformPoint(point)

解决方法

暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！

如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。

小编邮箱:dio#foxmail.com (将#修改为@）

2d3dgeometryitkitkplane