Code covered by the BSD License  

Highlights from
von Mises image entropy

von Mises image entropy

by

 

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)
%%
% [f,kappa,mu,A,phi,error]=vonmisesentropy(X,N,mask,curve,i0,j0,r) finds the 
% 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)
% r: radius of the circular mask (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.
% By Salvador Gabarda 2011
% salvador@optica.csic.es
%%


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

% mask image
[ro co]=size(X);

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

Area=sum(C(:));


X=orlain(X,N/2);

for k=1:8
    ang=(k-1)*dang;
    W=localwigner(X,N,or+ang,'radian');
    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







%%

Contact us