2次元 Poisson disc（点集合）

こんにちは。
2次元 Poisson disc 点集合を生成してみました。ランダムに分布する点の集合ですが、定めた距離以下には互いに接近しない制約付きです。（なお、ここ： Random Points on a Sphere Topology-Preserving では球面上の Poisson disc 点分布も扱っています。）

poissondisc2d.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from __future__ import print_function
import sys
import scipy
import numpy as np
import matplotlib.pyplot as plt

def addrandom(p, scale):
    theta = np.random.uniform(0, 2) * scipy.pi
    r = scipy.sqrt(np.random.uniform(1, 4)) * scale
    return [p[0] + r * scipy.cos(theta), p[1] + r * scipy.sin(theta)]

def dist2(x, y):
    return sum([(a - b)**2 for a, b in zip(x, y)])

def near(p, candid, scale):
    for q in candid:
        if dist2(p, q) < scale**2:
            return True
    return False

def idxpnt(p, scale):
    return [int(np.floor(x/scale/3.0)) for x in p]

def getproximity(points, idx):
    pnts = []
    for k in range(-1,2):
        for l in range(-1,2):
            pnts += points.get((idx[0]+k,idx[1]+l), [])
    return pnts

def pointsappend(points, p, scale):
    idx = tuple(idxpnt(p, scale))
    pnts = points.get(idx, [])
    pnts.append(p)
    points[idx] = pnts
    return 0

def poissondisc2d(scale, ntry, nmax):
    npnt = 1
    candid = [[0.0 for _ in range(2)]]
    points = {}
    points[tuple(idxpnt(candid[0], scale))] = list(candid)
    for _ in range(nmax*2):
        i = np.random.randint(0, len(candid))  
        created = False
        p0 = candid[i]
        pnts = getproximity(points, idxpnt(p0, scale))
        for _ in range(ntry):
            p = addrandom(p0, scale)
            if not near(p, pnts, scale):
                candid.append(p)
                pointsappend(points, p, scale)
                created = True
                break
        if created:
            npnt += 1
            if npnt >= nmax:
                break
        else:
            del candid[i]
    return sum(points.values(), [])  # flatten

def plotpoints(points):
    psiz = 8
    plt.scatter([x[0] for x in points], [x[1] for x in points], s=psiz, c='r', edgecolor='none')
    plt.grid()
    plt.axes().set_aspect(1)
    plt.show()
    return 0

def main():
    scale = 1.0
    ntry = 18
    nmax = int(sys.argv[1]) if len(sys.argv) >= 2 else 1000
    points = poissondisc2d(scale, ntry, nmax)
    plotpoints(points)
    exit(0)

if __name__ == '__main__':
    main()