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.

W=localwigner(X,N,theta,units,period,window_shape)
```function W=localwigner(X,N,theta,units,period,window_shape)
%%
% W=localwigner(X,N,theta,units,window_shape)
% Pixel-wise Pseudo-Wigner distribution for image X, calculated with a
% N-pixels-1D oriented window
% Inputs:
% X: input grayscale image. Gray-levels must be double precision
% non-negative values in the range 0 to 255 (or normalized to
% the [0, 1] interval)
% N: operational window size (MUST BE AN EVEN NUMBER). The algorithm takes
% an array of N pixels arranged in direction theta.
% theta: a number giving the direction in radians, degrees or slope
% units: according to input theta: 'radian', 'degree' or 'slope'
% period: (optional, when ommited the function is calculated as aperiodic)
%        'periodic'(takes the N+1 pixel value equal to the value in N=1)
%        'aperiodic' (takes the N+1 pixel value equal to the value
%                       determined by the image in position N+1)
% window_shape: (optional, whe ommited the window is considered as square)
%              'square', 'circular'
% Outputs:
% W: pixelwise Pseudo-Wigner Distribution of image X
% W(i,j,:) is the PWD of pixel X(i,j)

%%
X=double(X);
X=X+eps; % avoid zeros
[ro,co]=size(X);
if nargin<6
v=direction(N,theta,units,'square');
else % when 'circular' shape is determined in the input, values are
%calculated by bilinear interpolation and some smoothing effect
% will be expected
v=direction(N,theta,units,window_shape);
end

AI=zeros(ro,co,N);
A=zeros(ro,co,N);

for k=1:N
if k<=N/2+1
iz=v(:,:,k);
[IIx,IIy]=find(iz~=0);
izx=zeros(N+1,1);
izx(IIx,1)=1;
izy=zeros(1,N+1);
izy(1,IIy)=1;
IIIx=(N+2)-IIx;
IIIy=(N+2)-IIy;
dex=zeros(N+1,1);
dex(IIIx,1)=1;
dey=zeros(1,N+1);
dey(1,IIIy)=1;
AI(:,:,k)=conv2(conv2(X,izx,'same'),izy,'same');
if nargin>4
switch period
case 'periodic'
% this bound is required when the function is not real
if k==1
A(:,:,k)=AI(:,:,k).*conj(AI(:,:,k));
end
case 'aperiodic'
end
end
else
A(:,:,k)=A(:,:,N+2-k);
end
end

Y=fft(A,[],3);
W=fftshift(Y,3);
W=real(W);

end

function H=direction(N,theta,units,window_shape)
%%
% Position correction
switch units
case 'degree'
theta=theta-90;
case 'slope'
m=theta;
m=-1/(m+eps);
end

% Generatriz
if abs(m)<1
x=(-N/2:N/2);
y=round(m*x);
h=zeros(N+1,N+1);
x=x+N/2+1;
y=y+N/2+1;
for k=1:N+1
h(x(k),y(k))=1;
end
elseif abs(m)>1
y=(-N/2:N/2);
x=round(m^-1*y);
h=zeros(N+1,N+1);
x=x+N/2+1;
y=y+N/2+1;
for k=1:N+1
h(x(k),y(k))=1;
end
elseif abs(m)==1
h=diag(ones(1,N+1));
end

% Individualisation
H=zeros(N+1,N+1,N+1);
if abs(m)<=1
for k=1:N+1
H(:,:,k)=zeros(N+1);
H(k,:,k)=h(k,:);
end
elseif abs(m)>1
for k=1:N+1
H(:,:,k)=zeros(N+1);
H(:,k,k)=h(:,k);
end
end

switch window_shape
case 'circular'

% Interpolation
U=zeros(N+1,N+1,N+1);

% angle

%for k=1:N+1
for k=1:N+1
% interpolation position
x=(N+2)/2-(k-(N+2)/2)*sin(phi);
y=(N+2)/2+(k-(N+2)/2)*cos(phi);

% boundary positions
x0=fix(x);
x1=fix(x)+1;
y0=fix(y);
y1=fix(y)+1;

% weight factors
fx=x-fix(x);
fy=y-fix(y);
w00=(1-fx)*(1-fy);
w01=(1-fx)*fy;
w10=fx*(1-fy);
w11=fx*fy;

if and(x0<=N+1,y0<=N+1)
U(x0,y0,N+2-k)=w00;
end
if and(x0<=N+1,y1<=N+1)
U(x0,y1,N+2-k)=w01;
end
if and(x1<=N+1,y0<=N+1)
U(x1,y0,N+2-k)=w10;
end
if and(x1<=N+1,y1<=N+1)
U(x1,y1,N+2-k)=w11;
end

end

H=U;

case 'square'

end

end

```