Login  Register

Re: FFT along line in image

Posted by seb-7 on Oct 24, 2007; 10:10pm
URL: http://imagej.273.s1.nabble.com/FFT-along-line-in-image-tp3698130p3698131.html

Jon Harman wrote:

> Hi,
>
> Is there a plugin to ImageJ that can calculate the FFT along a line in
> the image and display the result (say for instance plot the power
> spectrum)?
>
> Jon
>
> .
>


Hi Jon,

Just for fun,
here is an infamous SlowFT macro, not even well tested.
A better idea would be to use an external FFT lib, like
jnt.fft or others you could use for 1D transform.
http://math.nist.gov/~BMiller/java/
http://www.google.com/codesearch?hl=en&lr=&q=fft+lang%3Ajava+&btnG=Search

and last but not least
http://rsb.info.nih.gov/ij/developer/source/ij/process/FHT.java.html
:-)

seb



SlowFT.txt
-----------------------------------------------------
//"Brute Force" Fourier Transform
//very naive implementation of usual DFTs
//the slowest DFT you can find
//neither symmetry/parity/memoize nor "butterfly" trick
//just raw and slow Discrete Fourier Transform
//unusable with large arrays!

//actually, large is quite small :-)
//direct+inverse transform took 13 sec for 800 pixels
//on a (old) Pentium4 1.6Ghz

var PI=3.1415926536;


macro "LineFTTest"
{
    if (selectionType != 5)
       exit("Line selection required!");

    x=getProfile();
    n=x.length;
    setBatchMode(true);
    T0=getTime();
    F=slowRFT1D(x);
    ift=slowIFT1D(F);
    T1=getTime();
    setBatchMode(false);
    name = "[SlowFT]";
    run("New... ", "name="+name+" type=Table");
    f = name;
 
print(f,"\\Headings:x\tReF=Re(FT(x))\tImF=Im(FT(x))\tRe(IFT(F))\tIm(IFT(F)");
    for(i=0;i<n;i++)
    {
       print(f,x[i]+"\t"+F[i]+"\t"+F[n+i]+"\t"+ift[i]+"\t"+ift[i+n]);
    }
    print("direct and inverse transform of "+n+" values took "+(T1-T0)+
"ms");
}
       I


//direct real FT
function slowRFT1D(real_in)
{
    N=real_in.length;
    S=newArray(2*N);

    for (m=0;m<N;m++)
    {

       RSm=0;
       ISm=0;
       for (k=0;k<N;k++)
       {
          RSm+=real_in[k]*cos(2*PI*k*m/N);
          ISm-=real_in[k]*sin(2*PI*k*m/N);
       }
       S[m]=RSm; //real part
       S[N+m]=ISm; //imaginary part
    }
    return S;
    //S:Complex: {Re[0],..Re[N-1],Im[0]..Im[N-1]}
}

//direct complex FT

function slowFT1D(z_in)
{//z_in complex: {Re[0]..Re[N-1],Im[0]..Im[N-1]}
    N=z_in.length/2;
    S=newArray(2*N);
    for(m=0;m<N;m++)
    {
       RSm=0;
       ISm=0;
       for (k=0;k<N;k++)
       {
          c=cos(2*PI*k*m/N);
          s=-sin(2*PI*k*m/N);
          RSm+=(z_in[k]*c-s*z_in[k+N]);
          ISm+=(z_in[k]*s+c*z_in[k+N]);
       }
       S[m]=RSm;
       S[m+N]=ISm;
    }
    return S;
    // output= complex (as z_in...)
}
//Inverse complex FT
function slowIFT1D(S_in)
{ //S_in complex: {Re[0]..Re[N-1],Im[0]..Im[N-1]}
    N=S_in.length/2;
    sout=newArray(2*N);
    for (k=0;k<N;k++)
    {
       Rsk=0;
       Isk=0;
       for (m=0;m<N;m++)
       {
          c=cos(2*PI*k*m/N);
          s=sin(2*PI*k*m/N);
          Rsk += (S_in[m]*c-s*S_in[m+N]);
          Isk += (S_in[m]*s+c*S_in[m+N]);
       }
       sout[k] = Rsk/N;
       sout[N+k] = Isk/N;
    }
    return sout;
}