quantize_stereo = {
const ROUND_TO_INT = 2**53; // to round to nearest integer
var precision=2**-prec;
var max_prec = 2 / prev_pow2(precision);
return function quantize_stereo(x, y) {
const sx = Math.sign(x), sy = Math.sign(y);
x *= sx, y *= sy;
const w = 0.5 * (1 + x*x + y*y);
if (w * (precision*precision) > 4) { // near the south pole
var xq = x / (x*x + y*y), yq = y / (x*x + y*y);
if ((xq + yq !== xq + yq) | (xq + yq < 0.5/max_prec))
x = Infinity, y = 0;
else {
x = sx * max_prec * ((yq - xq)*max_prec <= 0.5);
y = sy * max_prec * ((xq - yq)*max_prec <= 0.5);
}
}
else {
const
round = prev_pow2(ROUND_TO_INT * precision * w),
x2 = (x - round) + round,
y2 = (y - round) + round,
// potentially round the original number again using the
// precision from the latitude of the first rounding
w2 = Math.min(w, 0.5 * (1 + x2*x2 + y2*y2)),
round2 = prev_pow2(ROUND_TO_INT * precision * w2);
x = (x - round2) + round2;
y = (y - round2) + round2;
x *= sx, y *= sy; // restore signs
}
return [x, y];
}
}
cursor_to_stereo = ([px, py]) => {
const
[width, height] = size,
back = px >= width,
oy = -zoom * (2/height * py - 1),
ox = zoom * (2/width * (px - back * width) - 1);
if (ox*ox + oy*oy > 1) return [NaN, NaN];
let [x, y] = orthographic ? ortho_to_stereo(ox, oy) : [ox, oy];
if (back) {
[x, y] = [-x/(x*x+y*y), y/(x*x+y*y)];
if (x+y !== x+y) [x, y] = [Infinity, 0];
}
[x, y] = stereo_rotate(rotation, [x, y]);
return [x, y];
}
data1 = new Uint8ClampedArray(2*size[0] * 2*size[1] * 4);
// northern hemisphere
imageData = {
let [width, height] = size;
width *= 2; height *= 2;
for (let py = 0; py < height; py++) {
const oy = -zoom * (2/height * py - 1);
for (let px = 0; px < width; px++) {
const ox = zoom * (2/width * px - 1);
const index = 4 * (py*width + px);
let [x, y] = orthographic ? ortho_to_stereo(ox, oy) : [ox, oy];
[x, y] = stereo_rotate(rotation, [x, y]);
data1.set(colorxy(quantize_stereo(x, y)), index);
}
}
return new ImageData(data1, width, height);
}
data2 = new Uint8ClampedArray(2*size[0] * 2*size[1] * 4);
// southern hemisphere
imageData2 = {
let [width, height] = size;
width *= 2; height *= 2;
for (let py = 0; py < height; py++) {
const oy = -zoom * (2/height * py - 1);
for (let px = 0; px < width; px++) {
const ox = zoom * (2/width * px - 1);
const index = 4 * (py*width + px);
let [x, y] = orthographic ? ortho_to_stereo(ox, oy) : [ox, oy];
[x, y] = [-x/(x*x+y*y), y/(x*x+y*y)];
if (x+y !== x+y) [x, y] = [Infinity, 0];
[x, y] = stereo_rotate(rotation, [x, y]);
data2.set(colorxy(quantize_stereo(x, y)), index);
}
}
return new ImageData(data2, width, height);
}
colorxy = {
const
HI = little_endian, LO = 1 - HI,
EMASK = 0x7ff00000, // 0111 1111 1111 0000 0000 0000 0000 0000
float_array = new Float64Array([1.0, 1.0]),
int_array = new Int32Array(float_array.buffer);
return function colorxy([x, y]) {
if (x + y !== x + y) return [255, 255, 0, 255]; // NaN -> yellow
if ((x + y)*(x + y) === Infinity) return [0, 0, 0, 255]; // Inf -> black
if (x*x + y*y === 0) return [0, 0, 0, 255]; // 0 -> black
float_array[0] = x;
float_array[1] = y;
const
xlo = int_array[0+LO],
ylo = int_array[2+LO],
xhi = int_array[0+HI],
yhi = int_array[2+HI],
xtz = ctz32(xlo) + !xlo * ctz32(xhi | 0x100000),
ytz = ctz32(ylo) + !ylo * ctz32(yhi | 0x100000),
xexp = ((xhi & EMASK) >>> 20) - 1023,
yexp = ((yhi & EMASK) >>> 20) - 1023;
let
x_least_bit = xexp + xtz - 52,
y_least_bit = yexp + ytz - 52;
if (x_least_bit === -1023) x_least_bit = 1023;
if (y_least_bit === -1023) y_least_bit = 1023;
const
least_bit = Math.min(x_least_bit, y_least_bit),
x_tail = ((2**(-least_bit) * x) % 4) & 3,
y_tail = ((2**(-least_bit) * y) % 4) & 3;
// an index between 0 and 11, choosing a color from the list
const index = (
3 * (((x_tail & 2)) + 0.5 * (y_tail & 2)) +
(2 * (x_tail & 1) + (y_tail & 1) - 1));
return grouped_colors[least_bit & 3][index];
}
}
// d3 gives rotation angles anticlockwise; we are using reverse
// rotations here because we start with a position in image space
// and want to find out what place on the sphere it came from
rotation = rotate_d3geo_to_quaternion(-rot_xy, -rot_xz, -rot_yz)
// count trailing zeros
ctz32 = function ctz32(k) {
// fill all bits to the left of the rightmost 1 bit with 1s
k |= k << 16, k |= k << 8, k |= k << 4, k |= k << 2, k |= k << 1;
// subtract count of leading 1s from 32
return 32 - Math.clz32(~k);
}
prev_pow2_x = (x) =>
2**(Math.floor(Math.log2(x)))
// return the largest power of 2 less than or equal to x
prev_pow2 = {
const
HI = little_endian, LO = 1 - HI,
EMASK = 0x7ff00000, TOP = 0xfff00000,
float_array = new Float64Array([1.0]),
int_array = new Int32Array(float_array.buffer);
return function (x) {
float_array[0] = x;
if ((int_array[HI] & EMASK) === EMASK) return x; // inf or nan
if ((int_array[HI] & EMASK) === 0) // 0 or subnormal
return (x === 0) ? x : 2**(Math.floor(Math.log2(x)));
int_array[HI] &= TOP, int_array[LO] = 0; // zero out mantissa
return float_array[0];
}
}
ortho_to_stereo = (X, Y) => {
let Q = 1 + Math.sqrt(1 - X*X - Y*Y);
if (Q !== Q) Q = Math.sqrt(X*X + Y*Y); // clamp to circle
return [X/Q, Y/Q];
}
stereo_rotate = function stereo_rotate(quaternion, stereovector) {
const [x, y] = stereovector, [r1, rxy, ryz, rzx] = quaternion;
let px = r1*x + rxy*y - rzx, py = r1*y - rxy*x + ryz,
q1 = rzx*x - ryz*y + r1, qxy = rzx*y + ryz*x + rxy;
if (x*x + y*y === Infinity) px = r1, py = -rxy, q1 = rzx, qxy = ryz;
if (q1*q1 + qxy*qxy === 0) return [Infinity, 0];
return [(px*q1 + py*qxy) / (q1*q1 + qxy*qxy),
(py*q1 - px*qxy) / (q1*q1 + qxy*qxy)];
}
quatmul = function quatmul(quaternionR, quaternionS) {
const
[r1, rxy, ryz, rzx] = quaternionR,
[s1, sxy, syz, szx] = quaternionS;
// TODO: maybe consider edge cases with very large/small exponents or infinity
return [
r1*s1 - rxy*sxy - ryz*syz - rzx*szx,
r1*sxy + rxy*s1 - ryz*szx + rzx*syz,
r1*syz + rxy*szx + ryz*s1 - rzx*sxy,
r1*szx - rxy*syz + ryz*sxy + rzx*s1];
}
rotate_d3geo_to_quaternion = function rotate_d3geo_to_quaternion(xy_angle, xz_angle, yz_angle){
const halfrad = Math.PI / 360;
xy_angle *= halfrad, yz_angle *= halfrad, xz_angle *= halfrad;
const
R2 = [Math.cos(xy_angle), -Math.sin(xy_angle), 0, 0],
R1 = [Math.cos(xz_angle), 0, 0, Math.sin(xz_angle)],
R0 = [Math.cos(yz_angle), 0, -Math.sin(yz_angle), 0];
return quatmul(R0, quatmul(R1, R2));
}
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.
