problem with check_array in opencl sinofilter.py

[jcesardasilva]

I have a complex filter to be used in the Filtered Back Projection (the Hilbert filter), but with the "check_array" of sinofilter.py of the new Silx version, it only allows me to use numpy.float32 filter. Could you please either remove this checking or take complex filter into account? Something like:

if arr.dtype != np.float32 or arr.dtype != np.complex64:
    raise ValueError("Expected data type = numpy.float32 or numpy.complex64")

[pierrepaleo]

Hi @jcesardasilva The silx SinoFilter has changed its internal Fourier computation to Real-to-Complex for performances reasons. Ultimately, the sinogram is "convolved" by a real filter, so the filter we pass to Backprojection.sino_filter is always the FFT of some real filter.

In general, SinoFilter.set_filter() takes a complex array as an input ; this input is the the R2C FFT of some real filter. So the input of set_filter() is already in the Fourier domain.
The Hilbert filter 1/(2*j*pi) * sign(f) is a "valid" filter : since it is purely imaginary, it is the Fourier Transform of some real (and odd) function.

import numpy as np
from silx.test.utils import utilstest
from silx.opencl.backprojection import Backprojection
sino = np.load(utilstest.getfile("sino500.npz"))["data"]
sino_diff = np.diff(sino, axis=1) # simulate differential sinogram
B = Backprojection(sino_diff.shape)

Nf = B.sino_filter.sino_padded_shape[1]
freqs = np.fft.fftfreq(Nf)
hilbert_filter = 1./(2*pi*1j)*np.sign(freqs)
hilbert_filter_r2c = hilbert_filter[:Nf//2+1] # Hermitian symmetry
B.sino_filter.set_filter(hilbert_filter_r2c)
rec = B(sino_diff)

[jcesardasilva]

Hi @pierrepaleo Thank you very much for the reply and the sorry about the delay to come back to the discussion. Well, I understand the point, but I don't know why I still get the problem, but in a different way. This is when I don't use numpy backend. I do the backprojection once, then I do the projection to get again the sinogram to compare with my experimental data and, finally, when I want to do the reconstruction again, I get the warning:

/mntdirect/_data_id16a_inhouse4/visitor/ma4351/id16a/v97_v_nfptomo/v97_v_nfptomo_analysis/NearPtychoTomo/toupy/tomo/iradon.py in mod_iradonSilx(radon_image, theta, output_size, filter_type, derivatives, interpolation, circle, freqcutoff, use_numpy)
    318 
    319     # actual reconstruction
--> 320     recons = B(sinogram)
    321 
    322     if circle:

/mntdirect/_data_id16a_inhouse1/sware/pyvenvid16a/py3id16a/lib/python3.6/site-packages/silx/opencl/backprojection.py in filtered_backprojection(self, sino, output)
    378         """
    379         # Filter
--> 380         self.sino_filter(sino, output=self.d_sino)
    381         # Backproject
    382         res = self.backprojection(self.d_sino, output=output)

/mntdirect/_data_id16a_inhouse1/sware/pyvenvid16a/py3id16a/lib/python3.6/site-packages/silx/opencl/sinofilter.py in filter_sino(self, sino, output)
    368         """
    369         # Handle input sinogram
--> 370         self._prepare_input_sino(sino)
    371         # FFT
    372         self._do_fft()

/mntdirect/_data_id16a_inhouse1/sware/pyvenvid16a/py3id16a/lib/python3.6/site-packages/silx/opencl/sinofilter.py in _prepare_input_sino(self, sino)
    253         :param sino: sinogram
    254         """
--> 255         self.check_array(sino)
    256         self.d_sino_padded.fill(0)
    257         if self.fft_backend == "opencl":

/mntdirect/_data_id16a_inhouse1/sware/pyvenvid16a/py3id16a/lib/python3.6/site-packages/silx/opencl/sinofilter.py in check_array(self, arr)
    205     def check_array(self, arr):
    206         if arr.dtype != np.float32:
--> 207             raise ValueError("Expected data type = numpy.float32")
    208         if arr.shape != self.sino_shape:
    209             raise ValueError("Expected sinogram shape %s, got %s" %

ValueError: Expected data type = numpy.float32

[jcesardasilva]

Hi again, everything goes very smooth with Silx version 0.9.0 and not much difference in computing time compared to Silx version 0.11.0. At least from my side. I am wondering if by using numpy backend in Silx 0.11.0 as you suggested, I get exactly the same backprojector and projector as the one in version Silx 0.9.0 or not. I admit that I am tempted to keep Silx version 0.9.0 for the tomo reconstruction, but not for everything else. Could you keep the backprojector and projector from version 0.9.0 somehow inside the new versions, please? Or any other ideas are welcome.

[pierrepaleo]

It looks like this error is not related anymore to the filter side. I suspect that you are doing some iterative processing, and the sinogram changes its data type at some point. Schematically:

image = B(sino) # reconstruction. sino and image are float32
sino = P(image) # projection. sino is float32
sino = some_processing(sino) # float32 -> float64
image = B(sino) # error: sino is float64 here
# ...

The float64 dtype error was not raised in 0.9, as the sinogram data was copied anyway in another array. 0.10 tries to be more efficient, but is less easy to use.

Can you try casting the sinogram data to float32 before doing the backprojection ?

[jcesardasilva]

Hi @pierrepaleo, thank you for your answer. I indeed need the projector and backprojector in an iterative alignment method. I modify the sinogram in several parts of my code. Thus, what I decided to do is to initialize the projector and the backprojector at each time. The efficiency is still high, so, the time spent doing this is not important. At every call of the Projector as P and Backprojector as B, I initialize them with None and inside the function/method I call global P and global B, respectively.