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

% [ev1, ev2, ori] = localOri(Ix, Iy)
%
% Analyze local orientation.  
% [IY,IX] are two images providing the gradient measurements.  
% BFILT is the filter used to compute local averages (default='binom5') when 
% computing the covariance. 
% [EV1, EV2] give the max and min eigenvalues of the local orientation
% covariance.  
% ORI specifies 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.

if (exist('bfilt') ~= 1)
   bfilt = namedFilter('binom5');
elseif (ischar(bfilt))
   bfilt = namedFilter(bfilt);
end

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

% 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 filters (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];

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');
