2b:null 2c:["$","title",null,{"children":"Quantizing stereographic coordinates | Jacob Rus | Observable"}] 31:T46f,md`# Quantizing stereographic coordinates To save space when transmitting or storing data, we [quantize][] it (that’s a fancy word for [rounding][]). This notebook contains my proposal for a way to quantize points on the sphere – for instance describing geographical features such as roads or coastlines – which are represented using [stereographic][] coordinates. For more about the stereographic representation of the sphere, see [@jrus/planisphere](/@jrus/planisphere). The maximum quantization error stays roughly constant over the sphere, with differences between nearby points at different places taking about the same number of meaningful [bits][] to represent, leaving trailing irrelevantly precise bits all set to \`0\`. In combination with a “projective number” format (see [@jrus/pj](/@jrus/pj)) this should facilitate aggressive compression of spherical data. [quantize]: https://en.wikipedia.org/wiki/Quantization_(signal_processing) [rounding]: https://en.wikipedia.org/wiki/Rounding [stereographic]: https://en.wikipedia.org/wiki/Stereographic_projection [bits]: https://en.wikipedia.org/wiki/Bit `32:T427,viewof plot = { const [w, h] = size; const context = DOM.context2d(w*2, h, 2); context.canvas.style.cursor = 'crosshair'; d3.select(context.canvas) .on("touchmove", e => e.preventDefault()) // prevent scrolling .on("pointerenter pointerup", (e) => trackPointer(e, { move: (p) => mutable mouseover = cursor_to_stereo(p.point), out: (p) => mutable mouseover = [NaN, NaN] })); // add colored points context.putImageData(imageData, 0, 0); context.putImageData(imageData2, w*2, 0); context.beginPath(); // rectangle enclosing whole canvas context.moveTo(0,0); context.lineTo(2*w,0); context.lineTo(2*w,h); context.lineTo(0,h); context.lineTo(0,0); // subtract circle on the left context.moveTo(w/2 + w/2/zoom, w/2); context.arc(w/2, w/2, w/2/zoom, 0, 2 * Math.PI); // subtract circle on the right context.moveTo(3*w/2 + w/2/zoom, w/2); context.arc(3*w/2, w/2, w/2/zoom, 0, 2 * Math.PI); // fill with white context.fillStyle = 'white'; context.fill('evenodd'); return context.canvas; }33:T497,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]; } }34:T5b2,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]; } }35:[] 2f:{"id":"db3ed197ab069745","slug":"quantizing-stereographic-coordinates","title":"Quantizing stereographic coordinates","thumbnail":"c4c624432a150bf9728b7857c6f8aae634a85dc957131fdcf1e2d845efebcc7f","access_level":"public","create_time":"$D2021-03-23T00:02:48.929Z","update_time":"$D2021-06-21T00:44:43.639Z","view_time":null,"first_public_time":"$D2021-03-26T02:43:45.008Z","fork":false,"stars":13,"starred":false,"owner":{"type":"team","id":"9b8a8e177a49059a","login":"jrus","name":"Jacob Rus","avatar_url":"https://avatars-next.observableusercontent.com/c33dbeed76527c50e64208b0e005e0b6e54a3bef2a14ea6a823f1c7de4867662","tier":"starter_2024"},"version":1359,"publish_version":1359,"role":null,"authored":false,"authors":[{"type":"individual","id":"ddd7235babc79ce0","login":"jrus","name":"Jacob Rus","avatar_url":"https://avatars-next.observableusercontent.com/c33dbeed76527c50e64208b0e005e0b6e54a3bef2a14ea6a823f1c7de4867662","tier":null}],"members":[{"type":"individual","id":"ddd7235babc79ce0","login":"jrus","name":"Jacob Rus","avatar_url":"https://avatars-next.observableusercontent.com/c33dbeed76527c50e64208b0e005e0b6e54a3bef2a14ea6a823f1c7de4867662","tier":null,"role":"editor"}],"comments":[],"custom_thumbnail":true,"trashed":false,"subscribed":null,"forks":8,"collection":{"id":"0e11b15cb878449d","slug":"maps","title":"Maps","owner":{"id":"9b8a8e177a49059a","login":"jrus","name":"Jacob Rus"}},"type":"notebook","body":{"title":"Quantizing stereographic coordinates","cells":[{"id":0,"value":"$31","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":19,"value":"$32","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1122,"value":"html`\n

\nMouse coordinates:\n   (${mouseover.map(print_hexfloat).join(', ')})\n-> (${mouseover_quantized.map(print_hexfloat).join(', ')})\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":332,"value":"viewof orthographic = Toggle({label: \"Orthographic?\", value: true})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":165,"value":"viewof prec = Range([-1, 10], {label: \"Bits of Precision\", value: 4, width: 300})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1014,"value":"viewof rot_xy = Range([-180, 180], {label: \"λ\", value: 0, width: 300, step: 1})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1020,"value":"viewof rot_xz = Range([-180, 180], {label: \"ϕ\", value: 0, width: 300, step: 1})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1023,"value":"viewof rot_yz = Range([-180, 180], {label: \"γ\", value: 0, width: 300, step: 1})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":745,"value":"md`\n\nThe image above shows an orthographic view of the northern hemisphere (on the left) and southern hemisphere (on the right), with points represented via a stereographic projection of the sphere to a plane where the north pole maps to the origin and the south pole maps to a point at infinity.\n\nUntick the ‘orthographic’ box to see the hemispheres under stereographic projection instead of orthographic. Drag the ‘precision’ slider to round values more or less aggressively.\n\nFor practical use storing geographic data, the precision should generally be at least 15 bits (often 25 bits or more), but with such a dense grid it is hard to see what is going on. \n\nNext step: explore new compression algorithms for large contiguous arrays of geographical data quantized and then stored as projective numbers, as the main component of a binary file format used as an alternative to [geojson][] and [topojson][].\n\n  [geojson]: https://geojson.org\n  [topojson]: https://github.com/topojson/topojson-specification\n\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":100,"value":"$33","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1125,"value":"mutable mouseover = [NaN, NaN]","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1130,"value":"mouseover_quantized = quantize_stereo(...mouseover)","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1106,"value":"cursor_to_stereo = ([px, py]) => {\n  const\n    [width, height] = size,\n    back = px >= width,\n    oy = -zoom * (2/height * py - 1),\n    ox = zoom * (2/width * (px - back * width) - 1);\n  if (ox*ox + oy*oy > 1) return [NaN, NaN];\n  let [x, y] = orthographic ? ortho_to_stereo(ox, oy) : [ox, oy];\n  if (back) {\n    [x, y] = [-x/(x*x+y*y), y/(x*x+y*y)];\n    if (x+y !== x+y) [x, y] = [Infinity, 0];\n  }\n  [x, y] = stereo_rotate(rotation, [x, y]);\n  return [x, y];\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1225,"value":"data1 = new Uint8ClampedArray(2*size[0] * 2*size[1] * 4);","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":4,"value":"// northern hemisphere\nimageData = {\n  let [width, height] = size;\n  width *= 2; height *= 2;\n  \n  for (let py = 0; py < height; py++) {\n    const oy = -zoom * (2/height * py - 1);\n    for (let px = 0; px < width; px++) {  \n      const ox = zoom * (2/width * px - 1);\n      const index = 4 * (py*width + px);\n      let [x, y] = orthographic ? ortho_to_stereo(ox, oy) : [ox, oy];\n      [x, y] = stereo_rotate(rotation, [x, y]);\n      data1.set(colorxy(quantize_stereo(x, y)), index);\n    }\n  }\n  return new ImageData(data1, width, height);\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1231,"value":"data2 = new Uint8ClampedArray(2*size[0] * 2*size[1] * 4);","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":265,"value":"// southern hemisphere\nimageData2 = {\n  let [width, height] = size;\n  width *= 2; height *= 2;\n  \n  for (let py = 0; py < height; py++) {\n    const oy = -zoom * (2/height * py - 1);\n    for (let px = 0; px < width; px++) {  \n      const ox = zoom * (2/width * px - 1);\n      const index = 4 * (py*width + px);\n      let [x, y] = orthographic ? ortho_to_stereo(ox, oy) : [ox, oy];\n      [x, y] = [-x/(x*x+y*y), y/(x*x+y*y)];\n      if (x+y !== x+y) [x, y] = [Infinity, 0];\n      [x, y] = stereo_rotate(rotation, [x, y]);\n      data2.set(colorxy(quantize_stereo(x, y)), index);\n    }\n  }\n  return new ImageData(data2, width, height);\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":444,"value":"$34","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":514,"value":"viewof grouped_colors = {\n  const w = Math.min(60, (width / 12) | 0);\n  const groups = [\n    [ \"#00b9d4\",\"#009f90\",\"#42cdff\",\"#00b8b5\",\"#0185cb\",\"#00d1de\",\n      \"#00b69a\",\"#008aa1\",\"#00d0bf\",\"#009fd2\",\"#008988\",\"#4baeff\"],\n    [ \"#c863f2\",\"#7768de\",\"#f075ff\",\"#a372fb\",\"#b535ac\",\"#c98dff\",\n      \"#7384f8\",\"#9c3dbb\",\"#a499ff\",\"#c349cc\",\"#7d4bc2\",\"#eb57da\"],\n    [ \"#f85f6f\",\"#cf3482\",\"#ff9386\",\"#f14c89\",\"#c05635\",\"#ff8096\",\n      \"#e442a2\",\"#c03e4b\",\"#ff68a9\",\"#dc5d4f\",\"#ba2a60\",\"#f47d52\"],\n    [ \"#6ec24a\",\"#97880e\",\"#56e66e\",\"#9cab2d\",\"#00a059\",\"#99d04b\",\n      \"#b89222\",\"#45952d\",\"#c0b636\",\"#24b650\",\"#708204\",\"#00d37d\"]];\n  const view = md`## Colors:\n  ${\n    groups.map(list => html`${\n      list.map(c => html``)\n    }`)\n  }`;\n  view.value = groups.map(list => list.map(from_hex));\n  return view;\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1009,"value":"// d3 gives rotation angles anticlockwise; we are using reverse\n// rotations here because we start with a position in image space\n// and want to find out what place on the sphere it came from\nrotation = rotate_d3geo_to_quaternion(-rot_xy, -rot_xz, -rot_yz)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":399,"value":"// count trailing zeros\nctz32 = function ctz32(k) {\n  // fill all bits to the left of the rightmost 1 bit with 1s\n  k |= k << 16, k |= k << 8, k |= k << 4, k |= k << 2, k |= k << 1;\n\n  // subtract count of leading 1s from 32\n  return 32 - Math.clz32(~k);\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":102,"value":"prev_pow2_x = (x) =>\n  2**(Math.floor(Math.log2(x)))","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":927,"value":"// return the largest power of 2 less than or equal to x\nprev_pow2 = {\n  const\n    HI = little_endian, LO = 1 - HI,\n    EMASK = 0x7ff00000, TOP = 0xfff00000,\n    float_array = new Float64Array([1.0]),\n    int_array = new Int32Array(float_array.buffer);\n  return function (x) {\n    float_array[0] = x;\n    if ((int_array[HI] & EMASK) === EMASK) return x; // inf or nan\n    if ((int_array[HI] & EMASK) === 0) // 0 or subnormal\n      return (x === 0) ? x : 2**(Math.floor(Math.log2(x)));\n    int_array[HI] &= TOP, int_array[LO] = 0; // zero out mantissa\n    return float_array[0];\n  }\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":8,"value":"ortho_to_stereo = (X, Y) => {\n  let Q = 1 + Math.sqrt(1 - X*X - Y*Y);\n  if (Q !== Q) Q = Math.sqrt(X*X + Y*Y); // clamp to circle\n  return [X/Q, Y/Q];\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1003,"value":"stereo_rotate = function stereo_rotate(quaternion, stereovector) {\n  const [x, y] = stereovector, [r1, rxy, ryz, rzx] = quaternion;\n  let px = r1*x + rxy*y - rzx, py = r1*y - rxy*x + ryz,\n      q1 = rzx*x - ryz*y + r1, qxy = rzx*y + ryz*x + rxy;\n  if (x*x + y*y === Infinity) px = r1, py = -rxy, q1 = rzx, qxy = ryz;\n  if (q1*q1 + qxy*qxy === 0) return [Infinity, 0];\n  return [(px*q1 + py*qxy) / (q1*q1 + qxy*qxy),\n          (py*q1 - px*qxy) / (q1*q1 + qxy*qxy)];\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1005,"value":"quatmul = function quatmul(quaternionR, quaternionS) {\n  const\n    [r1, rxy, ryz, rzx] = quaternionR,\n    [s1, sxy, syz, szx] = quaternionS;\n\n  // TODO: maybe consider edge cases with very large/small exponents or infinity\n  return [\n    r1*s1 - rxy*sxy - ryz*syz - rzx*szx,\n    r1*sxy + rxy*s1 - ryz*szx + rzx*syz,\n    r1*syz + rxy*szx + ryz*s1 - rzx*sxy,\n    r1*szx - rxy*syz + ryz*sxy + rzx*s1];\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1007,"value":"rotate_d3geo_to_quaternion = function rotate_d3geo_to_quaternion(xy_angle, xz_angle, yz_angle){\n  const halfrad = Math.PI / 360;\n  xy_angle *= halfrad, yz_angle *= halfrad, xz_angle *= halfrad;\n  const\n    R2 = [Math.cos(xy_angle), -Math.sin(xy_angle), 0, 0],\n    R1 = [Math.cos(xz_angle), 0, 0, Math.sin(xz_angle)],\n    R0 = [Math.cos(yz_angle), 0, -Math.sin(yz_angle), 0];\n  return quatmul(R0, quatmul(R1, R2));\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":32,"value":"// a previous version used a coloring method with a scattered set of\n// distinct colors instead of separate precision bands of different hues\nold_colors = [\n  \"#00aef6\",\"#c72f8b\",\"#babf3b\",\"#00a7b8\",\"#ff87ea\",\"#cc962d\",\"#00988e\",\"#de6dff\",\n  \"#ce6e36\",\"#00ebbb\",\"#aa71fd\",\"#c44945\",\"#00db85\",\"#6e7aeb\",\"#ff97ab\",\"#5ebf4b\",\n  \"#0281ca\",\"#ff61b0\",\"#889b20\",\"#00d8fa\",\"#e149c0\",\"#997601\",\"#00cbc6\",\"#b846c8\",\n  \"#ffa057\",\"#00bb9b\",\"#8b4ac6\",\"#fd7560\",\"#00aa70\",\"#aaafff\",\"#ed5271\",\"#279736\"\n].map(from_hex)","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":5,"value":"size = [width/2, width/2]","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":43,"value":"zoom = 1.02","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":30,"value":"from_hex = function from_hex(hex) {\n  if (!/^#?[0-9a-fA-F]{6}$/.test(hex)) {\n    throw new Error('Bad Input: Must be of form \"666FAD\" or \"#DEFACE\"'); }\n  const RGB = parseInt(hex.substr(-6), 16);\n  return [\n    (RGB >> 16)      ,  // first byte  -> R\n    (RGB >> 8 & 0xff),  // second byte -> G\n    (RGB & 0xff)     ,  // third byte  -> B\n    255 ];              // fourth byte -> alpha\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":159,"value":"import {Range, Toggle} from '@observablehq/inputs'","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":446,"value":"import {little_endian} from '@jrus/misc'","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1116,"value":"import {print_hexfloat} from '@jrus/hexfloat'","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1160,"value":"import {trackPointer} from '@fil/pointer-events'","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1163,"value":"d3 = require(\"d3-selection@2\")","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"}],"files":"$W35","stdlib":"1","titling":"auto"}}
30:"$2f:body"