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
Dear @haesleinhuepf, is there a fft filter in clij2? I thought so, since there is this page:
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
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
Iāve also discovered now the clijx correct drift commands, will give those a try.
I guess the fft was just a way to get the backgroundā¦ not to calculate driftā¦ sorry to bug you!
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
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
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.
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.
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
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
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
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
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
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.
Yes, clij philosophy always was: There are no types, everything is an image
Let me know how it goes!
An image of what type?
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.
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
Iāve made some progress working towards better FFTs. See the latest iteration of FFT Convolve. I now have a 3D version, which calls the 3D FFT from a Utility.
The next step is to make it a plugin. One complication is the FFT is a different size than the input. How we deal with outputs of different size? From what I understand if we call an image ādestinationā it is generated with specific size? However is there a way to generate it with a size different than input, but based on input size?