Poisson Disk Sampler

UniformPoissonDiskSampler.cs

// Adapted 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
using System;
using System.Collections.Generic;
using UnityEngine;
using Random = System.Random;

public static class UniformPoissonDiskSampler
{
    public const int DefaultPointsPerIteration = 30;

    static readonly float SquareRootTwo = Mathf.Sqrt(2);

    struct Settings
    {
        public Vector2 BottomLeft, TopRight, Center;
        public Vector2 Dimensions;
        public float? RejectionSqDistance;
        public Func<Vector2, float> MinimumDistanceSampler;
        public float CellSize;
        public int GridWidth, GridHeight;
    }

    struct State
    {
        public Vector2?[,] Grid;
        public List<Vector2> ActivePoints, Points;
    }

    public static List<Vector2> SampleCircle(Vector2 center, float radius, float minimumDistanceInfBound, Func<Vector2, float> minimumDistanceSampler = null, int pointsPerIteration = DefaultPointsPerIteration)
    {
        return Sample(center - new Vector2(radius, radius), center + new Vector2(radius, radius), radius, minimumDistanceInfBound, minimumDistanceSampler, pointsPerIteration);
    }

    public static List<Vector2> SampleRectangle(Vector2 bottomLeft, Vector2 topRight, float minimumDistanceInfBound, Func<Vector2, float> minimumDistanceSampler = null, int pointsPerIteration = DefaultPointsPerIteration)
    {
        return Sample(bottomLeft, topRight, null, minimumDistanceInfBound, minimumDistanceSampler, pointsPerIteration);
    }

    static List<Vector2> Sample(Vector2 bottomLeft, Vector2 topRight, float? rejectionDistance, float minimumDistanceInfBound, Func<Vector2, float> minimumDistanceSampler, int pointsPerIteration)
    {
        var settings = new Settings
        {
            BottomLeft = bottomLeft,
            TopRight = topRight,
            Dimensions = topRight - bottomLeft,
            Center = (bottomLeft + topRight) / 2,
            CellSize = minimumDistanceInfBound / SquareRootTwo,
            MinimumDistanceSampler = minimumDistanceSampler ?? (v => minimumDistanceInfBound),
            RejectionSqDistance = rejectionDistance == null ? null : rejectionDistance * rejectionDistance
        };
        settings.GridWidth = (int)(settings.Dimensions.x / settings.CellSize) + 1;
        settings.GridHeight = (int)(settings.Dimensions.y / settings.CellSize) + 1;

        var state = new State
        {
            Grid = new Vector2?[settings.GridWidth, settings.GridHeight],
            ActivePoints = new List<Vector2>(),
            Points = new List<Vector2>()
        };

        AddFirstPoint(ref settings, ref state);

        while (state.ActivePoints.Count != 0)
        {
            var listIndex = RandomHelper.Random.Next(state.ActivePoints.Count);

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

            for (var k = 0; k < pointsPerIteration; k++)
                found |= AddNextPoint(point, ref settings, ref state);

            if (!found)
                state.ActivePoints.RemoveAt(listIndex);
        }

        return state.Points;
    }

    static void AddFirstPoint(ref Settings settings, ref State state)
    {
        var added = false;
        while (!added)
        {
            var d = RandomHelper.Random.NextDouble();
            var xr = settings.BottomLeft.x + settings.Dimensions.x * d;

            d = RandomHelper.Random.NextDouble();
            var yr = settings.BottomLeft.y + settings.Dimensions.y * d;

            var p = new Vector2((float)xr, (float)yr);
            if (settings.RejectionSqDistance != null && Vector2.SqrMagnitude(settings.Center - p) > settings.RejectionSqDistance)
                continue;
            added = true;

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

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

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

    static bool AddNextPoint(Vector2 point, ref Settings settings, ref State state)
    {
        var found = false;

        var q = GenerateRandomAround(point, settings.MinimumDistanceSampler(point));

        if (q.x >= settings.BottomLeft.x && q.x < settings.TopRight.x &&
            q.y > settings.BottomLeft.y && q.y < settings.TopRight.y &&
            (settings.RejectionSqDistance == null || Vector2.SqrMagnitude(settings.Center - q) <= settings.RejectionSqDistance))
        {
            var qIndex = Denormalize(q, settings.BottomLeft, settings.CellSize);
            var tooClose = false;

            var localMinDistance = settings.MinimumDistanceSampler(q);
            var nbCell = Mathf.CeilToInt(localMinDistance/settings.CellSize) + 1;
            for (var i = Mathf.Max(0, (int)qIndex.x - nbCell); i < Math.Min(settings.GridWidth, qIndex.x + nbCell + 1) && !tooClose; i++)
                for (var j = Mathf.Max(0, (int)qIndex.y - nbCell); j < Math.Min(settings.GridHeight, qIndex.y + nbCell + 1) && !tooClose; j++)
                    if (state.Grid[i, j].HasValue && Vector2.Distance(state.Grid[i, j].Value, q) < localMinDistance)
                        tooClose = true;

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

    static Vector2 GenerateRandomAround(Vector2 center, float minimumDistance)
    {
        var d = RandomHelper.Random.NextDouble();
        var radius = minimumDistance + minimumDistance * d;

        d = RandomHelper.Random.NextDouble();
        var angle = MathHelper.TwoPi * d;

        var newX = radius * Math.Sin(angle);
        var newY = radius * Math.Cos(angle);

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

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

public static class RandomHelper
{
    public static readonly Random Random = new Random();
}

public static class MathHelper
{
    public const float Pi = Mathf.PI;
    public const float HalfPi = Mathf.PI / 2;
    public const float TwoPi = Mathf.PI * 2;
}