Dear Ben and Rene,
I am also interested in the final version of the code. Could you post it or send it to me?
Cheers,
Bernie
"Ben " <gogo.xa@gmail.com> wrote in message <jhv1qs$hel$1@newscl01ah.mathworks.com>...
> Hi Rene,
>
> Would you post your final code that works for this problem? Thank you very much!
>
> Best
>
> "Rene Heideklang" wrote in message <jen2is$gjf$1@newscl01ah.mathworks.com>...
> > "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.
