Skip to content

Instantly share code, notes, and snippets.

@anissen
Created September 14, 2016 09:06
Show Gist options
  • Save anissen/e27e1a769be7a46d550fecf7523bc9e2 to your computer and use it in GitHub Desktop.
Save anissen/e27e1a769be7a46d550fecf7523bc9e2 to your computer and use it in GitHub Desktop.
Fast uniform Poisson sampling
// 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));
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment