import numpy as np
import itertools as it

m, n = 64, 64

def sample_mask_bernoulli(p):
    return np.random.random((m, n)) < p

shape = "rectangle"
norm = dict(rectangle=np.inf, diamond=1, sphere=2)[shape]

def sample_mask_chunked(p, k):
    didjs = [didj for didj in it.product(range(-k, k + 1), repeat=2)
             if np.linalg.norm(didj, ord=norm) < k]

    area = len(didjs)
    mask = sample_mask_bernoulli(p ** (1./area))

    return np.prod([np.roll(mask, (di, dj), (0, 1)) for di, dj in didjs], axis=0)

%matplotlib inline
import matplotlib.pyplot as plt

def plot_masks(sample_mask):
    k = 2
    fig, axes = plt.subplots(1, k)
    masks = [sample_mask() for _ in range(k)]
    print np.mean(masks)
    for ax, mask in zip(axes, masks):
        ax.imshow(mask, interpolation="nearest", cmap="bone")

plot_masks(lambda: sample_mask_chunked(0.5, 1))

0.49658203125

plot_masks(lambda: sample_mask_chunked(0.5, 2))

0.494018554688

plot_masks(lambda: sample_mask_chunked(0.5, 4))

0.577026367188

plot_masks(lambda: sample_mask_chunked(0.5, 8))

0.494750976562

plot_masks(lambda: sample_mask_chunked(0.5, 12))

0.500732421875