function [est,yh,estQ,estQr,eigss] = dti_fit( y,grid )
% [est,yh,estQ,estQr,eigss] = dti_fit( y,grid )
% Fits the DTI model to a signal y, measured using directions grid
%
% Inputs:
%  - y, n x 1 vector : signal
%  - grid, 3 x n matrix: measurement directions
%
% Outputs:
%  - est, 3x1: an estimated primary diffusion direction
%  - yh, n x1: predicted signal
%  - estQ, 3x3 : fitted diffusion tensor
%  - estQr, 3x3 : fitted diffusion tensor with constrained eigenvalues lambda2=lambda3
%  - eigss, 3x1 : eigenvalues of tensor
n=size(y,1);
dtiX = dti_design_matrix(grid,ones(n,1));
estimate=dtiX\log(y);
estQ = [estimate(1), estimate(4), estimate(5);...
    estimate(4), estimate(2), estimate(6);...
    estimate(5), estimate(6), estimate(3)];
yh = exp(-squeeze(diag(grid' * estQ * grid)));
[u,v]=eig(estQ);
eigss = sort(v);
[vv,o] = sort(diag(v));
mult = sqrt(vv(3) - (vv(1)+vv(2))/2);
est = u(:,o(3)).*mult;

uu=u(:,o);
vv2 = [(vv(1)+vv(2))/2,(vv(1)+vv(2))/2,vv(3)];
estQr = uu*diag(vv2)*uu';


end

