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.)


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}]] &;

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

Using CompiledNelderMead:

With[{epsilon = $MachineEpsilon, 
   minimizer = 
     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"];];
(*{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 = #}, 
     Block[{x = xvals, y = yvals}, 
      Norm[data - a Exp[(-(x - b)^2 - (y - c)^2)/(2 d^2)]]], {a, b, c,
       d}]] &;
(*{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: enter image description here

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!).

  • 1
    $\begingroup$ Dumb question: why not just take the log of the images, then fit a 2nd order polynomial to them? That way, it's a linear least squares problem, and you can use LeastSquares $\endgroup$ – Niki Estner Feb 1 '16 at 10:10
  • $\begingroup$ Also what are the reasons for avoiding plain NMinimize or FindMinimum? I'm guessing they may come down to speed and correct convergence respectively, but it would be best to have some understanding based on actual testing of the problem class one has in mind. $\endgroup$ – Daniel Lichtblau Feb 1 '16 at 18:23
  • $\begingroup$ @DanielLichtblau - Yes it mostly comes down to speed due to the 50 ms overhead of NMinimze 4700. 10^6 fits * 0.05s / 4cores = 3.5h versus ~20 min for 0.005s per fit . Still not great, but much more reasonable. $\endgroup$ – aaaron Feb 1 '16 at 21:03
  • $\begingroup$ @nikie that seems to be a good idea for data with good S/N (ref1,ref2). I will try some timing tests in the future. $\endgroup$ – aaaron Feb 2 '16 at 8:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.