Gaussian fit to FFT Peak

Hi All, I am writing an image analysis program to use with Interfermotey data. In the second step I have a peak as shown:

One of the issues I am running into is that the maximum vlaue of this peak is not always the centre of the peak, and want to try and fit a gausssian to the data, and use that to find the centre.

I have never done any 2-D gaussian fitting, and so am a little out of my depth.

The actual question is: Is there a funtion that can fit a Gassian to a matrix, and output the location of the centre of the peak.

Many thanks in advnaced

Answer 1

[Bjorn Gustavsson]

Sure, you can do something like this:

fcn_G2D = @(p,x,y) p(1)*exp(-(x-p(2)).^2/p(4)-(y-p(3)).^2/p(5));
err_fcn = @(p,x,y,I,fcn) sum(sum((I-fcn(p,x,y)).^2)); % is using fminsearch
res_fcn = @(p,x,y,I,fcn) (I-fcn(p,x,y)); % if using lsqnonlin
fftI = fftshift(fft2(I));
x = 1:size(fftI,2);
y = 1:size(fftI,1);
[x,y] = meshgrid(x,y);
par0 = [sum(I(:)),size(fftI,2)/2,size(fftI,1)/2,3,3];
parFMS = fminsearch(@(p) err_fcn(p,x,y,abs(fftI),fcn_G2D),par0);
parNLS = lsqnonlin(@(p) res_fcn(p,x,y,abs(fftI),fcn_G2D),par0);

Now tested, should be OK...

You will have to check the convergence of the optimizations...

HTH

[Hugo Fortescue]

Hi There,

Would you be able to explain very quicky how the code works? It has an error and I am tryignt to debugg. Many thanks

[Bjorn Gustavsson]

OK, the first row is the defintion of fcn_G2D simply a 2-D Gaussian, you should really make it such that it accepts a rotation, this variant is constrained to have the semi-major and semi-minor axes aligned with the horizontal and vertical axes. You can do that reasonably easy by explicitly including a rotation.

The second row should've been a proper definition of a sum-of-squares function taken between an image and a function on a x-y-grid with some parameters. This is what fminsearch expects to do its minimization

The third row should've been a proper definition of a residuals function taken between an image and a function on a x-y-grid with some parameters. This is what lsqnonlin expects to do its minimization.

Then we have the definition of our x-y-grid.

Followed by a simple guesstimate of suitable starting parameters for the optimization, and calls to fminsearch and lsqnonlin.

HTH