Error when fitting a sum of 2D gaussians: unable to vary standard deviation.

1 view (last 30 days)
I am trying to fit a sum of 2D gaussians to an image of two connected particles. The equation, at least in principle, is quite simple and involves fitting 11 parameters:
gaussEqn_xy = 'a*exp(-0.5*( ((x-b)/c)^2 + ((y-d)/e)^2) ) + f*exp(-0.5*( ((x-g)/h)^2 + ((y-i)/j)^2) ) + k';
The fit works when I set the standard deviation for the two gaussians to be equal (c==h and e==j, so fitting only 9 parameters in this case). However, if I allow the standard deviation to vary for both gaussians, I get a warning, followed by an error:
"Warning: Too many bounds. Length of upper and lower bounds is greater than the number of coefficients. Ignoring extra bounds."
"Error using fit>iFit (line 304). Too many start points. You need 9 start points for this model."
It's not clear why I am getting the warning or the error, given my number of start points and bounds are correct. Anyone got any suggestions why this is happening?
I have attached an example image and example code. I have copied the fit below, where the error arises.
Thanks.
%% Perform fit
% HIS DOES NOT WORK
xy_startpoints = [100 3 3 -1 3 100 -3 3 1 3 0];
lower_lim = [0, -10, 0, -10, 0, 0, -10, 0, -10, 0, -10];
upper_lim = [255, 10, 10, 10, 10, 255, 10, 10, 10, 10, 10];
gaussEqn_xy = 'a*exp(-0.5*( ((x-b)/c)^2 + ((y-d)/e)^2) ) + f*exp(-0.5*( ((x-g)/h)^2 + ((y-i)/j)^2) ) + k';
% BUT THIS WORKS
% xy_startpoints = [100 3 3 -1 3 100 -3 1 0];
% lower_lim = [0, -10, 0, -10, 0, 0, -10, -10, -10];
% upper_lim = [255, 10, 5, 10, 5, 255, 10, 10, 10];
% gaussEqn_xy = 'a*exp(-0.5*( ((x-b)/c)^2 + ((y-d)/e)^2) ) + f*exp(-0.5*( ((x-g)/c)^2 + ((y-h)/e)^2) ) + k';
fxy = fit([x,y],tmpf,gaussEqn_xy,'StartPoint',xy_startpoints, 'Lower', lower_lim, 'Upper', upper_lim);
fxy_coeffs = coeffvalues(fxy)
  8 Comments
Bjorn Gustavsson
Bjorn Gustavsson on 23 Apr 2020
Since you guys get different behaviour this might be a version-issue, so what versions are you using?

Sign in to comment.

Accepted Answer

Bjorn Gustavsson
Bjorn Gustavsson on 22 Apr 2020
Edited: Bjorn Gustavsson on 22 Apr 2020
You'll get my untested argument: Your free-widths parameters model include i and j and they might both be interpreted as (-1)^(1/2). I would argue that you should use fminsearch (or lsqnonlin) instead of fit. That would require you to explicitly write a function that calculates the sum-of-squared-residuals (or the residuals). That way you'll have explicit control over what goes on. I would also include a rotation-angle for your Gaussians. Something like this:
function err = towGaussianResiduals(pars,x,y,I)
Gaussians = pars(1);
for iP = 2:6:numel(pars)
Imax = pars(iP);
x0 = pars(iP+1);
sigmaX = pars(iP+2);
y0 = pars(iP+3);
sigmaY = pars(iP+4);
phi = pars(iP+5);
Gaussians = Gaussians + Imax*exp(-((+cos(phi)*(x-x0)+sin(phi)*(y-y0)).^2/sigmaX^2 + ...
(-sin(phi)*(x-x0)+cos(phi)*(y-y0)).^2/sigmaY^2));
end
err = sum(( I(:) - Gaussians(:)).^2);
end
That function you should be able to fit your parameters for with something like:
par0 = [0 I1 x1 sX1 y1 sY1 phi1 I2 x2 sX2 y2 sY2 phi2]; % you'll have to set these
[x2D,y2D] = meshgrid(x,y); % if you have these as 1-D arrays
pars = fminsearch(@(pars) twoGaussianResiduals(pars,x2D,y2D,I),par0);
You might have to check that the fitting has converged, check the help for fminsearch.
HTH
  2 Comments
Cameron Boggon
Cameron Boggon on 23 Apr 2020
Hi Bjorn, thanks for this! The immediate problem with the fitting was, as you said, that i and j are recognised as complex numbers rather than parameters. Thanks for the suggestions on the fit method too, I agree there is more control.

Sign in to comment.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!