[Rtk-users] ITK ERROR: Geometry is empty error when generating DRRs

Mon Dec 19 13:49:29 CET 2022

Hi Simon,

I created the geometry by adjusting the source position, detector position,
and detector u v coordinates. However, when I attempt to rotate it for
generating differnt DRRs, the rotation is not in the correct axis. I tried
setting different U V coodinates but had no luck. How can I fix this by
altering the U V coodinates ?

I'm asking because I'd like to change the geometry settings rather than
changing the physical coordinates of the CT image volume (that is, without
changing the origin or/and the orientation). The gantry angle zero results
were appeared as expected with this geomtry settings and were well aligned
with the actual kV xray configuration, i.e. with the field of view
alignment (both vertical and horizontal alignment). Now I wanted to make
DRRs by rotating the gantry in the correct axis.

Please assist me in resolving this problem.

On Thu, Dec 15, 2022 at 2:06 PM Suranga W <isurusuranga.wijesinghe at gmail.com>
wrote:

> Hi Simon,
>
> Thnaks for the response. I fixed the field of view problem issue which I
> mentioned earlier.
>
> Now I'm dealing with the gantry rotation problem for generating DRRs. The
> gantry angle of zero is exactly correct. However, when I changed the gantry
> angle, I discovered that the gantry was not rotating on the correct axis.
> Could you please guide me on resolving this issue please ?
>
> I herewith attached the code, drr at gantry angle 0 and 324 (in degrees)
> for your further reference.
>
> On Wed, Dec 14, 2022 at 8:41 AM Simon Rit <simon.rit at creatis.insa-lyon.fr>
> wrote:
>
>> Hi,
>> 1) AddProjection with source / detector positions + u/v decomposes the
>> provided information in a set of 9 parameters (SID, SDD, source / detector
>> offsets in the u / v directions and 3 angles). Maybe try to retrieve these
>> parameters and pass them changing the GantryAngle only? The getters are
>> available here
>> <https://github.com/SimonRit/RTK/blob/master/include/rtkThreeDCircularProjectionGeometry.h#L152-L212>.
>> The problem will likely be that the parameters are defined in a given
>> convention and what you call GantryAngle might not be equal to that of RTK.
>> 2) No, you'll have to play with small tests and find out by yourself.
>> Simon
>>
>> On Tue, Dec 13, 2022 at 11:19 AM Suranga W <
>> isurusuranga.wijesinghe at gmail.com> wrote:
>>
>>>  MCR.nii.gz
>>> <https://drive.google.com/file/d/1-fv97l5u0J_-P0CoO7gMVpqJpMHXUOHB/view?usp=drive_web>
>>> Hi Simon,
>>>
>>> 1.) So, if I want to generate different DRRs by changing the gantry
>>> angle, how can I do so with these configurations?
>>> For example, if I invoke the other overloaded AddProjection method which
>>> can pass the gantry angle, I am unable to set the machine parameters as
>>> expected. In that case, how can I do it?
>>>
>>> 2.) How can I change the origin and size of teh input CT image if I only
>>> have the field-of-view origin at detector ? Could you please elaborate with
>>> sample codes if possible for these two cases?
>>> I herewith attached a sample CT volume for your further consideration.
>>>
>>> Could you please assist me in resolving these two problems ? Thanks in
>>> advance for your time and consideration.
>>>
>>>
>>>
>>> On Mon, Dec 12, 2022 at 3:08 PM Simon Rit <
>>> simon.rit at creatis.insa-lyon.fr> wrote:
>>>
>>>> Hi,
>>>> 1) You can't. This method uses another parametrization.
>>>> 2) You can set the field-of-view parameters by adjusting the origin and
>>>> size of the input image.
>>>> Simon
>>>>
>>>> On Mon, Dec 12, 2022 at 3:59 PM Suranga W <
>>>> isurusuranga.wijesinghe at gmail.com> wrote:
>>>>
>>>>> Hi Simon,
>>>>>
>>>>> Many thanks. I have another few questions though.
>>>>>
>>>>> 1.) How can I pass different gantry angles if I use the aforementioned
>>>>> AddProjection method, which lacks the gantryAngle parameter? What I want is
>>>>> to generate DRR images for angles 36, 72, 108, and so on.
>>>>>
>>>>> 2.) If I have the field of view origin from a corresponding kV X-ray
>>>>> dicom image (let's say for gantry angle zero Field of View Origin at
>>>>> detector is [116.853836060, 4.566089153] in pixels see
>>>>> <https://dicom.innolitics.com/ciods/enhanced-xa-image/enhanced-xa-image-multi-frame-functional-groups/52009229/00189432/00187030>),
>>>>> how can I use them to produce DRRs with aforementioned machine
>>>>> configurations in AddProjection method ?
>>>>>
>>>>>
>>>>>
>>>>> On Mon, Dec 12, 2022 at 8:43 AM Simon Rit <
>>>>> simon.rit at creatis.insa-lyon.fr> wrote:
>>>>>
>>>>>> Hi,
>>>>>> You might have missed an earlier message
>>>>>> WARNING: In
>>>>>> /home/srit/src/itk/itk/Modules/Remote/RTK/src/rtkThreeDCircularProjectionGeometry.cxx,
>>>>>> line 210
>>>>>> DataObject (0x55d436d7ee10): Failed to AddProjection
>>>>>> The problem is that the two vectors describing the projection
>>>>>> orientation are not unit vector. You can correct this with:
>>>>>> geometry.AddProjection([-60.0, 59.0 + sid, 34.0],[-60.0 + 115.,
>>>>>> 59.0+sid-sdd, 34.0],[1./math.sqrt(2.),1./math.sqrt(2.),0.],[0.,0.,1.])
>>>>>> Simon
>>>>>>
>>>>>> On Fri, Dec 9, 2022 at 1:43 PM Suranga W <
>>>>>> isurusuranga.wijesinghe at gmail.com> wrote:
>>>>>>
>>>>>>> Hi,
>>>>>>>
>>>>>>> I attempted to generate DRRs by adding the geometry configuration to
>>>>>>> the ThreeDCircularProjectionGeometry object, as shown in the attached code
>>>>>>> (python based code). However, I encountered a runtime error as below
>>>>>>> (please look at the end of the email).
>>>>>>>
>>>>>>> On the other hand, how can I pass the relevant gantry angle to it
>>>>>>> for generating DRRs at different angles if I add geometry configuration
>>>>>>> using the addProjection method as mentioned in my code?
>>>>>>>
>>>>>>> What I mean is that if I use the code below to set machine
>>>>>>> parameters, I'm not able to add the corresponding gantry angle for
>>>>>>> generating DRRs.
>>>>>>>
>>>>>>> *geometry.AddProjection([-60.0, 59.0 + sid, 34.0],[-60.0 + 115.,
>>>>>>> 59.0+sid-sdd, 34.0],[1.,1.,0.],[0.,0.,1.])*
>>>>>>>
>>>>>>> where [-60.0, 59.0, 34.0] represents the patient's isocenter.
>>>>>>>
>>>>>>> However, if I use the code below, I can use the gantry angle but I
>>>>>>> cannot set the machine geometry configuration.
>>>>>>>
>>>>>>> *geometry.AddProjection(sid, sdd, gantryAngle, projOffsetX,
>>>>>>> projOffsetY, outOfPlaneAngle, inPlaneAngle, sourceOffsetX, sourceOffsetY)*
>>>>>>>
>>>>>>> I attached a CT volume for your further reference.
>>>>>>>
>>>>>>> Please assist me in resolving this issue.
>>>>>>>
>>>>>>> ---------------------------------------------------------------------------RuntimeError                              Traceback (most recent call last)
>>>>>>> Cell In [29], line 45     43 rei.SetInput(0, constantImageSource.GetOutput())     44 rei.SetInput(1, CT)---> 45 rei.Update()     47 Dimension = 3     48 OutputPixelType = itk.UC
>>>>>>> RuntimeError: D:\a\im\include\rtkProjectionsRegionConstIteratorRayBased.hxx:88:
>>>>>>> ITK ERROR: Geometry is empty, cannot determine iterator type.
>>>>>>>
>>>>>>>  drr_projection_with_diff_angles.txt
>>>>>>> <https://drive.google.com/file/d/1Ags89Lpczvrk_rdHpGSqDcBrashzNSZr/view?usp=drive_web>
>>>>>>>  sample_CT.nii.gz
>>>>>>> <https://drive.google.com/file/d/1SGh0cjspFwTj2v8qTLXtcCfbUvyUENS5/view?usp=drive_web>
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.creatis.insa-lyon.fr/pipermail/rtk-users/attachments/20221219/4d3c6524/attachment.htm>
-------------- next part --------------
import os
import math
import itk
from itk import RTK as rtk
import numpy as np
import math

# Loading 3D CT image
CT = itk.imread("./sample_CT.nii.gz", pixel_type=itk.F)

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

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

# Defines the RTK geometry object
geometry = rtk.ThreeDCircularProjectionGeometry.New()

sid=1000. # source to isocenter distance
sdd=1536. # source to detector distance

# isocenter of the patient is = [-60.0, 59.0, 34.0]
geometry.AddProjection([-60.0, 59.0 + sid, 34.0],[-60.0-115., 59.0+sid-sdd, 34.0],[1.,0.,0.],[0.,0.,1.])

inPlaneAngle = math.degrees(geometry.GetInPlaneAngles()[0])
outOfPlaneAngle = math.degrees(geometry.GetOutOfPlaneAngles()[0])
sourceOffsetX = geometry.GetSourceOffsetsX()[0]
sourceOffsetY = geometry.GetSourceOffsetsY()[0]
projOffsetX = geometry.GetProjectionOffsetsX()[0]
projOffsetY = geometry.GetProjectionOffsetsY()[0]
sourceToIsoDis = geometry.GetSourceToIsocenterDistances()[0]
sourceToDetDis = geometry.GetSourceToDetectorDistances()[0]

print('inPlaneAngle:', inPlaneAngle)
print('outOfPlaneAngle:', outOfPlaneAngle)
print('sourceOffsetX:', sourceOffsetX)
print('sourceOffsetY:', sourceOffsetY)
print('projOffsetX:', projOffsetX)
print('projOffsetY:', projOffsetY)
print('sourceToIsoDis:', sourceToIsoDis)
print('sourceToDetDis:', sourceToDetDis)
print('gantry:', geometry.GetGantryAngles()[0])
geometry.Clear()
geometry.AddProjection(sourceToIsoDis, sourceToDetDis, 324., 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.
# 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('./drr_324.png')
writer.SetInput(flipFilter.GetOutput())
writer.Update()