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