NelderMeadMinimize for 2D Gaussian fitting

Question

I would like to fit a large number (10^5-10^6) of small images to 2D gaussians by taking advantage of the NelderMeadMinimize package from Oleksandr: Shaving the last 50 ms off NMinimize. The application is for analyzing localization microscopy data (fitting lots of single molecule fluorescent emitters in large movies.) There is plenty of good/free software out there for doing this (ImageJ/FIJI plugins), but having an all Mathematica solution with the NedlerMead package would be nice (and bring calculation time to a reasonable amount.)

Using the compiled NelderMead Package for 2D gaussian fitting has been posed previously by blochwave with some suggestions in the comments, but with no resolution. I have been hesitant to re-ask this question since it is certainly a duplicate, but I believe there has been some interest in the past for fast fitting for image analysis (in the comments!) so it would nice to see a working example if possible! (My simpleton attempts are adapted from blochwaves' original question.)

Needs["NelderMeadMinimize`"]

Simple 1D gaussian test with a little noise:

data = Table[a Exp[-(x - b)^2/(2 c^2)] + RandomReal[], {x, 1, 9}] /. {a -> 10, b -> 5, c -> 1};
fitter = With[{x = Range[Length[#]]}, NelderMeadMinimize[Norm[#[[x]] - a Exp[-(x - b)^2/(2 c^2)]], {a, b, c}]] &;
fitter[data]

(*{0.623781, {a -> 10.8638, b -> 4.96825, c -> 1.01535}}*)

Using CompiledNelderMead:

With[{epsilon = $MachineEpsilon, 
   minimizer = 
    NelderMeadMinimize`Dump`CompiledNelderMead[
     Function[{a, b, c}, 
      Block[{x = Range@Length[data]}, 
       Norm[data[[x]] - a Exp[-(x - b)^2/(2 c^2)]]]], {a, b, c}, 
     "ReturnValues" -> "OptimizedParameters"]}, 
  serialFitter = 
    Compile[{{data, _Real, 1}}, 
     minimizer[RandomReal[{0, 1}, {3 + 1, 3}], epsilon, -1], 
     CompilationOptions -> {"InlineCompiledFunctions" -> True}, 
     RuntimeOptions -> {"Speed", "EvaluateSymbolically" -> False}, 
     CompilationTarget -> "C"];];
serialFitter[data]
(*{10.5594, 5.02918, 1.06511}*)

Ok great, compiled 1D gaussian fit works. Make some 2D test data. The previous suggestions were to make sure that all the lists that are fed into NelderMeadMinimize are 1D lists, hence the Flattens.

xyz = Table[{x, y, 
     a Exp[(-(x - b)^2 - (y - c)^2)/(2 d^2)] + RandomReal[]}, {x, 1, 
     9}, {y, 1, 9}] /. {a -> 10, b -> 5, c -> 5, d -> 1};
xvals = Flatten@xyz[[All, All, 1]];
yvals = Flatten@xyz[[All, All, 2]];
zvals = Flatten@xyz[[All, All, 3]];
ArrayPlot[xyz[[All, All, 3]]];

Test 2D gaussian fitting with NelderMeadMinimize:

fitter = Block[{data = #}, 
    NelderMeadMinimize[
     Block[{x = xvals, y = yvals}, 
      Norm[data - a Exp[(-(x - b)^2 - (y - c)^2)/(2 d^2)]]], {a, b, c,
       d}]] &;
fitter[zvals]
(*{4.83445, {a -> 10.0792, b -> 5.03425, c -> 5.00488, d -> 1.11341}}*)

Good, now try compiled 2D gaussian with CompiledNelderMead:

fitter = Compile[{{xvals, _Real, 1}, {yvals, _Real, 1}, {data, _Real, 
     1}}, NelderMeadMinimize`Dump`CompiledNelderMead[
     Function[{a, b, c, d}, 
      Block[{x = xvals, y = yvals}, 
       Norm[Flatten[data] - 
         a Exp[(-(x - b)^2 - (y - c)^2)/(2 d^2)]]]], {a, b, c, d}, 
     "ReturnValues" -> "OptimizedParameters"][
    RandomReal[{0, 1}, {4 + 1, 4}], $MachineEpsilon, -1], 
   CompilationOptions -> {"InlineCompiledFunctions" -> True}, 
   RuntimeOptions -> {"Speed", "EvaluateSymbolically" -> False}, 
   CompilationTarget -> "C"];
fitter[xvals, yvals, zvals]

Correct result, but uncompiled evaluation and the same error about rank 2 tensor that was observed previously:

Perhaps I didn't understand something in the previous suggestion, but I am at a lost for what else to try (and apologies for the duplicated question!).