HH
Posts:
454
Registered:
12/6/04


Re: 2D gaussian fit of image
Posted:
Jan 12, 2012 11:43 AM


"John D'Errico" <woodchips@rochester.rr.com> wrote in message <fqm7ao$k6v$1@fred.mathworks.com>... > 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
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.

