Code covered by the BSD License

# von Mises image entropy

22 Aug 2012 (Updated )

Algorithm to determine the von Mises distribution of the image Rényi entropy.

vonmisesentropy.m
function [f,kappa,mu,A,phi,error]=vonmisesentropy(X,N,mask,curve,i0,j0,r)
%%
% von Mises distribution of image X.
% Inputs:
% X: The input image as a double precission matrix
% N: number of pixel in the operational window
% mask: 'yes' or 'not' to indicate if a circular mask is used to localize
%       the measure in a desired area.
% curve: 'yes' to plot f, 'no' to do not plot f
% i0: row where the circular mask must be located (optional)
% j0: column where the circular mask must be located (optional)
% Outputs:
% f: the von Mises distribution of the entropy for image X
% kappa: the shape parameter for the von Mises distribution
% mu: the main direction in the von Mises distribution of the image
% A: anisotropy of the image as standard deviation of the image directional
%    entropy
% phi: fitness accuracy of the outcoming distribution. A value of phi close to
%      one indicates that the image fits accurately to the distribution.
%%

or=-pi+pi/8;
dang=pi/4;

[ro co]=size(X);

case 'yes'
if nargin>5
else
C=zeros(ro,co);
i0=fix(ro/2);
j0=fix(co/2);
r=min(i0,j0);
end

case 'no'
C=ones(ro,co);
end

Area=sum(C(:));

X=orlain(X,N/2);

for k=1:8
ang=(k-1)*dang;
R=renyientropy(W);
R=orlaoff(R,N/2);
save 'entropy' 'R'
% global-directional entropy
R=C.*R;
S(k)=sum(R(:))./Area;

% entropy vectors
[u(k),v(k)] = pol2cart(or+ang,S(k));
t(k,:)=[u(k) v(k)];
end

% image anisotropy
A=std(S);

% estimated initial values
% principal direction
[U,Sv,V] = svd(t);
mu=atan(V(2,1)/V(1,1));
if mu<0
mu=mu+pi;
end
emu=mu; % initial estimation of mu

% kappa stimation
tp=t(1:4,:);
T=sum(tp);
RT=norm(T);
RT=RT/4;
kappa=(1/2)/(1-RT);

% fitting
c=0.01;
[f1,kappa1,mu1,ab1,error1,phi1]=fitting(S,emu,kappa,c);
[f,kappa,mu,ab,error,phi]=fitting(S,emu,kappa,-c);

% select best option
if error1<error
f=f1;
kappa=kappa1;
mu=mu1;
ab=ab1;
error=error1;
phi=phi1;
end

switch curve
case 'yes'
figure
%plot(-pi:pi/128:pi-pi/128,f,'-k',-pi+pi/8:pi/4:pi-pi/8,S,'k^','LineWidth',2,...
plot(-180:180/128:180-180/128,f,'-b',-180+180/8:180/4:180-180/8,S,'k^','LineWidth',2,...
'MarkerEdgeColor','k',...
'MarkerFaceColor','r',...
'MarkerSize',10)
xlabel('angle (\theta)')
ylabel('VMD of directional entropy')
case 'no'
end

end

function R=renyientropy(W,alpha)
% normalize pseudo-Wigner distribution
[ro co N]=size(W);
W2=reshape(W,ro*co,N);
P=W2.*conj(W2);
S=sum(P,2);
SS=repmat(S,1,N);
P=P./(SS+eps);

if nargin==1
alpha=3;
end

if alpha==1
% Rnyi entropy = Shannon entropy
Pp=P.*log2(P+eps);
Q=-sum(Pp,2);
else
% Rnyi entropy
P=P.^alpha;
Q=(1/(1-alpha))*log2(sum(P,2)+eps);
end

% round-off error correction
I=find(Q<0);
Q(I)=0;
II=find(Q>log2(N));
Q(II)=0;
U=reshape(Q,ro,co);
R=U./log2(N);
end

function Y=orlain(X,p)
[ro co]=size(X);
Xup=X(1:p,:);
Xup=flipud(Xup);
Xbu=X(ro-p+1:ro,:);
Xbu=flipud(Xbu);
Xp=[Xup;X;Xbu];
Xle=Xp(:,1:p);
Xle=fliplr(Xle);
Xri=Xp(:,co-p+1:co);
Xri=fliplr(Xri);
Y=[Xle Xp Xri];
end

function Y=orlaoff(X,n)
[ro co]=size(X);
Y=X(n+1:ro-n,n+1:co-n);
end

function [f,kappa,mu,ab,error,phi]=fitting(S,mu,kappa,c)

f=zeros(1,256);
fmin=f;
kappamin=kappa;
mumin=mu;
ab=[1 0];
abmin=ab;
error=inf;
errormin=error;
phi=0;
i=0;
while errormin>0
i=i+1;
for k=1:256
ang=-pi+(k-1)*(pi/128);
f(k)=(1/(2*pi*besseli(0,kappa)))*cosh(kappa*cos(ang-mu));
end

% minimum square error fitting
for k=1:8
p=17+(k-1)*32;
g(k)=f(p);
end
G=ones(8,2);
G(:,2)=g;
if rank(G,10^(-3))<2
%disp('ill posed solution')
break

end
ab=polyfit(g,S,1);
f=ab(1,1)*f+ab(1,2);

error=norm(ab-[1 0]);

% minimum error
if error<errormin
errormin=error;
fmin=f;
kappamin=kappa;
mumin=mu;
abmin=ab;

% update mu
II=find(f==max(f(:)),1);
mu=-pi+(II-1)*(pi/128);
if mu>pi
mu=mu-pi;
end
else
break
end

kappa=(1+c)*kappa;

end
%disp(i)
% final result
error=errormin;
f=fmin;
kappa=kappamin;
mu=mumin;
ab=abmin;

% reliability of the measure
phi=exp(-error);

end

%%