[Rtk-users] Troubleshooting Conjugate Gradient Filter in Python

[Clément Mirabel]

Hi,

I created another topic since the content is different from the previous
one and people could land here with a different search in the archives.

While trying to reconstruct a set of images from 0 and 180° without having
to add the fan angle, Simon recommended to use the conjugate gradient
algorithm in RTK. So, I first tried to replicate in Python the conjugate
gradient application written in C++ (
https://github.com/SimonRit/RTK/tree/master/applications/rtkconjugategradient).
But whatever types of input I set, my python kernel dies when updating the
filter without any error or debug message. I spent some time thinking about
it but I can't find out what is wrong.
I have attached my code at the end of this message if you have suggestions
about what would be wrong.

Thanks for your help,

Clément

------------------------------------------------------------------------------------------------------------

# Defines the image type
ImageType = itk.Image[itk.F,3]

# Defines the RTK geometry object
geometry = rtk.ThreeDCircularProjectionGeometry.New()
numberOfProjections = 180
firstAngle = 90
sid = 600 # source to isocenter distance
sdd = 1200 # source to detector distance
angularArc = 180
for x in range(0,numberOfProjections):
  angle = firstAngle + x * angularArc / numberOfProjections
  geometry.AddProjection(sid,sdd,angle)

# Create a stack of empty projection images
ConstantImageSourceType = rtk.ConstantImageSource[ImageType]
constantImageSource = ConstantImageSourceType.New()
origin = [ -199, -199, 0. ]
sizeOutput = [ 200, 200, 180 ]
spacing = [ 2.0, 2.0, 2.0 ]
constantImageSource.SetOrigin( origin )
constantImageSource.SetSpacing( spacing )
constantImageSource.SetSize( sizeOutput )
constantImageSource.SetConstant(0.)

# Forward projection of input volume
proj = rtk.JosephForwardProjectionImageFilter[ImageType, ImageType].New()
proj.SetInput(constantImageSource)
proj.SetInput(1, input) # input is ImageType from a reader
proj.SetGeometry(geometry)

# Create weights
weightsSource = ConstantImageSourceType.New()
weightsSource.SetInformationFromImage(proj.GetOutput())
weightsSource.SetConstant(1.0)

# Create output like volume
outputImageSource = ConstantImageSourceType.New()
sizeOutput = [ 400, 400, 400]
origin = [ -200, -200, -200 ]
spacing = [ 1.0, 1.0, 1.0 ]
outputImageSource.SetOrigin(origin)
outputImageSource.SetSpacing(spacing)
outputImageSource.SetSize(sizeOutput)
outputImageSource.SetConstant(0.)

# Create reconstructed image
RCGType = rtk.ConjugateGradientConeBeamReconstructionFilter[ImageType]
conjugateGradient = RCGType.New()
conjugateGradient.SetInput(outputImageSource.GetOutput())
conjugateGradient.SetInput(1, proj.GetOutput())
conjugateGradient.SetInput(2, weightsSource.GetOutput())
conjugateGradient.SetGeometry(geometry)
conjugateGradient.SetGamma(2.0)
conjugateGradient.SetNumberOfIterations(30)
conjugateGradient.SetDisableDisplacedDetectorFilter(True)
conjugateGradient.Update()
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://public.kitware.com/pipermail/rtk-users/attachments/20201216/b9c2c8b2/attachment.htm>