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]]]

 J[i_, j_, ?_, ?_, xpos_, ypos_, clip_] := N[Clip[(? ?^2*?)/ 4*((Erf[(j - xpos)/? ] - Erf[(j + 1 - xpos)/? ])*(Erf[(i - ypos)/? ] - Erf[(i + 1 - ypos)/? ])), {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[(?1 ?1^2*?)/4*((Erf[(j - xpos1)/?1 ] - Erf[(j + 1 - xpos1)/?1 ])*(Erf[(i-ypos1)/?1 ] - Erf[(i + 1 - ypos1)/?1]))], {0,MyThreshold}];
 
 model2 = Clip[N[(?2 ?2^2*?)/4*((Erf[(j - xpos2)/?2 ] - Erf[(j + 1 - xpos2)/?2 ])*(Erf[(i -ypos2)/?2 ] - Erf[(i + 1 - ypos2)/?2]))], {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];

Data = Transpose[{data1[[All, 1]], data1[[All, 2]], Clip[data1[[All, 3]] + data2[[All, 3]], {0, MyThreshold}]}];

PixelImage[Round[Partition[Data[[All, 3]], 14]]]



nlm = NonlinearModelFit[Data, {Model,

    {0.4 <= ?1 <= 0.8, 0 <= ?1 <= 6000, 0 <= xpos1 <= 8,

     0 <= ypos1 <= 8, 0.4 <= ?2 <= 0.8, 0 <= ?2 <= 6000,

     0 <= xpos2 <= 8, 0 <= ypos2 <= 8}}, {{?1, 0.6}, {?1,

      300}, {xpos1, 6.5}, {ypos1, 6.5}, {?2, 0.6}, {?2,

     300}, {xpos2, 4.5}, {ypos2, 4}}, {i, j}, Method -> "NMinimize"];

nlm["BestFitParameters"]