You mentioned in a [related question](http://mathematica.stackexchange.com/q/27795/862) that you would like to do large numbers of these fits quickly.

If each data set has the same dimensions, you can write a fairly fast implementation of Rahul Narain's method by precomputing arrays of x and y coordinates for the data grid, and flattening the data and using `Dot` to calculate the mean and the elements of the covariance matrix:

    x = Table[i, {i, 9}, {j, 9}]//N;
    y = Transpose[x];
    x = Flatten[x]; y = Flatten[y];

    semiaxes[data_] := Module[{min, p, mx, my},
      min = Min[data];
      p = Flatten[data] - min;
      p /= Total[p];
      mx = x.p;
      my = y.p;
      With[{a = (x - mx)^2.p, b = ((x - mx) (y - my)).p, c = (y - my)^2.p},
       Sqrt @ Eigenvalues[{{a, b}, {b, c}}]]]

    semiaxes[data]
    (* {1.86325, 1.50567} *)

This runs in about 340 µs on my PC

Compiling can give you even more speed, but you need to replace `Eigenvalues` with the explicit symbolic expressions:

    semiaxesc = 
      With[{x = x, y = y}, 
       Compile[{{data, _Real, 2}}, Block[{min, p, mx, my, a, b, c},
         min = Min[data];
         p = Flatten[data] - min;
         p /= Total[p];
         mx = x.p;
         my = y.p;
         a = (x - mx)^2.p;
         b = ((x - mx) (y - my)).p;
         c = (y - my)^2.p;
          {Sqrt[1/2 (a + c - Sqrt[a^2 + 4 b^2 - 2 a c + c^2])],
           Sqrt[1/2 (a + c + Sqrt[a^2 + 4 b^2 - 2 a c + c^2])]}], 
        CompilationTarget -> "C", RuntimeOptions -> "Speed"]];
    
    semiaxesc[data]
    (* {1.50567, 1.86325} *)

This runs in about 5.7 µs on my PC.