I'm having trouble using Newton iterations to fit a 2D gaussian function. Does anyone have experience with this?

The overall goal of what I'm trying to do is locate the center of a gaussian intensity profile, like a dot in an image. (I expect the dots to have a roughly gaussian profile.) I have successfully done this using lsqcurvefit, and the profiles it gives me are pretty good, but it's much slower than the method I'd like to use and I'm having trouble getting my own code to work.

The method I'm following is explained in this paper: http://www.opticsinfobase.org/ao/abstract.cfm?uri=ao-47-36-6842

I think you should be able to see it if you're viewing it from a university, but I'm sorry if you can't access it. This one is open to everyone, but is only for 1D gaussians: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2464285/ This one also goes over the method itself. The short version is this (quoted from the paper):

Locating a solution via Newton’s method involves an iterative approach with the following procedure [8] (where a superscript indicates the iteration index):

    Obtain initial guesses for the parameters θ(0) = (A(0), x̄ (0), w(0))
    Calculate the gradient vector ∇ℓ(k)
    Calculate the Hessian matrix H(k)
    Check the algorithm stopping conditions
    Update the parameter estimate with the equation

θ(k+1)=θk−(H(k))^−1∇ℓ(k).

Except that for the 2D case, you have several more parameters in θ. The first paper I linked defines them all. I'm aiming for the most general case: a 2D gaussian that is wider in one direction and is rotated so it doesn't line up with the x and y axis. My stopping condition is that the summed squared errors needs to be smaller than a certain threshold. I have attached my code as it stands and an image to test it with. The image is a perfectly symmetrical gaussian, which I chose because I know my initial estimate will be almost exactly correct, making it (supposedly) quite easy to converge. I know it's a long shot, but I'm wondering if anyone else has ever done gaussian estimation like this, and if so, can you help me spot my error?

You can call my code by putting everything in the same folder and in the Matlab window, typing:

pic = double(imread('gaussian16.png'));
peak = gaussfit3(pic);

Of course, it doesn't work. You should see several lines of output get spit out before warnings that a matrix (my Hessian matrix) is nearly singular, followed by failure and some surface plots. At every iteration, it tells you how the sum of squared errors has changed. It goes from positive to negative back to positive again (so it's getting better, worse, better, etc.). If it converged correctly, it would all be negative (and small) until the change in the squared errors reached the error tolerance.

Thanks!

Nathan