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