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; } }