[Rtk-users] Error in the FFTRampImpageFilter

[Simon Rit]

Thanks. This was a multi-threading problem which is now fixed by processing
each projection at a time within a thread:
https://github.com/SimonRit/RTK/commit/a719d8913beb7bé6a39da787980ead720b353f545
<https://github.com/SimonRit/RTK/commit/a719d8913beb7b6a39da787980ead720b353f545>
I don't really understand the number you gave but they should depend on the
number of threads used by ITK. Normally, ITK splits the last dimension
(here the number of projections) in almost equal pieces for each thread. So
the last dimension of the processed region will depend on this...

On Sat, May 9, 2020 at 11:08 AM Jerome Lesaint <lesaint.jerome at gmail.com>
wrote:

> See the code below.
>
> I gave you a wrong information. Sorry. My stack has 720 projections and
> not 360. It actually works fine with 360. But fails with 720 !!
> After a few trials, the behavior is really strange. It works with 13, 17,
> 23, 300, 301, 360, 361, 480. But it fails with 400 (13), 401 (13), 540
> (17), 720 (23), 721 (23). In parentheses, the number which triggers the
> error (no 2s, 3s 5s in the prime decomp).
>
> Jérome
>
>
> -------------------------------------------------------------------
> sid = 500
> sdd = 1000
>
> spacing = 1
> nb_pixels = int(257/spacing)
> nb_pixels_u = nb_pixels
> nb_pixels_v = nb_pixels_u
>
> geo_0 = rtk.ThreeDCircularProjectionGeometry.New()
> for i in range(nprojs_0) :
>     geo_0.AddProjection(sid,sdd,i)
>
> ImageType = itk.Image[itk.F, 3]
> const = rtk.ConstantImageSource[ImageType].New()
> const.SetConstant(0.)
> const.SetSpacing([spacing,spacing,spacing])
> const.SetSize([nb_pixels_u,nb_pixels_v,len(geo_0.GetGantryAngles())])
> const.SetOrigin(-0.5*spacing*(np.array(const.GetSize())-1))
> sl = rtk.SheppLoganPhantomFilter[ImageType,ImageType].New()
> sl.SetInput(const.GetOutput())
> sl.SetGeometry(geo_0)
> sl.SetPhantomScale(60)
>
> sl.Update()
>
> #weights = rtk.FDKWeightProjectionFilter[ImageType].New()
> #weights.SetGeometry(geo_0)
> #weights.SetInput(sl.GetOutput())
> #weights.SetNumberOfWorkUnits(1)
> #weights.InPlaceOff()
>
> ramp = rtk.FFTRampImageFilter[ImageType,ImageType,itk.D].New()
> ramp.SetInput(sl.GetOutput())
> ramp.Update()
>
>
> Le ven. 8 mai 2020 à 22:15, Simon Rit <simon.rit at creatis.insa-lyon.fr> a
> écrit :
>
>> Hi,
>> I cannot reproduce the issue with this code
>>
>> import itk
>> from itk import RTK as rtk
>>
>> ImageType = itk.Image[itk.F,3]
>>
>> constantImageSource = rtk.ConstantImageSource[ImageType].New()
>> constantImageSource.SetOrigin( [ -127.875, -127.875, 0. ] )
>> constantImageSource.SetSpacing( [ .25 ] *3 )
>> constantImageSource.SetSize( [ 1024, 512, 360 ] )
>> constantImageSource.SetConstant(1.)
>>
>> ramp = rtk.FFTRampImageFilter[ImageType, ImageType, itk.D].New()
>> ramp.SetInput(constantImageSource.GetOutput())
>> ramp.Update()
>>
>> Can you tell me more on the pipeline you're using?
>> Simon
>>
>> On Fri, May 8, 2020 at 9:36 AM Jerome Lesaint <lesaint.jerome at gmail.com>
>> wrote:
>>
>>> Salut,
>>>
>>> thanks again for the answer. The streaming workaround seems to work.
>>> As for the second solution you mentioned : actually there are 360
>>> projections in my stack (not 23). 360 has only 2s, 3s and 5s in its prime
>>> decomp. Of course, 23 does not. But where this 23 comes form ???
>>>
>>> Anyway, it works fine now. Thank you.
>>> Jérome
>>>
>>>
>>> Le jeu. 7 mai 2020 à 20:48, Simon Rit <simon.rit at creatis.insa-lyon.fr>
>>> a écrit :
>>>
>>>> Hi,
>>>> Nope, sorry. We can't distribute the package with FFTW because its
>>>> license does not allow it. What you can do is put an
>>>> itkStreamingImageFilter after the fft ramp to process one projection at a
>>>> time. Or add one empty projection after you 23d.
>>>> Let me know if we can help with one of the two options.
>>>> Simon
>>>>
>>>> On Thu, May 7, 2020 at 7:24 PM Jerome Lesaint <lesaint.jerome at gmail.com>
>>>> wrote:
>>>>
>>>>> Hi,
>>>>>
>>>>> I am trying to use the ramp filter FFTRampImageFilter and it fails at
>>>>> Update time with the error message :
>>>>>
>>>>> itk::ERROR:
>>>>> VnlRealToHalfHermitianForwardFFTImageFilter(00000230EC986BA0): Cannot
>>>>> compute FFT of image with size [1024, 512, 23].
>>>>> VnlRealToHalfHermitianForwardFFTImageFilter operates only on images whose
>>>>> size in each dimension has a prime factorization consisting of only 2s, 3s,
>>>>> or 5s.
>>>>>
>>>>> I already read Simon's answer to a previous similar issue, saying that
>>>>> he recommends using FFTW by turning ITK_USE_FFTW when building ITK.
>>>>> But I did not compile ITK (on this computer). I simply installed the
>>>>> pre-compiled Python package.
>>>>> My question : can I turn any parameter in Python to force the ramp
>>>>> filter to use this library ?
>>>>>
>>>>> Thanks in advance for your help,
>>>>>
>>>>> Jérome
>>>>> _______________________________________________
>>>>> 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/20200513/ea89977f/attachment.html>