function [sn,p] = dsn(s,mu,lambda)
%DSN Compute a discrete spot noise texture.
%
% [sn,p] = dsn(s,mu,lambda) computes a piece of the
% normalized discrete spot noise associated to the kernel s,
% with intensity lambda, on which is reimposed the mean image mu.
% The output sn is of same size as mu.
%
% The optional output p allows to save the underlying Poisson point
% Process (more precisely the restriction to the image domain).
%
% The non-normalized spot noise on Z² consists in the superposition
% of copies of u which are randomly 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 the last autocorrelation is automatically non periodic.
%
% NB :
% - The mean value of s is not substracted.
% - The input s can have multiple channels. In that case, mu must
% be a vector containing the mean value of each channel.
%
% Author : Arthur Leclaire
% v. 1.0 (03/2014) : first public version
[m,n,~] = size(s);
[M,N,C] = size(mu);
ix = 1:m;
iy = 1:n;
sn = zeros(M,N,C);
p = zeros(M,N);
nb = poissrnd(lambda*(M+m-1)*(N+n-1));
for k=1:nb
a = randi([-m+1 M-1]);
b = randi([-n+1 N-1]);
indx = a+1:a+m;
condx = ((indx>0)&(indx<=M));
indx = indx(condx);
indy = b+1:b+n;
condy = ((indy>0)&(indy<=N));
indy = indy(condy);
% Increment the spot noise
sn(indx,indy,:) = sn(indx,indy,:) + s(ix(condx),iy(condy),:);
% Increment the Poisson point process
at = a + floor(m/2);
bt = b + floor(n/2);
if (at>0)&&(at<=M)&&(bt>0)&&(bt<=N)
p(at,bt) = p(at,bt)+1;
end
end
for c=1:C
sn(:,:,c)=mu(:,:,c)+(sn(:,:,c)-lambda*sum(sum(s(:,:,c))))/sqrt(lambda);
end
end