[Rtk-users] DRR with limited FOV from 3D-CT

Fri Nov 11 11:42:46 CET 2022

Hi Simon,

I cannot change the format of the image volume since I use another tool to
generate these nifti image volumes and corresponding DVFs in the same
format.
I want to use RTK to generate DRRs for these nifti instances with differnt
projection angles.

Therefore I first reorient the CT volme to have RAS orientation which is
now the direction matrix is preserved and the point (0,0,0) mm is in the
image .

However, the genreted DRRs do not appear correctly.

The previous example code you provided works correctly, but the origin is
reordered.

If I use this RAS-oriented new CT, the DRRs should be in the correct
versions. What could be the cause ?
 RAS_oriented_ref_CT.nii.gz
<https://drive.google.com/file/d/1ypdqeQQpbvjy7k4T_zxY3e8g3-fhmI5q/view?usp=drive_web>

On Thu, Nov 10, 2022 at 7:22 AM Simon Rit <simon.rit at creatis.insa-lyon.fr>
wrote:

> Hi,
> ref_CT_updated.nii.gz is still not oriented correctly, please check the
> RTK geometry doc more carefully. If I use the following lines
> 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)
> CTarray  += 1000
>
> and set the outOfPlaneAngle to 0, the result looks more as expected. I
> advise you to avoid the nifti file format which has many issues according
> to ITK discussions on github, prefer mha or mhd.
> Simon
>
> On Wed, Nov 9, 2022 at 5:42 PM Isuru Suranga Wijesinghe <
> isurusuranga.wijesinghe at gmail.com> wrote:
>
>> Hi Simon,
>>
>> Many thanks. I'm in need of some assistance once more please.
>>
>> First, I reoreinted the CT volume and used the field of view origin in
>> the corresponding real kV dicom file (i.e. with projection angle zero) to
>> adjust the correct field of view for projection angle zero (FOV origin
>> [116.368270874, 3.837730408]). I simply replace the default values of
>> projOffsetX and projOffsetY with the values shown below.
>> projOffsetY = -3.837730408
>>
>> I simply rotated around the x-axis by setting the outOfPlaneAngle
>> parameter to 90 degrees.
>>
>> However, If I project through a 324 (degrees) projection angle (with FOV
>> origin [117.085998535, 3.948999882]), it will rotate in the worng axis.
>> For this step, I simply replace the gantryAngle parameter with 324 and
>> set the new projOffsetX and projOffsetY.
>>
>> I've attached the updated code, the reoriented CT volume, the correct DRR
>> at projection angle zero, and the incorrect DRR at projection angle 324.
>>
>> Can you assist me in resolving that error ?
>>
>>  drr_0.png
>> <https://drive.google.com/file/d/1qOoMQ6jvOSY-7KivjmcLuAoIIMhBafmE/view?usp=drive_web>
>>  drr_324.png
>> <https://drive.google.com/file/d/1RjT6_hGF2vkog3YNUfv-03RAmTClrSbK/view?usp=drive_web>
>>  ref_CT_updated.nii.gz
>> <https://drive.google.com/file/d/1BntJvZ7JRF9kTxQsycB5JoipWVFESNQ_/view?usp=drive_web>
>>  rtk_drr_projection_code.txt
>> <https://drive.google.com/file/d/1mEaIbj9NjuU5BZB4NKlmv0S9CUN0FME3/view?usp=drive_web>
>>
>>
>>
>>
>> On Wed, Nov 9, 2022 at 2:26 PM Simon Rit <simon.rit at creatis.insa-lyon.fr>
>> wrote:
>>
>>> Hi,
>>> The first input of JosephForwardProjectionImageFilter is the image in
>>> which you'll calculate the DRR, typically an empty projection (filled with
>>> zeros). In your code, this is computed with constantImageSource but you're
>>> setting the volume origin, spacing size which is probably wrong. For an
>>> Elekta projection (you showed an Elekta drawing), this would be
>>> constantImageSource.SetOrigin( [-204.4,-204.4,0] )
>>> constantImageSource.SetSpacing( [0.8,0.8,0.8] )
>>> constantImageSource.SetSize( [512, 512,1]  )
>>> Besides, for the CT volume you sent, point (0,0,0) mm is not in the
>>> image and this is the default center of rotation in RTK (see doc
>>> <http://www.openrtk.org/Doxygen/DocGeo3D.html>). Moreover the
>>> orientation is wrong. So either you change the origin and the direction to
>>> use RTK convention or you build your geometry differently, e.g. by setting
>>> source center, detector origin position and detector u v coordinates:
>>> geometry.AddProjection([-57.,20.+sid,25.],[-57+115.,20.+sid-sdd
>>> ,25.],[1.,1.,0.],[0.,0.,1.])
>>> where the coordinates (-57,28,25) are those of a point somewhere in the
>>> center of the patient. This will give you line integrals, I don't know what
>>> is the content of your png file but you might have to take the exponential
>>> to have something close to real data. you also have to add 1000 to yourCT
>>> numbers to put air at 0.
>>> In short, read the geometry doc and don't expect us to be able to match
>>> exactly the DRR you sent... This will require some work on your side
>>> Cheers,
>>> Simon
>>>
>>> On Wed, Nov 9, 2022 at 11:32 AM Isuru Suranga Wijesinghe <
>>> isurusuranga.wijesinghe at gmail.com> wrote:
>>>
>>>> Can I ask any possible solution to my issue please (I just want to get
>>>> the DRR with the limited FOV based on the machine parameters attached
>>>> above) ???
>>>>
>>>> On Tue, Nov 8, 2022 at 2:12 PM Isuru Suranga Wijesinghe <
>>>> isurusuranga.wijesinghe at gmail.com> wrote:
>>>>
>>>>> In addition to that I want to generate DRRs with a limited field of
>>>>> view as attached in the sample real kV Xray image.
>>>>> I've included a sample CT volume that I'm currently using to generate
>>>>> such DRR images using the machine parameters.  Please assist in resolving
>>>>> this matter. Thank you in advance for your consideration and time.
>>>>>
>>>>>  ref_CT_volume.nii.gz
>>>>> <https://drive.google.com/file/d/1PilfGAhqWe-NFD4Up1CYwSHqDXns04CV/view?usp=drive_web>
>>>>>  sample_real_kv_xray.png
>>>>> <https://drive.google.com/file/d/1zX2HUEEL4-lGxS4daDrhXuhTy1uFBY89/view?usp=drive_web>
>>>>>
>>>>>
>>>>> On Tue, Nov 8, 2022 at 11:23 AM Isuru Suranga Wijesinghe <
>>>>> isurusuranga.wijesinghe at gmail.com> wrote:
>>>>>
>>>>>> Hi,
>>>>>>
>>>>>> I have a 3D-CT volume and want to obtain the DRR projection by
>>>>>> adjusting the machine parameters as shown in the image below.
>>>>>>
>>>>>> Here's my code, and I didn't get the expected DRRs with the limited
>>>>>> FOV. Is there anything missing from the parameter setup?
>>>>>>
>>>>>> Could you please assist me in resolving this issue? Any code level
>>>>>> assistance or advice would be greatly appreciated.
>>>>>>
>>>>>> # Loading 3D CT image
>>>>>> CT = itk.imread("./ct.nii.gz", pixel_type=itk.F)
>>>>>>
>>>>>> # Defines the image type
>>>>>> Dimension_CT = 3
>>>>>> PixelType = itk.F
>>>>>> ImageType = itk.Image[PixelType, Dimension_CT]
>>>>>>
>>>>>> origin_x, origin_y, origin_z = CT.GetOrigin()
>>>>>> space_x, space_y, space_z = CT.GetSpacing()
>>>>>>
>>>>>> # Define origin, sizeOutput and spacing (still need to change these)
>>>>>> origin = [origin_x, origin_y, origin_z]
>>>>>> spacing = [space_x, space_y, space_z]
>>>>>>
>>>>>> sizeOutput = [ 512, 512, 1]
>>>>>>
>>>>>> # Create a stack of empty projection images
>>>>>> ConstantImageSourceType = rtk.ConstantImageSource[ImageType]
>>>>>> constantImageSource = ConstantImageSourceType.New()
>>>>>>
>>>>>> constantImageSource.SetOrigin( origin )
>>>>>> constantImageSource.SetSpacing( spacing )
>>>>>> constantImageSource.SetSize( sizeOutput )
>>>>>> constantImageSource.SetConstant(0.)
>>>>>>
>>>>>> # Defines the RTK geometry object
>>>>>> geometry = rtk.ThreeDCircularProjectionGeometry.New()
>>>>>> firstAngle = 0.
>>>>>> angularArc = 360.
>>>>>> sid = 1000 # source to isocenter distance
>>>>>> sdd = 1536 # source to detector distance
>>>>>> gantryAngle = 0. # rot around y-axis
>>>>>> projOffsetX = 0.
>>>>>> projOffsetY = 0.
>>>>>> outOfPlaneAngle = 90. # rot around x-axis
>>>>>> inPlaneAngle = 0. # rot around z-axis
>>>>>> sourceOffsetX = 0.
>>>>>> sourceOffsetY = 0.
>>>>>> geometry.AddProjection(sid, sdd, gantryAngle)
>>>>>>
>>>>>> 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()
>>>>>>
>>>>>> WriteType = itk.ImageFileWriter[OutputImageType]
>>>>>> writer = WriteType.New()
>>>>>> writer.SetFileName('./drr_0.png')
>>>>>> writer.SetInput(rescaler.GetOutput())
>>>>>> writer.Update()
>>>>>>
>>>>>>
>>>>>> _______________________________________________
>>>> Rtk-users mailing list
>>>> Rtk-users at public.kitware.com
>>>> https://public.kitware.com/cgi-bin/mailman/listinfo/rtk-users
>>>>
>>>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.creatis.insa-lyon.fr/pipermail/rtk-users/attachments/20221111/a1b108b5/attachment.htm>
-------------- next part --------------
# Loading 3D CT image
CT = itk.imread("./RAS_oriented_ref_CT.nii.gz", pixel_type=itk.F)

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

constantImageSource.SetOrigin([-204.4,-204.4,0])
constantImageSource.SetSpacing([0.8,0.8,0.8])
constantImageSource.SetSize([512, 512,1])
constantImageSource.SetSize(sizeOutput)
constantImageSource.SetConstant(0.)

# Defines the RTK geometry object
geometry = rtk.ThreeDCircularProjectionGeometry.New()
firstAngle = 0.
angularArc = 360.
sid = 1000 # source to isocenter distance
sdd = 1536 # source to detector distance
gantryAngle = 324. # rot around y-axis
projOffsetX = 0.
projOffsetY = 0.
outOfPlaneAngle = 90. # rot around x-axis
inPlaneAngle = 0. # rot around z-axis
sourceOffsetX = 0.
sourceOffsetY = 0.

geometry.AddProjection(sid, sdd, gantryAngle, -117.085998535, -3.948999882, 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.
# take this part from the ITK drr projection code
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('./stack_324.png')
writer.SetInput(flipFilter.GetOutput())
writer.Update()