Dear all,

 

I’m quite new to Python (and the use of RTK and ITK), but I have a project for which I want to create digitally reconstructed radiographs (DRRs) of a CT. I therefore have successfully installed RTK for python and I have run the FirstReconstruction.py example. Now I want to edit this so that it is possible to create DRRs, thus forwardprojection.

 

I have read in the archives of RTK-users and people are mostly directed towards the C++ code https://github.com/SimonRit/RTK/blob/b32cffdc6f9d7a432c50023c370ed996a7985b69/applications/rtkforwardprojections/rtkforwardprojections.cxx for forward projections. However, I find this very hard to read since I don’t know C++ so I hope you can help. I have also read https://discourse.itk.org/t/itk-to-rtk-migration-forward-projection-issues-questions/2107 regarding this matter, but I cannot reproduce this.

 

What I do in my code:

  • Load 3D CT image and change it to float Image Type ( sitk.sitkFloat32)
  • Define the image type I want for DRR
  • Create a stack of empty projection images using ConstantImageSource, using ImageType of DRR
  • Define the RTK geometry object (ThreeDCircularProjectionGeometry)
  • Then create JosephForwardProjectionImageFilter

Here, the problems arise.

The filter only seems to accept 3D imageTypes, so I have given it the [itk.F,3], [itk.F,3] input.

Then, if I try to setInput(0) with constantImageSource.GetOutput(), it throws the error:

 
Expecting argument of type itkImageF3 or itkImageSourceIF3.
Additional information:
Wrong number or type of arguments for overloaded function 'itkImageToImageFilterIF3IF3_SetInput'.

 

Probably since I have used the ImageType of the DRR (2D) for the constantImageSource.

Additionally, it seems that the setInput(1) = CT also doesn’t work, maybe because I have used sitk.sitkFloat32) to get it to a float image?

 

Does anyone know what I have to do here to solve this? And/or is there some python script available which does forward projection so I can have a look how it is done in Python?

 

I also don’t know what I would have to do after this step. Do I just do an FDK reconstruction etc. just as in the FirstReconstruction.py script?

 

Thank you in advance for anyone who can answer (some) of my questions.

 

M

 

 

 

 

 

My code:

import sys

import itk

from itk import RTK as rtk

import SimpleITK as sitk

 

# Utility method that either downloads data from the Girder repository or

# if already downloaded returns the file name for reading from disk (cached data).

%run update_path_to_download_script

from downloaddata import fetch_data as fdata

 

# Always write output to a separate directory, we don't want to pollute the source directory.

import os

OUTPUT_DIR = 'Output_test'

 

# Loading 3D CT image

CT = sitk.ReadImage(fdata("training_001_ct.mha"), sitk.sitkFloat32)

 

# Defines the image type CT

Dimension_CT = 3

PixelType = itk.F

ImageType_CT = itk.Image[PixelType, Dimension_CT]

 

# Defines the image type DRR

Dimension_DRR = 2

ImageType_DRR = itk.Image[PixelType, Dimension_DRR]

 

# Create a stack of empty projection images

ConstantImageSourceType = rtk.ConstantImageSource[ImageType_DRR]

constantImageSource = ConstantImageSourceType.New()

# Define origin, sizeOutput and spacing (still need to change these)

origin = [ -127, -127]

sizeOutput = [ 512, 512]

spacing = [ 0.653595, 2.0]

constantImageSource.SetOrigin( origin )

constantImageSource.SetSpacing( spacing )

constantImageSource.SetSize( sizeOutput )

constantImageSource.SetConstant(0.)

 

# Defines the RTK geometry object

geometry = rtk.ThreeDCircularProjectionGeometry.New()

numberOfProjections = 360

firstAngle = 0.

angularArc = 360.

sid = 600 # source to isocenter distance

sdd = 1200 # source to detector distance

for x in range(0,numberOfProjections):

  angle = firstAngle + x * angularArc / numberOfProjections

  geometry.AddProjection(sid,sdd,angle)

   

# Writing the geometry to disk

xmlWriter = rtk.ThreeDCircularProjectionGeometryXMLFileWriter.New()

xmlWriter.SetFilename ( sys.argv[2] )

xmlWriter.SetObject ( geometry )

xmlWriter.WriteFile()

 

REIType = rtk.JosephForwardProjectionImageFilter[ImageType_CT, ImageType_CT]

rei = REIType.New()

semiprincipalaxis = [ 50, 50, 50]

center = [ 0, 0, 10]

# Set GrayScale value, axes, center...

#rei.SetDensity(2)

#rei.SetAngle(0)

#rei.SetCenter(center)

#rei.SetAxis(semiprincipalaxis)

rei.SetGeometry( geometry )

rei.SetInput(0, constantImageSource.GetOutput())

rei.SetInput(1, CT)

 

 

_______________________________________________
Rtk-users mailing list
[email protected]
https://public.kitware.com/mailman/listinfo/rtk-users

Reply via email to