2D gaussian fit of image
Steffen
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
Steffen
Thanks,
Steffen
ImageAnalyst
How about taking the log of the image and then fitting the resulting
image to a quadratic?
image = a*exp( b*(x-x0)^2)
log(image/a) = b*(x-x0)^2
0 = 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
John D'Errico
<e18c27e3-7cc1-4805-863c-
eae148...@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?
> 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
Steffen
thanks for ur advices. Looks a little more complicated than
i thought it would be, but I hope I can manage it to work.
Cheers,
Steffen
ImageAnalyst
> ImageAnalyst <imageanal...@mailinator.com> wrote in message
>
> <e18c27e3-7cc1-4805-863c-
> eae148e22...@34g2000hsz.googlegroups.com>...
>
> - 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
medical imaging Ibrahim
Hello
If you manage to make it work. Please send me the final code.
Thanks
John Kuehne
The key is to remove the noise floor, leaving only pixels that are appreciably signal. CCD noise is a complicated mixture of things, and a science unto itself, but in the end CCD noise - e.g. readout, shot, dark, and sky brightness - looks pretty much like a Gaussian sea - a plane with pixels tossed randomly above and below it, narrowly bounded. Starting near the top of that sea, simply detect singleton pixels with the convolution matrix
[1 1 1
1 -1 1
1 1 1]
and raise the noise floor until no more singleton pixels remain. You can also substitute -2 or -3 in the matrix to detect double and triple shot pixels, although these can be connected corner or edgewise to other pixels for this 3x3 detection matrix.
What you are left with is the just the data, floating on zeros. If your Gaussian is cut by a slit, for example, just mask this out to zeros. Now you can safely form the least squares matrix and solve for the centroid and sigma. You can also constrain the sigma as a parameter to get just the centroid.
Problems: if you have stray pixels in your image that didn't get snipped out, obviously these are going to bias you results. Remember, you're fitting a paraboloid, and stray points are going to pull hard. If A^T * A is not invertible, you'll know it. And if the sigma is not positive, obviously the result is bad. You can also sanity check your result to make sure it's within the data pixels - a centroid that lands in the zeros surrounding your image isn't correct.
A few dozen lines of Matlab.
There are non-linear least squares methods for Gaussian data. They are recursive, complicated, and sensitive to the initial guess. Have fun. Leave this to the nuts doing aperture photometry.
As the other John (I think) mentioned, logarithms turn the scale factor of your 2-D Gaussian into a vertical shift of the paraboloid, putting all pixels in your data on even footing. This can be a good thing - why should I trust bright pixels more than dim ones? Carefully removing the noise floor is the key to good results with least squares.
If you can't remove the noise floor, or your data don't fit this model, then go with quasi-statistical methods like 2-D convolution via 2-D FFT to find the centroid. TMW has some nice example of that in their image toolbox.
Sanjay
Did you finally implement this? I need to do exactly this. Also, I would like to create an image of higher resolution than the original using the parameters obtained from the gaussian fit.
Thanks,
Sanjay
"Steffen " <ril...@gmail.com> wrote in message <fqj32l$m4n$1...@fred.mathworks.com>...
Sanjay
Did you finally implement this? I need to do exactly this. Also, I would like to create an image of higher resolution than the original using the parameters obtained from the gaussian fit.
Thanks,
Sanjay
"Steffen " <ril...@gmail.com> wrote in message <fqj32l$m4n$1...@fred.mathworks.com>...
ImageAnalyst
you think Steffen is still monitoring this newsgroup?
Rene Heideklang
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.
Ben
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>...
Sven
Hey Ben,
thanks for the remarks, Ben. Can you say a few more words to the model process?
I'm a chemist and have wide-field fluorescence images containing a gaussian shaped background. So far I made the function running. It finds the correct center of the gaussian. But how I come now from the fitted params/covariance matrix to the 2D gaussian function??
Would be happy about every answer!!
Cheers
sven
Alistair
Let M be an image matrix with the a single gaussian spot (if you have lots of dots to fit gaussians in a single image you can pull them out by using bwlabel and regionprops Bounding box).
[h,w] = size(M);
[X,Y] = meshgrid(1:h,1:w);
X = X(:); Y=Y(:); Z = M(:);
figure(1); clf; scatter3(X,Y,Z);
% 2D gaussian fit object
gauss2 = fittype( @(a1, sigmax, sigmay, x0,y0, x, y) a1*exp(-(x-x0).^2/(2*sigmax^2)-(y-y0).^2/(2*sigmay^2)),...
'independent', {'x', 'y'},'dependent', 'z' );
a1 = max(M(:)); % height, determine from image. may want to subtract background
sigmax = 2; % guess width
sigmay = 2; % guess width
x0 = w/2; % guess position (center seems a good place to start)
y0 = h/2;
% compute fit
sf = fit([X,Y],double(Z),gauss2,'StartPoint',[a1, sigmax, sigmay, x0,y0]);
figure(6); clf; plot(sf,[X,Y],Z);
sf.x0 and sf.y0 is the center of gaussian.
sf.sigmax etc will get you the other parameters.
Bernhard
I am also interested in the final version of the code. Could you post it or send it to me?
Cheers,
Bernie
"Ben " <gog...@gmail.com> wrote in message <jhv1qs$hel$1...@newscl01ah.mathworks.com>...
Bernhard
params = lsqnonlin(@(P) objfun(P,px,py,img),params0,LB,UB,options)
type
P = lsqnonlin(@(params) objfun(params,px,py,img),params0,LB,UB,options)
:D
As I am pretty new to MatLab, could anyone give me a hint on how to visualize/plot the 2D Gaussian fit?
Cheers,
Bernie
Sozy Fuad
I have the same question about Matlab code for Gaussian PSF do u get any idea about it
thanks and best kinds