Published
Edited
Jun 6, 2022
Importers
19 stars
function SphereDiskSampleNPoints(n) {
return SphereDiskSample(2.949 / Math.sqrt(n));
}
function* SphereDiskSample(radius) {
const k = 11; // maximum number of samples before rejection
const m = 4; // a number mutually prime to k
const eps = 1e-6;
const rotx = Math.cos((2 * Math.PI * m) / k);
const roty = Math.sin((2 * Math.PI * m) / k);
const sradius = Math.tan(0.5 * radius);
const sradius2 = sradius * sradius;
const gridsize = Math.ceil((1 - eps) * Math.PI * (1 + Math.sqrt(2)) / 3 / radius);
const grid_padded_rowsize = 2 * gridsize + 8;
const grid = new Float64Array(8 * (gridsize + 4) * (gridsize + 4)).fill(NaN);
const queue = [];
// pick a random point near the north pole for the first sample
const [Px, Py] = random_northern_point();
yield sample(Px, Py, far(Px, Py));
// pick a random existing sample from the queue.
pick: while (queue.length) {
const i = (Math.random() * queue.length) | 0;
const parent = queue[i];
const t = tanpi_2(2 * Math.random() - 1);
const q = 1 / (1 + t * t);
let dw, dx, dy;
dx = q ? (1 - t * t) * q : -1; // [dx, dy] = random unit vector
dy = q ? 2 * t * q : 0;
// make a new candidate
for (let j = 0; j < k; ++j) {
dw = dx; // temp name for dx
dx = dw * rotx - dy * roty;
dy = dw * roty + dy * rotx;
const r = candidate_radius();
// stereo rotate a point near the north pole to be near the parent
const ax = parent[0], ay = parent[1];
const bx = r * dx, by = r * dy;
const px = ax + bx, py = ay + by;
const qr = 1 - ax*bx - ay*by, qi = ax*by - ay*bx;
const q2 = 1 / (qr*qr + qi*qi);
if (q2 === Infinity) continue; // skip the point at infinity
const x = q2 * (px * qr - py * qi);
const y = q2 * (py * qr + px * qi);
// accept candidates that are far enough from all existing samples
const indices = far(x, y)
if (indices) {
yield sample(x, y, indices);
continue pick;
}
}
// if none of k candidates were accepted, remove the point from the queue
const r = queue.pop();
if (i < queue.length) queue[i] = r;
}
function candidate_radius () {
const rand0 = Math.random();
const rand1 = Math.random();
return sradius * (1 + eps + 0.65 * rand0 * rand1);
}
// fast approximation of tan(πx/2)
function tanpi_2(a) {
let b = 1 - a * a;
return a * (-0.0187108 * b + 0.31583526 + 1.27365776 / b);
}
// random stereographic point near the north pole
function random_northern_point() {
const
t = tanpi_2(2*Math.random() - 1),
u = 0.25*Math.random(),
v = Math.sqrt(u/(1 - u)),
s = 1 / (1 + t*t);
return (s !== 0) ? [(1 - t*t)*(v*s), 2*t*(v*s)] : [-v, 0];
}
function stereo_compare (ax, ay, bx, by, rr) {
let dx = bx - ax, dy = by - ay;
let qr = 1 + ax*bx + ay*by, qi = ax*by - ay*bx;
return (dx*dx + dy*dy) <= rr * (qr*qr + qi*qi)
}
function far(x, y) {
const stereoX = x, stereoY = y;
// COMPUTE SAC QUINCUNCIAL PROJECTION GRID CELL
const sx = (x >= 0) - (x < 0);
const sy = (y >= 0) - (y < 0);
const s = sx * sy;
// reflect into positive quadrant
x *= sx, y *= sy;
const r_2 = 1 / (x*x + y*y);
const rbig = r_2 < 1;
if ((r_2 !== 0) | (r_2 !== Infinity)) {
// invert large values across unit circle
x *= r_2*rbig + !rbig;
y *= r_2*rbig + !rbig;
} else { x = 0, y = 0; }
// project fundamental octant
let xt = x;
x = (4 / Math.PI) * Math.atan2(xt, 1 + y);
y = (4 / Math.PI) * Math.atan2(y, 1 + xt);
// TODO: investigate fast approximations for atan2
// invert across x + y = 1
xt = x;
x = !rbig * x + rbig * (1 - y);
y = !rbig * y + rbig * (1 - xt);
// round to nearest grid cell
x *= gridsize, y *= gridsize;
let xg = x | 0, yg = y | 0; // grid index
const xr = x - xg, yr = y - yg; // remainder
xg -= s * (s*(2*xr + yr) < 1.5*s - 0.5) * (s * (yr - xr) > 0);
yg -= s * (s*(2*yr + xr) < 1.5*s - 0.5) * (s * (yr - xr) <= 0);
// reflect back into negative quadrants
xg = sx * xg - (sx < 0);
yg = sy * yg - (sy < 0);
// CHECK NEAREST 7x7 CELLS FROM THE EXISTING GRID
const j0 = yg + gridsize + 1; // yg - 3 translated and padded
const i0 = xg + gridsize + 1; // xg - 3 translated and padded
for (let j = j0; j < j0 + 7; ++j) {
const index0 = 2 * (j * grid_padded_rowsize + i0);
for (let i = index0; i < index0 + 14; i += 2) {
const gridx = grid[i];
if (gridx !== gridx) continue; // skip NaNs
if (stereo_compare(stereoX, stereoY, gridx, grid[i+1], sradius2)) {
return false; // too close to an existing point
}
}
}
return [xg, yg]; // far enough from every existing point
}
function populate_grid(i, j, x, y) {
const jpad = j + gridsize + 4, ipad = i + gridsize + 4
const index = 2 * (grid_padded_rowsize * jpad + ipad);
grid[index] = x;
grid[index + 1] = y;
}
// flip across the edge either gridsize or -gridsize
function flip(xg) {
const sx = (xg > 0) - (xg <= 0);
return ~(xg - sx*gridsize) + sx*gridsize;
}
function sample(x, y, [xg, yg]) {
const w = gridsize;
const sample_point = [x, y];
populate_grid(xg, yg, x, y);
// also add the sample to grid cells in the padded area
// representing the same part of the sphere, if applicable
const xg_near_side = (xg > w - 4) | (xg <= -w + 4);
const yg_near_side = (yg > w - 4) | (yg <= -w + 4);
if (xg_near_side) {
populate_grid(flip(xg), ~yg, x, y);
if (yg_near_side) populate_grid(~flip(xg), ~flip(yg), x, y);
}
if (yg_near_side) populate_grid(~xg, flip(yg), x, y);
queue.push(sample_point);
return sample_point;
}
}
Purpose-built for displays of data
Observable is your go-to platform for exploring data and creating expressive data visualizations. Use reactive JavaScript notebooks for prototyping and a collaborative canvas for visual data exploration and dashboard creation.
