RTKConsortium / RTK

Reconstruction Toolkit
Apache License 2.0
240 stars 143 forks source link

Python FDK does not filter data #483

Closed lesaintjerome closed 2 years ago

lesaintjerome commented 2 years ago

Using RTK 2.3.0, ITK 5.3.0rc3 wrappings with python 3.9 on MacOs. FDK reconstruction seems to not filter the data prior to back-projection. This is the case both in Spyder or by direct execution of the script in command line (python testfdk.py). Anaconda python. Pb occurs when script launched within an env. or not.

Steps to reproduce :

import itk as itk
from itk import RTK as rtk
import matplotlib.pyplot as plt

geo = rtk.ThreeDCircularProjectionGeometry.New()
for i in range(360) : 
    geo.AddProjection(1000,1500,i)

imtype = itk.Image[itk.F,3]

cst = rtk.ConstantImageSource[imtype].New(Constant=0.,Size=[256,3,360],Origin=[-127.5,-1.,0.])
proj = rtk.SheppLoganPhantomFilter[imtype,imtype].New()
proj.SetGeometry(geo)
proj.SetInput(cst)
proj.SetPhantomScale(80)

csttomo = rtk.ConstantImageSource[imtype].New(Constant=0.,Size=[256,1,256],Spacing=1.,Origin=[-127.5,0,-127.5])

fdk = rtk.FDKConeBeamReconstructionFilter[imtype].New()
fdk.SetGeometry(geo)
fdk.SetInput(0,csttomo.GetOutput())
fdk.SetInput(1,proj.GetOutput())

fdk.Update()
recons = fdk.GetOutput()

itk.imwrite(recons,'recons.mha')
SimonRit commented 2 years ago

I don't seem to be able to reproduce on Linux. Which RTK package did you use, the one on pypi?

lesaintjerome commented 2 years ago

Nope. I downloaded the Wheels from last action on master here.

SimonRit commented 2 years ago

I can reproduce it on linux and it's an FFT issue

import itk as itk
from itk import RTK as rtk

geo = rtk.ThreeDCircularProjectionGeometry.New()
for i in range(360) : 
    geo.AddProjection(1000,1500,i)
imtype = itk.Image[itk.F,3]
cst = rtk.ConstantImageSource[imtype].New(Constant=0.,Size=[256,3,360],Origin=[-127.5,-1.,0.])
proj = rtk.SheppLoganPhantomFilter[imtype,imtype].New()
proj.SetGeometry(geo)
proj.SetInput(cst)
proj.SetPhantomScale(80)

ramp = rtk.FFTRampImageFilter[imtype, imtype, itk.D].New()
ramp.SetInput(0,proj.GetOutput())

ramp.UpdateLargestPossibleRegion()
recons = ramp.GetOutput()

itk.imwrite(recons,'ramp.mha')

I'm guessing that compiling new RTK packages with new ITK code but using an "old" ITK release causes a problem due to, e.g., InsightSoftwareConsortium/ITK#3156. I'll ask if we can get new ITK wheels...

SimonRit commented 2 years ago

Locally recompiling the wrappings from ITK and RTK fixes the issue. If I add after the imports the line

itk.VnlFFTImageFilterInitFactory.RegisterFactories()

it seems to fix the issue... So there is nothing else to fix the issue than compiling locally, using this workaround or wait for a new ITK release.

SimonRit commented 2 years ago

New packages provided in InsightSoftwareConsortium/ITK#3156 show that the issue is not solved.

SimonRit commented 2 years ago
itk.VnlFFTImageFilterInitFactory.RegisterFactories()

does not help.