function [fitresult, gof, fout] = gauss2DRotFit(im,varargin)
%gauss2DRotFit(im,fop,title)
%  Create a 2D Gaussian  fit (distribution might be squeezed and rotated).
%  Input :
%   
%      im  : input image of 2D gaussian
%      fop : fitoptions - a structure with 3 fields : Lower,StartPoint and
%            Upper. Each of the fields is an array of 7 values : 
%             offset/background
%             amplitude of the 2D gaussian
%             centroid X of the 2D gaussian
%             centroid Y of the 2D gaussian
%             angle of rotation for the 2D gaussian
%             width X of the 2D gaussian
%             width Y of the 2D gaussian
%
%  Output:
%      fitresult : a fit object representing the fit.
%      gof : structure with goodness-of fit info.
%
%  Created by:  Steffen B. Petersen (12 Dec 2012)
%  Modified by: Eugeniy E. Mikhailov (Nov 2014)

 

if nargin==0
    errordlg('gauss2DRotFit has to be called with at least 1 argument');
    return;
end
[sx,sy,sz]=size(im);
if sz==1
    im=double(im);
end
if sz==3
    im=rgb2gray(im);
end
[X,Y]=meshgrid(1:sy,1:sx);
 

mystr='Untitled fit';
 
[xData, yData, zData] = prepareSurfaceData( X,Y, im );

% Set up fittype and options.
% Note we will fit the expression for intensity (not field) of the Gaussian beam,
% thus extra 2 in exponents.
ft = fittype( 'a + b*exp(-2*(((x-c1)*cosd(t1)+(y-c2)*sind(t1))/w1)^2 - 2*((-(x-c1)*sind(t1)+(y-c2)*cosd(t1))/w2)^2)', 'independent', {'x', 'y'}, 'dependent', 'z' );


opts = fitoptions( ft );
flagMakeStartPointGuess = true;
if numel(varargin)>0
    s=varargin{1};
    if isfield(s,'Lower')
        opts.Lower=s.Lower;
    end
    if isfield(s,'StartPoint')
        opts.StartPoint=s.StartPoint;
        flagMakeStartPointGuess = false;
    end
    if isfield(s,'Upper')
        opts.Upper=s.Upper;
    end
end
if ( flagMakeStartPointGuess )
    % let's attempt to guess starting point
    a0 = min(zData);
    b0 = max(zData) - a0;
    cx0=mean(xData.*zData)/mean(zData);
    cy0=mean(yData.*zData)/mean(zData);
    t=0;
    % see above note about intensity of Gaussian beam
    w10=sqrt(2* mean((xData-cx0).^2.*zData)/mean(zData) );
    w20=sqrt(2* mean((yData-cy0).^2.*zData)/mean(zData) );
    opts.StartPoint= [a0, b0, cx0, cy0, t, w10, w20]
end
if numel(varargin)==2
    mystr=varargin{2};
end
    
opts.Display = 'Off';
 

% Fit model to data.
[fitresult, gof,fout] = fit( [xData, yData], zData, ft, opts );


opts
gof
fout

imgfit = fitresult(xData, yData);
imgfit = reshape(imgfit, sx, sy);


% Create a figure for the plots.
figure( 1);
imagesc(1:sy, 1:sx, imgfit);
title('Fitted image');
xlabel( 'X' );
ylabel( 'Y' );
zlabel( 'image fitted' );
colormap(jet);

% Plot residuals.
figure(2);
imagesc(1:sy, 1:sx, imgfit-im);
title('Residuals of fitted image');
xlabel( 'X' );
ylabel( 'Y' );
zlabel( 'residuals' );
colormap(jet);