Drexel dragonThe Math ForumDonate to the Math Forum

Search All of the Math Forum:

Views expressed in these public forums are not endorsed by NCTM or The Math Forum.

Math Forum » Discussions » Software » comp.soft-sys.matlab

Topic: 2D gaussian fit of image
Replies: 14   Last Post: Apr 19, 2012 4:56 PM

Advanced Search

Back to Topic List Back to Topic List Jump to Tree View Jump to Tree View   Messages: [ Previous | Next ]
John D'Errico

Posts: 9,127
Registered: 12/7/04
Re: 2D gaussian fit of image
Posted: Mar 5, 2008 8:31 AM
  Click to see the message monospaced in plain text Plain Text   Click to reply to this topic Reply

ImageAnalyst <imageanalyst@mailinator.com> wrote in message
> 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

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.


Point your RSS reader here for a feed of the latest messages in this topic.

[Privacy Policy] [Terms of Use]

© The Math Forum at NCTM 1994-2016. All Rights Reserved.