function [ev1, ev2, ori] = localOri(Ix, Iy)

% [ev1, ev2, ori] = localOri(Ix, Iy)
%
% Analyze local orientation, by computing the "orientation tensor"
% (2x2 covariance matrix of the local gradient vectors).  
% [Iy,Ix] are two images providing the gradient measurements.  
%
% [EV1, EV2] are images containing the max and min eigenvalues of the local orientation
% covariance.  
% ORI is an image of local orientations, in radians, clockwise
% from the X-axis.
%
% A perfectly oriented region should produce an EV2 of zero.
%
% A measure of local energy: EV1+EV2.
% A measure of orientedness: (EV1-EV2)/(EV1+EV2) 

% Eero Simoncelli, 1997.  Revised 2003.

% blur filter
bfilt = [1 2 1]/4;

Mxx = conv2(conv2(Ix.*Ix, bfilt, 'valid'), bfilt', 'valid');
Myy = conv2(conv2(Iy.*Iy, bfilt, 'valid'), bfilt', 'valid');
Mxy = conv2(conv2(Ix.*Iy, bfilt, 'valid'), bfilt', 'valid');

% compute eigenvalues of matrix [Mxx Mxy; Mxy Myy]
term1 = (Mxx + Myy)/2;
term2 = sqrt(term1.^2 - (Mxx.*Myy - Mxy.^2));
ev1 = term1 + term2;
ev2 = term1 - term2;

% compute angle of first eigenvector
ori = atan2(Mxx-ev2, Mxy)-pi/2;

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% TEST:

% deriv kernels (from Farid&Simoncelli 97)
gp = [0.036420 0.248972 0.429217 0.248972 0.036420];
gd = [-0.108415 -0.280353 0 0.280353 0.108415];

%% kernels from Farid & Simoncelli, IEEE Trans Im Proc, 13(4):496-508, 2004.
gp = [0.037659 0.249153 0.426375 0.249153 0.037659]';
gd = [-0.109604 -0.276691 0.000000 0.276691 0.109604]';

im = pgmRead('einstein.pgm');
Ix = corrDn(corrDn(im,gd), gp');
Iy = corrDn(corrDn(im,gp), gd');

%% Here, evaluate lines in function body above

%% Check that these determinants are zero everywhere
det = (Mxx - ev1).*(Myy-ev1) - Mxy.*Mxy;  imStats(det);
det = (Mxx - ev2).*(Myy-ev2) - Mxy.*Mxy;  imStats(det);

%% Now try on some examples
o = rand*pi-pi/2;
im = mkSine(64,10,o);
180*[-o, mean2(ori)]/pi

%% more examples:
im=mkImpulse(64);
im=mkDisc(64);
im = zeros(64); im(1:32,1:32)=1;
im = zeros(64); im(1:32,32) = 1;
im = rand(64);

Ix = corrDn(corrDn(im,gd), gp');
Iy = corrDn(corrDn(im,gp), gd');
[ev1,ev2,ori] = localOri(Ix,Iy);
energy = ev1+ev2;
orientedness = power((ev1-ev2)./(ev1+ev2+eps),2);
subplot(2,2,1);  showIm('im');
subplot(2,2,2); showIm('energy');
subplot(2,2,3);  showIm('orientedness');
subplot(2,2,4); showIm(orientedness.*sqrt(energy)); title('product');