[Rtk-users] Troubleshooting Conjugate Gradient Filter in Python

[Simon Rit]

Thanks for getting back on this. I really need to update the pypi packages
but I'm fighting a bit with the CI, it's no longer working (see #395
<https://github.com/SimonRit/RTK/pull/395>). I'll keep everyone posted on
the mailing list of course.

On Fri, Dec 18, 2020 at 12:17 PM Clément Mirabel <clement.mirabel at gmail.com>
wrote:

> 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/62cdd806/attachment.htm>