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

Thu Dec 15 15:06:48 CET 2022

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/20221215/7d2ec014/attachment.htm>
-------------- next part --------------
# 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()
-------------- next part --------------
A non-text attachment was scrubbed...
Name: drr_angle_0.png
Type: image/png
Size: 36413 bytes
Desc: not available
URL: <http://www.creatis.insa-lyon.fr/pipermail/rtk-users/attachments/20221215/7d2ec014/attachment.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: drr_angle_324.png
Type: image/png
Size: 32232 bytes
Desc: not available
URL: <http://www.creatis.insa-lyon.fr/pipermail/rtk-users/attachments/20221215/7d2ec014/attachment-0001.png>