Backprojection: GpyFFT_Error: 'Invalid plan.

[jcesardasilva]

I saw the PR #2373 which is closed. So, I decided to open an issue.

I am having problems with the gpyfft for filtering projections on GPU on the version 0.10.0. This does not happen in version 0.9.0. When I try to use backprojections, I get the error: GpyFFT_Error: 'Invalid plan.' See the traceback below:

GpyFFT_Error                              Traceback (most recent call last)
/mntdirect/_data_id16a_inhouse5/visitor/ma3495/id16a/O2_LSCF_CGO/O2_LSCF_CGO_25nm_analysis/analysis_pynx3/NearPtychoTomoID16A/iradon.py in mod_iradon2(radon_image, theta, output_size, filter_type, derivative, interpolation, circle, freqcutoff)
--> 232     B = Backprojection(radon_image.T.shape, angles=np.pi*(theta)/180.)
    233     print("Initialized OpenCL backprojector on %s" % B.device)
    234     cust_filter = compute_filter(radon_image.shape[0], filter_type=filter_type,derivative=derivative, freqcutoff=freqcutoff)

/mntdirect/_data_id16a_inhouse1/sware/pyvenvid16a/py3id16a/lib/python3.6/site-packages/silx/opencl/backprojection.py in __init__(self, sino_shape, slice_shape, axis_position, angles, filter_name, ctx, devicetype, platformid, deviceid, profile, extra_options)
    111         self._compute_angles()
    112         self._init_kernels()
--> 113         self._init_filter(filter_name)
    114 
    115     def _init_geometry(self, sino_shape, slice_shape, angles, axis_position,

/mntdirect/_data_id16a_inhouse1/sware/pyvenvid16a/py3id16a/lib/python3.6/site-packages/silx/opencl/backprojection.py in _init_filter(self, filter_name)
    259             ctx=self.ctx,
    260             filter_name=self.filter_name,
--> 261             extra_options=self.extra_options,
    262         )
    263 

/mntdirect/_data_id16a_inhouse1/sware/pyvenvid16a/py3id16a/lib/python3.6/site-packages/silx/opencl/sinofilter.py in __init__(self, sino_shape, filter_name, ctx, devicetype, platformid, deviceid, profile, extra_options)
    173 
    174         self._calculate_shapes(sino_shape)
--> 175         self._init_fft()
    176         self._allocate_memory()
    177         self._compute_filter(filter_name, extra_options)

/mntdirect/_data_id16a_inhouse1/sware/pyvenvid16a/py3id16a/lib/python3.6/site-packages/silx/opencl/sinofilter.py in _init_fft(self)
    218                 axes=(-1,),
    219                 backend="opencl",
--> 220                 ctx=self.ctx,
    221             )
    222         else:

/mntdirect/_data_id16a_inhouse1/sware/pyvenvid16a/py3id16a/lib/python3.6/site-packages/silx/math/fft/fft.py in FFT(shape, dtype, template, shape_out, axes, normalize, backend, **kwargs)
     92         axes=axes,
     93         normalize=normalize,
---> 94         **kwargs
     95     )
     96     return F

/mntdirect/_data_id16a_inhouse1/sware/pyvenvid16a/py3id16a/lib/python3.6/site-packages/silx/math/fft/clfft.py in __init__(self, shape, dtype, template, shape_out, axes, normalize, ctx, fast_math, choose_best_device)
     92         self.allocate_arrays()
     93         self.real_transform = np.isrealobj(self.data_in)
---> 94         self.compute_forward_plan()
     95         self.compute_inverse_plan()
     96         self.refs = {

/mntdirect/_data_id16a_inhouse1/sware/pyvenvid16a/py3id16a/lib/python3.6/site-packages/silx/math/fft/clfft.py in compute_forward_plan(self)
    186             axes=self.axes,
    187             fast_math=self.fast_math,
--> 188             real=self.real_transform,
    189         )
    190 

/mntdirect/_data_id16a_inhouse1/sware/pyvenvid16a/py3id16a/lib/python3.6/site-packages/gpyfft/fft.py in __init__(self, context, queue, in_array, out_array, axes, fast_math, real, callbacks)
     92 
     93         plan = GFFT.create_plan(context, t_shape)
---> 94         plan.inplace = t_inplace
     95         plan.strides_in = t_strides_in
     96         plan.strides_out = t_strides_out

gpyfft/gpyfftlib.pyx in gpyfft.gpyfftlib.Plan.inplace.__set__()

gpyfft/gpyfftlib.pyx in gpyfft.gpyfftlib.errcheck()

GpyFFT_Error: 'Invalid plan.'

[pierrepaleo]

Hi @jcesardasilva thanks for the report.

The error is due to the size of the transform. The underlying GPU FFT backend (clFFT) only supports sizes which are powers of primes [2, 3, 5, 7]. The data in "O2_LSCF_CGO" has 2392 pixels columns, and 2392 = 2**3 * 13 * 23. In the SinoFilter implementation, the sinogram width is padded to 2*width, which leaves "big primes" in the integer factorization of width.

The error appears in 0.10.0 because clFFT became the default backend for doing the sinogram filtering. Other backends (ex numpy) are less regarding on the transform size.

The PR #2490 attempts at fixing this. You can also directly use built-in filters:

B = Backprojection(radon_image.T.shape, angles=np.deg2rad(theta), filter_name="hann")

It is also possible to force the usage of the numpy backend:

B = Backprojection(radon_image.T.shape, angles=np.deg2rad(theta), filter_name="hann", extra_options={"use_numpy_fft": True})