Thread Subject:

2D gaussian fit of image

HI there,

I?d like to fit a gaussian to a bead image representing the
PSF of my microscope. The ultimate goal is to define the
'exact' position of the peak in x and y. I need to fit a
gaussian to the image rather than just calculating the
center of mass via centroids.
Could someone give me a hint how to do so with an image, are
there eventually example codes available (haven?t yet found
a proper one)?

Any pointers are much appreciated!

Thanks and best regards,
Steffen

No one with an idea??

Thanks,
Steffen

> HI there,
>
> I?d like to fit a gaussian to a bead image representing the
> PSF of my microscope. The ultimate goal is to define the
> 'exact' position of the peak in x and y. I need to fit a
> gaussian to the image rather than just calculating the
> center of mass via centroids.
> Could someone give me a hint how to do so with an image, are
> there eventually example codes available (haven?t yet found
> a proper one)?
>
> Any pointers are much appreciated!
>
> Thanks and best regards,
> Steffen

On Mar 4, 4:00=A0am, "Steffen " <rile...@gmail.com> wrote:
> HI there,
>
> I?d like to fit a gaussian to a bead image representing the
> PSF of my microscope. The ultimate goal is to define the
> 'exact' position of the peak in x and y. I need to fit a
> gaussian to the image rather than just calculating the
> center of mass via centroids.
> Could someone give me a hint how to do so with an image, are
> there eventually example codes available (haven?t yet found
> a proper one)?
>
> Any pointers are much appreciated!
>
> Thanks and best regards,
> Steffen

How about taking the log of the image and then fitting the resulting
image to a quadratic?
image =3D a*exp( b*(x-x0)^2)
log(image/a) =3D b*(x-x0)^2
0 =3D c*x2 + d*x + e - log(image)
Do least square to solve for c, d, e, then plug back into your
original equation. Do in two dimensions. I'll leave the math details
to you but this is a common way of doing it.
Regards,
ImageAnalyst

ImageAnalyst <imageanalyst@mailinator.com> wrote in message
<e18c27e3-7cc1-4805-863c-
eae148e22e6f@34g2000hsz.googlegroups.com>...
> On Mar 4, 4:00=A0am, "Steffen " <rile...@gmail.com> wrote:
> > HI there,
> >
> > I?d like to fit a gaussian to a bead image representing the
> > PSF of my microscope. The ultimate goal is to define the
> > 'exact' position of the peak in x and y. I need to fit a
> > gaussian to the image rather than just calculating the
> > center of mass via centroids.
> > Could someone give me a hint how to do so with an image, are
> > there eventually example codes available (haven?t yet found
> > a proper one)?
> >
> > Any pointers are much appreciated!
> >
> > Thanks and best regards,
> > Steffen
>
> How about taking the log of the image and then fitting the resulting
> image to a quadratic?
> image =3D a*exp( b*(x-x0)^2)
> log(image/a) =3D b*(x-x0)^2
> 0 =3D c*x2 + d*x + e - log(image)
> Do least square to solve for c, d, e, then plug back into your
> original equation. Do in two dimensions. I'll leave the math details
> to you but this is a common way of doing it.
> Regards,
> ImageAnalyst

Multiple problems with logging the image data
to fit a quadratic model:

1. This will tend to accentuate the fit in exactly
the WRONG areas, i.e., where the peak is not.
Logging the data will bring up the small numbers.
In effect, logging the image implicitly assume a
lognormal error structure in the regression model.

2. The model for a general gaussian in 2-d
requires a positive definite covariance matrix,
something that a linear regression will not be
able to constrain. In the logged domain, I would
not be amazed to find that you have just fit a
hyperbolic surface to your data. Remember that
once you log the image, it brings up those tails.
Crap in the corners will suddenly become
important.

3. Logging an image will be problematic if you
have zero pixels.

4. If you really wanted to fit a gaussian model
PLUS a constant offset in z, this is impossible
to do so via linear regression on the logged
image data.
 
There are probably some other reasons I've
skipped over as to why logging the image is
not the right way to do this.

How would I do it?

I'd use lsqnonlin or lsqcurvefit to fit a model.
Generate the pixels coordinates for each pixel
with meshgrid first.

Formulate the 2-d gaussian model using
a covariance matrix built from a Cholesky
factorization. You use three parameters to
define the cholesky factor, not the actual
covariance matrix. This assures you that
the covariance matrix is ALWAYS positive
definite. (Cute, huh?)

Constrain the gaussian parameters in
logical ways. For example, an amplitude
parameter must be positive. The gaussian
mode must lie inside the image boundaries.

Thus, if img is the image array, then your
code might look vaguely like this, for a
gaussian peak with a constant offset in z.


[ny,nx] = size(img);
[px,py] = meshgrid(1:nx,1:ny);

% starting values:
% 1 - constant offset
% 2 - mode_x
% 3 - mode_y
% 4 - amplitude
% 5,6,7 - cholesky parameters for cov matrix
params0 = [mean(img(:)),nx/2,ny/2,1,5,5,0];
LB = [-inf,1,1,0,-inf,-inf-inf];
UB = [inf,nx,ny,inf,inf,inf,inf];

options = optimset('lsqnonlin');
options.Display = 'iter';
params = lsqnonlin(@(P) objfun(P,px,py,img),params0,LB,UB,options);

cov = [params(5),0;params(6),params(7)];
cov = cov*cov';

% ==============================
% The objective function will be vaguely like this:
function resids = objfun(params,px,py,img);
covchol = [params(5),0;params(6),params(7)];
temp = [px(:) - params(2),py(:)-params(3)]*covchol;
pred = params(1) + params(4)*exp(-sum(temp.*temp,2)/2);
resids = pred - img(:);


The above code should be reasonably close,
although I predict roughly one typo since
I've not tested the code.

I you don't have the optimization toolbox
to run lsqnonlin, then get it.

John

Hi John,

thanks for ur advices. Looks a little more complicated than
i thought it would be, but I hope I can manage it to work.

Cheers,
Steffen

On Mar 5, 8:31=A0am, "John D'Errico" <woodch...@rochester.rr.com> wrote:
> ImageAnalyst <imageanal...@mailinator.com> wrote in message
>
> <e18c27e3-7cc1-4805-863c-
> eae148e22...@34g2000hsz.googlegroups.com>...
>
>
>
>
>
> > On Mar 4, 4:00=3DA0am, "Steffen " <rile...@gmail.com> wrote:
> > > HI there,
>
> > > I?d like to fit a gaussian to a bead image representing the
> > > PSF of my microscope. The ultimate goal is to define the
> > > 'exact' position of the peak in x and y. I need to fit a
> > > gaussian to the image rather than just calculating the
> > > center of mass via centroids.
> > > Could someone give me a hint how to do so with an image, are
> > > there eventually example codes available (haven?t yet found
> > > a proper one)?
>
> > > Any pointers are much appreciated!
>
> > > Thanks and best regards,
> > > Steffen
>
> > How about taking the log of the image and then fitting the resulting
> > image to a quadratic?
> > image =3D3D a*exp( b*(x-x0)^2)
> > log(image/a) =3D3D b*(x-x0)^2
> > 0 =3D3D c*x2 + d*x + e - log(image)
> > Do least square to solve for c, d, e, then plug back into your
> > original equation. =A0Do in two dimensions. =A0I'll leave the math detai=
ls
> > to you but this is a common way of doing it.
> > Regards,
> > ImageAnalyst
>
> Multiple problems with logging the image data
> to fit a quadratic model:
>
> 1. This will tend to accentuate the fit in exactly
> the WRONG areas, i.e., where the peak is not.
> Logging the data will bring up the small numbers.
> In effect, logging the image implicitly assume a
> lognormal error structure in the regression model.
>
> 2. The model for a general gaussian in 2-d
> requires a positive definite covariance matrix,
> something that a linear regression will not be
> able to constrain. In the logged domain, I would
> not be amazed to find that you have just fit a
> hyperbolic surface to your data. Remember that
> once you log the image, it brings up those tails.
> Crap in the corners will suddenly become
> important.
>
> 3. Logging an image will be problematic if you
> have zero pixels.
>
> 4. If you really wanted to fit a gaussian model
> PLUS a constant offset in z, this is impossible
> to do so via linear regression on the logged
> image data.
>
> There are probably some other reasons I've
> skipped over as to why logging the image is
> not the right way to do this.
>
> How would I do it?
>
> I'd use lsqnonlin or lsqcurvefit to fit a model.
> Generate the pixels coordinates for each pixel
> with meshgrid first.
>
> Formulate the 2-d gaussian model using
> a covariance matrix built from a Cholesky
> factorization. You use three parameters to
> define the cholesky factor, not the actual
> covariance matrix. This assures you that
> the covariance matrix is ALWAYS positive
> definite. (Cute, huh?)
>
> Constrain the gaussian parameters in
> logical ways. For example, an amplitude
> parameter must be positive. The gaussian
> mode must lie inside the image boundaries.
>
> Thus, if img is the image array, then your
> code might look vaguely like this, for a
> gaussian peak with a constant offset in z.
>
> [ny,nx] =3D size(img);
> [px,py] =3D meshgrid(1:nx,1:ny);
>
> % starting values:
> % 1 - constant offset
> % 2 - mode_x
> % 3 - mode_y
> % 4 - amplitude
> % 5,6,7 - cholesky parameters for cov matrix
> params0 =3D [mean(img(:)),nx/2,ny/2,1,5,5,0];
> LB =3D [-inf,1,1,0,-inf,-inf-inf];
> UB =3D [inf,nx,ny,inf,inf,inf,inf];
>
> options =3D optimset('lsqnonlin');
> options.Display =3D 'iter';
> params =3D lsqnonlin(@(P) objfun(P,px,py,img),params0,LB,UB,options);
>
> cov =3D [params(5),0;params(6),params(7)];
> cov =3D cov*cov';
>
> % =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=
=3D=3D=3D=3D=3D=3D
> % The objective function will be vaguely like this:
> function resids =3D objfun(params,px,py,img);
> covchol =3D [params(5),0;params(6),params(7)];
> temp =3D [px(:) - params(2),py(:)-params(3)]*covchol;
> pred =3D params(1) + params(4)*exp(-sum(temp.*temp,2)/2);
> resids =3D pred - img(:);
>
> The above code should be reasonably close,
> although I predict roughly one typo since
> I've not tested the code.
>
> I you don't have the optimization toolbox
> to run lsqnonlin, then get it.
>
> John- Hide quoted text -
>
> - Show quoted text -

John:
You are correct. The method I gave does not give the absolutely best
and optimal fit. Your way is better. But it does seem more
complicated and the simple method does, in practice, often work fairly
well. While not optimal it's often OK for something quick and dirty.
But your way would be better for utmost accuracy.
Regards,
ImageAnalyst

"Steffen " <rileksn@gmail.com> wrote in message <fqj32l$m4n$1@fred.mathworks.com>...
> HI there,
>
> I?d like to fit a gaussian to a bead image representing the
> PSF of my microscope. The ultimate goal is to define the
> 'exact' position of the peak in x and y. I need to fit a
> gaussian to the image rather than just calculating the
> center of mass via centroids.
> Could someone give me a hint how to do so with an image, are
> there eventually example codes available (haven?t yet found
> a proper one)?
>
> Any pointers are much appreciated!
>
> Thanks and best regards,
> Steffen

Hello
If you manage to make it work. Please send me the final code.

Thanks

Logarithm-to-linear least squares methods can work quite well for determining the centroid and sigma of a 2-D Gaussian, including images where part of the Gaussian is masked out, when the Gaussian is floating in a sea of noise, e.g. an image of a star.

The key is to remove the noise floor, leaving only pixels that are appreciably signal. CCD noise is a complicated mixture of things, and a science unto itself, but in the end CCD noise - e.g. readout, shot, dark, and sky brightness - looks pretty much like a Gaussian sea - a plane with pixels tossed randomly above and below it, narrowly bounded. Starting near the top of that sea, simply detect singleton pixels with the convolution matrix

[1 1 1
 1 -1 1
 1 1 1]

and raise the noise floor until no more singleton pixels remain. You can also substitute -2 or -3 in the matrix to detect double and triple shot pixels, although these can be connected corner or edgewise to other pixels for this 3x3 detection matrix.

What you are left with is the just the data, floating on zeros. If your Gaussian is cut by a slit, for example, just mask this out to zeros. Now you can safely form the least squares matrix and solve for the centroid and sigma. You can also constrain the sigma as a parameter to get just the centroid.

Problems: if you have stray pixels in your image that didn't get snipped out, obviously these are going to bias you results. Remember, you're fitting a paraboloid, and stray points are going to pull hard. If A^T * A is not invertible, you'll know it. And if the sigma is not positive, obviously the result is bad. You can also sanity check your result to make sure it's within the data pixels - a centroid that lands in the zeros surrounding your image isn't correct.

A few dozen lines of Matlab.

There are non-linear least squares methods for Gaussian data. They are recursive, complicated, and sensitive to the initial guess. Have fun. Leave this to the nuts doing aperture photometry.

As the other John (I think) mentioned, logarithms turn the scale factor of your 2-D Gaussian into a vertical shift of the paraboloid, putting all pixels in your data on even footing. This can be a good thing - why should I trust bright pixels more than dim ones? Carefully removing the noise floor is the key to good results with least squares.

If you can't remove the noise floor, or your data don't fit this model, then go with quasi-statistical methods like 2-D convolution via 2-D FFT to find the centroid. TMW has some nice example of that in their image toolbox.

Steffen,

Did you finally implement this? I need to do exactly this. Also, I would like to create an image of higher resolution than the original using the parameters obtained from the gaussian fit.

Thanks,

Sanjay

"Steffen " <rileksn@gmail.com> wrote in message <fqj32l$m4n$1@fred.mathworks.com>...
> HI there,
>
> I?d like to fit a gaussian to a bead image representing the
> PSF of my microscope. The ultimate goal is to define the
> 'exact' position of the peak in x and y. I need to fit a
> gaussian to the image rather than just calculating the
> center of mass via centroids.
> Could someone give me a hint how to do so with an image, are
> there eventually example codes available (haven?t yet found
> a proper one)?
>
> Any pointers are much appreciated!
>
> Thanks and best regards,
> Steffen

Steffen,

Did you finally implement this? I need to do exactly this. Also, I would like to create an image of higher resolution than the original using the parameters obtained from the gaussian fit.

Thanks,

Sanjay

"Steffen " <rileksn@gmail.com> wrote in message <fqj32l$m4n$1@fred.mathworks.com>...
> HI there,
>
> I?d like to fit a gaussian to a bead image representing the
> PSF of my microscope. The ultimate goal is to define the
> 'exact' position of the peak in x and y. I need to fit a
> gaussian to the image rather than just calculating the
> center of mass via centroids.
> Could someone give me a hint how to do so with an image, are
> there eventually example codes available (haven?t yet found
> a proper one)?
>
> Any pointers are much appreciated!
>
> Thanks and best regards,
> Steffen

Sanjay, you are aware that this was almost 3 years ago don't you? Do
you think Steffen is still monitoring this newsgroup?

"John D'Errico" <woodchips@rochester.rr.com> wrote in message <fqm7ao$k6v$1@fred.mathworks.com>...
> ImageAnalyst <imageanalyst@mailinator.com> wrote in message
> <e18c27e3-7cc1-4805-863c-
> eae148e22e6f@34g2000hsz.googlegroups.com>...
> > On Mar 4, 4:00=A0am, "Steffen " <rile...@gmail.com> wrote:
> > > HI there,
> > >
> > > I?d like to fit a gaussian to a bead image representing the
> > > PSF of my microscope. The ultimate goal is to define the
> > > 'exact' position of the peak in x and y. I need to fit a
> > > gaussian to the image rather than just calculating the
> > > center of mass via centroids.
> > > Could someone give me a hint how to do so with an image, are
> > > there eventually example codes available (haven?t yet found
> > > a proper one)?
> > >
> > > Any pointers are much appreciated!
> > >
> > > Thanks and best regards,
> > > Steffen
> >
> > How about taking the log of the image and then fitting the resulting
> > image to a quadratic?
> > image =3D a*exp( b*(x-x0)^2)
> > log(image/a) =3D b*(x-x0)^2
> > 0 =3D c*x2 + d*x + e - log(image)
> > Do least square to solve for c, d, e, then plug back into your
> > original equation. Do in two dimensions. I'll leave the math details
> > to you but this is a common way of doing it.
> > Regards,
> > ImageAnalyst
>
> Multiple problems with logging the image data
> to fit a quadratic model:
>
> 1. This will tend to accentuate the fit in exactly
> the WRONG areas, i.e., where the peak is not.
> Logging the data will bring up the small numbers.
> In effect, logging the image implicitly assume a
> lognormal error structure in the regression model.
>
> 2. The model for a general gaussian in 2-d
> requires a positive definite covariance matrix,
> something that a linear regression will not be
> able to constrain. In the logged domain, I would
> not be amazed to find that you have just fit a
> hyperbolic surface to your data. Remember that
> once you log the image, it brings up those tails.
> Crap in the corners will suddenly become
> important.
>
> 3. Logging an image will be problematic if you
> have zero pixels.
>
> 4. If you really wanted to fit a gaussian model
> PLUS a constant offset in z, this is impossible
> to do so via linear regression on the logged
> image data.
>
> There are probably some other reasons I've
> skipped over as to why logging the image is
> not the right way to do this.
>
> How would I do it?
>
> I'd use lsqnonlin or lsqcurvefit to fit a model.
> Generate the pixels coordinates for each pixel
> with meshgrid first.
>
> Formulate the 2-d gaussian model using
> a covariance matrix built from a Cholesky
> factorization. You use three parameters to
> define the cholesky factor, not the actual
> covariance matrix. This assures you that
> the covariance matrix is ALWAYS positive
> definite. (Cute, huh?)
>
> Constrain the gaussian parameters in
> logical ways. For example, an amplitude
> parameter must be positive. The gaussian
> mode must lie inside the image boundaries.
>
> Thus, if img is the image array, then your
> code might look vaguely like this, for a
> gaussian peak with a constant offset in z.
>
>
> [ny,nx] = size(img);
> [px,py] = meshgrid(1:nx,1:ny);
>
> % starting values:
> % 1 - constant offset
> % 2 - mode_x
> % 3 - mode_y
> % 4 - amplitude
> % 5,6,7 - cholesky parameters for cov matrix
> params0 = [mean(img(:)),nx/2,ny/2,1,5,5,0];
> LB = [-inf,1,1,0,-inf,-inf-inf];
> UB = [inf,nx,ny,inf,inf,inf,inf];
>
> options = optimset('lsqnonlin');
> options.Display = 'iter';
> params = lsqnonlin(@(P) objfun(P,px,py,img),params0,LB,UB,options);
>
> cov = [params(5),0;params(6),params(7)];
> cov = cov*cov';
>
> % ==============================
> % The objective function will be vaguely like this:
> function resids = objfun(params,px,py,img);
> covchol = [params(5),0;params(6),params(7)];
> temp = [px(:) - params(2),py(:)-params(3)]*covchol;
> pred = params(1) + params(4)*exp(-sum(temp.*temp,2)/2);
> resids = pred - img(:);
>
>
> The above code should be reasonably close,
> although I predict roughly one typo since
> I've not tested the code.
>
> I you don't have the optimization toolbox
> to run lsqnonlin, then get it.
>
> John

Thank you John for this code, it helped me a lot!
Some remarks for future users:
1) You model the inverse of the covariance matrix (which is also symmetric positive definit), not the covariance matrix itself
2) The cholesky decomposition matrix, which you model with the last 3 parameters, is defined to have positive values on its diagonal. So LB should be 0 for parameter 5 and 7.

Hi Rene,

Would you post your final code that works for this problem? Thank you very much!

Best

"Rene Heideklang" wrote in message <jen2is$gjf$1@newscl01ah.mathworks.com>...
> "John D'Errico" <woodchips@rochester.rr.com> wrote in message <fqm7ao$k6v$1@fred.mathworks.com>...
> > ImageAnalyst <imageanalyst@mailinator.com> wrote in message
> > <e18c27e3-7cc1-4805-863c-
> > eae148e22e6f@34g2000hsz.googlegroups.com>...
> > > On Mar 4, 4:00=A0am, "Steffen " <rile...@gmail.com> wrote:
> > > > HI there,
> > > >
> > > > I?d like to fit a gaussian to a bead image representing the
> > > > PSF of my microscope. The ultimate goal is to define the
> > > > 'exact' position of the peak in x and y. I need to fit a
> > > > gaussian to the image rather than just calculating the
> > > > center of mass via centroids.
> > > > Could someone give me a hint how to do so with an image, are
> > > > there eventually example codes available (haven?t yet found
> > > > a proper one)?
> > > >
> > > > Any pointers are much appreciated!
> > > >
> > > > Thanks and best regards,
> > > > Steffen
> > >
> > > How about taking the log of the image and then fitting the resulting
> > > image to a quadratic?
> > > image =3D a*exp( b*(x-x0)^2)
> > > log(image/a) =3D b*(x-x0)^2
> > > 0 =3D c*x2 + d*x + e - log(image)
> > > Do least square to solve for c, d, e, then plug back into your
> > > original equation. Do in two dimensions. I'll leave the math details
> > > to you but this is a common way of doing it.
> > > Regards,
> > > ImageAnalyst
> >
> > Multiple problems with logging the image data
> > to fit a quadratic model:
> >
> > 1. This will tend to accentuate the fit in exactly
> > the WRONG areas, i.e., where the peak is not.
> > Logging the data will bring up the small numbers.
> > In effect, logging the image implicitly assume a
> > lognormal error structure in the regression model.
> >
> > 2. The model for a general gaussian in 2-d
> > requires a positive definite covariance matrix,
> > something that a linear regression will not be
> > able to constrain. In the logged domain, I would
> > not be amazed to find that you have just fit a
> > hyperbolic surface to your data. Remember that
> > once you log the image, it brings up those tails.
> > Crap in the corners will suddenly become
> > important.
> >
> > 3. Logging an image will be problematic if you
> > have zero pixels.
> >
> > 4. If you really wanted to fit a gaussian model
> > PLUS a constant offset in z, this is impossible
> > to do so via linear regression on the logged
> > image data.
> >
> > There are probably some other reasons I've
> > skipped over as to why logging the image is
> > not the right way to do this.
> >
> > How would I do it?
> >
> > I'd use lsqnonlin or lsqcurvefit to fit a model.
> > Generate the pixels coordinates for each pixel
> > with meshgrid first.
> >
> > Formulate the 2-d gaussian model using
> > a covariance matrix built from a Cholesky
> > factorization. You use three parameters to
> > define the cholesky factor, not the actual
> > covariance matrix. This assures you that
> > the covariance matrix is ALWAYS positive
> > definite. (Cute, huh?)
> >
> > Constrain the gaussian parameters in
> > logical ways. For example, an amplitude
> > parameter must be positive. The gaussian
> > mode must lie inside the image boundaries.
> >
> > Thus, if img is the image array, then your
> > code might look vaguely like this, for a
> > gaussian peak with a constant offset in z.
> >
> >
> > [ny,nx] = size(img);
> > [px,py] = meshgrid(1:nx,1:ny);
> >
> > % starting values:
> > % 1 - constant offset
> > % 2 - mode_x
> > % 3 - mode_y
> > % 4 - amplitude
> > % 5,6,7 - cholesky parameters for cov matrix
> > params0 = [mean(img(:)),nx/2,ny/2,1,5,5,0];
> > LB = [-inf,1,1,0,-inf,-inf-inf];
> > UB = [inf,nx,ny,inf,inf,inf,inf];
> >
> > options = optimset('lsqnonlin');
> > options.Display = 'iter';
> > params = lsqnonlin(@(P) objfun(P,px,py,img),params0,LB,UB,options);
> >
> > cov = [params(5),0;params(6),params(7)];
> > cov = cov*cov';
> >
> > % ==============================
> > % The objective function will be vaguely like this:
> > function resids = objfun(params,px,py,img);
> > covchol = [params(5),0;params(6),params(7)];
> > temp = [px(:) - params(2),py(:)-params(3)]*covchol;
> > pred = params(1) + params(4)*exp(-sum(temp.*temp,2)/2);
> > resids = pred - img(:);
> >
> >
> > The above code should be reasonably close,
> > although I predict roughly one typo since
> > I've not tested the code.
> >
> > I you don't have the optimization toolbox
> > to run lsqnonlin, then get it.
> >
> > John
>
> Thank you John for this code, it helped me a lot!
> Some remarks for future users:
> 1) You model the inverse of the covariance matrix (which is also symmetric positive definit), not the covariance matrix itself
> 2) The cholesky decomposition matrix, which you model with the last 3 parameters, is defined to have positive values on its diagonal. So LB should be 0 for parameter 5 and 7.

Hey together,
Hey Ben,

thanks for the remarks, Ben. Can you say a few more words to the model process?
I'm a chemist and have wide-field fluorescence images containing a gaussian shaped background. So far I made the function running. It finds the correct center of the gaussian. But how I come now from the fitted params/covariance matrix to the 2D gaussian function??
Would be happy about every answer!!

Cheers
sven

Alterantive approach (requires curve fitting toolbox).

Let M be an image matrix with the a single gaussian spot (if you have lots of dots to fit gaussians in a single image you can pull them out by using bwlabel and regionprops Bounding box).

[h,w] = size(M);
    
    [X,Y] = meshgrid(1:h,1:w);
    X = X(:); Y=Y(:); Z = M(:);
    figure(1); clf; scatter3(X,Y,Z);

% 2D gaussian fit object
    gauss2 = fittype( @(a1, sigmax, sigmay, x0,y0, x, y) a1*exp(-(x-x0).^2/(2*sigmax^2)-(y-y0).^2/(2*sigmay^2)),...
        'independent', {'x', 'y'},'dependent', 'z' );

    a1 = max(M(:)); % height, determine from image. may want to subtract background
    sigmax = 2; % guess width
    sigmay = 2; % guess width
    x0 = w/2; % guess position (center seems a good place to start)
    y0 = h/2;
    
% compute fit
    sf = fit([X,Y],double(Z),gauss2,'StartPoint',[a1, sigmax, sigmay, x0,y0]);
    figure(6); clf; plot(sf,[X,Y],Z);

sf.x0 and sf.y0 is the center of gaussian.
sf.sigmax etc will get you the other parameters.


"Steffen " <rileksn@gmail.com> wrote in message <fqj32l$m4n$1@fred.mathworks.com>...
> HI there,
>
> I?d like to fit a gaussian to a bead image representing the
> PSF of my microscope. The ultimate goal is to define the
> 'exact' position of the peak in x and y. I need to fit a
> gaussian to the image rather than just calculating the
> center of mass via centroids.
> Could someone give me a hint how to do so with an image, are
> there eventually example codes available (haven?t yet found
> a proper one)?
>
> Any pointers are much appreciated!
>
> Thanks and best regards,
> Steffen

Dear Ben and Rene,

I am also interested in the final version of the code. Could you post it or send it to me?

Cheers,
Bernie



"Ben " <gogo.xa@gmail.com> wrote in message <jhv1qs$hel$1@newscl01ah.mathworks.com>...
> Hi Rene,
>
> Would you post your final code that works for this problem? Thank you very much!
>
> Best
>
> "Rene Heideklang" wrote in message <jen2is$gjf$1@newscl01ah.mathworks.com>...
> > "John D'Errico" <woodchips@rochester.rr.com> wrote in message <fqm7ao$k6v$1@fred.mathworks.com>...
> > > ImageAnalyst <imageanalyst@mailinator.com> wrote in message
> > > <e18c27e3-7cc1-4805-863c-
> > > eae148e22e6f@34g2000hsz.googlegroups.com>...
> > > > On Mar 4, 4:00=A0am, "Steffen " <rile...@gmail.com> wrote:
> > > > > HI there,
> > > > >
> > > > > I?d like to fit a gaussian to a bead image representing the
> > > > > PSF of my microscope. The ultimate goal is to define the
> > > > > 'exact' position of the peak in x and y. I need to fit a
> > > > > gaussian to the image rather than just calculating the
> > > > > center of mass via centroids.
> > > > > Could someone give me a hint how to do so with an image, are
> > > > > there eventually example codes available (haven?t yet found
> > > > > a proper one)?
> > > > >
> > > > > Any pointers are much appreciated!
> > > > >
> > > > > Thanks and best regards,
> > > > > Steffen
> > > >
> > > > How about taking the log of the image and then fitting the resulting
> > > > image to a quadratic?
> > > > image =3D a*exp( b*(x-x0)^2)
> > > > log(image/a) =3D b*(x-x0)^2
> > > > 0 =3D c*x2 + d*x + e - log(image)
> > > > Do least square to solve for c, d, e, then plug back into your
> > > > original equation. Do in two dimensions. I'll leave the math details
> > > > to you but this is a common way of doing it.
> > > > Regards,
> > > > ImageAnalyst
> > >
> > > Multiple problems with logging the image data
> > > to fit a quadratic model:
> > >
> > > 1. This will tend to accentuate the fit in exactly
> > > the WRONG areas, i.e., where the peak is not.
> > > Logging the data will bring up the small numbers.
> > > In effect, logging the image implicitly assume a
> > > lognormal error structure in the regression model.
> > >
> > > 2. The model for a general gaussian in 2-d
> > > requires a positive definite covariance matrix,
> > > something that a linear regression will not be
> > > able to constrain. In the logged domain, I would
> > > not be amazed to find that you have just fit a
> > > hyperbolic surface to your data. Remember that
> > > once you log the image, it brings up those tails.
> > > Crap in the corners will suddenly become
> > > important.
> > >
> > > 3. Logging an image will be problematic if you
> > > have zero pixels.
> > >
> > > 4. If you really wanted to fit a gaussian model
> > > PLUS a constant offset in z, this is impossible
> > > to do so via linear regression on the logged
> > > image data.
> > >
> > > There are probably some other reasons I've
> > > skipped over as to why logging the image is
> > > not the right way to do this.
> > >
> > > How would I do it?
> > >
> > > I'd use lsqnonlin or lsqcurvefit to fit a model.
> > > Generate the pixels coordinates for each pixel
> > > with meshgrid first.
> > >
> > > Formulate the 2-d gaussian model using
> > > a covariance matrix built from a Cholesky
> > > factorization. You use three parameters to
> > > define the cholesky factor, not the actual
> > > covariance matrix. This assures you that
> > > the covariance matrix is ALWAYS positive
> > > definite. (Cute, huh?)
> > >
> > > Constrain the gaussian parameters in
> > > logical ways. For example, an amplitude
> > > parameter must be positive. The gaussian
> > > mode must lie inside the image boundaries.
> > >
> > > Thus, if img is the image array, then your
> > > code might look vaguely like this, for a
> > > gaussian peak with a constant offset in z.
> > >
> > >
> > > [ny,nx] = size(img);
> > > [px,py] = meshgrid(1:nx,1:ny);
> > >
> > > % starting values:
> > > % 1 - constant offset
> > > % 2 - mode_x
> > > % 3 - mode_y
> > > % 4 - amplitude
> > > % 5,6,7 - cholesky parameters for cov matrix
> > > params0 = [mean(img(:)),nx/2,ny/2,1,5,5,0];
> > > LB = [-inf,1,1,0,-inf,-inf-inf];
> > > UB = [inf,nx,ny,inf,inf,inf,inf];
> > >
> > > options = optimset('lsqnonlin');
> > > options.Display = 'iter';
> > > params = lsqnonlin(@(P) objfun(P,px,py,img),params0,LB,UB,options);
> > >
> > > cov = [params(5),0;params(6),params(7)];
> > > cov = cov*cov';
> > >
> > > % ==============================
> > > % The objective function will be vaguely like this:
> > > function resids = objfun(params,px,py,img);
> > > covchol = [params(5),0;params(6),params(7)];
> > > temp = [px(:) - params(2),py(:)-params(3)]*covchol;
> > > pred = params(1) + params(4)*exp(-sum(temp.*temp,2)/2);
> > > resids = pred - img(:);
> > >
> > >
> > > The above code should be reasonably close,
> > > although I predict roughly one typo since
> > > I've not tested the code.
> > >
> > > I you don't have the optimization toolbox
> > > to run lsqnonlin, then get it.
> > >
> > > John
> >
> > Thank you John for this code, it helped me a lot!
> > Some remarks for future users:
> > 1) You model the inverse of the covariance matrix (which is also symmetric positive definit), not the covariance matrix itself
> > 2) The cholesky decomposition matrix, which you model with the last 3 parameters, is defined to have positive values on its diagonal. So LB should be 0 for parameter 5 and 7.

Found the typo. Instead of

      params = lsqnonlin(@(P) objfun(P,px,py,img),params0,LB,UB,options)

type

      P = lsqnonlin(@(params) objfun(params,px,py,img),params0,LB,UB,options)
:D

As I am pretty new to MatLab, could anyone give me a hint on how to visualize/plot the 2D Gaussian fit?

Cheers,
Bernie