Help needed for FFT

fft
Tags: #<Tag:0x00007fd5421744a8>

#1

Hello,

I need to compute a 3D FFT and am thus aiming to use:
Plugins>Process>Fast FFT (2D/3D) of @StephanPreibisch.

In order to practice I used this image sine_pattern.zip, where the peak to peak distance is 15 pixels. My problem is that I do not understand how to compute the number 15 from the computed FFT.

Any help is greatly appreciated!

Best, Christian

PS: I am having the same issue using this plugin: https://sites.google.com/site/piotrwendykier/software/parallelfftj


FFT understanding
#2

For an alternative approach using ImgLib2 data structures, you could look at the fft ops of ImageJ Ops. The API is still undergoing some changes and simplifications, but IMO it is worth a shot—@bnorthan has done a lot of work on it.


#3

@Christian_Tischer

Hi Christian

Below is an example script that runs an fft using ops. It wraps imglib2 so the result should be the same as the imglib2 result.

The FFT result should have a peak at 0, 0 (DC frequency) and a second peak corresponding to the frequency of your signal. The result may not be zero every where else (as you would expect with a pure frequency). It is tricky to construct a “pure” frequency with a finite signal because the signal ends up getting padded to a compatible FFT size. So even if you construct a pattern with an integer number of cycles, the padding might mess things up.

# @OpService ops
# @UIService ui
# @Dataset image

from net.imglib2.img.display.imagej import ImageJFunctions
from net.imglib2.type.numeric.complex import ComplexFloatType
from net.imglib2.outofbounds import OutOfBoundsMirrorExpWindowingFactory
from net.imglib2.outofbounds import OutOfBoundsMirrorFactory

from net.imglib2.converter import ComplexImaginaryFloatConverter
from net.imglib2.converter import ComplexPhaseFloatConverter
from net.imglib2.converter import ComplexRealFloatConverter

# perform fft of the templatehttps://www.facebook.com/

# basic fft call with no parameters
#FFT=ops.filter().fft(template.getImgPlus())

# alternatively to pass an outofbounds factory we have to pass every parameter.  We want:
# output='None', input=template, borderSize=10 by 10, fast='True', outOfBoundsFactor=OutOfBoundsMirrorExpWindowingFactory
FFT=ops.filter().fft(image.getImgPlus(), [0, 0], True, OutOfBoundsMirrorFactory(OutOfBoundsMirrorFactory.Boundary.SINGLE))
# display fft (by default in generalized log power spectrum)
ImageJFunctions.show(FFT).setTitle("fft power spectrum")

# display fft phase spectrum
ImageJFunctions.show( FFT,ComplexPhaseFloatConverter() ).setTitle( "fft phase spectrum" )

# display fft real values
ImageJFunctions.show( FFT,ComplexRealFloatConverter() ).setTitle( "fft real values" )

#4

@Christian_Tischer

Here is a groovy script that uses ops to calculate fft size, generates a sin function, and then performs fft and displays the power spectrum. If you use size 90, 90 no padding is done, and you get an FFT with an impulse at the expected frequency. If you use size 91, 91 padding is done, and the signal is no longer a “pure” frequency and you get a more complicated FFT. Let me know if you have any questions.

I am planning on polishing this and putting it in the templates. If you find parts of it confusing let me know, and I can work on explaining it better.

// @OpService ops
// @UIService ui
// @ConvertService convert
// @DatasetService data

import net.imglib2.FinalDimensions
import net.imglib2.type.numeric.real.FloatType
import net.imglib2.img.display.imagej.ImageJFunctions

// define an xSize and ySize, feel free to try different sizes, most sizes will be padded, however 90, 90 is not padded 
// TRY THIS:  After running script change size to 91, 91 and notice the differences in 'paddedSize' and final output. 
xSize=90
ySize=90

// convert original dimensions to the FinalDimensions class
origDims=new FinalDimensions([xSize, ySize] as long[])

// define arrays for the padded size and fft size
paddedSize=new long[2]
fftSize=new long[2]

// use ops to get the padded size and fft size for the original size
ops.filter().fftSize(origDims,paddedSize, fftSize, true, true)

// print out the results and check if extra padding will be done inside the fft
println "Original size: "+xSize+" "+ySize
println "Padded size: "+paddedSize[0]+" "+paddedSize[1]
// note the fft size will be n/2+1 the original size in the x dimension... 
// this is due to the fact the fft is an array of complex numbers, 
// and known to be symmetrical for a real signal, thus only n/2+1 of the coefficients
// are needed (the others are redundant). 
println "FFT size: "+fftSize[0]+" "+fftSize[1]

// now use ops to create an empty image
blank=ops.create().img(new FinalDimensions(xSize,ySize),new FloatType())

// use ops to create a sin wave pattern
formula = "Math.sin(2*Math.PI*0.1*p[0])+1"
sinImg = ops.image().equation(blank, formula)
ui.show(sinImg)

// perform fft and display using ImageJFunctions 
fft=ops.filter().fft(sinImg)
// ImageJFunctions displays the power spectrum of our complex signal
ImageJFunctions.show(fft).setTitle("fft power spectrum")

#5

@bnorthan:
Thanks so much for the answers and example code!

I updated my Fiji and tried running your examples, however I got errors for both scripts.
For the first one (running as jython script):

TypeError: fft(): 2nd arg can't be coerced to net.imglib2.RandomAccessibleInterval, net.imglib2.img.Img

And for the second one (running as groovy script):

groovy.lang.MissingMethodException: No signature of method: net.imagej.ops.filter.FilterNamespace.fftSize() is applicable for argument types: (net.imglib2.FinalDimensions, [J, [J, java.lang.Boolean, java.lang.Boolean) values: [net.imglib2.FinalDimensions@314ebd81, [0, 0], [0, 0], true, ...]
Possible solutions: fftSize([J, [J, [J, java.lang.Boolean, java.lang.Boolean)

#6

However, I got a simplified version of your code running:

// @OpService ops
// @UIService ui
// @ConvertService convert
// @DatasetService data

import net.imglib2.FinalDimensions
import net.imglib2.type.numeric.real.FloatType
import net.imglib2.img.display.imagej.ImageJFunctions

// define an xSize and ySize, feel free to try different sizes, most sizes will be padded, however 90, 90 is not padded 
// TRY THIS:  After running script change size to 91, 91 and notice the differences in 'paddedSize' and final output. 
xSize=91
ySize=91

// now use ops to create an empty image
blank=ops.create().img(new FinalDimensions(xSize,ySize),new FloatType())

// use ops to create a sin wave pattern
formula = "Math.sin(2*Math.PI*0.1*p[0])+1"
sinImg = ops.image().equation(blank, formula)
ui.show(sinImg)

// perform fft and display using ImageJFunctions 
fft=ops.filter().fft(sinImg)
// ImageJFunctions displays the power spectrum of our complex signal
ImageJFunctions.show(fft).setTitle("fft power spectrum")
'

My question is actually simply how to compute the frequencies from the resulting FFT image, i.e. what is the mapping from the x,y positions in the FFT image to the frequencies in the input image.

For instance in the example above the FFT image shows a quite strong peak at [x=39, y=0], while the frequencies in the input data are 0 in y-direction and 0.1 in x direction, right? Basically, I simply lack the mathematical formula to compute f_x=0.1 from the FFT image


#7

Hi @Christian_Tischer

I’m glad you got the code running. Yes you are correct, the frequency is 0 in y, and 0.1 cycles per sample in z. The highest frequency you can capture (nyquist) is 0.5 cycles per sample.

I changed the formula to use cos and removed the ‘+1’ dc frequency

// use ops to create a cos pattern
formula = "Math.cos((0.1)*2*Math.PI*p[0])"
cosImg = ops.image().equation(blank, formula)

I also added a line plot of the fourier transform. Here are my plots for an image size of 90,90 (a nice peak at 9) and 91,91 (peak still at 10, but many other frequencies are present). Do you get similar results??


The complete script is below.

// @OpService ops
// @UIService ui

import net.imglib2.FinalDimensions
import net.imglib2.type.numeric.real.FloatType
import net.imglib2.img.display.imagej.ImageJFunctions
import ij.IJ

// define an xSize and ySize, feel free to try different sizes, most sizes will be padded, however 90, 90 is not padded 
// TRY THIS:  After running script change size to 91, 91 and notice the differences in 'paddedSize' and final output. 
xSize=90
ySize=90

// define arrays for the padded size and fft size
origSize=[xSize, ySize] as long[]
paddedSize=new long[2]
fftSize=new long[2]

// use ops to get the padded size and fft size for the original size
ops.filter().fftSize(origSize,paddedSize, fftSize, true, true)

// print out the resuts and check if extra padding will be done inside the fft
println "Original size: "+xSize+" "+ySize
println "Padded size: "+paddedSize[0]+" "+paddedSize[1]
// note the fft size will be n/2+1 the original size in the x dimension... 
// this is due to the fact the fft is an array of complex numbers, 
// and known to be symmetrical for a real signal, thus only n/2+1 of the coefficients
// are needed (the others are redundant). 
println "FFT size: "+fftSize[0]+" "+fftSize[1]

// now use ops to create an empty image
blank=ops.create().img(new FinalDimensions(xSize,ySize),new FloatType())

// use ops to create a cos pattern
formula = "Math.cos((0.1)*2*Math.PI*p[0])"
cosImg = ops.image().equation(blank, formula)
ui.show(cosImg)

// perform fft and display using ImageJFunctions 
fft=ops.filter().fft(cosImg)
// ImageJFunctions displays the power spectrum of our complex signal
ImageJFunctions.show(fft).setTitle("fft power spectrum")

IJ.resetMinAndMax();
IJ.makeLine(0, 0, xSize/2+1, 0);
IJ.run("Plot Profile");

#8

Here is the plot I get for a frequency of 0.5 cycles per sample. There is a peak at the very end. So the x-scale actually goes from 0 to 0.5 cycles/sample or equivalently 0 to PI.


#9

Hi @bnorthan,

Thanks a lot again! Yes, I could reproduce your results!

I think I got it now: if the FFT image has Nx pixels in the x-dimension and there is a peak at xPeak, the corresponding frequency is fx = xPeak/Nx*0.5


#10

Hi @bnorthan,

The differences of the FFT images between 90x90 and 91x91 pixels are a bit shocking to me, as the frequencies in the image do not ‘really’ change.

I also tried this plugin to compute the FFT: http://imagej.net/plugins/fftj.html

And, interestingly, the code used in this plugin is unimpressed by the 91x91 issue and still shows the same nice and clean FFT as for the 90x90 image. However the code runs much slower than yours. Do you happen to know the reason for this?

Maybe if we run your code with an OutOfBoundsFactory that is not mirroring but rather 0-padding, we would see less or no difference between the 90 vs 91 image size?! (as said this part of the code crashed for me)


#11

Hi Christian

I suspect in that plugin the image is not extended, which produces a nicer result for a simple cosine. However the reason the image is often extended is because there are certain sizes that an FFT will run much faster at (+ edge artifacts occur for realistic signals (it only works out nicely for periodic signals with an integer number of cycles, such that they “wrap” smoothly to the other side)).

When you extend the image you can introduce a sharp edge. The FFT code does default to “zero-pad”, which means, in this case, the signal is a multiplication of a cosine and a square in the spatial domain, which makes it an impulse convolved with a sync function in frequency, which is the result you see.

It seems counter intuitive at first but if you look at tables of fourier transforms you see that simple signals (like deltas, and square waves) have complicated looking fourier transforms (ie see this google image search for ‘common fourier pairs’.

I am currently working on some changes to the fft and filter interfaces, it should be integrated into Fiji soon and I can make some examples showing the effects of different padding strategies on FFT and Convolution.


#12

Hi Brian,
Thanks for the explanations!
Talking also to other people it seems that some way of producing smooth edges at the (padded) image boundaries is the way to go. However I am not sure about the best way to do this. Do you have ideas?


#13

Hi Christian

Ops is built on imglib2 which has several OutOfBounds factories. You can pass different OutOfBounds factories to the various fft and filtering ops. Ops also exposes the padding operations. So you can display the padded signals for verification. Using these ops you can experiment with different out of bounds strategies and see how they effect the results.

The only problem is that the interfaces are going to change slightly in the next week or two (which was why my scripts didn’t work for you, they were written for an un-released version of ops).


#14

Hi Brian,

This sounds great! Maybe you could let me know when I can try it out?!

I have two more questions:

To me, it looks like as if the Fourier-spectrum that your function shows is the logarithm of the actual spectrum. Is that true? If so, could you add an option to turn the logarithm on and off?

You are using image.equation() to generate images; is it somewhere documented with string expressions, such as “Math.sin()”, one can use there?


#15

I released imagej-ops 0.30.0 and uploaded it to the Java-8 update site, along with several other core components. So if you update your Fiji, you should have the latest version now, including the FFT refactoring that @bnorthan mentioned.

The image.equation op uses JavaScript to evaluate the expression, so you can do (nearly) anything you can do in JS. @hinerm is currently investigating potentially beefing up that op to support mutation of existing images, a la the Process :arrow_forward: Math :arrow_forward: Macro command.


#16

Hi @Christian_Tischer

Here is a script that pads the image two ways (mirror and zero) and then performs the fft on each. You could pass the BoundaryFactory directly to the FFT, but I explicitly padded the images instead, so we can see exactly what the padded image looks like. If you let me know more details on how you are trying to filter your image I might be able to come up with other examples that are relevant to you.


// @OpService ops
// @UIService ui

import net.imglib2.FinalDimensions
import net.imglib2.type.numeric.real.FloatType
import net.imglib2.img.display.imagej.ImageJFunctions
import net.imglib2.outofbounds.OutOfBoundsMirrorFactory
import net.imglib2.outofbounds.OutOfBoundsMirrorFactory.Boundary
import net.imglib2.outofbounds.OutOfBoundsConstantValueFactory
import ij.IJ

// define an xSize and ySize, feel free to try different sizes, most sizes will be padded (if doing fast fft), however 90, 90 is not padded 
// 91, 91 will be padded  
xSize=91
ySize=91

// the ops need the size in the 'Dimensions' format
origDims=new FinalDimensions([xSize, ySize] as long[])

// use ops to create an empty image
blank=ops.create().img(new FinalDimensions(xSize,ySize),new FloatType())

// use ops to create a cos pattern
formula = "Math.cos((0.5)*2*Math.PI*p[0])"
cosImg = ops.image().equation(blank, formula)
ui.show(cosImg)

// pad the input using two different strategies
padded1=ops.filter().padFFTInput(cosImg, origDims, true, new OutOfBoundsMirrorFactory(Boundary.SINGLE))
padded2=ops.filter().padFFTInput(cosImg, origDims, true, new OutOfBoundsConstantValueFactory(new FloatType()))

// create buffers for output
output1=ops.filter().createFFTOutput(origDims, new FloatType())
output2=ops.filter().createFFTOutput(origDims, new FloatType())

// use ops to perform fft on padded and place in the output
ops.filter().fft(output1, padded1)
ops.filter().fft(output2, padded2)

// show the padded signals
ui.show("padded1", padded1)
ui.show("padded2", padded2)

// now show the outputs and their line profiles. 
ImageJFunctions.show(output1).setTitle("power spectrum output1")
ui.showDialog("1st show power spectrum mirror")
IJ.resetMinAndMax();
IJ.makeLine(0, 0, output1.dimension(0)-1, 0);
IJ.run("Plot Profile");
ImageJFunctions.show(output2).setTitle("power spectrum output2")

ui.showDialog("Now power spectrum of zero pad")
IJ.resetMinAndMax();
IJ.makeLine(0, 0, output1.dimension(0)-1, 0);
IJ.run("Plot Profile");


#17

Hi Brian,
Thanks a lot for the educative examples!