function [kbest, finalKernel,finalFit, resnorm,residual] = fit_harmonicGaussian2D_Kernel_k7_new_up2(init_guess, xdata, ydata, varargin)

%fit_gauss2D_PD(init_guess, xdata, ydata, fixed, lowerbound,upperboud)
% X. Tian, UCLA, 2017, modify to find the kernel of conv/deconv
% Y. Yang, UCLA, 2017
% This function is just after slight modification from
% original version by M. Bartels, UCLA, 2014
% To make the Gaussian always positively defined

%set standard values for optional parameters
fixed=zeros(size(init_guess));
lowerbound= zeros(size(init_guess)) - Inf;
upperbound= zeros(size(init_guess)) + Inf;

%check optional parameters
if nargin>3
   fixed=cell2mat(varargin(1));
end

if nargin>4
   lowerbound=cell2mat(varargin(2));
end

if nargin>5
   upperbound=cell2mat(varargin(3));
end

%convert to double
init_guess    = double(init_guess);
ydata = double(ydata);

%store all initial parameters
init_guess_all = init_guess;

%only fit the variable parameters
init_guess = init_guess_all(fixed == 0);
lowerbound = lowerbound(fixed == 0);
upperbound = upperbound(fixed == 0);


% addpadding
dimx=size(ydata,1);
dimy=size(ydata,2);

pad = 30;
halfPad = pad/2;
upSample = 2;

% xNew = zeros(dimx+pad,dimy+pad)*min(xdata(:)); %zeropadded
% xNew(halfPad+1:halfPad+dimx,halfPad+1:halfPad+dimy)=xdata;

yNew = zeros(dimx+pad,dimy+pad)*min(ydata(:)); %zeropadded
yNew(halfPad+1:halfPad+dimx,halfPad+1:halfPad+dimy)=ydata;

simuFinal = xdata(:,:,1);
mask = My_paddzero(xdata(:,:,2),size(simuFinal)+pad);

%do the fit
opt=optimset('TolFun',1e-12);
opt=optimset(opt,'Display','off');
[kbest, resnorm,residual] = lsqcurvefit(@Fhelp, init_guess, simuFinal, yNew, lowerbound, upperbound, opt);

% %add fixed parameters to the final result again
% init_guess_all(fixed == 0) = kbest;
% kbest = init_guess_all;


PPP = fspecial_XTHarmonicGaussian_new_copy('motion',kbest(1),kbest(2),kbest(4),upSample);
finalKernel = My_paddzero(PPP,[31 31]);
finalFit = Fhelp(kbest,simuFinal);

function y = Fhelp(k, simuFinal)       
    %merge with fixed parameters
    k_current = k;
    k = init_guess_all;
    k(fixed == 0) = k_current;

    %calculate the acutal function, k(2)=0:90,????;
    PSF = fspecial_XTHarmonicGaussian_new_copy('motion',k(1),k(2),k(4),upSample); %length,angle,stdGaussian

    kernel = My_paddzero(PSF,[31 31]);
    y = (k(7)*conv2(simuFinal,kernel)+k(3)).*mask;
%     y = abs(My_FourierShift(ytemp,k(5) ,k(6)));

end

end


