Speed up Gaussian Fitting of many points

I have an image where the location of the peaks are identified with red spots

I have a vector of each x and y coordinate of the red spots. I then apply a Gaussia fit to each one (with a span of 6 pixels). My issue is whilst it all works great, I am looping through the coordinates and performing the fit. Its incredibly slow and I was wondering if there was a faster way.

e.g. xf= yf=

%Perform fit to all objects
 for indx=1:l
    ....
    %Now do Gaussian fit      
    [a(indx),b(indx),c(indx),d(indx),xpeak(indx),ypeak(indx)]=myGaussianFit(xdata,ydata, b0,c0);
    % %Increase resolution of x data (by 30)
    xdataFine=(linspace(xdata(1),xdata(end),100))';  
    % Create high res Gauss function
    fitGaus =  a(indx)*exp(-0.5*((xdataFine-b(indx))/c(indx)).^2)+d(indx);
    fwhm(indx) = c(indx) * sqrt(log(256));
    I(indx)=indx;
end

[Thorsten]

Before optimizing, it is helpful to use profile to identify where most of the time is spent.

help profile

And it seems that fitGaus is overwritten in each cycle and not used.

And if would be helpful to post your myGaussianFit and sample data to run the code.

[Jason]

Thanks for your advice. I have run the profiler and its result is indicating its the actual fitting thats taking time:

This is my fitting routine that I call evertime in the loop.

function [a,b,c,d,xpeak,ypeak]=myGaussianFit(xdata,ydata, b0,c0)

%--------------------Gaussian Fit----------------------------------------
%Input is xdata, ydata and must be in column  vector format
%only two initial guesses required as the following are determine by ydata

a0 = max(ydata(:))-min(ydata(:));
d0 = min(ydata(:));

%b0 is the guess at x-location of peak
%c0 is guess at width of peak

%Output gives the fitted parameters as well as xpeak location and its 
%yvalue (i.e. xpeak & ypeak
%-------------------------------------------------------------------------

%Define Gauss Equation (remember the . notation
gaussEqn ='a*exp(-0.5*((x-b)/(c^2)).^2)+d';

%Use startpoint to define guess vlaues (otherwise fit doesn't perform well)
f = fit(xdata,ydata,gaussEqn,'Normalize','off', 'StartPoint',[a0,b0,sqrt(c0),d0]);  %use c^2 instrad of c to enforce postive sigma/fwhm

%Get x location of max value
coeffs=coeffvalues(f);
a=coeffs(1);
b=coeffs(2);
c=coeffs(3);     %need to square it as used c^2 in fitting equation to enforce +ve values
c=c^2;      
d=coeffs(4);

% %Increase resolution of x data (by 30)
xdataFine=(linspace(xdata(1),xdata(end),30))';  
% Create high res Gauss function
fitGaus =  a*exp(-0.5*((xdataFine-b)/c).^2)+d;

%Find max value and its x location
ypeak=max(fitGaus);
xpeak=b;
xpeak=mean(xpeak(:));   %Take mean incase more than one peak found

If I dig deeper, the particular line that is taking up 99.8% of the time is

f = fit(xdata,ydata,gaussEqn,'Normalize','off', 'StartPoint',[a0,b0,sqrt(c0),d0]);  %use c^2 instrad of c to enforce postive sigma/fwhm

Answer 1

[Thorsten]

So the problem is the fitting. There are some alternative methods if you google "matlab fit 2D gaussian", like http://www.mathworks.com/matlabcentral/fileexchange/37087-fit-2d-gaussian-function-to-data . But think about it: Maybe you spend two hours to find a better function, that gives a 10s speedup. In this case you would have to call your function more than 720 times to compensate for the time it has taken to optimize the function.

[Jason]

I was also wondering if my approach using a for loop was the only way to do this i.e. sequentially.

BTW, I only fit along 1 direction.

[Thorsten]

The for loop is not the problem in this case.

[Jason]

Is there an alternative fit routine to use rather than

f =

fit(xdata,ydata,gaussEqn,'Normalize','off', 'StartPoint',[a0,b0,sqrt(c0),d0]);

[Thorsten]

Yes. See my answer above.