[Rtk-users] Troubleshooting Conjugate Gradient Filter in Python

[Clément Mirabel]

Hi,

Thank you for your answer. I still had the issue with your updated code but
I noticed I was still using itk 5.0.1 because of an issue
<https://tinyurl.com/y3dm3h86> with itk 5.1.2 due to numpy 1.9.4 with
Windows 10 recent update. I downgraded numpy to 1.9.3 so itk could install
5.1.2 and then installed the artifact you pointed for RTK.
Everything works fine now, time to play with the parameters!

Clément



Le ven. 18 déc. 2020 à 08:00, Simon Rit <simon.rit at creatis.insa-lyon.fr> a
écrit :

> 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/39ba7ec3/attachment-0001.htm>