Fitting a two-dimensional Gaussian to a set of 2D pixels

Sjoerd's answer applies the power of Mathematica's very general model fitting tools. Here's a more low-tech solution.

If you want to fit a Gaussian distribution to a dataset, you can just find its mean and covariance matrix, and the Gaussian you want is the one with the same parameters. For Gaussians this is actually the optimal fit in the sense of being the maximum likelihood estimator -- for other distributions this may not work equally well, so there you will want NonlinearModelFit.

One wrinkle is that your data doesn't fall off to zero but to something like Min[data] $\approx 0.0387$, but we can easily get rid of that by just subtracting it off. Next, I normalize the array to sum to $1$, so that I can treat it like a discrete probability distribution. (All this really does is allow me to avoid dividing by Total[data, Infinity] every time in the following.)

min = Min[data];
sum = Total[data - min, Infinity];
p = (data - min)/sum;

Now we find the mean and covariance. Mathematica's functions don't seem to work with weighted data, so we'll just roll our own.

{m, n} = Dimensions[p];
mean = Sum[{i, j} p[[i, j]], {i, 1, m}, {j, 1, n}];
cov = Sum[Outer[Times, {i, j} - mean, {i, j} - mean] p[[i, j]], {i, 1, m}, {j, 1, m}];

We can easily get the probability distribution for the Gaussian with this mean and covariance. Of course, if we want to match the original data, we have to rescale and shift it back.

f[i_, j_] := PDF[MultinormalDistribution[mean, cov], {i, j}] // Evaluate;
g[i_, j_] := f[i, j] sum + min;

The match is not too bad, although the data has two humps where the Gaussian has one. You can also compare the fit through a contour plot. (Use the tooltips to compare contour levels.)

Show[ListPlot3D[data, PlotStyle -> None], 
 Plot3D[g[i, j], {j, 1, 9}, {i, 1, 9}, MeshStyle -> None, 
  PlotStyle -> Opacity[0.8]], PlotRange -> All]
Show[ListContourPlot[data, Contours -> Range[0.05, 0.25, 0.025], 
  ContourShading -> None, ContourStyle -> ColorData[1, 1], InterpolationOrder -> 3], 
 ContourPlot[g[i, j], {j, 1, 9}, {i, 1, 9}, Contours -> Range[0.05, 0.25, 0.025], 
  ContourShading -> None, ContourStyle -> ColorData[1, 2]]]

The variances along the principal axes are the eigenvalues of the covariance matrix, and the standard deviations (which I guess you're calling the semi-axes) are their square roots.

Sqrt@Eigenvalues[cov]
(* {1.86325, 1.50567} *)

data3D = Flatten[MapIndexed[{#2[[1]], #2[[2]], #1} &, data, {2}], 1];

fm = NonlinearModelFit[data3D, 
        a E^(-(((-my + y) Cos[b] - (-mx + x) Sin[b])^2/(2 sy^2)) - 
               ((-mx + x) Cos[b] + (-my + y) Sin[b])^2/(2 sx^2)), 
        {{a, 0.1}, {b, 0}, {mx, 4.5}, {my, 4.5}, {sx, 3}, {sy, 3}}, 
        {x, y}
     ]

Show[
   ListPlot3D[data3D, PlotStyle -> None, MeshStyle -> Red],
   Plot3D[fm["BestFit"], {x, 1, 9}, {y, 1, 9}, PlotRange -> All, 
          PlotStyle -> Opacity[0.9], Mesh -> None]
]

fm["BestFitParameters"]

{a -> 0.1871830545, b -> -0.4853901689, mx -> 5.152549499, my -> 5.092511036, sx -> 2.756524919, sy -> 2.072163433}

Starting with your data, I would map out the data something like this:

data3D = Flatten[MapIndexed[{#2[[1]], #2[[2]], #1} &, data, {2}], 1]

and then count it like:

data2 = Flatten[Table[{#[[1]], #[[2]]}, {100 #[[3]]}] & /@ data3D, 1];

and then use this function:

FindDistributionParameters[data2, 
MultinormalDistribution[{a, b}, {{c, d}, {e, f}}]]

to find the means and coveriance of the data. And then I might do other stuff with the distribution, such as create a distribution object using:

edist = EstimatedDistribution[data2, 
MultinormalDistribution[{a, b}, {{c, d}, {e, f}}]]

and then check the distribution against the data, using:

dtest = DistributionFitTest[data2, edist, "HypothesisTestData"]

and create a table of the test results, using:

N[dtest["TestDataTable", All]]

or individually

AndersonDarlingTest[data2, edist, "TestConclusion"]
CramerVonMisesTest[data2, edist, "TestConclusion"]
JarqueBeraALMTest[data2, "TestConclusion"]
MardiaKurtosisTest[data2, "TestConclusion"]
MardiaSkewnessTest[data2, "TestConclusion"]
PearsonChiSquareTest[data2, edist, "TestConclusion"]
ShapiroWilkTest[data2, "TestConclusion"]

Now, I'm not a statistician, and I probably only know enough statistics to be dangerous, and I'm not familar with all the tests listed in the table, but I suspect that one of those small numbers is significant in measuring how close the actual data fits the distribution, i.e. how "Gaussian" the set of points is.

Also I might just plot the data as follows and visually judge the gaussiness of the data using these functions (as per VLC's answer to question 2984)

Show[DensityHistogram[d, {.2}, ColorFunction -> (Opacity[0] &), 
Method -> {"DistributionAxes" -> "Histogram"}], ListPlot[d]]

or these:

GraphicsColumn[
Table[SmoothDensityHistogram[d, ColorFunction -> "DarkRainbow", 
Method -> {"DistributionAxes" -> p}, ImageSize -> 500, 
BaseStyle -> {FontFamily -> "Helvetica"}, 
LabelStyle -> Bold], {p, {True, "SmoothHistogram", "BoxWhisker"}}]]