HH
Posts:
453
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
> <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.
|
|
