[Rtk-users] RTK first reconstruction

[MG Vallet]

> Dear RTK users
>
> I need to perform a simple cone beam reconstruction.
> Starting from the SimpleRTK example, I split it into two python scripts:
>  - a test case generator that produces data files
>  - and a reconstruction script reading these files
> Both scripts are attached.
>
> I tried several file format, including tif, for the data files. But I
> always get a segmentation fault in the feldkamp execution step.
>
> Is it my error (probably) or does the write/read process drop some data ?
>
> MG Vallet
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/rtk-users/attachments/20151214/efd7c833/attachment-0009.html>
-------------- next part --------------
__author__ = 'ORS'
# !/usr/bin/env python
import Wrapping.SimpleRTK as srtk
import sys
import os
import matplotlib.pyplot as plt
import matplotlib.cm as cm

srtk.Image(10,10,10,srtk.srtkFloat32)
# Defines the RTK geometry object
geometry = srtk.ThreeDCircularProjectionGeometry()  # Defines the projections
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)

writer = srtk.ThreeDCircularProjectionGeometryXMLFileWriter()
print("Writing file d:/ellipse_geometry.xml")
writer.SetFileName('d:/ellipse_geometry.xml')
writer.Execute(geometry)

constantImageSource = srtk.ConstantImageSource()  # Initializes empty 2D images
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()  # Compute the projections
# Defines the object
semiprincipalaxis = [50, 30, 20]
center = [-10, 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)

# Output the initial projections
writer = srtk.ImageFileWriter()
print("Writing file d:/ellipse_projection.mhd")
writer.SetFileName('d:/ellipse_projection.mhd')
writer.Execute(reiImage)
-------------- next part --------------
__author__ = 'ORS'
#!/usr/bin/env python
import Wrapping.SimpleRTK as srtk
import sys
import os
import matplotlib.pyplot as plt
import matplotlib.cm as cm


# Input the RTK geometry object
reader = srtk.ThreeDCircularProjectionGeometryXMLFileReader()
print("Reading file d:/ellipse_geometry.xml")
reader.SetFileName ('d:/ellipse_geometry.xml' )
geometry = reader.Execute()

# Input the projections
reader = srtk.ImageFileReader()
print("Reading file d:/ellipse_projection.mhd")
reader.SetFileName ('d:/ellipse_projection.mhd' )
reiImage = reader.Execute()
#the origin is not set correctly
reiImage.SetOrigin([-127.5, -127.5, 0])

# Create (an empty) reconstructed image
constantImageRecons = srtk.ConstantImageSource()
origin = [ -63.5, -63.5, -63.5 ]
spacing = [ 1.0, 1.0, 1.0 ]
sizeOutput = [ 128, 128, 128 ]
constantImageRecons.SetOrigin( origin )
constantImageRecons.SetSpacing( spacing )
constantImageRecons.SetSize( sizeOutput )
constantImageRecons.SetConstant(0.0)
recons = constantImageRecons.Execute()

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

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

# Output the 3D image
writer = srtk.ImageFileWriter()
writer.SetFileName ('d:/testrtkwithdata.tif' )
writer.Execute ( image )