Now, decomposition of 3 particles (including the results and timing...)
Doron
J[i_, j_, \[Alpha]_, \[Beta]_, xpos_, ypos_, clip_] :=
N[Clip[(\[Beta] \[Alpha]^2*\[Pi])/
4*((Erf[(j - xpos)/\[Alpha] ] -
Erf[(j + 1 - xpos)/\[Alpha] ])*(Erf[(i - ypos)/\[Alpha] ] -
Erf[(i + 1 - ypos)/\[Alpha] ])), {0, clip}]];
PixelImage[mat_]:=ArrayPlot[Rescale[mat,{0,Max[mat]},{Max[mat],0}],Epilog->{Red,MapIndexed[Text[#1,Reverse[#2-1/2]]&,Reverse[mat],{2}]},Mesh->True];
MyThreshold = 4043;
model1 = Clip[
N[(\[Beta]1 \[Alpha]1^2*\[Pi])/
4*((Erf[(j - xpos1)/\[Alpha]1 ] -
Erf[(j + 1 - xpos1)/\[Alpha]1 ])*(Erf[(i -
ypos1)/\[Alpha]1 ] -
Erf[(i + 1 - ypos1)/\[Alpha]1]))], {0, MyThreshold}];
model2 = Clip[
N[(\[Beta]2 \[Alpha]2^2*\[Pi])/
4*((Erf[(j - xpos2)/\[Alpha]2 ] -
Erf[(j + 1 - xpos2)/\[Alpha]2 ])*(Erf[(i -
ypos2)/\[Alpha]2 ] -
Erf[(i + 1 - ypos2)/\[Alpha]2]))], {0, MyThreshold}];
model3 = Clip[
N[(\[Beta]3 \[Alpha]3^2*\[Pi])/
4*((Erf[(j - xpos3)/\[Alpha]3 ] -
Erf[(j + 1 - xpos3)/\[Alpha]3 ])*(Erf[(i -
ypos3)/\[Alpha]3 ] -
Erf[(i + 1 - ypos3)/\[Alpha]3]))], {0, MyThreshold}];
Model = Clip[(model1 + model2), {0, MyThreshold}];
data1 = Flatten[
Table[{i, j, J[i, j, 0.7, 500, 6.5, 6.5, MyThreshold]}, {i, 1,
14}, {j, 1, 14}], 1];
data2 = Flatten[
Table[{i, j, J[i, j, 0.6, 300, 4.5, 4.2, MyThreshold]}, {i, 1,
14}, {j, 1, 14}], 1];
data3 = Flatten[
Table[{i, j, J[i, j, 0.6, 300, 7.5, 4.2, MyThreshold]}, {i, 1, 14},
{j, 1, 14}], 1];
Data = Transpose[{data1[[All, 1]], data1[[All, 2]],
Clip[data1[[All, 3]] + data2[[All, 3]] + data3[[All, 3]], {0,
MyThreshold}]}];
PixelImage[Round[Partition[Data[[All, 3]], 14]]]
nlm = NonlinearModelFit[Data, {Model,
{0.4 <= \[Alpha]1 <= 0.8, 0 <= \[Beta]1 <= 500, 4 <= xpos1 <= 8,
4 <= ypos1 <= 8,
0.4 <= \[Alpha]2 <= 0.8, 0 <= \[Beta]2 <= 500, 6 <= xpos2 <= 10,
4 <= ypos2 <= 8,
0.4 <= \[Alpha]3 <= 0.8, 0 <= \[Beta]3 <= 500, 4 <= xpos3 <= 10,
4 <= ypos3 <= 8}}, {{\[Alpha]1, 0.6}, {\[Beta]1, 300}, {xpos1,
6.5}, {ypos1, 6.5}, {\[Alpha]2, 0.6}, {\[Beta]2, 300}, {xpos2,
7.5}, {ypos2, 4}, {\[Alpha]3, 0.6}, {\[Beta]3, 300}, {xpos3,
7.5}, {ypos3, 4}}, {i, j}, Method -> "NMinimize"]; // Timing
nlm["BestFitParameters"]
{950.124090, Null}
{\[Alpha]1 -> 0.600768, \[Beta]1 -> 299.514, xpos1 -> 7.49939,
ypos1 -> 4.20022, \[Alpha]2 -> 0.700811, \[Beta]2 -> 499.503,
xpos2 -> 6.49934,
ypos2 -> 6.49934, \[Alpha]3 -> 0.644292, \[Beta]3 -> 115.722,
xpos3 -> 10., ypos3 -> 4.}
PixelImage[Round[Partition[nlm["FitResiduals"], 14]]]