That’s my current implementation, working in 1D/2D/3D. It supports multiple types of norm (argument p_norm).
It does not directly create a WHD image, but just a list of point coordinates (and a list of their corresponding parent points).
My plans is to have something generic that can be used and extended for many filters, like extending it for non-rectangular domains (e.g. masked images), or having the possibility to guide the location of the newly sampled points (e.g. with the contour direction given by another image).
I’m not really sure it works properly in 3D now, I have to check.
Feel free to play with this implementation and fix it if needed 
#@cli sample_poissondisk : width>0,_height>0,_depth>0,_radius>1,_max_sample_attempts>0,_p_norm>0
+sample_poissondisk : check "isint($1) && $1>0 && isint(${2=$1}) && $2>0 && isint(${3=1}) && $3>0 && "\
"${4=20}>1 && isint(${5=30}) && $5>0 && ${6=2}>0"
# Determine optimal list of neighbors to be tested.
W,H,D,R,K=${1-5}
C={$R/sqrt(sum([$W,$H,$D]>1))} # Cell size
0 eval "
dmin(dx,dy,dz,C) = (
_dx = dx*C; _dy = dy*C; _dz = dz*C;
norm$6(min(abs(_dx),abs(_dx + C),abs(_dx - C)),
min(abs(_dy),abs(_dy + C),abs(_dy - C)),
min(abs(_dz),abs(_dz + C),abs(_dz - C)));
);
for (n = 1, n<4, ++n,
is_none_added = 1;
for (k = -n, k<=n, ++k,
for (j = -n, j<=n, ++j,
for (i = -n, i<=n, ++i,
abs(maxabs(i,j,k))==n && dmin(i,j,k,$C)<$R?(da_push([i,j,k]); is_none_added = 0)
)
)
);
is_none_added?break();
);
resize(#-1,1,da_size(),1,3,0)"
f. "($W<2 && i0) || ($H<2 && i1) || ($D<2 && i2)?[8,8,8]:I" discard. 8 s. y,3 a[-3--1] x # Discard useless dimensions
neighbors={^} rm.
# Start point sampling.
{ceil([$W,$H,$D]/$C)},1,-1 nm. grid # Grid used for quick testing of neighboring points
1,32,1,3 nm. points # List of points
1,32,1,1 nm. parents # List of parents
1,32,1,1 nm. active # List of active indices
eval "
insert_active_point(X,ind_parent) = (
i(#$grid,int(X/$C)) = new_ind = da_size(#$points);
da_push(#$points,X);
da_push(#$parents,ind_parent);
da_push(#$active,new_ind);
);
neighbors = ["$neighbors"];
# Initialize list of active points.
X = int(u([$W,$H,$D] - 1e-4));
insert_active_point(X,0);
# Propagate new points.
while (siz = da_size(#$active),
indX = int(u(siz) - 1e-4);
X = I[#$points,i[#$active,indX]];
for (k = 0, k<$K, ++k,
U = [ $W>1?u(-1,1):0, $H>1?u(-1,1):0, $D>1?u(-1,1):0 ];
U*=u($R,1.5*$R)/(1e-9 + norm(U));
Y = round(X + U);
inrange(Y[0],0,$W - 1) && inrange(Y[1],0,$H - 1) && inrange(Y[2],0,$D - 1)?(
P = int(Y/$C);
indY = i(#$grid,P);
indY<0?(
is_sample_ok = 1;
repeat(size(neighbors)/3,n, # Test neighboring points
ind = i(#$grid,P + neighbors[3*n,3]);
ind>=0?(nY = I[#$points,ind]; norm$6(nY - Y)<=$R?(is_sample_ok = 0; break())); # Invalid point
);
is_sample_ok?(insert_active_point(Y,i[#$active,indX]); break()); # Add new point
);
);
);
k>=$K?da_remove(#$active,indX); # Remove point from active list
);
resize(#$points,1,da_size(#$points),1,3,0);
resize(#$parents,1,da_size(#$parents),1,1,0)
"
rm[grid,active]
it can be used like this:
$ gmic tic sample_poissondisk 400,400,1,15 toc 400,400 eval... "i(#-1,i0,i1)=1;I"
or, if you use the parent points to draw connected lines:
$ gmic tic sample_poissondisk 400,400,1,10 toc 400,400 eval... "X=[i0,i1];Y=I[i(#-2)][0,2];polygon(#-1,2,X,Y,1,1);I"
I’m pretty sure a generic version of this Poisson Disk Sampling algorithm could have lot of uses.
