Fast uniform Poisson sampling

UniformPoissonSampler.hx

// Adapated from java source by Herman Tulleken
// http://www.luma.co.za/labs/2008/02/27/poisson-disk-sampling/

// The algorithm is from the "Fast Poisson Disk Sampling in Arbitrary Dimensions" paper by Robert Bridson
// http://www.cs.ubc.ca/~rbridson/docs/bridson-siggraph07-poissondisk.pdf

// Code adapted from http://theinstructionlimit.com/fast-uniform-poisson-disk-sampling-in-c

typedef Vector2 = { x :Float, y :Float }

typedef Settings = {
    var TopLeft :Vector2;
    var LowerRight :Vector2;
    var Center :Vector2;
    var Dimensions :Vector2;
    var RejectionSqDistance :Float;
    var MinimumDistance :Float;
    var CellSize :Float;
    var GridWidth :Int;
    var GridHeight :Int;
}

typedef State = {
    var Grid :Array<Array<Vector2>>;
    var ActivePoints :Array<Vector2>;
    var Points :Array<Vector2>;
}
    
    class Test {
        public static function main() {
            //var res = Test.UniformPoissonSampler.SampleCircle({ x: 10, y: 10 }, 10, 1, 30);
            var res = Test.UniformPoissonSampler.SampleRectangle({ x: 0, y: 0 }, { x: 50, y: 50 }, 1, 30);
            //for (r in res) trace('${r.x},${r.y}');
            
            var canvas = js.Browser.document.createCanvasElement();
            canvas.width = 600;
            canvas.height = 600;
            js.Browser.document.body.appendChild(canvas);

            var ctx = canvas.getContext2d();
            for (r in res) {
                ctx.beginPath();
                ctx.fillStyle = '#0000ff';
                ctx.arc(r.x * 25, r.y * 25, 2.0, 0, Math.PI * 2);
                ctx.fill();
            }
        }
        
    }

class UniformPoissonSampler {
    static public var DefaultPointsPerIteration :Int = 30;
    static var SquareRootTwo :Float = Math.sqrt(2);

    public static function SampleCircle(center :Vector2, radius :Float, minimumDistance :Float, pointsPerIteration :Int) :Array<Vector2> {
        return Sample({ x: center.x - radius, y: center.y - radius }, { x: center.x + radius, y: center.y + radius }, radius, minimumDistance, pointsPerIteration);
    }

    public static function SampleRectangle(topLeft :Vector2, lowerRight :Vector2, minimumDistance :Float, pointsPerIteration :Int) :Array<Vector2> {
        return Sample(topLeft, lowerRight, null, minimumDistance, pointsPerIteration);
    }

    static function Sample(topLeft :Vector2, lowerRight :Vector2, rejectionDistance :Float, minimumDistance :Float, pointsPerIteration :Int) :Array<Vector2> {
        var dimensions = { x: lowerRight.x - topLeft.x, y: lowerRight.y - topLeft.y };
        var cellSize = minimumDistance / SquareRootTwo;
        var settings :Settings = {
            TopLeft: topLeft,
            LowerRight: lowerRight,
            Dimensions: dimensions,
            Center: { x: (topLeft.x + lowerRight.x) / 2, y: (topLeft.y + lowerRight.y) / 2 },
            CellSize: cellSize,
            MinimumDistance: minimumDistance,
            RejectionSqDistance: rejectionDistance == null ? null : rejectionDistance * rejectionDistance,
            GridWidth: Std.int((dimensions.x / cellSize) + 1),
            GridHeight: Std.int((dimensions.y / cellSize) + 1)
        };

        var state :State = {
            Grid: [ for (y in 0 ... settings.GridHeight) [ for (x in 0 ... settings.GridWidth) null ] ], // new Vector2?[settings.GridWidth, settings.GridHeight],
            ActivePoints: [],
            Points: []
        };

        AddFirstPoint(settings, state);

        while (state.ActivePoints.length != 0) {
            var ArrayIndex :Int = Math.floor(Math.random() * state.ActivePoints.length);

            var point = state.ActivePoints[ArrayIndex];
            var found = false;

            for (k in 0 ... pointsPerIteration) {
                if (AddNextPoint(point, settings, state)) found = true;
            }

            if (!found) state.ActivePoints.splice(ArrayIndex, 1);
        }

        return state.Points;
    }

    static function DistanceSquared(p1 :Vector2, p2 :Vector2) {
        var dx = p1.x - p2.x;
        var dy = p1.y - p2.y;
        return (dx * dx) + (dy * dy);
    }

    static function Distance(p1 :Vector2, p2 :Vector2) {
        return Math.sqrt(DistanceSquared(p1, p2));
    }

    static function AddFirstPoint(settings :Settings, state :State) {
        var added = false;
        while (!added) {
            var d = Math.random();
            var xr = settings.TopLeft.x + settings.Dimensions.x * d;

            d = Math.random();
            var yr = settings.TopLeft.y + settings.Dimensions.y * d;

            var p = { x: xr, y: yr }; // new Vector2((float) xr, (float) yr);
            if (settings.RejectionSqDistance != null && DistanceSquared(settings.Center, p) > settings.RejectionSqDistance)
                continue;
            added = true;

            var index = Denormalize(p, settings.TopLeft, settings.CellSize);

            state.Grid[Std.int(index.y)][Std.int(index.x)] = p;

            state.ActivePoints.push(p);
            state.Points.push(p);
        } 
    }

    static function AddNextPoint(point :Vector2, settings :Settings, state :State) :Bool {
        var found = false;
        var q = GenerateRandomAround(point, settings.MinimumDistance);

        if (q.x >= settings.TopLeft.x && q.x < settings.LowerRight.x && q.y > settings.TopLeft.y && q.y < settings.LowerRight.y &&
            (settings.RejectionSqDistance == null || DistanceSquared(settings.Center, q) <= settings.RejectionSqDistance)) {
            var qIndex = Denormalize(q, settings.TopLeft, settings.CellSize);
            var tooClose = false;

            for (i in Std.int(Math.max(0, qIndex.x - 2)) ... Std.int(Math.min(settings.GridWidth, qIndex.x + 3))){
                for (j in Std.int(Math.max(0, qIndex.y - 2)) ... Std.int(Math.min(settings.GridHeight, qIndex.y + 3))) {
                    if (state.Grid[j][i] != null && Distance(state.Grid[j][i], q) < settings.MinimumDistance) {
                        tooClose = true;
                        break;
                    }
                }
                if (tooClose) break;
            }

            if (!tooClose) {
                found = true;
                state.ActivePoints.push(q);
                state.Points.push(q);
                state.Grid[Std.int(qIndex.y)][Std.int(qIndex.x)] = q;
            }
        }
        return found;
  }

    static function GenerateRandomAround(center :Vector2, minimumDistance :Float) :Vector2 {
        var d = Math.random();
        var radius = minimumDistance + minimumDistance * d;

        d = Math.random();
        var angle = Math.PI * 2 * d;

        var newX = radius * Math.sin(angle);
        var newY = radius * Math.cos(angle);

        return { x: center.x + newX, y: center.y + newY }; // new Vector2((float) (center.x + newX), (float) (center.y + newY));
    }

    static function Denormalize(point :Vector2, origin :Vector2, cellSize :Float) :Vector2 {
        return { x: (point.x - origin.x) / cellSize, y: (point.y - origin.y) / cellSize }; //new Vector2((int) ((point.x - origin.x) / cellSize), (int) ((point.y - origin.y) / cellSize));
    }
}