octb

function C = octb(n,l)
if (l<1)
    C = ones(n,n,n);
    return;
end
[x y z] = meshgrid(0:n-1,0:n-1,0:n-1);
d1 = x+y+z;
d2 = n-1-x+y+z;
d3 = x+n-1-y+z;
d4 = n-1-x+n-1-y+z;
seps = (n-1)/2/l:(n-1)/l:(n-1);
seps(1:ceil(l/2)) = floor(seps(1:ceil(l/2)));
seps(ceil(l/2+1):l) = ceil(seps(ceil(l/2+1):l));
seps = [seps (2*n-2-seps) (seps+2*n-2)];
b = zeros(n,n,n);
for ii = 1:(2*l)
    s = seps(ii);
    b((d1==s)|(d2==s)|(d3==s)|(d4==s)) = 1;
end
C = b+1;
end