function out = rpn( s , mu )
%RPN Compute a Random Phase Noise texture.
%
% out = rpn(s,mu) computes a realization of the random phase noise
% mean mu and whose covariance function is the autocorrelation of s.
%
% Notice that we consider a periodic version of the autocorrelation.
%
% NB :
% - The mean value of s is not substracted.
% - If size(s,1)>M or size(s,2)>N , then s is cropped.
% - The input s can have multiple channels.
% - The RPN associated to a texture u of size [M N] is obtained
% by calling rpn with s = (u-mu)/sqrt(M*N) where mu is the empirical
% mean of u.
%
% This texture model is presented in the paper
% "Random Phase Textures: Theory and Synthesis",
% (B. Galerne, Y. Gousseau, J.-M. Morel),
% IEEE Transactions on Image Processing, 2011.
%
% Author : Arthur Leclaire, Lionel Moisan
% v. 1.0 (03/2014) : first public version
if (size(s,3)~=size(mu,3))
error('s and mu must have the same number of channels')
end
[M,N,C] = size(mu);
if size(s,1)>M
s = s(1:M,:,:);
end
if size(s,2)>N
s = s(:,1:N,:);
end
m = size(s,1);
n = size(s,2);
out = zeros(M,N,C);
% build a random phase function
x = 1:N; y = 1:M;
X = repmat(x,M,1); Y = repmat(y',1,N);
Xp = mod(N+1-X,N)+1; Yp = mod(M+1-Y,M)+1;
i0 = (X-1)*M+Y; ip = (Xp-1)*M+Yp;
A = find(ip>i0);
B = find(ip==i0);
phi = zeros(M,N);
phi(A) = rand(size(A))*2*pi-pi;
phi(ip(A)) = -phi(A);
phi(B) = pi*randi([0,1],size(B));
% apply the random phase shift to each channel of the FFT
for c=1:C
out(:,:,c) = [s(:,:,c) zeros(m,N-n); zeros(M-m,N)];
out(:,:,c) = sqrt(M*N)*real(ifft2(fft2(out(:,:,c)).*exp(1i*phi)));
end
out = mu + out;
end