[Rtk-users] 3D mesh points projection onto the DRR image plane
Suranga W
isurusuranga.wijesinghe at gmail.com
Sun Oct 1 12:30:03 CEST 2023
Hi Simon,
I attempted to implement point projection based on the information you
mentioned above. However, the points projected onto the DRR do not appear
to be accurate. Could you kindly assist me in identifying if I missed any
steps or how I should incorporate parameters such as origin, pixel spacing,
and size (which I used when implementing point projection in ITK - code is
attached above) to ensure correct point projection? For your reference, I
have attached a DRR image along with the projected point cloud.
Thanks,
Surang
On Thu, Sep 28, 2023 at 8:50 AM Simon Rit <simon.rit at creatis.insa-lyon.fr>
wrote:
> Hi,
> The easiest is to use the projection matrices from the geometry object
> with GetMatrix
> <http://www.openrtk.org/Doxygen/classrtk_1_1ProjectionGeometry.html#ac848531d266395624ee8f1f328ae67d3>
> which takes the projection number as a paramter or GetMatrices
> <http://www.openrtk.org/Doxygen/classrtk_1_1ProjectionGeometry.html#aa00734aa99a710a163a0eb30334cbac7>
> to get the full list. You can then multiply your 3D homogeneous coordinate
> in image space (x,y,z,1) by multiplying it with the matrix to get the
> homogenous projection coordinates (u,v,w) (so you get the actual 2D
> coordinates (u/w, v/w)).
> Simon
>
> On Mon, Sep 25, 2023 at 8:29 AM Suranga W <
> isurusuranga.wijesinghe at gmail.com> wrote:
>
>> 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
>> _______________________________________________
>> Rtk-users mailing list
>> rtk-users at openrtk.org
>> https://www.creatis.insa-lyon.fr/mailman/listinfo/rtk-users
>>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.creatis.insa-lyon.fr/pipermail/rtk-users/attachments/20231001/d7b1f966/attachment-0001.htm>
-------------- next part --------------
# pip3 install itk-rtk
import os
import itk
from itk import RTK as rtk
import numpy as np
import pydicom as dicom
import glob
from PIL import Image, ImageDraw
import meshio
def generateDRRsRTK(input_file=None, points=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]
# Defines the RTK geometry object
geometry = rtk.ThreeDCircularProjectionGeometry.New()
geometry.AddProjection(sid, sdd, gantryAngle, -projOffsetX, projOffsetY, outOfPlaneAngle, inPlaneAngle, sourceOffsetX, sourceOffsetY)
matrix = geometry.GetMatrix(0)
print(matrix)
# this entire bottom part should be iterated over all the points
point_list = points.tolist()
# converted (u,v,w) => (u,v,w,1)
for nested_point_list in point_list:
nested_point_list.append(1)
u = np.array([])
v = np.array([])
for point in point_list:
# coordinate to the coordinate system defined by the detector coordinate system
drrPixelWorld = np.matmul(matrix, point)
#print(drrPixelWorld)
drrPixelWorld[0] = drrPixelWorld[0]/drrPixelWorld[2]
drrPixelWorld[1] = drrPixelWorld[1]/drrPixelWorld[2]
u = np.append(u, drrPixelWorld[0])
v = np.append(v, drrPixelWorld[1])
print('u:', u)
print('v:', v)
v = (dy - 1.) - v
return u, v
def from_meshio(mesh):
r"""Converts a :.msh file to a
:class:`torch_geometric.data.Data` instance.
Args:
mesh (meshio.read): A :obj:`meshio` mesh.
"""
if meshio is None:
raise ImportError('Package `meshio` could not be found.')
pos = torch.from_numpy(mesh.points).to(torch.float)
tetra = torch.from_numpy(mesh.cells_dict['tetra']).to(torch.long).t().contiguous()
return Data(pos=pos, tetra=tetra)
def getAndSaveDRR(ctPath, rawReaKVPath, refMeshPath, drr_img_path, output_img_path):
drr_img_name = 'ref_'
ref_mesh = meshio.read(refMeshPath)
ds = dicom.dcmread(rawReaKVPath)
angle = float(ds.get("PositionerPrimaryAngle"))
print('angle:', angle)
# read FOV origin
fovOriginX = float(ds.get("FieldOfViewOrigin")[0])
fovOriginY = float(ds.get("FieldOfViewOrigin")[1])
gtX, gtY = generateDRRsRTK(input_file=ctPath, points=ref_mesh.points, gantryAngle=angle, projOffsetX=fovOriginX, projOffsetY=fovOriginY)
gtX = gtX.tolist()
gtY = gtY.tolist()
i = Image.open(drr_img_path).convert("RGBA")
draw = ImageDraw.Draw(i)
w, h = i.size
# plot U, V points
for gt_x, gt_y in zip(gtX, gtY):
draw.point((gt_x, gt_y), fill="red")
i.save(output_img_path)
rawReaKVPath = './real_kv/IDAP775756_20191028_142329.719211_Image_75.dcm'
ctPath = './ref_CT/MCR.nii.gz'
ref_mesh_path = './REF_MESH_GEN/mesh/MCR.vtk'
drr_img_path = './ref_drr/MCR_290.922225952.png'
output_img_path = './out/MCR_out_290.922225952.png'
getAndSaveDRR(ctPath, rawReaKVPath, ref_mesh_path, drr_img_path, output_img_path)
-------------- next part --------------
A non-text attachment was scrubbed...
Name: MCR_out_290.922225952.png
Type: image/png
Size: 54379 bytes
Desc: not available
URL: <http://www.creatis.insa-lyon.fr/pipermail/rtk-users/attachments/20231001/d7b1f966/attachment-0001.png>
More information about the Rtk-users
mailing list