[Rtk-users] Ray-tracing through a binary mask image volume

[Suranga W]

Hi,

I have a CT volume and its corresponding binary mask image volume of the
liver.

I just want to ray-trace through a binary mask image volume to generate
mask-projection DRRs.

The DRRs can be sucessfully genrated from the CT volume, but when I passed
the binary mask volume, I did not get the correct result.

Could anyone please assist me in resolving this problem ?  What kind of
change should I make to generate these DRRs ?

The code, resulted DRR versions for CT and mask, CT volume and the binary
mask volume attached herewith for your further reference.

 reference_CT.nii.gz
<https://drive.google.com/file/d/1p6ODn6gK11v2pU1OmRd-8SyXG5nwzZf6/view?usp=drive_web>
 reference_liver_binary_mask.nii.gz
<https://drive.google.com/file/d/1I3J3zUtdgT6ge6kY8EspiPjIqhQi8-eh/view?usp=drive_web>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.creatis.insa-lyon.fr/pipermail/rtk-users/attachments/20221128/5b33a3ff/attachment.htm>
-------------- next part --------------
import sys
import itk
from itk import RTK as rtk
import itk
import numpy as np

def generateDRRs(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()


ct_file = './reference_CT.nii.gz'
binary_mask_file = './reference_liver_binary_mask.nii.gz'

out_ct_drr_file = './ct_drr_324.png'
out_mask_drr_file = './mask_drr_324.png'

offset_X = 117.085998535
offset_Y = 3.948999882

# generate DRRs fro CT volume
generateDRRs(input_file=ct_file, output_file=out_ct_drr_file, projOffsetX=offset_X, projOffsetY=offset_Y, gantryAngle=324)

# generate DRRs from binary mask volume
generateDRRs(input_file=binary_mask_file, output_file=out_mask_drr_file, projOffsetX=offset_X, projOffsetY=offset_Y, gantryAngle=324)
-------------- next part --------------
A non-text attachment was scrubbed...
Name: mask_drr_324.png
Type: image/png
Size: 6889 bytes
Desc: not available
URL: <http://www.creatis.insa-lyon.fr/pipermail/rtk-users/attachments/20221128/5b33a3ff/attachment.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: ct_drr_324.png
Type: image/png
Size: 35943 bytes
Desc: not available
URL: <http://www.creatis.insa-lyon.fr/pipermail/rtk-users/attachments/20221128/5b33a3ff/attachment-0001.png>