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

2d gaussian fittingcurve fittingparticle tracking

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)

Best Answer

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
Related Question