Clij2 fft filter?

acrevenna · October 30, 2020, 11:42am

Dear @haesleinhuepf, is there a fft filter in clij2? I thought so, since there is this page:

but then if I type:
Ext.CLIJ2_
on my macro, there is no FFT. How could I call the FFT from clij2?
thanks and congrats again on the new position!
Alvaro

haesleinhuepf · October 30, 2020, 12:51pm

Hi Alvaro @acrevenna,

sorry to tell you that there is no FFT function yet in CLIJ2. As you can read on the repository you linked, it contains FFT-based alogrithms, e.g. as convolution and deconvolution.

May I ask what you would like to do with an FFT-function? Maybe we can guide you in how to get it done on the GPU.

Cheers,
Robert

acrevenna · October 30, 2020, 1:28pm

Dear Robert @haesleinhuepf, thanks for the fast reply. We want to do motion correction of data from a miniscope.
We were using a FFT as a bandpass filter:
run(“Bandpass Filter…”, “filter_large=10000 filter_small=100 suppress=None tolerance=5 process”);
fft=getTitle();
and once having that, would do: original - fft and then that result would go to turboreg. I have your nice example of motion_corr in clij2 but I would like to use the FFT and not an intensity threshold. What would you the clij-ean way now?
thanks,
Alvaro

acrevenna · October 30, 2020, 1:42pm

I’ve also discovered now the clijx correct drift commands, will give those a try.

acrevenna · October 30, 2020, 2:01pm

I guess the fft was just a way to get the background… not to calculate drift… sorry to bug you!

acrevenna · October 30, 2020, 4:05pm

Ok, after trying a few things a bit more…
Dear @haesleinhuepf ,
I’d like to use the fft to get the cross-correlation between two images and calculate a spatial shift between them, which then I pass to the affinetransform2D. Is there a way to do that in clij2?
thanks

haesleinhuepf · October 30, 2020, 6:29pm

Hi @acrevenna,

you can compute the cross correlation between two images. If you know that you are only looking for a translation (not a full affine transform), you can also try the translation registration.

Furthermore, I put the FFT on the wish list for the future CLIJ. It might not be too complicated though. Would you be interested in helping? If so, how are your Java coding skills?

Cheers,
Robert

bnorthan · November 1, 2020, 1:37pm

There is one flavor of FFT currently exposed here.

That is a 2d, 32 float type FFT that uses ‘long pointers’ to opencl memory as the input. This should be callable from a CLIJ program. Problem is we currently haven’t exposed the inverse FFT, so it probably isn’t that useful yet.

The above is a wrapper to clfft. One option to get FFT into CLIJ is to simply wrap more functionality from clfft. This would require some skills in both java and c.

schmid · November 2, 2020, 10:38am

Hi everyone,
if you do not use the ‘suppress stripes’ of the bandpass filter, in principle you can replace it by Gaussian blurs. You have to multiply the numbers by 0.45, so a ’ filter_small=100’ should become a Gaussian Blur with a sigma value of 45.
For the ‘filter large’ you have to do a high-pass filter, i.e. subtract the result of the Gaussian Blur from the original filter.
The only difference between the Gaussian filter and the FFT-based bandpass is the handling of near-edge pixels: Smoothing needs an assumption for what the out-of image pixels would be. The ImageJ Gaussian blur uses replication of the edge pixels. When using large radii, this gives improperly high weight to the pixels at the very edge of the image, and even more so for the corner pixels. The FFT method mirrors the image at the borders (at least for some distance). For objects close to or touching the edge, this can lead to smearing out towards the edge. The image below shows these differences (red arrows; the blue contour line is at equal height for both). For Gaussian blurs, smoothing of near-edge pixels better suppresses the noise of the edge pixels when performing two (or more) blur operations after each other, starting with smaller sigma, and such that the sum of the squared sigma values is the square of the desired sigma.

haesleinhuepf · November 2, 2020, 11:15am

Hey Michael @schmid,

that’s so cool. I knew about that relationship but I can’t wrap my mind around the 0.45. Can you tell why it’s this factor?

Thanks!

Cheers,
Robert

schmid · November 2, 2020, 12:43pm

Hi Robert,

good question! I have determined the factor of 0.45 experimentally, by blurring a line and fitting a Gaussian. I did not follow all the code, but I think that the exact value would be √2/π = 0.450158, with the factor of √2 because the Gaussian Blur is defined in terms of standard deviation, which is 1/√2 of the length of decay to 1/e, and π from the Fourier transform.
In any case, the math behind blurring via multiplying with a Gaussian in Fourier space is

Fourier transform of (exp(-ax²)) = √(π/a) × exp(-π²k²/a) ,

where k = 2π/wavelength and c is a constant.
The convolution theorem of the Fourier transform tells that instead of convolution with a Gaussian one may multiply the Fourier components of the image with the Fourier transform of the Gaussian, and then do the back transformation.

–Michael

bnorthan · November 2, 2020, 4:32pm

Hi @haesleinhuepf and @acrevenna

I exposed the c++ float 2d inverse FFT in java and added an example showing how to use it here.

The example takes an FFT of the image applies a filter and then performs inverse FFT,. The example is messy but potentially could be used as a template for more sophisticated operations.

A couple of caveats to consider

I don’t think CLIJ has a formal way to handle complex numbers (correct me if I am wrong @haesleinhuepf). Keep in mind even though the FFT is put in a “float” buffer it actually represents a complex float with real and imaginary parts interleaved.
In the current code the FFT “plans” are created every time an FFT is called. In a time critical application you would want to refactor plan and creation out of the FFT call to improve performance.

Brian

bnorthan · February 13, 2021, 7:35pm

This has come up on Twitter again so it a good time to revisit it

The output of CLFFT is a contiguous gpu array of complex interleaved format so it’s just a matter of designing CLIJ kernels that can handle complex interleaved operations. Thus far I’ve written the complex kernels as strings for example here. I am sure there is a nicer way to do it long term. A good starting point as @axtimwalde mentioned was to study the imglib2 api then figure out how to implement that with GPU buffers.

The second point isn’t trivial. In clfft every different dimensionality, size, data type, and direction (forward or inverse) of FFT needs a plan. For example this code initiates an FFT plan.

So as I mentioned above in time critical applications, where the same size FFT is used over and over again you may want to expose plan creation at a high level, so the user can keep the plan, and re-use it. Or maybe you wouldn’t want to do that, at least not right away. You may want to just expose FFTs in a simple way, then advanced users could dive into the clfft c-api. There’d be some busy work either way.

Brian

haesleinhuepf · February 13, 2021, 8:49pm

Alright! Challenge accepted. I just added a MultiplyComplexImages operation with a clij-opencl-kernel in the back made from your C++ code. Additionally, I added two functions for making complex images and for splitting them into a real and an imaginary part. You find all that together with a test as branch on the clij2-fft repository. It multiplies a complex image with itself and outputs this:

Complex input:
1.0 5.0 2.0 6.0
3.0 7.0 4.0 8.0
Computing input squared...
Complex output:
-24.0 10.0 -32.0 24.0
-40.0 42.0 -48.0 64.0
Real output
-24.0 -32.0
-40.0 -48.0
Imaginary output
10.0 24.0
42.0 64.0

Does that make sense? Apolgies for my ignorance, I don’t do complex algebra on a daily basis

Btw. you find a description of CLIJs OpenCL-dialect here.

Cheers,
Robert

bnorthan · February 14, 2021, 2:55pm

Hi @haesleinhuepf

That looks good. I haven’t had a chance to test it yet but will this week.

So it seems for CLIJ we would need to specify complex operations by name, (ie multiply_complex_images instead of just multiply).

I think that is OK. Imglib2 and Matlab have overloaded versions of most operations, which handle complex numbers automatically, however that is probably more work then it is worth for CLIJ. At least for the early iterations of this.

haesleinhuepf · February 14, 2021, 4:27pm

Yes, clij philosophy always was: There are no types, everything is an image

Let me know how it goes!

axtimwalde · February 14, 2021, 8:48pm

An image of what type?

haesleinhuepf · February 14, 2021, 9:18pm

Depends on what your GPU supports. Typically, (u)int 8/16/32, float16/32[/64] are supported. Something like a complex type doesn’t exist in OpenCL afaik.

bnorthan · February 22, 2021, 1:35pm

Hi @haesleinhuepf

I tried writing an FFT Convolution using your new complex mulitply kernel. The good news is it works. The bad news is that there are still a bunch of annoying complications that Imglib2 and Ops handle under the hood. My example is a bit messy but I’ve added a bunch of “TODOs” which hopefully help convey what needs to be done to get a “easy to use” FFT in CLIJ.

Add FFT Convolution example · clij/clij2-fft@ac938d6 · GitHub