function [F_dir,dir]=txsim_dir(f,t,n)
% TXSIM_COM  Texture directionality measure
%	F_crs=txsim_dir(f,t,n) computes the directionality measure suggested by
%	H. Tamura, S. Mori, T. Yamawaki, ``Textural Features Corresponding
%	to Visual Perception'', IEEE Trans. Sys., Man, Cyb., June 1978
%	with threshold t and n histogram bins. The default value of t is 12
%	and of n is 16.
%	
%	[F_dir,dir]=txsim_dir(f,t,n) also returns the dominating directions.
%	
%	The measure is modified to be low when little of the histogram
%	"integral" is used for computing the moments.

if exist('t')~=1
  t=12;
end
if exist('n')~=1
  n=16;
end

% Compute angle histogram
h_v=[-1 0 1;-1 0 1;-1 0 1];
h_h=[1 1 1;0 0 0;-1 -1 -1];
G=filter2(h_v,f) + sqrt(-1)*filter2(h_h,f);
theta=angle(G) .* (abs(G)>t);
theta=theta .* (theta>0) + (pi+theta) .* (theta<0);
theta=theta(find(theta));
angles=pi*(0:n-1)/(n-1);
H_D=hist(theta,angles)/length(theta);

% Find peaks and valleys
Extrema=filter([-1 1],1,[H_D H_D H_D]);
Extrema=Extrema(n+1:2*n);
Extrema=(Extrema.*hshift(Extrema,-1))<0;
% Find peaks
[dummy,phi]=sort(-H_D .* Extrema);
phi=sort(phi(1:2));
% Find neighboring valleys to the right
[dummy,v(1,2)]=min([Inf*(1:phi(1)) H_D(phi(1)+1:phi(2)-1) Inf*(phi(2):n)]);
[dummy,v(2,1)]=min([H_D(1:phi(1)-1) Inf*(phi(1):phi(2)) H_D(phi(2)+1:n)]);

% Determine whether to use one or two peaks
if    H_D(v(1,2)) / max(H_D(phi)) < 0.5 ...
    & H_D(v(2,1)) / max(H_D(phi)) < 0.5 ...
    & min(H_D(phi)) / max(H_D(phi)) > 0.2
  n_p=2;
else
  n_p=1;
end

% Determine the directionality measure
F_dir=0;
H_D_used=0;
[dummy,p]=sort(-H_D(phi));
for k=1:n_p
  if p(k)==1
    if v(2,1)>phi(2)
      w_p=[v(2,1):n 1:v(1,2)];
      phi_phi=angles([phi(1)+n-v(2,1)+1:-1:1   2:v(1,2)-phi(1)+1]);
    else
      w_p=[v(2,1):v(1,2)];
      phi_phi=angles([phi(1)-v(2,1)+1:-1:1   2:v(1,2)-phi(1)+1]);
    end
  else
    if v(2,1)>phi(2)
      w_p=[v(1,2):v(2,1)];
      phi_phi=angles([phi(2)-v(1,2)+1:-1:1   2:v(2,1)-phi(2)+1]);
    else
      w_p=[v(1,2):n 1:v(2,1)];
      phi_phi=angles([phi(2)-v(1,2)+1:-1:1   2:n-phi(2)+v(2,1)+1]);
    end
  end
  % [w_p;phi_phi.^2 ; H_D(w_p)]
  F_dir=F_dir+sum(phi_phi.^2 .* H_D(w_p));
  H_D_used=H_D_used+sum(H_D(w_p));
end
dir=angles(phi(p(1:n_p)));
F_dir=1 - n_p*F_dir/H_D_used - (sum(H_D)-H_D_used)/sum(H_D);
