function z=gaussiennes(x,np);


[n, K] = size(x);

ni = 1 + round((0:np-1)*(n-1)/(np-1));
interp = zeros(n,np);
for i = 1:np-1
  k = ni(i+1) - ni(i);
  interp(ni(i):ni(i+1)-1,i) = (k:-1:1)'/k;
  interp(ni(i)+1:ni(i+1),i+1) = (1:k)'/k;
end



%[rx, p] = sort(x);                      % p temporally stores the indexes
%for k=1:K
%  rx(p(:,k),k) = (1:n)';
%end

p0 = exp(-erfinv((2*ni-1)'/n - 1).^2)/sqrt(2*pi);
%p0 = ones(np,1);
p = interp*p0;

%compute the quantile z from p, the auxilary array lambda and the entropy
z = diff(p);
lamda = z./p(1:n-1);
non0 = abs(lamda) > 1e-5;
z(non0) = log(1+lamda(non0))./(z(non0)*n);
z(~non0) = (1/n)./p(~non0);
entz = sum(log(z))/(n-1);
lamda(non0) = 1./log(lamda(non0)+1) - 1./lamda(non0);
lamda(~non0) = .5;
z = [0; cumsum(z)];
z = z - sum(z)/n;

