[Rtk-users] Troubleshooting Conjugate Gradient Filter in Python

[Simon Rit]

Hi,
I slightly modified your script and it worked for me. What version are you
using? I used the latest ITK v5.1.2 with, for RTK, an artifact published
here <https://github.com/SimonRit/RTK/actions/runs/344449575> (visible only
if you're logged in github).
Simon

PS: I don't think you need the displaced detector, it is centered on the
central ray in your simulation.

import itk
from itk import RTK as rtk

# 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 = [ 100, 100, 180 ]
spacing = [ 4.0, 4.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, itk.imread('test.mha')) # input is ImageType from a reader
proj.SetGeometry(geometry)
proj.Update()

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

# Create output like volume
outputImageSource = ConstantImageSourceType.New()
sizeOutput = [ 100, 100, 100]
origin = [ -200, -200, -200 ]
spacing = [ 4.0, 4.0, 4.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(2)
conjugateGradient.SetDisableDisplacedDetectorFilter(False)
conjugateGradient.Update()

itk.imwrite(conjugateGradient.GetOutput(), 'output.mha')

On Wed, Dec 16, 2020 at 10:36 AM Clément Mirabel <clement.mirabel at gmail.com>
wrote:

> 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()
> _______________________________________________
> Rtk-users mailing list
> Rtk-users at public.kitware.com
> https://public.kitware.com/mailman/listinfo/rtk-users
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://public.kitware.com/pipermail/rtk-users/attachments/20201218/d108cbb3/attachment.htm>