function m = vdgrid2(N,n,sigma,cent)

% function m = vdgrid2(N,n,sigma,[cent])
% N:       Matrix size
% n:       number of samples (must be <= N^2
% sigma:   sigma in DK units
% cent:    radius of center patch to guarantee appear
%
% Output:  sampling pattern as a mask of 0s and 1s
%          desired density is Gaussian with width sigma
%          sampled repeattedly until you get n non-overlapping samples

if n>N^2
    error('n is greater than N^2');
end

if nargin<4 | cent<0
	cent = 0;
end;

if (cent^2 * pi) > n
	error('center patch includes more than n points');
end;

m = zeros(N);
[X Y] = meshgrid(-N/2:N/2-1);
R = X.^2 + Y.^2;
m(R<cent^2) = 1;

dmin = -0.5*(N+1);
xmin = dmin/sigma;
xmax = 0.5*(N-1)/sigma;
dx = xmax-xmin;

amin = erf(xmin);
amax = erf(xmax);
da = amax-amin;

tot = sum(m(:));
while tot<n
    samp = erfinv(amin + rand(1,n-tot)*da);
    x = N/2 + 1 + round(dmin + (samp-xmin)*N/dx);
    samp = erfinv(amin + rand(1,n-tot)*da);
    y = N/2 + 1 + round(dmin + (samp-xmin)*N/dx);
    m = m + sparse(x,y,1,N,N,n-tot);
    m = min(m,1);
    tot = sum(m(:));
end

return;

