SimpleRTK: Difference between revisions

From Openrtk
Jump to navigation Jump to search
No edit summary
No edit summary
Line 26: Line 26:


<source lang="python">
<source lang="python">
  #!/usr/bin/env python
#!/usr/bin/env python
  from __future__ import print_function
from __future__ import print_function
  import SimpleRTK as srtk
import SimpleRTK as srtk
  import sys
import sys
  import os
import os
import matplotlib.pyplot as plt
import matplotlib.cm as cm


  if len ( sys.argv ) < 2:
if len ( sys.argv ) < 2:
     print( "Usage: RTKFirstReconstruction <output>" )
     print( "Usage: RTKFirstReconstruction <output>" )
     sys.exit ( 1 )
     sys.exit ( 1 )


  # Defines the RTK geometry object
# Defines the RTK geometry object
  geometry = srtk.ThreeDimCircularProjectionGeometry()
geometry = srtk.ThreeDimCircularProjectionGeometry()
  numberOfProjections = 360
numberOfProjections = 360
  firstAngle = 0
firstAngle = 0
  angularArc = 360
angularArc = 360
  sid = 600 # source to isocenter distance in mm
sid = 600 # source to isocenter distance in mm
  sdd = 1200 # source to detector distance in mm
sdd = 1200 # source to detector distance in mm
  isox = 0 # X coordinate on the projection image of isocenter
isox = 0 # X coordinate on the projection image of isocenter
  isoy = 0 # Y coordinate on the projection image of isocenter
isoy = 0 # Y coordinate on the projection image of isocenter
  for x in range(0,numberOfProjections):
for x in range(0,numberOfProjections):
    angle = firstAngle + x * angularArc / numberOfProjections
  angle = firstAngle + x * angularArc / numberOfProjections
    geometry.AddProjection(sid,sdd,angle,isox,isoy)
  geometry.AddProjection(sid,sdd,angle,isox,isoy)


  constantImageSource = srtk.ConstantImageSource()
constantImageSource = srtk.ConstantImageSource()
  origin = [ -127.0, -127.0, -127.0 ]
origin = [ -127.5, -127.5, 0. ]
  sizeOutput = [ 256, 256,  numberOfProjections ]
sizeOutput = [ 256, 256,  numberOfProjections ]
  spacing = [ 1.0, 1.0, 1.0 ]
spacing = [ 1.0, 1.0, 1.0 ]
  constantImageSource.SetOrigin( origin )
constantImageSource.SetOrigin( origin )
  constantImageSource.SetSpacing( spacing )
constantImageSource.SetSpacing( spacing )
  constantImageSource.SetSize( sizeOutput )
constantImageSource.SetSize( sizeOutput )
  constantImageSource.SetConstant(0.0)
constantImageSource.SetConstant(0.0)
  source = constantImageSource.Execute()
source = constantImageSource.Execute()


  rei = srtk.RayEllipsoidIntersectionImageFilter()
rei = srtk.RayEllipsoidIntersectionImageFilter()
  semiprincipalaxis = [ 50, 50, 50]
semiprincipalaxis = [ 50, 50, 50]
  center = [ 0, 0, 0]
center = [ 0, 0, 0]
  # Set GrayScale value, axes, center...
# Set GrayScale value, axes, center...
  rei.SetDensity(20)
rei.SetDensity(20)
  rei.SetAngle(0)
rei.SetAngle(0)
  rei.SetCenter(center)
rei.SetCenter(center)
  rei.SetAxis(semiprincipalaxis)
rei.SetAxis(semiprincipalaxis)
  rei.SetGeometry( geometry )
rei.SetGeometry( geometry )
  reiImage = rei.Execute(source)
reiImage = rei.Execute(source)
    
    
   # Create reconstructed image
    
  constantImageSource2 = srtk.ConstantImageSource()
# Create reconstructed image
  sizeOutput = [ 256, 256, 256 ]
constantImageSource2 = srtk.ConstantImageSource()
  constantImageSource2.SetOrigin( origin )
origin = [ -63.5, -63.5, -63.5 ]
  constantImageSource2.SetSpacing( spacing )
sizeOutput = [ 128, 128, 128 ]
  constantImageSource2.SetSize( sizeOutput )
constantImageSource2.SetOrigin( origin )
  constantImageSource2.SetConstant(0.0)
constantImageSource2.SetSpacing( spacing )
  source2 = constantImageSource2.Execute()
constantImageSource2.SetSize( sizeOutput )
constantImageSource2.SetConstant(0.0)
source2 = constantImageSource2.Execute()


  print("Performing reconstruction")
print("Performing reconstruction")
  feldkamp = srtk.FDKConeBeamReconstructionFilter()
feldkamp = srtk.FDKConeBeamReconstructionFilter()
  feldkamp.SetGeometry( geometry );
feldkamp.SetGeometry( geometry );
  feldkamp.SetTruncationCorrection(0.0);
feldkamp.SetTruncationCorrection(0.0);
  feldkamp.SetHannCutFrequency(0.0);
feldkamp.SetHannCutFrequency(0.0);
  image = feldkamp.Execute(source2,reiImage)  
image = feldkamp.Execute(source2,reiImage)  
    
    
  pixelID = image.GetPixelIDValue()
pixelID = image.GetPixelIDValue()
 
plt.imshow(srtk.GetArrayFromImage(image[:,64,:]), cmap = cm.Greys_r)


  caster = srtk.CastImageFilter()
caster = srtk.CastImageFilter()
  caster.SetOutputPixelType( pixelID )
caster.SetOutputPixelType( pixelID )
  image = caster.Execute( image )
image = caster.Execute( image )


  writer = srtk.ImageFileWriter()
writer = srtk.ImageFileWriter()
  writer.SetFileName ( sys.argv[1] )
writer.SetFileName ( sys.argv[1] )
  writer.Execute ( image );
writer.Execute ( image );
</source>
</source>



Revision as of 09:50, 26 September 2014

Introduction

SimpleRTK provides a simple interface to RTK in a variety of languages. For now only Python and C# binding have been tested and are supported. SimpleRTK will be provided as part of RTK making it easier to configure, build and use. Note that SimpleRTK is derived from the SimpleITK project for the Insight Toolkit, however it's under its own project.

Building RTK with SimpleRTK

SimpleRTK uses Swig to generate the wrapping for different languages. As the version of Swig should be controlled, SimpleRTK will automatically download the right version of Swig and its dependencies for the different platforms. However you can always point to the installed version of Swig if you'd like to, but this is not recommended.

Here are the simple steps to build and configure wrapping for RTK:

  1. Run CMake on RTK as you normally do
  2. Check the option BUILD_SIMPLERTK (default is OFF)
  3. Ename WRAP_PYTHON and/or WRAP_CSHARP (On Windows)
  4. Configure
  5. Build as usual. Note that building will download the require libraries so you should have an internet connection enabled.

Python installation

To install a built python package into the system Python, as root run:

cd SimpleRTK-build/Wrapping
python PythonPackage/setup.py install

Alternatively, a Python virtual environment can be created and the distribution installed there. If you build the "dist" target a Python egg will be created in the "Wrapping/dist" directory. Building Python wheels can be enabled with a CMake flag.

Testing SimpleRTK

A simple example is located in Utilities/SimpleRTK/Examples/RTKFirstReconstruction.py and shows how to use SimpleRTK

#!/usr/bin/env python
from __future__ import print_function
import SimpleRTK as srtk
import sys
import os
import matplotlib.pyplot as plt
import matplotlib.cm as cm

if len ( sys.argv ) < 2:
    print( "Usage: RTKFirstReconstruction <output>" )
    sys.exit ( 1 )

# Defines the RTK geometry object
geometry = srtk.ThreeDimCircularProjectionGeometry()
numberOfProjections = 360
firstAngle = 0
angularArc = 360
sid = 600 # source to isocenter distance in mm
sdd = 1200 # source to detector distance in mm
isox = 0 # X coordinate on the projection image of isocenter
isoy = 0 # Y coordinate on the projection image of isocenter
for x in range(0,numberOfProjections):
  angle = firstAngle + x * angularArc / numberOfProjections
  geometry.AddProjection(sid,sdd,angle,isox,isoy)

constantImageSource = srtk.ConstantImageSource()
origin = [ -127.5, -127.5, 0. ]
sizeOutput = [ 256, 256,  numberOfProjections ]
spacing = [ 1.0, 1.0, 1.0 ]
constantImageSource.SetOrigin( origin )
constantImageSource.SetSpacing( spacing )
constantImageSource.SetSize( sizeOutput )
constantImageSource.SetConstant(0.0)
source = constantImageSource.Execute()

rei = srtk.RayEllipsoidIntersectionImageFilter()
semiprincipalaxis = [ 50, 50, 50]
center = [ 0, 0, 0]
# Set GrayScale value, axes, center...
rei.SetDensity(20)
rei.SetAngle(0)
rei.SetCenter(center)
rei.SetAxis(semiprincipalaxis)
rei.SetGeometry( geometry )
reiImage = rei.Execute(source)
  
  
# Create reconstructed image
constantImageSource2 = srtk.ConstantImageSource()
origin = [ -63.5, -63.5, -63.5 ]
sizeOutput = [ 128, 128, 128 ]
constantImageSource2.SetOrigin( origin )
constantImageSource2.SetSpacing( spacing )
constantImageSource2.SetSize( sizeOutput )
constantImageSource2.SetConstant(0.0)
source2 = constantImageSource2.Execute()

print("Performing reconstruction")
feldkamp = srtk.FDKConeBeamReconstructionFilter()
feldkamp.SetGeometry( geometry );
feldkamp.SetTruncationCorrection(0.0);
feldkamp.SetHannCutFrequency(0.0);
image = feldkamp.Execute(source2,reiImage) 
   
pixelID = image.GetPixelIDValue()

plt.imshow(srtk.GetArrayFromImage(image[:,64,:]), cmap = cm.Greys_r)

caster = srtk.CastImageFilter()
caster.SetOutputPixelType( pixelID )
image = caster.Execute( image )

writer = srtk.ImageFileWriter()
writer.SetFileName ( sys.argv[1] )
writer.Execute ( image );

Extending SimpleRTK

Adding new classes to SimpleRTK should be fairly straightforward. In this section we cover how to add common types as well as filters. Note that common types should only be added when necessary

Common Type

One example is for the ThreeDCircularProjectionGeometry object. The code is located in utilities/SimpleRTK/Code/Common. Note that the class is named srtkThreeDimCircularProjectionGeometry to avoid some conflicting type resolution by Swig (this could be tricky). This class in the namespace rtk::simple provides a PIMPL class implementation of the object as well as a conversion from SimpleRTK to RTK itself to expose the necessary functions.

Filters

Filters are implemented using a json description for simplicity. The best way is to look into already wrapped filters in the: utilities/SimpleRTK/Code/BasicFilters/json directory.