function [fittedImage, vesicles] = fGauss2DFit(img, vesicles, parameters)
% V3. This function fits the peaks whose coordinates are in vesicles
% structure to the 2D-gaussian function. The fit is done in a defined MxN
% area (13x13 usually). It is assumed to be called from parallel pool if
% there are more than 2 arguments, then the third one replaces global
% PARAM. R2 squared of the fit calculated

global PARAM;

% check if called with third argument and assigns parameters accordingly
if nargin > 2
    fitM=parameters.fitM;  
    fitN=parameters.fitN;
else
    fitM=PARAM.fitM;  
    fitN=PARAM.fitN;%Defining the area (a 2FitM*2FitN rectangle) that will be taken out of the image for 2D Gaussian fitting
end

options = optimoptions('lsqcurvefit','Display','none');

fittedImage = img;

for i = 1:length(vesicles)
   rowStart = vesicles(i).m - fitM;
   if (rowStart < 1); rowStart = 1; end
   rowEnd = vesicles(i).m + fitM;
   if (rowEnd > size(img, 1)); rowEnd = size(img, 1); end
   colStart = vesicles(i).n - fitN;
   if (colStart < 1); colStart = 1; end
   colEnd = vesicles(i).n + fitN;
   if (colEnd > size(img, 2)); colEnd = size(img, 2); end
   
   % unfold the image as [x, y, z] matrix for the fitting
   xyzData = []; 
   for iCol = colStart:colEnd
       for iRow = rowStart:rowEnd
           xyzData(end+1,:) = [iRow, iCol, img(iRow, iCol)];
       end
   end
   
   % defining the bounds for fitting [x0, y0, sx, sy, a, b]
   lowerBounds = [rowStart, colStart, 0, 0, 0, 0];
   upperBounds = [rowEnd, colEnd, 3*fitM, 3*fitN, double(2*max(max(img))), double(max(max(img)))];
   
   % defining starting parameters for fitting
   params0 = [vesicles(i).m, vesicles(i).n, 5, 5, 100, 130];
   
   area = img(rowStart:rowEnd, colStart:colEnd);
   
   % fitting
   [fittedParam, resNorm, residual] = lsqcurvefit(@fGauss2D, params0, xyzData(:,1:2), xyzData(:,3), lowerBounds, upperBounds, options);
   
   % writing the resulting fitting parameters into vesicles structure
   vesicles(i).x0 = fittedParam(1);
   vesicles(i).y0 = fittedParam(2);
   vesicles(i).sx = fittedParam(3);
   vesicles(i).sy = fittedParam(4);
   vesicles(i).a = fittedParam(5);
   vesicles(i).b = fittedParam(6);
   vesicles(i).r2 = 1 - resNorm/sum(sum((int32(area)-mean2(area)).^2));
   vesicles(i).pointerror = reshape(abs(residual)./xyzData(:,3), rowEnd-rowStart+1, colEnd-colStart+1);
   
   % reconstruting the portion of the fitted image
   fittedImage(rowStart:rowEnd, colStart:colEnd) = reshape(fGauss2D(fittedParam, xyzData(:,1:2)), rowEnd-rowStart+1, colEnd-colStart+1);
end



end
