% Compute color gradients, use p order binomial filter, and threshold t

function [Gx,Gy,m] = ColorGradients(I,p,t);

Ir = BiSmooth(I(:,:,1),p);
Ig = BiSmooth(I(:,:,2),p);
Ib = BiSmooth(I(:,:,3),p);
[rx,ry] = Gradient(Ir, p);
[gx,gy] = Gradient(Ig, p);
[bx,by] = Gradient(Ib, p);

cx = rx+gx+bx;
cy = ry+gy+by;

aa = rx.^2 + gx.^2 + bx.^2;
cc = ry.^2 + gy.^2 + by.^2;

tr = aa + cc;

% eigenvalues cannot be bigger than a+c (trace)
ind = find(tr >= t);

Gx = zeros(size(rx));
Gy = zeros(size(rx));
m = zeros(size(rx));

% NEW Vectorized implementation: 2x2 eigenvalues and eigenvectors computed
% directly

%eigenvalues
a = aa(ind);
c = cc(ind);
b = rx(ind).*ry(ind) + gx(ind).*gy(ind) + bx(ind).*by(ind);

D = sqrt((a+c).^2 - 4*(a.*c-b.^2));
lmax = 0.5*((a+c) + D);
% lmin = 0.5*((a(ind)+c(ind)) - D);

% eigenvectors
ex = zeros(length(ind),1);
ey = zeros(length(ind),1);

% zero out low threshold  edges

% ind1 = find(lmax<t);
% lmax(ind1) = 0;      % ignore those points

ind1 = find(lmax>0);
ex(ind1) = 1;
xx = a - lmax;
ind1 = find((lmax>0) & (abs(xx)>=0.00001));

% compute normalized eigenvectors
x = -b(ind1)./xx(ind1);
xd = 1./sqrt(x.^2+1);
ex(ind1) = x.*xd;
ey(ind1) = xd;
ind1 = find((ex.*cx(ind)+ey.*cy(ind))<0);
ex(ind1) = -ex(ind1);
ey(ind1) = -ey(ind1);
m(ind) = sqrt(lmax);
Gx(ind) = ex.*m(ind);
Gy(ind) = ey.*m(ind);



% 
% OLD NON-VECTORIZED VERSION: eigenvalues and eigenvectors computed using
% eig; may be a bit slow due to the sequential nature of the implementation
% 
% Gx1 = zeros(size(rx));
% Gy1 = zeros(size(rx));
% m1 = zeros(size(rx));
% a = rx.^2 + gx.^2 + bx.^2;
% c = ry.^2 + gy.^2 + by.^2;
% b = rx.*ry + gx.*gy + bx.*by;
% length(ind)
% 
% for i=1:length(ind),
%     S = [a(ind(i)) b(ind(i)); b(ind(i)) c(ind(i))];
%     [V,D] = eig(S);
%     if (D(2,2) > D(1,1))
%         m11 = D(1,1);
%         D(1,1) = D(2,2);
%         D(2,2) = m11;
%         e = V(:,1);
%         V(:,1) = V(:,2);
%         V(:,2) = e;
%     end;
%     m1(ind(i)) = sqrt(D(1,1));
%     e = V(:,1)*m1(ind(i));
%     Gx1(ind(i)) = e(1);
%     Gy1(ind(i)) = e(2);
%     if (cx(ind(i))*Gx1(ind(i))+cy(ind(i))*Gy1(ind(i))<0)
%         Gx1(ind(i)) = -Gx1(ind(i));
%         Gy1(ind(i)) = -Gy1(ind(i));
%     end;
% end;
% 
% Gx = Gx1;
% Gy = Gy1;
% m = m1;