function [sn,p] = dsn_periodic( s,mu,lambda )
%DSN_PERIODIC Compute a periodic discrete spot noise texture.
%
%   [sn,p] = dsn_periodic(s,mu,lambda) computes the normalized
%   circular discrete spot noise associated to the kernel s, 
%   with intensity lambda, on which we reimpose the mean image mu.
%   The output is of same size as image mu.
%
%   The optional output p allows to save the underlying 
%       Poisson point Process.
%
%   The non-normalized circular spot noise consists in the superposition 
%   of copies of u which are randomly and circularly translated 
%   according to a Poisson process of intensity lambda.
%   The normalization allows to set the expectation of the output random
%   field to mu, and its covariance to the autocorrelation of s.
%
%   Notice that we consider a periodic version of the autocorrelation, but
%   one can obtain the nonperiodic version by embedding s in a twice larger
%   domain by zero-padding.
%
%   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.
%
%   Author : Arthur Leclaire
%   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);

sn = zeros(M,N,C);
p = zeros(M,N);

nb = poissrnd(lambda*M*N);

for k=1:nb
    a = randi(M);
    b = randi(N);
    indx = 1+mod(a:a+m-1,M);
    indy = 1+mod(b:b+n-1,N);
    sn(indx,indy,:) = sn(indx,indy,:) + s;
    % Increment the Poisson point process
    at = a + floor(m/2);
    bt = b + floor(n/2);
    if at>M
        at = at-M;
    end
    if bt>N
        bt = bt-N;
    end
    p(at,bt) = p(at,bt)+1;
end

for c=1:C
    sn(:,:,c)=mu(:,:,c)+(sn(:,:,c)-lambda*sum(sum(s(:,:,c))))/sqrt(lambda);
end

end

