"John D'Errico" <firstname.lastname@example.org> wrote in message <email@example.com>... > ImageAnalyst <firstname.lastname@example.org> wrote in message > <e18c27e3-7cc1-4805-863c- > email@example.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.