Autocorrelation using Fourier

[vischer]

Hello all,

I created a test macro that performs 1D autocorrelation.
It runs well in a few seconds, but I wonder if the speed can be improved by Fourier techniques for larger data sets.

The macro first creates a test image that is randomly covered with white horizontal lines with a mean length of 20 pixels. The image is then periodically multiplied with a horizontally shifted copy, so I get the autocorrelation plot from which I can derive the mean line length. The output looks like this:

  Simulated Length = 20.17
  Measured Length = 20.06

Later in reality, the white lines will represent the fluorescent duration of molecules, and the data set will be much larger.
I doubt that a native plugin will do the same thing much faster, but has anyone a suggestion how to use Fourier techniques?

My demo macro can be downloaded from:
http://simon.bio.uva.nl/objectj/GenericMacros/AutoCorrelation/AutoCorrelation-2.txt
 

Norbert Vischer

[Michael Schmid]

Hi Norbert,

here is a rough equivalent to the core of your correlation macro  
using FFT:

k100 = 100; //how far out we need the autocorrelation
run("FD Math...", "image1=AutoCorr.tif operation=Correlate  
image2=AutoCorr.tif result=Result do");
width = getWidth();
run("Divide...", "value="+(width*width)); //FD Math gives sum over  
all pixels; convert to mean
half = width/2;  //origin is at the center
makeLine(half, half, half+k100, half);
run("Plot Profile");

It won't be 100% exact because the FFT version extends the image  
assuming periodic boundary conditions, not using a value of zero as  
the 'translate' command in your macro.
This won't matter for your test image which has a rather large border  
with zero value.
If periodic boundary conditions are a problem, simply enlarge the  
image to the next power of 2 with zero background.

Hope this helps.

Michael
________________________________________________________________


On 13 Dec 2010, at 23:19, Norbert Vischer wrote:

> Hello all,
>
> I created a test macro that performs 1D autocorrelation.
> It runs well in a few seconds, but I wonder if the speed can be  
> improved by Fourier techniques for larger data sets.
>
> The macro first creates a test image that is randomly covered with  
> white horizontal lines with a mean length of 20 pixels. The image  
> is then periodically multiplied with a horizontally shifted copy,  
> so I get the autocorrelation plot from which I can derive the mean  
> line length. The output looks like this:
>
>   Simulated Length = 20.17
>   Measured Length = 20.06
>
> Later in reality, the white lines will represent the fluorescent  
> duration of molecules, and the data set will be much larger.
> I doubt that a native plugin will do the same thing much faster,  
> but has anyone a suggestion how to use Fourier techniques?
>
> My demo macro can be downloaded from:
> http://simon.bio.uva.nl/objectj/GenericMacros/AutoCorrelation/ 
> AutoCorrelation-2.txt
>
>
> Norbert Vischer

[Norbert Vischer-2]

Hi Michael,

thanks for this simple and elegant solution. I updated my macro, so the user
can now choose between an FFT and non-FFT macro:

http://simon.bio.uva.nl/objectj/GenericMacros/AutoCorrelation/


Norbert

[dpoburko]

Michael,

    This little gem looks like it should be very helpful with automated
registration of paired, offset images. I am probably missing something
simple, but the correlation for partially overlapping images seems to
fail when the alignment offset of two images  is more than half the
width or height of the image. Could you suggest a method to work around
this?

Thanks,
Damon

On 12/14/2010 1:39 AM, Michael Schmid wrote:

> Hi Norbert,
>
> here is a rough equivalent to the core of your correlation macro using
> FFT:
>
> k100 = 100; //how far out we need the autocorrelation
> run("FD Math...", "image1=AutoCorr.tif operation=Correlate
> image2=AutoCorr.tif result=Result do");
> width = getWidth();
> run("Divide...", "value="+(width*width)); //FD Math gives sum over all
> pixels; convert to mean
> half = width/2;  //origin is at the center
> makeLine(half, half, half+k100, half);
> run("Plot Profile");
>
> It won't be 100% exact because the FFT version extends the image
> assuming periodic boundary conditions, not using a value of zero as
> the 'translate' command in your macro.
> This won't matter for your test image which has a rather large border
> with zero value.
> If periodic boundary conditions are a problem, simply enlarge the
> image to the next power of 2 with zero background.
>
> Hope this helps.
>
> Michael
> ________________________________________________________________
>
>
> On 13 Dec 2010, at 23:19, Norbert Vischer wrote:
>
>> Hello all,
>>
>> I created a test macro that performs 1D autocorrelation.
>> It runs well in a few seconds, but I wonder if the speed can be
>> improved by Fourier techniques for larger data sets.
>>
>> The macro first creates a test image that is randomly covered with
>> white horizontal lines with a mean length of 20 pixels. The image is
>> then periodically multiplied with a horizontally shifted copy, so I
>> get the autocorrelation plot from which I can derive the mean line
>> length. The output looks like this:
>>
>>   Simulated Length = 20.17
>>   Measured Length = 20.06
>>
>> Later in reality, the white lines will represent the fluorescent
>> duration of molecules, and the data set will be much larger.
>> I doubt that a native plugin will do the same thing much faster, but
>> has anyone a suggestion how to use Fourier techniques?
>>
>> My demo macro can be downloaded from:
>> http://simon.bio.uva.nl/objectj/GenericMacros/AutoCorrelation/AutoCorrelation-2.txt 
>>
>>
>>
>> Norbert Vischer

[Michael Schmid]

Hi Damon,

you are right, that's a fundamental property of the Fourier transform:
You have periodic boundary conditions, so, e.g., a shift right by 0.6  
image widths is equivalent to a shift left by 0.4 image widths.

What you could try:
Enlarge the images, with a background equal to the average brightness  
of the image, and preferably some smearing from the image border into  
the background to avoid artifacts caused by the shap edge.

For the correlation image, you will typically need some highpass  
filter or background subtraction to find the well-defined maximum.
Here is a quick example, using my 'Fast Filters' for highpass filtering
   http://imagejdocu.tudor.lu/doku.php?
id=plugin:filter:fast_filters:start

run("Boats (356K)");

makeRectangle(117, 78, 333, 333);
run("Duplicate...", "title=boats-1.gif");
run("Canvas Size...", "width=512 height=512 position=Center zero");
makeRectangle(89, 89, 333, 333);  //the previous size
run("Make Inverse");
run("Set...", "value=140");       //background: average of image
run("Gaussian Blur...", "sigma=5");
run("Gaussian Blur...", "sigma=4");
run("Gaussian Blur...", "sigma=3");

selectWindow("boats.gif");
makeRectangle(337, 226, 333, 333); //x shift 220, y shift 148
run("Duplicate...", "title=boats-2.gif");
run("Canvas Size...", "width=512 height=512 position=Center zero");
makeRectangle(89, 89, 333, 333);
run("Make Inverse");
run("Set...", "value=90");
run("Gaussian Blur...", "sigma=5");
run("Gaussian Blur...", "sigma=4");
run("Gaussian Blur...", "sigma=3");

run("FD Math...", "image1=boats-1.gif operation=Correlate  
image2=boats-2.gif result=Result do");
run("Fast Filters", "link filter=[background from minima] x=20 y=20  
preprocessing=none subtract offset=0");
run("Enhance Contrast", "saturated=0.001");

(remove any linebreaks caused by the mailer)

Michael
________________________________________________________________

On 14 Dec 2010, at 21:04, Damon Poburko wrote:

> Michael,
>
>    This little gem looks like it should be very helpful with  
> automated registration of paired, offset images. I am probably  
> missing something simple, but the correlation for partially  
> overlapping images seems to fail when the alignment offset of two  
> images  is more than half the width or height of the image. Could  
> you suggest a method to work around this?
>
> Thanks,
> Damon
>
> On 12/14/2010 1:39 AM, Michael Schmid wrote:
>> Hi Norbert,
>>
>> here is a rough equivalent to the core of your correlation macro  
>> using FFT:
>>
>> k100 = 100; //how far out we need the autocorrelation
>> run("FD Math...", "image1=AutoCorr.tif operation=Correlate  
>> image2=AutoCorr.tif result=Result do");
>> width = getWidth();
>> run("Divide...", "value="+(width*width)); //FD Math gives sum over  
>> all pixels; convert to mean
>> half = width/2;  //origin is at the center
>> makeLine(half, half, half+k100, half);
>> run("Plot Profile");
>>
>> It won't be 100% exact because the FFT version extends the image  
>> assuming periodic boundary conditions, not using a value of zero  
>> as the 'translate' command in your macro.
>> This won't matter for your test image which has a rather large  
>> border with zero value.
>> If periodic boundary conditions are a problem, simply enlarge the  
>> image to the next power of 2 with zero background.
>>
>> Hope this helps.
>>
>> Michael
>> ________________________________________________________________
>>
>>
>> On 13 Dec 2010, at 23:19, Norbert Vischer wrote:
>>
>>> Hello all,
>>>
>>> I created a test macro that performs 1D autocorrelation.
>>> It runs well in a few seconds, but I wonder if the speed can be  
>>> improved by Fourier techniques for larger data sets.
>>>
>>> The macro first creates a test image that is randomly covered  
>>> with white horizontal lines with a mean length of 20 pixels. The  
>>> image is then periodically multiplied with a horizontally shifted  
>>> copy, so I get the autocorrelation plot from which I can derive  
>>> the mean line length. The output looks like this:
>>>
>>>   Simulated Length = 20.17
>>>   Measured Length = 20.06
>>>
>>> Later in reality, the white lines will represent the fluorescent  
>>> duration of molecules, and the data set will be much larger.
>>> I doubt that a native plugin will do the same thing much faster,  
>>> but has anyone a suggestion how to use Fourier techniques?
>>>
>>> My demo macro can be downloaded from:
>>> http://simon.bio.uva.nl/objectj/GenericMacros/AutoCorrelation/ 
>>> AutoCorrelation-2.txt
>>>
>>>
>>> Norbert Vischer

[dpoburko]

Hello again,

   Thanks for the tip. Padding the images up to the next power of two
worked nicely. I didn't try smoothing the transition from the image to
the padded region. Rather, I added a little bit of noise, then ran the
correlation image through a band pass filter. Seems like it works just
fine. Many thanks!!! I've been scheming for months over how to implement
this, not knowing it was under my nose the whole time.

Damon

On 12/14/2010 12:44 PM, Michael Schmid wrote:

> Hi Damon,
>
> you are right, that's a fundamental property of the Fourier transform:
> You have periodic boundary conditions, so, e.g., a shift right by 0.6
> image widths is equivalent to a shift left by 0.4 image widths.
>
> What you could try:
> Enlarge the images, with a background equal to the average brightness
> of the image, and preferably some smearing from the image border into
> the background to avoid artifacts caused by the shap edge.
>
> For the correlation image, you will typically need some highpass
> filter or background subtraction to find the well-defined maximum.
> Here is a quick example, using my 'Fast Filters' for highpass filtering
>   http://imagejdocu.tudor.lu/doku.php?id=plugin:filter:fast_filters:start
>
> run("Boats (356K)");
>
> makeRectangle(117, 78, 333, 333);
> run("Duplicate...", "title=boats-1.gif");
> run("Canvas Size...", "width=512 height=512 position=Center zero");
> makeRectangle(89, 89, 333, 333);  //the previous size
> run("Make Inverse");
> run("Set...", "value=140");       //background: average of image
> run("Gaussian Blur...", "sigma=5");
> run("Gaussian Blur...", "sigma=4");
> run("Gaussian Blur...", "sigma=3");
>
> selectWindow("boats.gif");
> makeRectangle(337, 226, 333, 333); //x shift 220, y shift 148
> run("Duplicate...", "title=boats-2.gif");
> run("Canvas Size...", "width=512 height=512 position=Center zero");
> makeRectangle(89, 89, 333, 333);
> run("Make Inverse");
> run("Set...", "value=90");
> run("Gaussian Blur...", "sigma=5");
> run("Gaussian Blur...", "sigma=4");
> run("Gaussian Blur...", "sigma=3");
>
> run("FD Math...", "image1=boats-1.gif operation=Correlate
> image2=boats-2.gif result=Result do");
> run("Fast Filters", "link filter=[background from minima] x=20 y=20
> preprocessing=none subtract offset=0");
> run("Enhance Contrast", "saturated=0.001");
>
> (remove any linebreaks caused by the mailer)
>
> Michael
> ________________________________________________________________
>
> On 14 Dec 2010, at 21:04, Damon Poburko wrote:
>
>> Michael,
>>
>>    This little gem looks like it should be very helpful with
>> automated registration of paired, offset images. I am probably
>> missing something simple, but the correlation for partially
>> overlapping images seems to fail when the alignment offset of two
>> images  is more than half the width or height of the image. Could you
>> suggest a method to work around this?
>>
>> Thanks,
>> Damon
>>
>> On 12/14/2010 1:39 AM, Michael Schmid wrote:
>>> Hi Norbert,
>>>
>>> here is a rough equivalent to the core of your correlation macro
>>> using FFT:
>>>
>>> k100 = 100; //how far out we need the autocorrelation
>>> run("FD Math...", "image1=AutoCorr.tif operation=Correlate
>>> image2=AutoCorr.tif result=Result do");
>>> width = getWidth();
>>> run("Divide...", "value="+(width*width)); //FD Math gives sum over
>>> all pixels; convert to mean
>>> half = width/2;  //origin is at the center
>>> makeLine(half, half, half+k100, half);
>>> run("Plot Profile");
>>>
>>> It won't be 100% exact because the FFT version extends the image
>>> assuming periodic boundary conditions, not using a value of zero as
>>> the 'translate' command in your macro.
>>> This won't matter for your test image which has a rather large
>>> border with zero value.
>>> If periodic boundary conditions are a problem, simply enlarge the
>>> image to the next power of 2 with zero background.
>>>
>>> Hope this helps.
>>>
>>> Michael
>>> ________________________________________________________________
>>>
>>>
>>> On 13 Dec 2010, at 23:19, Norbert Vischer wrote:
>>>
>>>> Hello all,
>>>>
>>>> I created a test macro that performs 1D autocorrelation.
>>>> It runs well in a few seconds, but I wonder if the speed can be
>>>> improved by Fourier techniques for larger data sets.
>>>>
>>>> The macro first creates a test image that is randomly covered with
>>>> white horizontal lines with a mean length of 20 pixels. The image
>>>> is then periodically multiplied with a horizontally shifted copy,
>>>> so I get the autocorrelation plot from which I can derive the mean
>>>> line length. The output looks like this:
>>>>
>>>>   Simulated Length = 20.17
>>>>   Measured Length = 20.06
>>>>
>>>> Later in reality, the white lines will represent the fluorescent
>>>> duration of molecules, and the data set will be much larger.
>>>> I doubt that a native plugin will do the same thing much faster,
>>>> but has anyone a suggestion how to use Fourier techniques?
>>>>
>>>> My demo macro can be downloaded from:
>>>> http://simon.bio.uva.nl/objectj/GenericMacros/AutoCorrelation/AutoCorrelation-2.txt 
>>>>
>>>>
>>>>
>>>> Norbert Vischer

[Albert Cardona-2]

2010/12/14 Damon Poburko <[hidden email]>:
> Hello again,
>
>  Thanks for the tip. Padding the images up to the next power of two worked
> nicely. I didn't try smoothing the transition from the image to the padded
> region. Rather, I added a little bit of noise, then ran the correlation
> image through a band pass filter. Seems like it works just fine. Many
> thanks!!! I've been scheming for months over how to implement this, not
> knowing it was under my nose the whole time.


Dear Damon:

http://pacific.mpi-cbg.de/wiki/index.php/Stitching_2D/3D

It's even macro recordable. Uses a multithreaded FFT implementation,
and will work on 2D images and on 3D stacks. The padding is automatic
and done smartly with mirroring to prevent image edge effects.

Using it's internals, we put up a script to correct for 3D drift in 4D
image acquisitions. See "Plugins > Registration > Correct 3D drift"

Best,

Albert
--
http://albert.rierol.net