
Re: 2D gaussian fit of image
Posted:
Mar 5, 2008 8:31 AM


ImageAnalyst <imageanalyst@mailinator.com> wrote in message <e18c27e37cc14805863c 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*(xx0)^2) > log(image/a) =3D b*(xx0)^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 2d 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 2d 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,infinf]; 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

