|
|
Re: 2D gaussian fit of image
Posted:
Mar 5, 2008 11:18 PM
|
|
On Mar 5, 8:31 am, "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=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- 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
|
|
